Gauss-Jordan法實(shí)矩陣求逆

1、gauss-jordan法實(shí)矩陣求逆 下面是實(shí)現(xiàn)gauss-jordan法實(shí)矩陣求逆。 #include stdlib.h #include math.h #include stdio.h 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(al);

2、if (pd) d=p; isk=i; jsk=j; if (d+1.0=1.0) free(is); free(js); printf(err*not invn); return(0); if (isk!=k) for (j=0; j=n-1; j+) u=k*n+j; v=isk*n+j; p=au; au=av; av=p; if (jsk!=k) for (i=0; i=n-1; i+) u=i*n+k; v=i*n+jsk; p=au; au=av; av=p; l=k*n+k; al=1.0/al; for (j=0; j=n-1; j+) if (j!=k) u=k*n+j; a

3、u=au*al; for (i=0; i=n-1; i+) if (i!=k) for (j=0; j=n-1; j+) if (j!=k) u=i*n+j; au=au-ai*n+k*ak*n+j; for (i=0; i=n-1; i+) if (i!=k) u=i*n+k; au=-au*al; for (k=n-1; k=0; k-) if (jsk!=k) for (j=0; j=n-1; j+) u=k*n+j; v=jsk*n+j; p=au; au=av; av=p; if (isk!=k) for (i=0; i=n-1; i+) u=i*n+k; v=i*n+isk; p=

4、au; au=av; av=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; cu=0.0; for (l=0; l=n-1; l+) cu=cu+ai*n+l*bl*k+j; return; int main() int i,j; static double a44= 0.2368,0.2471,0.2568,1.2671, 1.11

5、61,0.1254,0.1397,0.1490, 0.1582,1.1675,0.1768,0.1871, 0.1968,0.2071,1.2168,0.2271; static double b44,c44; for (i=0; i=3; i+) for (j=0; j=3; j+) bij=aij; 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 ,bij); printf(n); printf(n); printf(mat a- is:n); f

6、or (i=0; i=3; i+) for (j=0; j=3; j+) printf(%13.7e ,aij); 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 ,cij); printf(n); 矩陣求逆的快速算法 算法介紹 矩陣求逆在3d程序中很常見(jiàn)，主要應(yīng)用于求billboard矩陣。根據(jù)定義的計(jì)算方法乘法運(yùn)算，嚴(yán)峻影響了性能。在需要大量billboard矩陣運(yùn)算時(shí)，矩陣求逆的優(yōu)化能極大提高性能。這里要介

7、紹的矩陣求逆算法稱為全選主元高斯-約旦法。 高斯-約旦法（全選主元）求逆的步驟如下： 首先，對(duì)于 k 從 0 到 n - 1 作如下幾步： 從第 k 行、第 k 列開頭的右下角子陣中選取肯定值最大的元素，并記住次元素所在的行號(hào)和列號(hào)，在通過(guò)行交換和列交換將它交換到主元素位置上。這一步稱為全選主元。 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)

8、= -m(i, k) * m(k, k)，i = 0, 1, ., n-1；i != k 最終，依據(jù)在全選主元過(guò)程中所記錄的行、列交換的信息進(jìn)行恢復(fù)，恢復(fù)的原則如下：在全選主元過(guò)程中，先交換的行（列）后進(jìn)行恢復(fù)；原來(lái)的行（列）交換用列（行）交換來(lái)恢復(fù)。 實(shí)現(xiàn)(4階矩陣) float inverse(claymatrix mout, const claymatrix rhs) claymatrix m(rhs); dword is4; dword js4; float fdet = 1.0f; int f = 1; for (int k = 0; k 4; k +) / 第一步，全選主元 flo

9、at 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; isk = i; jsk = j; if (abs(fmax) 0.0001f) return 0; if (isk != k) f = -f; swap(m(k, 0), m(isk, 0); swap(m(k, 1), m(isk, 1); swap(m(k, 2), m(isk, 2); swap(m(k, 3), m(isk, 3); if (

10、jsk != k) f = -f; swap(m(0, k), m(0, jsk); swap(m(1, k), m(1, jsk); swap(m(2, k), m(2, jsk); swap(m(3, k), m(3, jsk); / 計(jì)算行列值 fdet *= m(k, k); / 計(jì)算逆矩陣 / 其次步 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 (jsk != k) swap(m(k, 0), m(jsk, 0); swap(m(k, 1), m(jsk, 1); swap(m(k, 2), m(jsk, 2); swap