高斯列选主元算法.cpp

这是计算方法上的所有算法

#include <math.h>
#include "stdio.h"
#include "conio.h"
#define N 3
static double aa[N][N]={{1,2,-1},{1,-1,5},{4,1,-2}};
static double bb[N+1]={3,0,2};

void main()
{
	int i,j;
	double a[N+1][N+1],b[N+1],x[N+1],det;
	double gaussl(double a[][N+1],double b[],double x[]);

	for(i=1;i<=N;i++)
	{
		for(j=1;j<=N;j++) a[i][j]=aa[i-1][j-1];
		b[i]=bb[i-1];
	}
		
	det=gaussl(a,b,x);
	if(det!=0)
	{	printf("\n方程组的解为:");
		for(i=1;i<=N;i++)printf(" x[%d]=%f",i,x[i]);
		printf("\n\n系数矩阵的行列式值＝%f",det);
	}
	else printf("\n\n系数矩阵奇异阵，高斯方法失败 !:");
	getch();
		
}

double gaussl(double a[][N+1],double b[],double x[])
//a传入增广矩阵(有效元素为a[1,1]...a[n,n+1]),x负责传入迭代初值并返回解向量
//(有效值为x[1]...x[n]);返回值：系数矩阵行列式值detA
{double det=1.0,c0,F,m,temp;
 int i,j,k,r;
 void disp(double a[][N+1],double b[]);
 printf("消元前增广矩阵 :\n");
 disp(a,b);
 for(k=1;k<N;k++)
 { temp=a[k][k];
   r=k;
   for(i=k;i<=N;i++)
       if(temp<a[i][k])
	{temp=a[i][k];
	 r=i;
	}//按列选主元，即确定ik
 if(a[r][k]==0)
	  return 0;//如果aik,k=0,则A为奇异矩阵，停止计算
  if(r!=k)
    {for(j=k;j<=N;j++)
      {c0=a[k][j];
       a[k][j]=a[r][j];
       a[r][j]=c0;
       }
    c0=b[k];
    b[k]=b[r];
    b[r]=c0;
    det=-det;//如果ik!=k,则交换[A,b]第ik行与第k行元素
    printf("交换第 %d, %d行:\n",k,r);
    disp(a,b);
    }
  for(i=k+1;i<=N;i++)
  {m=a[i][k]/a[k][k];
    for(j=1;j<=k;j++)
		a[i][j]=0.0;
    for(j=k+1;j<=N;j++)
       a[i][j]-=m*a[k][j];
	b[i]-=m*b[k];
	//消元计算
	//aik<-mik=aik/akk (i=k+1,...,n)
	//aij<-aij-mik*akj (i,j=k+1,...,n)
	//bi<-bi-mik*bk (i=k+1,...,n)
  }
    printf("第%d次消元:\n",k);
    disp(a,b);
    det*=a[k][k]; 
 }
  x[N]=b[N]/a[N][N];
  for(i=N-1;i>0;i--)
    {F=0.0;
     for(j=i+1;j<=N;j++)
	    F+=a[i][j]*x[j];
     x[i]=(b[i]-F)/a[i][i];
     }
  //回代计算
  //bn<-bn/ann (ann!=0)
  //bi<-(bi-@@aijbj)/aii (i=n-1,...,2,1)

  det*=a[N][N];
 
  return det;
	
}

void disp(double a[][N+1],double b[])
//显示选主元及消元运算中间增广矩阵结果
{int i,j;
 for(i=1;i<=N;i++)
   {for(j=1;j<=N;j++)
      printf("%f    ",a[i][j]);
    printf("%f\n",b[i]);
    }
	
}