gauss列主元消去法1.c

来自「高斯列主元消去法」· C语言 代码 · 共 83 行

#include<stdio.h>
#include<math.h>
#define MAX_N 20
main()
{
 int i,j,k;
 int n;
 int i_k;                              //列主元所在行的指标
 double m;                             //列主元的值
 double a,b;                           //换行时做终结制值
 double D=1;                           //行列式的值
 double A[MAX_N][MAX_N];
 double B[MAX_N];
// double M[MAX_N][MAX_N];
 printf("请输入系数矩阵的阶数n:\n");
 scanf("%d",&n);
 printf("请输入系数矩阵A:\n");
 for(i=0;i<n;i++)
   for(j=0;j<n;j++)
   scanf("%lf",&A[i][j]);
 printf("请输入右端项b：\n");
 for(i=0;i<n;i++)
 scanf("%lf",&B[i]);
 for(k=0;k<n-1;k++)         //最后一行就不用选了
 {
  m=0;                      //开始按列选主元
  i_k=k;
  for(i=k;i<n;i++)
  {
   if(fabs(A[i][k])>m)
   {
    m=fabs(A[i][k]);
    i_k=i;                  //记录主元所在行的指标
   } //if
  }  //for i
  if(m==0)
  {
   printf("系数矩阵奇异！\ndetA=0\n");
   break;
  } //if
  else if(i_k!=k)
  {
   for(j=k;j<n;j++)       //换行     书上的算法这里错了！！！
   {
    a=A[i_k][j];
    A[i_k][j]=A[k][j];
    A[k][j]=a;
   } //for j
   b=B[i_k];
   B[i_k]=B[k];
   B[k]=b;
   D=(-1)*D;                   //换行，行列式变号
  } //else if
  D=D*A[k][k];
  for(i=k+1;i<n;i++)           //计算乘子
  {
   A[i][k]=A[i][k]/A[k][k];
  } //for i
  for(i=k+1;i<n;i++)           //消元
  {
   for(j=k+1;j<n;j++)
   {
    A[i][j]=A[i][j]-A[i][k]*A[k][j];
   } //for j
   B[i]=B[i]-A[i][k]*B[k];
  }  //for i
 } //for k
  for(i=n-1;i>=0;i--)          //回代；临时用b来存储这里的和
  {
   b=0;
   for(j=i+1;j<n;j++)
   {b+=A[i][j]*B[j];}
   B[i]=(B[i]-b)/A[i][i];
  } //for i
 D=D*A[n-1][n-1];             //在k进行循环的时候只循环到了n-2   （也可以在上面让k循环到n-1）
 printf("系数矩阵行列式的值为%lf。\n",D);
 printf("解为x_i=\n");
 for(i=0;i<n;i++)
 printf("%lf\n",B[i]);
 return 0;
}