求解线性方程组之全选主元高斯消元法

//全选主元Gauss消去法
#include <iostream>
#include <fstream>
#include <math.h>
#include <process.h>

using namespace std;

class gauss
{
private:
 int i, j, k, n, *pivrow, **pivcol;
 double eps, pivot, sum, aik, al, *x, **a;

public:
 void gauss_input();
 void gauss_elimination();
 void gauss_output();
 ~gauss()
 {
  delete[] x;
  delete[] pivrow;
  for (i = 0; i < n; i++)
  {
   delete[] a[i];
  }
  delete[] a;
  for (i = 0; i < n; i++)
  {
   delete[] pivcol[i];
  }
  delete[] pivcol;
 }
};

void main()
{
 gauss solution;
 solution.gauss_input();
 solution.gauss_elimination();
 solution.gauss_output();
}

void gauss::gauss_input()
{
 ifstream fin("gauss.txt");
 fin >> n;
 x = new double[n];
 pivrow = new int[n];
 pivcol = new int*[n];
 for (i = 0; i < n; i++)
 {
  pivcol[i] = new int[2];
 }
 a = new double*[n];
 for (i = 0; i < n; i++)
 {
  a[i] = new double[n+1];
 }
 for (i = 0; i < n; i++)
  for (j = 0; j < n; j++)
  {
   fin >> a[i][j];
  }
 for(i = 0; i < n; i++)
 {
  fin >> a[i][n];
 }
 fin >> eps;
 fin.close();
}

void gauss::gauss_elimination()
{
 for (k = 0; k < (n-1); k++)
 {
  pivot = fabs(a[k][k]);
  pivrow[k] = k;
  pivcol[k][0] = k;
  pivcol[k][1] = k;
  for (i = k; i < n; i++)
   for (j = k; j < n; j++)
   {
    if (pivot < fabs(a[i][j]))
    {
     pivot = fabs(a[i][j]);
     pivrow[k] = i;
     pivcol[k][1] = j;
    }
   }
  if (pivot < eps)
  {
   cout << "\n矩阵可能奇异，失败..." << endl;
   exit(0);
  }
  if (pivrow[k] != k)
  {
   for (j = k; j < (n+1); j++)
   {
    al = a[pivrow[k]][j];
    a[pivrow[k]][j] = a[k][j];
    a[k][j] = al;
   }
  }
  if (pivcol[k][1] != k)
  {
   for (i = 0; i < n; i++)
   {
    al = a[i][pivcol[k][1]];
    a[i][pivcol[k][1]] = a[i][k];
    a[i][k] = al;
   }
  }
  if (k != (n-1))
  {
   for (i = (k+1); i < n; i++)
   {
    aik = a[i][k];
    for (j = k; j < (n+1); j++)
    {
     a[i][j] -= aik * a[k][j] / a[k][k];
    }
   }
  }
 }
 x[n-1] = a[n-1][n] / a[n-1][n-1];
 for (i = (n-2); i >= 0; i--)
 {
  sum = 0.0;
  for (j = (i+1); j < n; j++)
  {
   sum += a[i][j] * x[j];
  }
  x[i] = (a[i][n] - sum) / a[i][i];
 }
 for (k = (n-2); k >= 0; k--)
 {
  al = x[pivcol[k][1]];
  x[pivcol[k][1]] = x[pivcol[k][0]];
  x[pivcol[k][0]] = al;
 }
}

void gauss::gauss_output()
{
 ofstream fout("gauss.out");
 fout << "\n结果是：" << endl;
 for (i = 0; i < n; i++)
 {
  fout << "\nx[" << i << "] = " << x[i] << endl;
 }
 fout.close();
}