Matrix& Matrix::getInverse()
{
if(len==0) throw Excp(BOp);
double* pm;
int row,col;
if((pm=new double [len*len])==NULL) throw Excp(OMem);
this->toPrimitive(pm,row,col);
Matrix* m=new Matrix(len);
int k,i,j,i1;
double temp;
for(i=0;i<len;i++)
(*m)[i][i]=1;
for(k=0;k<row-1;k++){
i1=k;
for(i=k+1;i<row;i++){
if(fabs(*(pm+i1*col+k))<fabs(*(pm+i*col+k))){
i1=i;
}
}
if(*(pm+i1*col+k)==0.0){
delete [] pm;
throw Excp(BOp);
}
if(i1!=k){
for(j=0;j<col;j++){
temp=*(pm+k*col+j);
*(pm+k*col+j)=*(pm+i1*col+j);
*(pm+i1*col+j)=temp;
temp=(*m)[k][j];
(*m)[k][j]=(*m)[i1][j];
(*m)[i1][j]=temp;
}
}
temp=*(pm+k*col+k);
for(j=col-1;j>=0;j--){
*(pm+k*col+j)/=temp;
(*m)[k][j]/=temp;
}
for(i=k+1;i<row;i++){
temp=*(pm+i*col+k);
for(j=0;j<col;j++){
*(pm+i*col+j)-=temp*(*(pm+k*col+j));
(*m)[i][j]-=temp*(*m)[k][j];
}
}
}
//k=row-1;
for(j=0;j<col;j++) (*m)[k][j]/=*(pm+k*col+k);
*(pm+k*col+k)/=*(pm+k*col+k);
for(k=row-1;k>0;k--){
for(i=k-1;i>=0;i--){
temp=*(pm+i*col+k);
for(j=col-1;j>=0;j--){
*(pm+i*col+j)-=temp*(*(pm+k*col+j));
(*m)[i][j]-=temp*(*m)[k][j];
}
}
}
delete [] pm;
return *m;
}
呵呵，原来以前我也会C＋＋的*_*
评论0