矩阵的求逆

#include   
#include   
#include   
  
#define MAX 20  
#define E 0.000000001  
  
/** 
 * 计算矩阵src的模 
 */  
double calculate_A( double src[][MAX], int n )  
{  
    int i,j,k,x,y;  
    double tmp[MAX][MAX], t;  
    double result = 0.0;  
      
    if( n == 1 )  
    {  
        return src[0][0];  
    }  
      
    for( i = 0; i < n; ++i )  
    {  
        for( j = 0; j < n - 1; ++j )  
        {  
            for( k = 0; k < n - 1; ++k )  
            {  
                x = j + 1;  
                y = k >= i   k + 1 : k;  
                  
                tmp[j][k] = src[x][y];  
            }  
        }  
          
        t = calculate_A( tmp, n - 1 );  
          
        if( i % 2 == 0 )  
        {  
            result += src[0][i] * t;  
        }  
        else  
        {  
            result -= src[0][i] * t;  
        }  
    }  
  
    return result;  
}  
  
/** 
 * 计算伴随矩阵 
 */  
void calculate_A_adjoint( double src[MAX][MAX], double dst[MAX][MAX], int n )  
{  
    int i, j, k, t, x, y;  
    double tmp[MAX][MAX];  
  
    if( n == 1 )  
    {  
        dst[0][0] = 1;  
        return;  
    }  
      
    for( i = 0; i < n; ++i )  
    {  
        for( j = 0; j < n; ++j )  
        {  
            for( k = 0; k < n - 1; ++k )  
            {  
                for( t = 0; t < n - 1; ++t )  
                {  
                    x = k >= i   k + 1 : k ;  
                    y = t >= j   t + 1 : t;  
                      
                    tmp[k][t] = src[x][y];  
                }  
            }  
              
            dst[j][i]  =  calculate_A( tmp, n - 1 );  
              
            if( ( i + j ) % 2 == 1 )  
            {  
                dst[j][i] = -1*dst[j][i];  
            }  
        }  
    }  
}  
  
/** 
 * 得到逆矩阵 
 */  
int calculate_A_inverse( double src[MAX][MAX], double dst[MAX][MAX], int n )  
{  
    double A = calculate_A( src, n );  
    double tmp[MAX][MAX];  
    int i, j;  
  
    if ( fabs( A - 0 ) <= E )  
    {  
        printf("不可能有逆矩阵！\n");  
        return 0;  
    }  
  
    calculate_A_adjoint( src, tmp, n );    
  
    for( i = 0; i < n; ++i )    
    {    
        for( j = 0; j < n; ++j )    
        {    
            dst[i][j] = (double)( tmp[i][j] / A );  
        }    
    }  
  
    return 1;  
}  
  
/** 
 * 输出矩形查看 
 */  
void print_M( double M[][MAX], int n )  
{  
    int i, j;  
      
    for ( i = 0; i < n; ++i )  
    {  
        for ( j = 0; j < n; ++j )  
        {  
            printf("%lf ", M[i][j]);  
        }  
          
        printf("\n");  
    }  
}  
  
/** 
 * main 
 */  
int main()  
{  
    /** 
     * 测试矩阵 
     */  
    double test[MAX][MAX], dst[MAX][MAX];  
    int n = 3;  
    int is_exist;  
      
    /** 
     * 构造最简单的： 
     *  1, 0, 0, 
     *  0, 2, 0, 
     *  0, 0, 5 
     */  
    memset( test, 0, sizeof( test ) );  
      
    test[0][0] = 1.0;  
    test[1][1] = 2.0;  
    test[2][2] = 5.0;  
      
      
    is_exist = calculate_A_inverse( test, dst, n );  
  
    if ( is_exist )  
    {  
        print_M(dst, n);  
    }  
    else  
    {  
        printf("不存在!\n");  
    }  
      
    return 0;  
}