Gauss-Jordan法实矩阵求逆
下面是实现Gauss-Jordan法实矩阵求逆。
#include
#include
#include
int brinv(double a[], int n)
{ int *is,*js,i,j,k,l,u,v;
double d,p;
is=malloc(n*sizeof(int));
js=malloc(n*sizeof(int));
for (k=0; k<=n-1; k++)
{ d=0.0;
for (i=k; i<=n-1; i++)
for (j=k; j<=n-1; j++)
{ l=i*n+j; p=fabs(a[l]);
if (p>d) { d=p; is[k]=i; js[k]=j;}
}
if (d+1.0==1.0)
{ free(is); free(js); printf("err**not inv\n");
return(0);
}
if (is[k]!=k)
for (j=0; j<=n-1; j++)
{ u=k*n+j; v=is[k]*n+j;
p=a[u]; a[u]=a[v]; a[v]=p;
}
if (js[k]!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+js[k];
p=a[u]; a[u]=a[v]; a[v]=p;
}
l=k*n+k;
a[l]=1.0/a[l];
for (j=0; j<=n-1; j++)
if (j!=k)
{ u=k*n+j; a[u]=a[u]*a[l];}
for (i=0; i<=n-1; i++)
if (i!=k)
for (j=0; j<=n-1; j++)
if (j!=k)
{ u=i*n+j;
a[u]=a[u]-a[i*n+k]*a[k*n+j];
}
for (i=0; i<=n-1; i++)
if (i!=k)
{ u=i*n+k; a[u]=-a[u]*a[l];}
}
for (k=n-1; k>=0; k--)
{ if (js[k]!=k)
for (j=0; j<=n-1; j++)
{ u=k*n+j; v=js[k]*n+j;
p=a[u]; a[u]=a[v]; a[v]=p;
}
if (is[k]!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+is[k];
p=a[u]; a[u]=a[v]; a[v]=p;
}
}
free(is); free(js);
return(1);
}
void brmul(double a[], double b[],int m,int n,int k,double c[])
{ int i,j,l,u;
for (i=0; i<=m-1; i++)
for (j=0; j<=k-1; j++)
{ u=i*k+j; c[u]=0.0;
for (l=0; l<=n-1; l++)
c[u]=c[u]+a[i*n+l]*b[l*k+j];
}
return;
}
int main()
{ int i,j;
static double a[4][4]={ {0.2368,0.2471,0.2568,1.2671},
{1.1161,0.1254,0.1397,0.1490},
{0.1582,1.1675,0.1768,0.1871},
{0.1968,0.2071,1.2168,0.2271}};
static double b[4][4],c[4][4];
for (i=0; i<=3; i++)
for (j=0; j<=3; j++)
b[i][j]=a[i][j];
i=brinv(a,4);
if (i!=0)
{ printf("MAT A IS:\n");
for (i=0; i<=3; i++)
{ for (j=0; j<=3; j++)
printf("%13.7e ",b[i][j]);
printf("\n");
}
printf("\n");
printf("MAT A- IS:\n");
for (i=0; i<=3; i++)
{ for (j=0; j<=3; j++)
printf("%13.7e ",a[i][j]);
printf("\n");
}
printf("\n");
printf("MAT AA- IS:\n");
brmul(b,a,4,4,4,c);
for (i=0; i<=3; i++)
{ for (j=0; j<=3; j++)
printf("%13.7e ",c[i][j]);
printf("\n");
}
}
}
矩阵求逆的快速算法
算法介绍
矩阵求逆在3D程序中很常见，主要应用于求Billboard矩阵。按照定义的计算
方法
快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载
乘法运算，严重影响了性能。在需要大量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(k, k)，i = 0, 1, ..., n-1;i != k m(i, k) = -
最后，根据在全选主元过程中所记录的行、列交换的信息进行恢复，恢复的原则如下:在全选主元过程中，先交换的行(列)后进行恢复;原来的行(列)交换用列(行)交换来恢复。
实现(4阶矩阵)
float Inverse(CLAYMATRIX& mOut, const CLAYMATRIX& 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;
}
浙公网安备 33021202002483