逆矩阵.txt

高斯赛德尔逆矩阵算法类

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＋＋的*_*