矩阵求你.docVIP

矩阵求你 前段时间写到一个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;kn;k++) { double pivot = Matrix[k][k]; for (i =k+1;in;i++) { Matrix[i][k] /= pivot; for (j = k+1;jn;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) cerrPai Malloc ERROR; for (i=0;in;i++) Pai[i] = i; for (k = 0;kn;k++) { double p = 0; for (i = k;in;i++) if (fabs(Matrix[i][k])p) { p = fabs(Matrix[i][k]); kp = i; } if (0==p) { cerrsingular Matrix; PAUSE } swap(Pai[k],Pai[kp]); for (i = 0;in;i++) swap(Matrix[k][i],Matrix[kp][i]); for (i =k+1;in;i++) { Matrix[i][k] /= Matrix[k][k]; for (j = k+1;jn;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)) cerrLUP_SOLVE xMalloc ERROR; for (i =0;irow;i++) { x[i] = 0; y[i] = 0; } for (i = 0;in;i++) { double S1 = 0; for (j = 0;ji;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;jn-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