matrixinverse.txt

通过CSharp写的关于矩阵求逆的程序

矩阵求逆的快速算法   龚敏敏　算法介绍
   矩阵求逆在3D程序中很常见，主要应用于求Billboard矩阵。按照定义的计算方法乘法运算，严重影响了性能。在需要大量Billboard矩阵运算时，矩阵求逆的优化能极大提高性能。这里要介绍的矩阵求逆算法称为全选主元高斯-约旦法。
高斯-约旦法（全选主元）求逆的步骤如下：
首先，对于 k 从 0 到 n - 1 作如下几步：
   1. 从第 k 行、第 k 列开始的右下角子阵中选取绝对值最大的元素，并记住次元素所在的行号和列号，在通过行交换和列交换将它交换到主元素位置上。这一步称为全选主元。
   2. m(k, k) = 1 / m(k, k)
   3. m(k, j) = m(k, j) * m(k, k)，j = 0, 1, ..., n-1；j != k
   4. m(i, j) = m(i, j) - m(i, k) * m(k, j)，i, j = 0, 1, ..., n-1；i, j != k
   5. m(i, k) = -m(i, k) * m(k, k)，i = 0, 1, ..., n-1；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;
}
比较　                 原算法                     原算法(经过高度优化)            新算法
加法次数     103                            61                            39
乘法次数     170                            116                     69
需要额外空间   16 * sizeof(float)         34 * sizeof(float)         25 * sizeof(float)
结果不言而喻吧。


//***************************
//求任何一个矩阵的逆矩阵
//***************************
#include <stdio.h>
#include <malloc.h>
void main( void )
{    float *buffer,*p;   //定义数组首地址指针变量
     short int row,num; //定义矩阵行数row及矩阵元素个数
     short int i,j;
     float determ;      //定义矩阵的行列式
     float comput_D(float *p,short int n);      //求矩阵的行列式
     float Creat_M(float *p, short int m,short int n,short int k); //求代数余子式
     void Print( float *p,short int n);     //打印n×n的矩阵
     printf("\nPlease input the number of rows: ");
     scanf("%d",&row);    
     num=2 * row * row;
     buffer = (float *)calloc(num, sizeof(float));     //分配内存单元
     p=buffer;
     if(p != NULL)
     { for(i=0;i<row;i++)                   //输入各单元值
         {
             printf("Input the number of %d row ",i+1);
             for(j=0;j<row;j++)
             {scanf("%f",p++);  }
         }  
     }
     else
         printf( "Can't allocate memory\n" );
     printf("\nThe original matrix is:\n");
     Print(buffer,row);     //打印该矩阵
     determ=comput_D(buffer,row);     //求整个矩阵的行列式
     p=buffer + row * row;
     if (determ != 0)
     {
         for (i=0;i<row; i++)       //求逆矩阵
             for (j=0; j<row; j++)
                    *(p+j*row+i)=   Creat_M(buffer,i,j,row)/determ;               
         printf("The determinant is %G\n",determ);
         p=buffer + row * row;
         printf("\nThe inverse matrix is:\n"); 
         Print(p,row);     //打印该矩阵
     }
     else
         printf("The determnant is 0, and there is no inverse matrix !\n");
     free( buffer );
}
//--------------------------------------------------------
//功能：求矩阵 n X n 的行列式
//入口参数：矩阵首地址 p；矩阵行数 n
//返回值：矩阵的行列式值
//--------------------------------------------------------
float comput_D(float *p,short int n)  
{    short int i,j,m;         //i--row; j--column
     short int lop=0;
     float result=0;
     float mid=1;    
     if (n!=1)
     {
         lop=(n==2)?1:n;     //控制求和循环次数，若为2阶，则循环1次，否则为n次
         for(m=0;m<lop;m++)
         {    mid=1;          //顺序求和
             for(i=0,j=m;i<n;i++,j++)
                 mid = mid * ( *(p+i*n+j%n) );
             result+=mid;
         }
         for(m=0;m<lop;m++)
         {                       
             mid=1;          //逆序相减
             for(i=0,j=n-1-m+n; i<n; i++,j--)
                 mid=mid * ( *(p+i*n+j%n));
             result-=mid;
         }
        }
     else result=*p;
     return(result);
}
//----------------------------------------------------
//功能：求k×k矩阵中元素A(mn)的代数余子式
//入口参数：k×k矩阵首地址；元素A的下标m,n; 矩阵行数 k
//返回值： k×k矩阵中元素A(mn)的代数余子式
//----------------------------------------------------
float Creat_M(float *p, short int m,short int n,short int k)
{
     short int len;
     short int i,j;
     float mid_result=0;
     short int quo=1;
     float *p_creat,*p_mid;
     len=(k-1)*(k-1);
     p_creat = (float *)calloc(len, sizeof(float));     //分配内存单元
     p_mid=p_creat;
     for(i=0;i<k;i++)
         for(j=0;j<k;j++)
         {
             if (i!=m && j!=n)
                 *p_mid++ =* (p+i*k+j);            
         } 
     //     Print(p_creat,k-1);
     quo = (m + n) %2==0 ? 1:-1; 
     mid_result = (float ) quo * comput_D(p_creat,k-1);
     free(p_creat);
     return(mid_result);
}
//-------------------------------------------
//功能：打印n×n的矩阵
//入口参数：n×n矩阵的首地址；该矩阵的行数 n
//返回值： 无
//-------------------------------------------
void Print( float *p,short int n)    
{
     int i,j;
     for (i=0;i<n;i++)
     {
         for (j=0; j<n;j++)
             printf("%10G ",*p++);
         printf("\n");
     }
     printf("--------------\n");
}