guangyinijuzhen.cpp

基于vc6.0的一个计算非方阵矩阵的广义逆矩阵的算法

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

//矩阵输出
void pr(double *a,int n,int m)
{
	int i,j;
	for(i=0;i<n;i++)
	{
		for(j=0;j<m;j++)
		{
			printf("%.6lf ",a[i*m+j]);
		}
		printf("\n");
	}
	printf("\n");
}
//矩阵的转置
void tr(double *a,double *t,int n,int m)
{
	int i,j;
	for(i=0;i<n;i++)
	{
		for(j=0;j<m;j++)
		{
			t[j*n+i]=a[i*m+j];
		}
	}
}

//矩阵乘法
void mult(double *a,double *b,double *c,int n,int m,int x)
{
	int i,j,k;
	for(i=0;i<n;i++)
	{
		for(j=0;j<m;j++)
		{
			c[i*m+j]=0;
			for(k=0;k<x;k++)
			{
				c[i*m+j]+=a[i*x+k]*b[k*m+j];
			}
		}
	}
}

//逆矩阵求解
double *ni(double *a,int n)
{
	int i,j;
	double b,temp,d,*s,*c;
	int k;
	s=(double*)malloc(n*2*n*sizeof(double));
	c=(double*)malloc(n*n*sizeof(double));

	//建立增广矩阵
	for(i=0;i<n;i++)	
	{
		for(k=0;k<2*n;k++)
			if(k<n)
				s[i*2*n+k]=a[i*n+k];
			else
			{
				if(i==k-n)
					s[i*2*n+k]=1;
				else
					s[i*2*n+k]=0;
			}
	}
	//完成建立增广矩阵
	
	//i从0变化到n
	for(i=0;i<n;i++)
	{
		b=s[i*2*n+i];			
		j=i;
		//2.1查找第i列绝对值最大的元素
		for(k=i;k<n;k++)		
		{	
			if(fabs(s[k*2*n+i])>fabs(b))
			{
				b=s[k*2*n+i];
				j=k;
			}
		
		
		}
		//2.1查找完成

		//2.3交换第i行和第j行
		if(b!=0)			
		{
			if(j!=i)
			{
				for(k=0;k<2*n;k++)
				{
					temp=s[j*2*n+k];
					s[j*2*n+k]=s[i*2*n+k];
					s[i*2*n+k]=temp;
				}
			}
		}
		//2.3交换完成

		//2.4 b!=1,i行元素除以b
		if(b!=1)
		{
			
			for(k=0;k<2*n;k++)
			{
				s[i*2*n+k]=s[i*2*n+k]/b;	
			}
		}
		//2.4

		//2.5 第i列k行为非0元素s[k][i],则第k行第j列元素为s[k][j]=s[k][j]-s[k][i]*s[k][i]
		for(k=0;k<2*n;k++)
		{
			if(k!=i&&s[k*2*n+i]!=0)
			{
				d=s[k*2*n+i];
				for(j=0;j<2*n;j++)
				{
					s[k*2*n+j]=s[k*2*n+j]-s[i*2*n+j]*d;
				}
			}
		}
		//2.5
	}
	for(k=0;k<n;k++)
	{
		for(j=0;j<n;j++)
		{
			c[k*n+j]=s[k*2*n+n+j];
		}
	}
	return c;
}
//逆矩阵求解完成

//
int max(int n,int m)
{
	int i;
	if(n>m)
		i=n;
	else
		i=m;
	return i;

}

//广义逆矩阵
double *q(double *A,int m,int n)
{

	
	double *An,*R,*R1,*R2,*R3,*C,*C1,*C3,*B,b,d,temp,*X,*Y,*Z,*P,*Q,*Pt,*Qt,*Bn,*E;
	int  i,j,bn,k,x,y;
	An=(double *)malloc(n*m*(sizeof(double)));	
	R=(double *)malloc(m*m*(sizeof(double)));
	C=(double *)malloc(n*n*(sizeof(double)));
	B=(double *)malloc(n*m*sizeof(double));
	Bn=(double *)malloc(n*m*sizeof(double));

	R1=(double *)malloc(m*m*sizeof(double));
	R2=(double *)malloc(m*m*sizeof(double));
	R3=(double *)malloc(m*m*sizeof(double));
	C1=(double *)malloc(n*n*sizeof(double));
	C3=(double *)malloc(n*n*sizeof(double));

	X=(double *)malloc(max(n,m)*max(n,m)*sizeof(double));
	Y=(double *)malloc(max(n,m)*max(n,m)*sizeof(double));
	Z=(double *)malloc(max(n,m)*max(n,m)*sizeof(double));
	//定义单位阵R
	for(x=0;x<m;x++)
	{
		for(k=0;k<m;k++)
		{
			if(x==k)
			{
				R[x*m+k]=1;
			}
			else{
				R[x*m+k]=0;
			}
		}
	}

//	pr(R,m,m);
	//定义单位阵C
	for(x=0;x<n;x++)
	{
		for(k=0;k<n;k++)
		{
			if(k==x)
			{
				C[x*n+k]=1;
			}
			else
			{
				C[x*n+k]=0;
			}
		}
	}
//	pr(C,n,n);
	//定义矩阵b=a
	for(i=0;i<m;i++)
	{
		for(j=0;j<n;j++)
			B[i*n+j]=A[i*n+j];
	}

//	pr(B,m,n);
	i=0;
	bn=n;


	//2
	while(i<bn)
	{
		//初始化
		b=B[i*n+i];			
		j=i;

//		printf("执行操作前的矩阵B\n");
//		pr(B,m,n);


		//2.1查找第i列绝对值最大的元素
		for(k=i;k<m;k++)		
		{	
			if(fabs(B[k*n+i])>fabs(b))
			{
				b=B[k*n+i];
				j=k;
			}
			
		}
		//2.1查找完成
		
//		printf("i=%d j=%d bn=%d b=%.6lf\n",i,j,bn,b);


		//2.2开始
		if(b==0&&i==bn-1)
		{
			break;
		}

		//2.2结束

		//2.3开始
		else if(b==0&&i!=bn-1)
		{
			for(k=0;k<m;k++)
			{
				temp=B[k*n+i];
				B[k*n+i]=B[k*n+bn-1];
				B[k*n+bn-1]=temp;
			}
			for(x=0;x<n;x++)
			{
				for(k=0;k<n;k++)
				{
					if(x==k)
						C1[x*n+k]=1;
					else
						C1[x*n+k]=0;
				}
			}
			C1[i*n+(bn-1)]=1;
			C1[(bn-1)*n+i]=1;
			C1[i*n+i]=0;
			C1[(bn-1)*n+(bn-1)]=0;

			mult(C,C1,X,n,n,n);

			//C[]=X[]
			for(x=0;x<n;x++)
			{
				for(k=0;k<n;k++)
				{
					C[x*n+k]=X[x*n+k];
				}
			}
/*
			printf("执行2.3后的C\n");
			pr(C,n,n);
			pr(B,m,n);
*/			

			bn--;
		}

		//2.3结束
		

		//2.4开始
		else if(j!=i)
		{
			for(k=0;k<n;k++)
			{
				temp=B[j*n+k];
				B[j*n+k]=B[i*n+k]; 
				B[i*n+k]=temp;
			}

			//R1(m,i,j)
			for(x=0;x<m;x++)
			{
				for(y=0;y<m;y++)
				{
					if(x==y)
						R1[x*m+y]=1;
					else
						R1[x*m+y]=0;
				}
			}

			R1[i*m+j]=1;
			R1[j*m+i]=1;
			R1[i*m+i]=0;
			R1[j*m+j]=0;
			//R1配置完成		

/*			printf("2.4R1\n");
			pr(R1,m,m);
			printf("2.4前R\n");
			pr(R,m,m);
*/
			mult(R1,R,X,m,m,m);

			//R[]=X[]
			for(x=0;x<m;x++)
			{
				for(k=0;k<m;k++)
				{
					R[x*m+k]=X[x*m+k];
				}
			}

//			printf("2.4后R\n");
//			pr(R,m,m);
			
		}
		//2.4结束

		//2.5开始
		if(b!=1&&b!=0)
		{
			for(k=0;k<n;k++){
				B[i*n+k]=B[i*n+k]/b;
			}

//			printf("2.5后的B\n");
//			pr(B,m,n);
			//R2(m,i,1/b)
			for(x=0;x<m;x++)
			{
				for(y=0;y<m;y++)
				{
					if(x==y)
						R2[x*m+y]=1;
					else
						R2[x*m+y]=0;
				}
			}
			R2[i*m+i]=1/b;

			//R2配置完成

//			printf("b=%.6lf,i=%d\n",b,i);
			//R2
//			printf("2.5R2\n");
//			pr(R2,m,m);	
//			printf("2.5前R\n");
//			pr(R,m,m);
			mult(R2,R,X,m,m,m);
			
			//R[]=X[]
			for(x=0;x<m;x++)
			{
				for(y=0;y<m;y++)
				{
					R[x*m+y]=X[x*m+y];
				}
			}

			//R
//			printf("2.5后的R\n");
//			pr(R,m,m);

			
		}
		//2.5结束

//		printf("2.6前B");
//		pr(B,m,n);
		//2.6开始
		for(k=0;k<m;k++)
		{
			if(k!=i&&B[k*n+i]!=0)
			{
				d=B[k*n+i];
				for(j=0;j<n;j++)
				{
					B[k*n+j]=B[k*n+j]-B[i*n+j]*d;
				}

				//配置R3(m,k,i,-d)
				for(x=0;x<m;x++)
				{
					for(y=0;y<m;y++)
					{
						if(x==y)
							R3[x*m+y]=1;
						else
							R3[x*m+y]=0;
					}
				}
				R3[k*m+i]=-d;
				//配置R3完成
				
//				printf("2.6R3\n");
//				printf("k=%d,i=%d,d=%.6lf\n",k,i,-d);
//				pr(R3,m,m);
//				printf("2.6前的R\n");
//				pr(R,m,m);

				mult(R3,R,X,m,m,m);	
				
				//R[]=X[]
				for(x=0;x<m;x++)
				{
					for(y=0;y<m;y++)
					{
						R[x*m+y]=X[x*m+y];
					}
				}

				//R
//				printf("2.6后的R\n");
//				pr(R,m,m);
			
			}		
		}
		//2.6结束

		//2.7开始
		for(k=0;k<n;k++)
		{
			if(k!=i&&B[i*n+k]!=0)
			{
				d=B[i*n+k];

				for(j=0;j<m;j++)
				{
					B[j*n+k]=B[j*n+k]-B[j*n+i]*d;
				}
				//配置C3
				for(x=0;x<n;x++)
				{
					for(y=0;y<n;y++)
					{
						if(x==y)
							C3[x*n+y]=1;
						else
							C3[x*n+y]=0;
					}
				}
				C3[i*n+k]=-d;
				//配置C3完成

//				printf("2.7前C\n");
//				pr(C,n,n);

//				printf("2.7C3\n");
//				pr(C3,n,n);

				mult(C,C3,X,n,n,n);	
				
				//C[]=X[]
				for(x=0;x<n;x++)
				{
					for(y=0;y<n;y++)
					{
						C[x*n+y]=X[x*n+y];
					}
				}

//				printf("2.7后的C d=%.6lf\n",-d);
//				pr(C,n,n);
//				printf("\n");

			}
		}
/*		
		printf("执行操作后的矩阵B\n");
		pr(B,m,n);
		printf("\n");
*/
		
		//2.8开始
		i++;
		//2.8结束
	}
	//2结束

	//3开始
	int r=0;
	for(r=0;B[r*n+r]!=0&&r<m&&r<n;r++){}
	printf("%d\n",r);
	//3结束	
	pr(C,n,n);
	pr(R,m,m);


	pr(A,m,n);
	tr(B,Bn,m,n);
	pr(Bn,n,m);
	mult(C,Bn,X,n,m,n);
	mult(X,R,An,n,m,m);
	pr(An,n,m);
	//Y为A+
/*验证
	mult(A,An,Z,m,m,n);
	mult(Z,A,X,m,n,m);
	pr(X,m,n);
*/
	return An;
}
void main()
{
	double a[3][2]={{1,-1},{2,3},{3,2}};
	double *A=&a[0][0],*B,*C;
	C=(double*)malloc(3*2*sizeof(double));
	B=q(A,3,2);
	mult(A,B,C,3,3,2);
	pr(B,2,3);

}