maxconv.c

一个关于处理矩阵求逆的代码

#include "dos.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "conio.h"


void Matrix_rev(double * pa, double * pb, int n);
void Matrix_mul( double * pa, double * pb, double * pc, int n, int p, int m );
void Matrix_turn(double * pa, double * pb, int n, int m);
void Matrix_plus(pa,pb,pc,n,m,p);

main()
{
	double pa[2][2] = {
		{1,2}, {3,4}
	};
	//double pb[2][2] ;
	double pc[2][2] ;

	//double *pap = pa ;
	//double *pbp = pb ;
	
	//Matrix_rev(pa, pb, 2);
	//printf("%f",pb[1][1]);

	//Matrix_mul( pa, pa, pc, 2, 2, 2 );
	//printf("%f",pc[0][0]);

	//Matrix_turn( pa, pb, 2, 2); 
	//printf("%f",pb[0][1]);

	Matrix_plus(pa,pa,pc,2,2,0);
	printf("%f",pc[0][0]);

}
  

void Matrix_rev(double * pa, double * pb, int n)
 // double *pa,*pb;
  //int n;
{
	int i, j, k, max ;
	double t, a[400] ;

	for( i=0; i<n; i++ )

		for(j=0; j<n; j++)
			a[i*n+j] = * (pa+i*n+j) ;

	for( i=0; i<n; i++ )
	{
		for( j=0; j<n; j++ )
		{
			if( i == j )   
				* (pb+i*n+j) = 1;
			else * (pb+i*n+j) = 0;
		}
	}

	for( i=0; i<n; i++ )
	{
		for( max=i, k=i; k<n; k++ )
			if( a[i*n+i] < a[k*n+i] )
				max = k ;

		if(k!=i)
			for(j=0;j<n;j++)
			{
				t = a[i*n+j];
				a[i*n+j] = a[max*n+j];
				a[max*n+j] = t;
				t = *(pb+i*n+j);
				*(pb+i*n+j) = *(pb+max*n+j);
				*(pb+max*n+j) = t;
			}

		if(a[i*n+i]==0)/*========================本列所有元素都等于零*/
		{
			printf("\nCAN NOT REVERSE THE MAXI.\n");
			exit(0);
		}


		for( j=0; j<n; j++ )
			*(pb+i*n+j) = *(pb+i*n+j) / a[i*n+i];

		for( j=i+1; j<n; j++ )
			a[i*n+j] = a[i*n+j] / a[i*n+i];

		for( k=i+1; k<n; k++ )
		{
			for( j=0; j<n; j++ )
				*(pb+k*n+j) = *(pb+k*n+j) - *(pb+i*n+j) * a[k*n+i];

			for( j=i+1; j<n; j++ )
				a[k*n+j] = a[k*n+j] - a[i*n+j] * a[k*n+i];
		}
	}/*=========================================至此,a[]已换成对角线元素为1的上三角阵*/

	for( i=n-1; i>=0; i-- )
		for( k=0; k<i; k++ )
			for( j=0; j<n; j++ )
				*(pb+k*n+j) = *(pb+k*n+j) - *(pb+i*n+j) * a[k*n+i];

}


void Matrix_mul( double * pa, double * pb, double * pc, int n, int p, int m )
  //double *pa,*pb,*pc;
  //int n,p,m;
{
	int i,j,k;
	double sum;
	for(i=0;i<n;i++)
	{
		for(j=0;j<m;j++)
		{
			for ( k=0,sum=0; k<p; k++ )
			sum = sum + *(pa+i*p+k) * ( *(pb+k*m+j) );
			*(pc+i*m+j)=sum;
		}
	 }
  }



void Matrix_turn(double * pa, double * pb, int n, int m)
  //double *pa,*pb;/*=======================这种声明方式不会出警告*/
  //int m,n;
 
{
	int i,j;
	for(i=0;i<n;i++)
	{
		for(j=0;j<m;j++)
		*(pb+j*n+i)=*(pa+i*m+j);
	}
}


/*======p=1为加法,否则为减法*/
void Matrix_plus(pa,pb,pc,n,m,p)
	double *pa,*pb,*pc;
	int n,m,p;
{
	int i,j;
	double ta,tb;
	for(i=0;i<n;i++)
	{
		for(j=0;j<m;j++)
		{
			ta=*(pa+i*m+j);
			tb=*(pb+i*m+j);
			if(p==1) 
				*(pc+i*m+j)=ta+tb;
			else  
				*(pc+i*m+j)=ta-tb;
	   
		}
	}
}


#if 0
void Matrix_cmul(pa,pb,pc,n1,m1,n2,m2)
	double *pa,*pb,*pc;
	int n1,n2,m1,m2;
{
	int i,j;
	double t1,t2;
	for(i=0;i<n1*n2;i++)
		for(j=0;j<m1*m2;j++)
		{
			t1 = *(pa+(int)fmod(i,n1)*m1+(int)fmod(j,m1));
			t2 = *(pb+(int)i/n1*m2+(int)(j/m1));
			*(pc+i*m1*m2+j)=t1*t2;
		}
}
#endif