convertmatrix.txt

矩阵求逆是数学运算中常常遇到的。本源码采用全主元矩阵求逆的方法求得矩阵的逆。

int ConvertMatrix(double matrix[MAX_DIM][MAX_DIM], int n, double 
answer[MAX_DIM][MAX_DIM]) 
{ //全主元矩阵求逆 
   double A[MAX_DIM * MAX_DIM]; 
   int i, j, k, *nx, *ny, *temp1, *temp2; 
   k = 0; 
   for (i = 0; i < n; i ++ ) 
       for (j = 0; j < n; j ++ ) A[k ++] = matrix[j]; 

   double temp,mid; 
   temp1 = nx = new int [n];  //存放行交换后的排列顺序 
   temp2 = ny = new int [n];  //存放列交换后的排列顺序 
   for(i = 0; i < n; i ++ ) { 
       nx[ i ] = i;  //初始行顺序 
       ny[ i ] = i;  //初始列顺序 
   } 
   for(k = 0; k < n; k ++ ) {//找出主元所在行列 
       double max = 0; 
       for(i = k; i < n; i ++) { 
           for(j = k; j < n; j ++) { 
               mid = fabs(*(A + i * n + j)); 
               if(mid > max) { max = mid; nx[k] = i; ny[k] = j; } 
           } 
       } 
       if (max == 0) { //主元所在行列与第k行k列交换 
           printf("\n the main element be zero!\n"); 
           return 0; 
       } 
       if (nx[k] != k) { 
           for(j = 0; j < n; j ++ ) {  //行交换 
               temp = *(A + k * n + j); 
               *(A + k * n + j) = *(A + nx[k] * n + j); 
               *(A + nx[k] * n + j) = temp; 
           } 
       } 
       if (ny[k] != k) { 
           for(i = 0;i < n; i ++ ){ //列交换 
               temp = *(A+i*n+k); 
               *(A+i*n+k) = *(A+i*n+ny[k]); 
               *(A+i*n+ny[k]) = temp; 
           } 
       } 
       *(A + k * n + k) = 1 / (*(A + k * n + k));//k行所有元素除以A[k][k] 
       for (j = 0; j < n; j ++) 
           if (j != k) 
           *(A + k * n + j) = (*(A + k * n + j))*(*(A + k * n + k)); 
           //对其余各行进行消去运算 
           for(i = 0; i < n; i ++ ) 
               if(i != k) 
                   for(j = 0; j < n; j ++ ) 
                       if(j != k) 
                           *(A + i * n + j) = *(A + i * n + j) - (*(A + i 
* n + k)) * (*(A + k * n + j)); 
           //填入系数 
           for (i = 0; i < n; i++) 
               if (i != k) 
                   *(A + i * n + k) = -(*(A + i * n + k))*(*(A + k * n + 
k)); 
   } 
   //在选主元时若使用了行交换，要用列交换进行恢复。 
   //若使用了列交换，要用行交换来恢复。 
   //并且恢复的顺序刚好与选主元时交换的次序相反。 
   for (k = n - 1; k >= 0; k -- ) { 
       for (j = 0; j < n; j ++) { 
           if(ny[k] != k) {  //行交换 
               temp = *(A + k * n + j); 
               *(A + k * n + j) = *(A + ny[k] * n + j); 
               *(A + ny[k] * n + j) = temp; 
           } 
       } 
       for (i = 0; i < n; i ++ ) { 
           if (nx[k] != k) {  //列交换 
               temp = *(A + i * n + k); 
               *(A + i * n + k) = *(A + i * n + nx[k]); 
               *(A + i * n + nx[k]) = temp; 
           } 
       } 
   }; 
   for (i = 0; i < n; i ++ ) 
       for(j = 0; j < n; j ++ ) 
           answer [ i ][j] = *(A + i * n + j); 

   delete temp1; 
   delete temp2; 
   return 1; 
} 


double MatrixValue(double tmp[MAX_DIM][MAX_DIM], int nn) 
{ 
   int i, j, k; 
   double t1, t2; 
   double matrix[MAX_DIM][MAX_DIM]; 
   for (i = 0; i < nn; i ++ ) 
       for (j = 0; j < nn; j ++ ) matrix[ i ][j] = tmp[ i ][j]; 
   for (i = 0; i < nn - 1; i ++ ) { 
       t1 = matrix[ i ] [ i ]; 
       for (j = i + 1; j < nn; j ++ ) { 
           t2 = matrix[j][ i ]; 
           matrix[j][ i ] = 0; 
           for (k = i + 1; k < nn; k ++ ) 
               matrix[j][k] -= matrix[ i ][k] / t1 * t2; 
       } 
   } 
   double answer = 1; 
   for (i = 0; i < nn; i ++ ) answer *= matrix[ i ][ i ]; 
   return answer; 
}