inv.h

实现矩阵加、减、乘和求逆的C++程序

/*****************C=inv(A)******************/
#include <math.h>

inline void swap(double &a,double &b){double c=a;a=b;b=c;};
void inv(double *A,double *C,int N)
{
	int i,j,k;
	double d;
	int JS[100],IS[100];
	for (k=0;k<N;k++)
	{
		d=0;
		for (i=k;i<N;i++)
			for (j=k;j<N;j++)
			{
				if (fabs(A[i*N+j])>d)
				{
					d=fabs(A[i*N+j]);
					IS[k]=i;
					JS[k]=j;
				}
			}
		if (d+1.0==1.0) return;
		if (IS[k]!=k)
			for (j=0;j<N;j++)
				swap(A[k*N+j],A[IS[k]*N+j]);
			if (JS[k]!=k)
				for (i=0;i<N;i++)
					swap(A[i*N+k],A[i*N+JS[k]]);
				A[k*N+k]=1/A[k*N+k];
				for (j=0;j<N;j++)
					if (j!=k) A[k*N+j]=A[k*N+j]*A[k*N+k];
				for (i=0;i<N;i++)
					if (i!=k)
						for (j=0;j<N;j++)
							if (j!=k) A[i*N+j]=A[i*N+j]-A[i*N+k]*A[k*N+j];
				for (i=0;i<N;i++)
					if (i!=k) A[i*N+k]=-A[i*N+k]*A[k*N+k];
	}

	for (k=N-1;k>=0;k--)
	{
		for (j=0;j<N;j++)
			if (JS[k]!=k) swap(A[k*N+j],A[JS[k]*N+j]);
			for (i=0;i<N;i++)
				if (IS[k]!=k) swap(A[i*N+k],A[i*N+IS[k]]);
	}

	for (i=0;i<N;i++)
		for (j=0;j<N;j++)
			C[i*N+j]=A[i*N+j];	
}