i.cpp

用最大主元法求矩阵的逆

/////////////矩阵相乘
void CalMatProduct(double* pDbSrc1, double *pDbSrc2, double *pDbDest, int ny, int nx, int nz)
{
	
	for(int i=0;i<ny;i++)
		for(int j=0;j<nx;j++)
		{
			pDbDest[i*nx + j] = 0;
			for(int m=0;m<nz;m++)
				pDbDest[i*nx + j] += pDbSrc1[i*nz + m]*pDbSrc2[m*nx + j];
		}
}



///////////////求矩阵的转置
void   CalTMatrix(double* pDbSrc,int y,int x,double *pDbDest)
{
	
	for(int i=0;i<x;i++)
		for(int j=0;j<y;j++)
		{
			pDbDest[i*y +j] = pDbSrc[j*x+i];
			
		}
}


/////求矩阵的逆
bool CalInvMatrix(double *pDbSrc, int nLen)
{
	int *is,*js,i,j,k;
	double d,p;
	is = new int[nLen];
	js = new int[nLen];
	for(k=0;k<nLen;k++)
	{
		d=0.0;
		for(i=k;i<nLen;i++)
			for(j=k;j<nLen;j++)
			{
				p=fabs(pDbSrc[i*nLen + j]);
				if(p>d)
				{
					d     = p; 
					is[k] = i;
					js[k] = j;
				}
			}
			if(d+1.0==1.0)
			{
				delete is;
				delete js;
				return  false;
			}
			if(is[k] != k)
				for(j=0;j<nLen;j++)
				{
					p = pDbSrc[k*nLen + j];
					pDbSrc[k*nLen + j] = pDbSrc[(is[k]*nLen) + j];
					pDbSrc[(is[k])*nLen + j] = p;
				}
				if(js[k] != k)
					for(i=0; i<nLen; i++)
					{
						p = pDbSrc[i*nLen + k];
						pDbSrc[i*nLen + k] = pDbSrc[i*nLen + (js[k])];
						pDbSrc[i*nLen + (js[k])] = p;
					}
					
					pDbSrc[k*nLen + k]=1.0/pDbSrc[k*nLen + k];
					for(j=0; j<nLen; j++)
						if(j != k)
						{
							pDbSrc[k*nLen + j]*=pDbSrc[k*nLen + k];
						}
						for(i=0; i<nLen; i++)
							if(i != k)
								for(j=0; j<nLen; j++)
									if(j!=k)
									{
										pDbSrc[i*nLen + j] -= pDbSrc[i*nLen + k]*pDbSrc[k*nLen + j];
									}
									for(i=0; i<nLen; i++)
										if(i != k)
										{
											pDbSrc[i*nLen + k] *= -pDbSrc[k*nLen + k];
										}
	}
	for(k=nLen-1; k>=0; k--)
	{
		if(js[k] != k)
			for(j=0; j<nLen; j++)
			{
				p = pDbSrc[k*nLen + j];
				pDbSrc[k*nLen + j] = pDbSrc[(js[k])*nLen + j];
				pDbSrc[(js[k])*nLen + j] = p;
			}
			if(is[k] != k)
				for(i=0; i<nLen; i++)
				{
					p = pDbSrc[i*nLen + k];
					pDbSrc[i*nLen + k] = pDbSrc[i*nLen +(is[k])];
					pDbSrc[i*nLen + (is[k])] = p;
				}
	}
	delete is;	
	return true;	
}