gauss.txt

全主元高斯消元法求解复系数线性方程组的解

// ****全主元高斯消去法解复系数方程组。线性方程组a.x=b求解x；形参b返回的即是x
//***a为n*n方阵；b为已知的n列阵；n为方阵行和列数
void agaus(complex **a,complex *b,int n)
{ int *js,k,i,j,is;
  complex d,t;
  js=new int [n];
  for (k=0;k<=n-2;k++)
  { d=0.0;
    for (i=k;i<=n-1;i++)
    for (j=k;j<=n-1;j++)
    { t=(a[i][j]);
      if (abs(t)>abs(d)) { d=t; js[k]=j; is=i;}
    }
    if (d+1.0==1.0)
    { delete []js; printf("fail--\n"); return;}
    else
    { if (js[k]!=k)
      for (i=0;i<=n-1;i++) { t=a[i][k]; a[i][k]=a[i][js[k]]; a[i][js[k]]=t;}
      if (is!=k)
      { for (j=k;j<=n-1;j++) { t=a[k][j]; a[k][j]=a[is][j]; a[is][j]=t;}
        t=b[k]; b[k]=b[is]; b[is]=t;
      }
    }
    d=a[k][k];
    for (j=k+1;j<=n-1;j++)  a[k][j]=a[k][j]/d;
    b[k]=b[k]/d;
    for (i=k+1;i<=n-1;i++)
    { for (j=k+1;j<=n-1;j++) a[i][j]= a[i][j]-a[i][k]*a[k][j];
      b[i]=b[i]-a[i][k]*b[k];
    }
  }
  d=a[n-1][n-1];
  if (abs(d)+1.0==1.0) { delete []js; printf("fail++\n"); getchar(); return; }
  b[n-1]=b[n-1]/d;
  for (i=n-2;i>=0;i--)
  { t=0.0;
    for (j=i+1;j<=n-1;j++) t=t+a[i][j]*b[j];
    b[i]=b[i]-t;
  }
  js[n-1]=n-1;
  for (k=n-1;k>=0;k--)
    if (js[k]!=k) { t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;}
  delete []js;
  return;
}