gs3.cpp

本程序实现了矩阵计算中的列主元消去法

// GS3.cpp : Defines the entry point for the console application.


#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
int GS1(int ,int,double **, double **,double);
double **TwoArrayAlloc(int,int);
void TwoArrayFree(double **);
int main(int argc, char* argv[])
{
	//printf("Hello World!\n");
	int i,j,m,n;
	double ep,**a,**b,c;
	m=2;
	n=3;
	ep=1e-12;
	a=TwoArrayAlloc(n,n);
	b=TwoArrayAlloc(n,m);
	//b=(double*)calloc(n,sizeof(double));
	if(b==NULL)
		exit(1);

	a[0][0]=1;a[0][1]=-1;a[0][2]=1;
	a[1][0]=5;a[1][1]=-4;a[1][2]=3;
	a[2][0]=2;a[2][1]=1;a[2][2]=1;
	b[0][0]=-4;b[0][1]=-4;b[1][0]=-12;
		b[1][1]=2;b[2][0]=11;b[2][1]=-2;



	/*

	a[0][0]=10;a[0][1]=-7;a[0][2]=0;
	a[1][0]=5;a[1][1]=-1;a[1][2]=5;
	a[2][0]=-3;a[2][1]=2;a[2][2]=6;
	b[0]=7;b[1]=6;b[2]=4;

  */

/*
	a[0][0]=10;a[0][1]=-7;a[0][2]=0;
	a[1][0]=-3;a[1][1]=2.099;a[1][2]=6;
	a[2][0]=5;a[2][1]=-1;a[2][2]=5;
	b[0]=7;b[1]=3.901;b[2]=6;
	*/
	/*
	a[0][0]=1;a[0][1]=-1;a[0][2]=1;
	a[1][0]=5;a[1][1]=-4;a[1][2]=3;
	a[2][0]=2;a[2][1]=1;a[2][2]=1;
	b[0]=-4;b[1]=-12;b[2]=11;
	*/






	if(!GS1(n,m,a,b,ep))
	{
		printf("failed\n");
		exit(0);
	}
	for(j=1;j<=m;j++)
	for(i=1;i<=n;i++)//改变一下显示顺序，先列后行
		//for(j=1;j<=m;j++)
		printf("%f\n",b[i-1][j-1]);
	TwoArrayFree(a);
	TwoArrayFree(b);
//	free(b);
	
}

int GS1(int n,int m,double **a,double **b,double ep)//输出结果为什么没交换？？？？ 按行交换
{
	int i,j,k,l;
	double c, t;
	for(k=1;k<=n-1;k++)
	{
		c=0.0;
		for(i=k;i<=n;i++)
		{
			if(fabs(a[i-1][k-1])>fabs(c))
			{
				c=a[i-1][k-1];
				l=i;
			}
		}
		if (fabs(c)<=ep) return(0);

		
				
				if(l!=k)
				{
					for(j=k;j<=n;j++)
					{
						t=a[k-1][j-1];
						a[k-1][j-1]=a[l-1][j-1];
						a[l-1][j-1]=t;
					}
					for(i=1;i<=m;i++)
					{
					t=b[k-1][i-1];
					b[k-1][i-1]=b[l-1][i-1];
					b[l-1][i-1]=t;
					}
				}
			
				t=1/c;
				
				for(j=k+1;j<=n;j++)//????
					a[k-1][j-1]=t*a[k-1][j-1];
				for(i=1;i<=m;i++)
				b[k-1][i-1]=b[k-1][i-1]*t;
			
			//	if(k==n)
				//	break;
				for(i=k+1;i<=n;i++)
				{
					for(j=k+1;j<=n;j++)//????
						a[i-1][j-1]=a[i-1][j-1]-a[i-1][k-1]*a[k-1][j-1];
				
				}
					
						for(i=k+1;i<=n;i++)
							for(j=1;j<=m;j++)
					b[i-1][j-1]=b[i-1][j-1]-a[i-1][k-1]*b[k-1][j-1];

				
	}

	b[n-1][m-1]=b[n-1][m-1]/a[n-1][n-1];
		b[n-1][m-2]=b[n-1][m-2]/a[n-1][n-1];
	for(i=n-1;i>=1;i--)//????
		for(j=i+1;j<=n;j++)
			for(k=1;k<=m;k++)
			b[i-1][k-1]=b[i-1][k-1]-a[i-1][j-1]*b[j-1][k-1];

		return(1);
}

double **TwoArrayAlloc(int r,int c)
{
	double *x,**y;
	int n;
	x=(double *)calloc(r*c,sizeof(double));
	y=(double **)calloc(r,sizeof(double *));
	for(n=0;n<=r-1;n++)
		y[n]=&x[c*n];
	return(y);
}


void TwoArrayFree(double **x)
{
	free(x[0]);
	free(x);
}