矩阵计算问题(50分)

矩阵计算问题
1、计算矩阵A的转置矩阵；
2、计算矩阵A和其转置矩阵的积；
3、计算第2步的逆矩阵；
4、计算矩阵A和矩阵B的积；
5、计算第3步的逆矩阵与第4步积的积；
即：
[a u]=（AT * A）-1*AT*B
AT:矩阵A的转置矩阵
-1:逆矩阵

DELPHI EXE/DLL 原码

矩阵A和矩阵B，形式如下：

| 159 1 |
| 234 1 |
矩阵A= | 83 1 |
| 6328.5 1 |
| -753.6 1 |



| 8386 |
| 13654 |
矩阵B= | 21202 |
| 37863 |
| 74761 |

查算法书

考，这是个反演问题。
下面是我编的FORTRAN代码，矩阵的转置和矩阵的相乘以及矩阵与向量的相乘，你可以改成Delphi的
关于求逆实际上可以转化为消元法求解一个线性方程组，我想B可能是个向量，如果是个矩阵则意味着求解多个线性方程组。

!矩阵的转置
Subroutine At(A,Ans,Ir,Ic,M,N)
Integer,Intent(in)::Ir,Ic,M,N
Real(8),Intent(In),Dimension(Ir:M,Ic:N)::A
Real(8),Intent(Out),Dimension(Ic:N,Ir:M)::Ans
Integer i,j
Do i=Ic,N
Do j=Ir,M
Ans(i,j)=A(j,i)
Enddo
Enddo
EndSubroutine At
!矩阵与矩阵的相乘
Subroutine AMultB(A,B,C,Ir,Is,Ic,m,n,l)
Integer,Intent(in)::m,n,l,Ir,Is,Ic
Real(8),Dimension(Ir:m,Is:n),Intent(in)::A
Real(8),Dimension(Is:n,Ic:l),Intent(in)::b
Real(8),Dimension(Ir:m,Ic:l),Intent(Out)::C
Integer i,j,k
Do i=Ir,m
do j=Ic,l
c(i,j)=0
Do k=Is,n
c(i,j)=c(i,j)+a(i,k)*b(k,j)
Enddo
Enddo
Enddo
EndSUbroutine AmultB

!矩阵和向量相乘
Subroutine ABMult(A,B,C,Ir,Ic,m,n)
Integer,Intent(in)::m,n,Ir,Ic
Real(8),Dimension(Ir:m,Ic:n),Intent(in)::A
Real(8),Dimension(Ic:n),Intent(in)::b
Real(8),Dimension(Ir:m),Intent(Out)::C
Integer i,j,k
Do i=Ir,m
c(i)=0
Do k=Ic,n
c(i)=c(i)+a(i,k)*b(k)
Enddo
Enddo
EndSUbroutine ABmult

这里的矩阵A是5*2矩阵,矩阵B是5*1矩阵,我在网上找到了C的矩阵算法，请高手变为delphi的好吗？

double * MatrixOpp(double A[],int m,int n) /*矩阵求逆*/
{
int i,j,x,y,k

double *SP=NULL,*AB=NULL,*B=NULL,X,*C

SP=(double *)malloc(m*n*sizeof(double))

AB=(double *)malloc(m*n*sizeof(double))

B=(double *)malloc(m*n*sizeof(double))

X=Surplus(A,m,n)

X=1/X

for(i=0;i<m;i++)
for(j=0;j<n;j++)
{for(k=0;k<m*n;k++)
B[k]=A[k]

{for(x=0;x<n;x++)
B[i*n+x]=0

for(y=0;y<m;y++)
B[m*y+j]=0

B[i*n+j]=1

SP[i*n+j]=Surplus(B,m,n)

AB[i*n+j]=X*SP[i*n+j]

}
}
C=MatrixInver(AB,m,n)

return C

}
double * MatrixInver(double A[],int m,int n) /*矩阵转置*/
{
int i,j

double *B=NULL

B=(double *)malloc(m*n*sizeof(double))

for(i=0;i<n;i++)
for(j=0;j<m;j++)
B[i*m+j]=A[j*n+i]

return B

}


矩阵相乘的快速算法
算法介绍
矩阵相乘在进行3D变换的时候是经常用到的。在应用中常用矩阵相乘的定义算法对其进行计算。这个算法用到了大量的循环和相乘运算，这使得算法效率不高。而矩阵相乘的计算效率很大程度上的影响了整个程序的运行速度，所以对矩阵相乘算法进行一些改进是必要的。

这里要介绍的矩阵算法称为斯特拉森方法，它是由v.斯特拉森在1969年提出的一个方法。

我们先讨论二阶矩阵的计算方法。

对于二阶矩阵

a11 a12 b11 b12
A = a21 a22 B = b21 b22

先计算下面7个量(1)

x1 = (a11 + a22) * (b11 + b22)

x2 = (a21 + a22) * b11

x3 = a11 * (b12 - b22)

x4 = a22 * (b21 - b11)

x5 = (a11 + a12) * b22

x6 = (a21 - a11) * (b11 + b12)

x7 = (a12 - a22) * (b21 + b22)


再设C = AB。根据矩阵相乘的规则，C的各元素为(2)

c11 = a11 * b11 + a12 * b21
c12 = a11 * b12 + a12 * b22
c21 = a21 * b11 + a22 * b21
c22 = a21 * b12 + a22 * b22

比较(1)(2)，C的各元素可以表示为(3)

c11 = x1 + x4 - x5 + x7
c12 = x3 + x5
c21 = x2 + x4
c22 = x1 + x3 - x2 + x6

根据以上的方法，我们就可以计算4阶矩阵了，先将4阶矩阵A和B划分成四块2阶矩阵，分别利用公式计算它们的乘积，再使用(1)(3)来计算出最后结果。

ma11 ma12 mb11 mb12
A4 = ma21 ma22 B4 = mb21 mb22

其中

a11 a12 a13 a14 b11 b12 b13 b14
ma11 = a21 a22 ma12 = a23 a24 mb11 = b21 b22 mb12 = b23 b24

a31 a32 a33 a34 b31 b32 b33 b34
ma21 = a41 a42 ma22 = a43 a44 mb21 = b41 b42 mb22 = b43 b44

实现
// 计算2X2矩阵
void Multiply2X2(float&amp
fOut_11, float&amp
fOut_12, float&amp
fOut_21, float&amp
fOut_22,
float f1_11, float f1_12, float f1_21, float f1_22,
float f2_11, float f2_12, float f2_21, float f2_22)
{
const float x1((f1_11 + f1_22) * (f2_11 + f2_22))

const float x2((f1_21 + f1_22) * f2_11)

const float x3(f1_11 * (f2_12 - f2_22))

const float x4(f1_22 * (f2_21 - f2_11))

const float x5((f1_11 + f1_12) * f2_22)

const float x6((f1_21 - f1_11) * (f2_11 + f2_12))

const float x7((f1_12 - f1_22) * (f2_21 + f2_22))


fOut_11 = x1 + x4 - x5 + x7

fOut_12 = x3 + x5

fOut_21 = x2 + x4

fOut_22 = x1 - x2 + x3 + x6

}

// 计算4X4矩阵
void Multiply(CLAYMATRIX&amp
mOut, const CLAYMATRIX&amp
m1, const CLAYMATRIX&amp
m2)
{
float fTmp[7][4]


// (ma11 + ma22) * (mb11 + mb22)
Multiply2X2(fTmp[0][0], fTmp[0][1], fTmp[0][2], fTmp[0][3],
m1._11 + m1._33, m1._12 + m1._34, m1._21 + m1._43, m1._22 + m1._44,
m2._11 + m2._33, m2._12 + m2._34, m2._21 + m2._43, m2._22 + m2._44)


// (ma21 + ma22) * mb11
Multiply2X2(fTmp[1][0], fTmp[1][1], fTmp[1][2], fTmp[1][3],
m1._31 + m1._33, m1._32 + m1._34, m1._41 + m1._43, m1._42 + m1._44,
m2._11, m2._12, m2._21, m2._22)


// ma11 * (mb12 - mb22)
Multiply2X2(fTmp[2][0], fTmp[2][1], fTmp[2][2], fTmp[2][3],
m1._11, m1._12, m1._21, m1._22,
m2._13 - m2._33, m2._14 - m2._34, m2._23 - m2._43, m2._24 - m2._44)


// ma22 * (mb21 - mb11)
Multiply2X2(fTmp[3][0], fTmp[3][1], fTmp[3][2], fTmp[3][3],
m1._33, m1._34, m1._43, m1._44,
m2._31 - m2._11, m2._32 - m2._12, m2._41 - m2._21, m2._42 - m2._22)


// (ma11 + ma12) * mb22
Multiply2X2(fTmp[4][0], fTmp[4][1], fTmp[4][2], fTmp[4][3],
m1._11 + m1._13, m1._12 + m1._14, m1._21 + m1._23, m1._22 + m1._24,
m2._33, m2._34, m2._43, m2._44)


// (ma21 - ma11) * (mb11 + mb12)
Multiply2X2(fTmp[5][0], fTmp[5][1], fTmp[5][2], fTmp[5][3],
m1._31 - m1._11, m1._32 - m1._12, m1._41 - m1._21, m1._42 - m1._22,
m2._11 + m2._13, m2._12 + m2._14, m2._21 + m2._23, m2._22 + m2._24)


// (ma12 - ma22) * (mb21 + mb22)
Multiply2X2(fTmp[6][0], fTmp[6][1], fTmp[6][2], fTmp[6][3],
m1._13 - m1._33, m1._14 - m1._34, m1._23 - m1._43, m1._24 - m1._44,
m2._31 + m2._33, m2._32 + m2._34, m2._41 + m2._43, m2._42 + m2._44)


// 第一块
mOut._11 = fTmp[0][0] + fTmp[3][0] - fTmp[4][0] + fTmp[6][0]

mOut._12 = fTmp[0][1] + fTmp[3][1] - fTmp[4][1] + fTmp[6][1]

mOut._21 = fTmp[0][2] + fTmp[3][2] - fTmp[4][2] + fTmp[6][2]

mOut._22 = fTmp[0][3] + fTmp[3][3] - fTmp[4][3] + fTmp[6][3]


// 第二块
mOut._13 = fTmp[2][0] + fTmp[4][0]

mOut._14 = fTmp[2][1] + fTmp[4][1]

mOut._23 = fTmp[2][2] + fTmp[4][2]

mOut._24 = fTmp[2][3] + fTmp[4][3]


// 第三块
mOut._31 = fTmp[1][0] + fTmp[3][0]

mOut._32 = fTmp[1][1] + fTmp[3][1]

mOut._41 = fTmp[1][2] + fTmp[3][2]

mOut._42 = fTmp[1][3] + fTmp[3][3]


// 第四块
mOut._33 = fTmp[0][0] - fTmp[1][0] + fTmp[2][0] + fTmp[5][0]

mOut._34 = fTmp[0][1] - fTmp[1][1] + fTmp[2][1] + fTmp[5][1]

mOut._43 = fTmp[0][2] - fTmp[1][2] + fTmp[2][2] + fTmp[5][2]

mOut._44 = fTmp[0][3] - fTmp[1][3] + fTmp[2][3] + fTmp[5][3]

}

比较
在标准的定义算法中我们需要进行n * n * n次乘法运算，新算法中我们需要进行7log2n次乘法，对于最常用的4阶矩阵： 　 原算法 新算法
加法次数 48 72(48次加法，24次减法)
乘法次数 64 49
需要额外空间 16 * sizeof(float) 28 * sizeof(float)

新算法要比原算法多了24次减法运算，少了15次乘法。但因为浮点乘法的运算速度要远远慢于加/减法运算，所以新算法的整体速度有所提高。


矩阵求逆的快速算法



算法介绍
矩阵求逆在3D程序中很常见，主要应用于求Billboard矩阵。按照定义的计算方法乘法运算，严重影响了性能。在需要大量Billboard矩阵运算时，矩阵求逆的优化能极大提高性能。这里要介绍的矩阵求逆算法称为全选主元高斯-约旦法。

高斯-约旦法（全选主元）求逆的步骤如下：

首先，对于 k 从 0 到 n - 1 作如下几步：

从第 k 行、第 k 列开始的右下角子阵中选取绝对值最大的元素，并记住次元素所在的行号和列号，在通过行交换和列交换将它交换到主元素位置上。这一步称为全选主元。
m(k, k) = 1 / m(k, k)
m(k, j) = m(k, j) * m(k, k)，j = 0, 1, ..., n-1；j != k
m(i, j) = m(i, j) - m(i, k) * m(k, j)，i, j = 0, 1, ..., n-1；i, j != k
m(i, k) = -m(i, k) * m(k, k)，i = 0, 1, ..., n-1；i != k
最后，根据在全选主元过程中所记录的行、列交换的信息进行恢复，恢复的原则如下：在全选主元过程中，先交换的行（列）后进行恢复；原来的行（列）交换用列（行）交换来恢复。

实现(4阶矩阵)
float Inverse(CLAYMATRIX&amp
mOut, const CLAYMATRIX&amp
rhs)
{
CLAYMATRIX m(rhs)

DWORD is[4]

DWORD js[4]

float fDet = 1.0f

int f = 1


for (int k = 0
k < 4
k ++)
{
// 第一步，全选主元
float fMax = 0.0f

for (DWORD i = k
i < 4
i ++)
{
for (DWORD j = k
j < 4
j ++)
{
const float f = Abs(m(i, j))

if (f > fMax)
{
fMax = f

is[k] = i

js[k] = j

}
}
}
if (Abs(fMax) < 0.0001f)
return 0


if (is[k] != k)
{
f = -f

swap(m(k, 0), m(is[k], 0))

swap(m(k, 1), m(is[k], 1))

swap(m(k, 2), m(is[k], 2))

swap(m(k, 3), m(is[k], 3))

}
if (js[k] != k)
{
f = -f

swap(m(0, k), m(0, js[k]))

swap(m(1, k), m(1, js[k]))

swap(m(2, k), m(2, js[k]))

swap(m(3, k), m(3, js[k]))

}

// 计算行列值
fDet *= m(k, k)


// 计算逆矩阵

// 第二步
m(k, k) = 1.0f / m(k, k)

// 第三步
for (DWORD j = 0
j < 4
j ++)
{
if (j != k)
m(k, j) *= m(k, k)

}
// 第四步
for (DWORD i = 0
i < 4
i ++)
{
if (i != k)
{
for (j = 0
j < 4
j ++)
{
if (j != k)
m(i, j) = m(i, j) - m(i, k) * m(k, j)

}
}
}
// 第五步
for (i = 0
i < 4
i ++)
{
if (i != k)
m(i, k) *= -m(k, k)

}
}

for (k = 3
k >= 0
k --)
{
if (js[k] != k)
{
swap(m(k, 0), m(js[k], 0))

swap(m(k, 1), m(js[k], 1))

swap(m(k, 2), m(js[k], 2))

swap(m(k, 3), m(js[k], 3))

}
if (is[k] != k)
{
swap(m(0, k), m(0, is[k]))

swap(m(1, k), m(1, is[k]))

swap(m(2, k), m(2, is[k]))

swap(m(3, k), m(3, is[k]))

}
}

mOut = m

return fDet * f

}


比较
　 原算法 原算法(经过高度优化) 新算法
加法次数 103 61 39
乘法次数 170 116 69
需要额外空间 16 * sizeof(float) 34 * sizeof(float) 25 * sizeof(float)

结果不言而喻吧。

你要搞得那么麻烦只能靠你自己，别人不可能给你花时间的。

矩阵计算中只要知道乘法和求逆的算法就可以解决这个问题
乘法的原则是：
用前面的矩阵的行元素乘以后面矩阵的列元素将他们的各个乘积相加，以前矩阵的行号，后矩阵的列号为行列号，相加的结果就放在这个行列号码所在矩阵的行列号上
求逆的方法很多，但都很麻烦：
可以求原方阵的n－1阶子式的行列式的值来构成过渡阵，但是非常复杂
如果用初等矩阵变换，算法难以实现，我还没有好办法～