前段时间写到一个LDA算法，要用到求解矩阵的逆。之前写的代码都是在网上Copy过来的，今天终于把算法导论上的求解逆矩阵看完了，这里代码贴出来分享。。。
用的方法是LUP分解求解矩阵<注：代码在VS2008中编译通过且结果正确>

#include <iostream>
#include <math.h>
#include <stdio.h>
#include <malloc.h>
#include <assert.h>
using namespace std;

#define PAUSE system("PAUSE");
#define Malloc(type,n) (type *)malloc((n)*sizeof(type))

//LU分解
void LU_DECOMPOSITION(const int row,double **Matrix)
{
assert(NULL!=Matrix);
int n = row;
int k,i,j;
for (k=0;k<n;k++)
{
   double pivot = Matrix[k][k];
   for (i =k+1;i<n;i++)
   {
    Matrix[i][k] /= pivot;
    for (j = k+1;j<n;j++)
     Matrix[i][j] -= Matrix[i][k]*Matrix[k][j];
   }
}
}

//LUP分解,返回一个置换矩阵
int *LUP_DECOMPOSITION(const int row,double **Matrix)
{
assert(NULL!=Matrix);
int kp;
int n = row;
int k,i,j;

int *Pai = Malloc(int,row);
if (NULL ==Pai) cerr<<"Pai Malloc ERROR";

for (i=0;i<n;i++)
   Pai[i] = i;

for (k = 0;k<n;k++)
{
   double p = 0;
   for (i = k;i<n;i++)
    if (fabs(Matrix[i][k])>p)
    {
     p = fabs(Matrix[i][k]);
     kp = i;
    }
   if (0==p)
   {
    cerr<<"singular Matrix";
    PAUSE
   }

   swap(Pai[k],Pai[kp]);

   for (i = 0;i<n;i++)
    swap(Matrix[k][i],Matrix[kp][i]);

   for (i =k+1;i<n;i++)
   {
    Matrix[i][k] /= Matrix[k][k];
    for (j = k+1;j<n;j++)
     Matrix[i][j] -= Matrix[i][k]*Matrix[k][j];
   }
}
return Pai;
}

//LUP计算返回x解
double *LUP_SOLVE(const int row,double **Matrix,double *b,int *Pai)
{
assert((NULL!=Matrix)||(NULL!=b)||(NULL!=Pai));
int n = row;
int i,j;

double *y = Malloc(double,row);
double *x = Malloc(double,row);
if ((NULL==x)||(NULL==y)) cerr<<"LUP_SOLVE xMalloc ERROR";

for (i =0;i<row;i++)
{
   x[i] = 0;
   y[i] = 0;
}

for (i = 0;i<n;i++)
{
   double S1 = 0;
   for (j = 0;j<i;j++)
    S1 +=Matrix[i][j]*y[j];

   y[i] = b[Pai[i]] - S1;
}

for (i=n-1;i>=0;i--)
{
   double S2 = 0;
   for (j = i;j<n-1;j++)
    S2 += Matrix[i][j+1]*x[j+1];

   x[i] = (y[i] - S2)/Matrix[i][i];
}
return x;

}

//用LUP分解算法来求解逆矩阵
void main()
{
const int row = 4;
double Matrix[4][4]={{0,0,0,0},{6,13,5,19},{2,19,10,23},{4,10,11,31}};//{2,3,1,5}
double INV[4][4]={0};
//double Matrix[4][4]={{2,0,2,0.6},{3,3,4,-2},{5,5,4,2},{-1,-2,3.4,-1}};
//double Matrix[3][3]={{1,2,0},{3,4,4},{5,6,3}};
//double b[3] = {0,0,1};
double **y = Malloc(double *,row);
double *x;
int i,j;
for (i=0;i<row;i++)
{
   y[i] = &Matrix[i][0];
}//绑定地址

int *Pai = LUP_DECOMPOSITION(row,y);
for (i=0;i<row;i++)
{
   double *b =Malloc(double,row);
   for (j=0;j<row;j++)
    b[j] = 0;
   b[i] = 1;

   x = LUP_SOLVE(row,y,b,Pai);
   for (j = 0;j<row;j++)
    INV[j][i] = x[j];
}

double *ptr = &INV[0][0];
for (i=0;i<row*row;i++)
{
   cout<<*ptr++<<endl;
}

//cout<<"\n";
//for (i=0;i<row;i++)
//{
// cout<<x[i]<<endl;
//}

PAUSE
}