矩阵求你

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

#include
#include
#include
#include
#include
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{
double pivot = Matrix[k][k];
for (i =k+1;i{
Matrix[i][k] /= pivot;
for (j = k+1;jMatrix[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;iPai[i] = i;

for (k = 0;k{
double p = 0;
for (i = k;iif (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;iswap(Matrix[k][i],Matrix[kp][i]);

for (i =k+1;i{
Matrix[i][k] /= Matrix[k][k];
for (j = k+1;jMatrix[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{
x[i] = 0;
y[i] = 0;
}

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

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

for (i=n-1;i>=0;i--)
{
double S2 = 0;
for (j = i;jS2 += 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{
y[i] = &Matrix[i][0];
}//绑定地址

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

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


double *ptr = &INV[0][0];
for (i=0;i

{
cout<<*ptr++<}

//cout<<"\n";
//for (i=0;i//{
// cout<//}

PAUSE
}

///////////////////////////////////
/* ==================================
*
* Copyright (C) Bill Hsu
* /probill
* 2009-12-11
********************************** */
#include < iostream >
#include < vector >
#include < math.h >
using namespace std;

typedef vector < float > s_line; // 用来表示一行
s_line line;

typedef vector < s_line > s_matrix; // 用来表示一个矩阵
s_matrix matrix;
s_matrix mat;
int nSize; // 矩阵维数
int nSign; // 标记行列式正负
void outprint(s_matrix & _mat);
void printstep(s_matrix & _mat);
int step = 0 ;
void line_add(s_matrix & _mat, int a, int b, float k = 1.0f ) // 第b行乘k加到第a行
{
int size = _mat[ 0 ].size();

for ( int i = 0 ;i < size; ++ i)
{
_mat[a][i] += _mat[b][i] * k;

} // end for
}



void work1(s_matrix & _mat) // 主计算函数
{

for ( int i = 1 ;i < nSize; ++ i)
{

if (fabs(_mat[i - 1 ][i - 1 ]) < 0.000001 )
{
int mm;
for (mm = i;mm < nSize; ++ mm)
{
if (fabs(_mat[mm - 1 ][i - 1 ]) > 0.000001 ) break ;
} // end for
line_add(_mat,i - 1 ,mm - 1 );
} // end if

for ( int j = i;j < nSize; ++ j)
{
line_add(_mat,j,i - 1 , - _mat[j][i - 1 ] / _mat[i - 1 ][i - 1 ]);

} // end for j
printstep(_mat);
} // end for i

}


void work2(s_matrix & _mat) // 第二部计算
{
for ( int i = nSize - 2 ;i >= 0 ; -- i)
{
for ( int j = i;j >= 0 ; -- j)
{
line_add(_mat,j,i + 1 , - _mat[j][i + 1 ] / _mat[i + 1 ][i + 1 ]);
}
printstep(_mat);
}

}


void makeunit(s_matrix & _mat) // 单位化
{

mat.clear();

for ( int i = 0 ;i < nSize; ++ i)
{
line.clear();
for ( int j = 0 ;j < nSize * 2 ; ++ j)
{
float tmp = _mat[i][j] / _mat[i][i];
if (fabs(tmp) < 0.000001 ) tmp = 0 ;
line.push_back(tmp);
}
mat.push_back(line);
// cout<}
_mat = mat;
}

void printstep(s_matrix & _mat) // 显示求的过程
{
cout << " 第 " <<++ step << " 步 " << endl;
for ( int i = 0 ;i < nSize; ++ i)
{

for ( int j = 0 ;j < 2 * nSize; ++ j)
{
if (fabs(_mat[i][j]) < 0.000001 ) _mat[i][j] = 0 ;
cout << _mat[i][j] << " " ;
if (j == nSize - 1 )cout << " | " ;
}
cout << endl;
}
cout << endl;

}

void outprint(s_matrix & _mat) // 输出函数
{
for ( int i = 0 ;i < nSize; ++ i)
{

for ( int j = nSize;j < 2 * nSize; ++ j)
{
cout << _mat[i][j] << " " ;
}
cout << endl;
}


}

int main()
{
step = 0 ;
ma

trix.clear();
line.clear();
cout << " *********矩阵 求逆********* " << endl;
cout << " *********Bill Hsu********* " << endl;
cout << " /probill " << endl << endl;

cout << " 请输入矩阵维数(输入0退出)： " ;
cin >> nSize;
if (nSize <= 0 ) return 0 ;
for ( int i = 0 ;i < nSize; ++ i)
{
line.clear();
cout << " 输入第 " << i + 1 << " 行: " << endl;
for ( int j = 0 ;j < nSize; ++ j)
{
float tmp;
cin >> tmp;
line.push_back(tmp); // 压入一个数到某行
}

for ( int j = 0 ;j < nSize; ++ j)
{
if (i == j) line.push_back( 1.0f );
else line.push_back( 0.0f );
}


matrix.push_back(line); // 压入一行到矩阵
}

cout << endl;
work1(matrix);
work2(matrix);
makeunit(matrix);
cout << endl << " ######################## " << endl
<< " 求逆结果： " << endl;
outprint(matrix);
cout << " ######################## " << endl;

main();


return 0 ;
}