//求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;
}