juzhenqiuni.txt

求满秩矩阵pMatrix的逆矩阵

//求pMatrix的逆矩阵，并存结果于矩阵_pMatrix中
void ContraryMatrix(double *const pMatrix, double *const _pMatrix, const int &dim)
{ 
     double *tMatrix = new double[2*dim*dim];
     for (int i=0; i<dim; i++)
	 {
         for (int j=0; j<dim; j++)
             tMatrix[i*dim*2+j] = pMatrix[i*dim+j];        
     }
     for (i=0; i<dim; i++)
	 {
         for (int j=dim; j<dim*2; j++)
             tMatrix[i*dim*2+j] = 0.0;
         tMatrix[i*dim*2+dim+i] = 1.0;        
     }
     //Initialization over!
     for (i=0; i<dim; i++)//Process Cols
     {
         double base = tMatrix[i*dim*2+i];
         if (fabs(base) < 1E-300)
		 {
             cout << endl << "求逆矩阵过程中被零除，无法求解!" << endl;
             exit(0);
         }
         for (int j=0; j<dim; j++)//row
         {
             if (j == i) continue;
             double times = tMatrix[j*dim*2+i]/base;
             for (int k=0; k<dim*2; k++)//col
             {        
                 tMatrix[j*dim*2+k] = tMatrix[j*dim*2+k] - times*tMatrix[i*dim*2+k];
             }
         }        
         for (int k=0; k<dim*2; k++)
		 {
             tMatrix[i*dim*2+k] /= base;
         }
     }
     for (i=0; i<dim; i++)
     {
         for (int j=0; j<dim; j++)
             _pMatrix[i*dim+j] = tMatrix[i*dim*2+j+dim];        
     }    
     delete[] tMatrix;
}