Solving Gaussian Elimination with Zero Pivot

[philosafari]

I found this code on my computer a while ago, and I want to adjust the algorithm so that I can have a zero pivot, whereas now the program will try to divide by zero and gives the #undef message that we all know and love. Is there another type of Gaussian Elimination that I can look into? I only include the source code so that you can see what algorithm I'm currently using.

#include <stdlib.h>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <new>
using namespace std;

double** matrixAlloc(int, int);
double* vectorAlloc(int);
void matrixFree(double**, int);
void vectorFree(double*);
void matrixPrint(double**, int);
void vectorPrint(double*, int);

int main()
{
  double** a;
  double* b;
  double x, sum, r;
  int n;

  cout <<  "******* gaussian elimination *******" << endl << endl;
  cout << "Input number of equations n=" << endl;
  cin >> n;

  a = matrixAlloc(n, n);
  b = vectorAlloc(n);

  cout << "Input matrix coefficients a(i,j)=" << endl;
  for ( int i = 0; i < n; i++)
    for ( int j = 0; j < n; j++)
      cin >> a[i][j];

  cout << "Input right-hand side vector b(i)" << endl;
  for (int i = 0; i < n; i++)
    cin >> b[i];

  cout << endl << endl;
  cout << "Coefficient matrix A:" << endl;
  matrixPrint(a, n);
  cout << "Right-hand-side vector b:" << endl;
  vectorPrint(b, n);

  // convert to upper triangular form
  for (int k=0;k<n-1;k++)
  {
    if ( fabs(a[k][k])>=1.e-6)
    {
      for (int i=k+1;i<n;i++)
      {
        x = a[i][k]/a[k][k];
        for (int j=k+1;j<n;j++) a[i][j] = a[i][j] -a[k][j]*x;
          b[i] = b[i] - b[k]*x;
      }
    }
  }

  // back substituti n
  b[n-1]=b[n-1]/a[n-1][n-1];
  for ( int i = n-2; i >= 0; i--)
  {
    sum = b[i];
      for (int j = i+1; j < n; j++)
        sum = sum - a[i][j]*b[j];
    b[i] = sum/a[i][i];
  }

  cout << "Solution:" << endl;
  vectorPrint(b, n);
  
  matrixFree(a, n);
  vectorFree(b);

  cout << "Input \'q\' to quit ...\n";
  cin >> r;
  return 0;
}

double** matrixAlloc(int n, int m)
{
  double** p; 
  
  try {
    p = new double* [n];
    for ( int i = 0; i < n; i++ )
        p[i] = new double[m];
  }
  catch (bad_alloc e) {
    cout << "Exception occurred: "
         << e.what() << endl;
  }

  return p;
} 

double* vectorAlloc(int n)
{
  double* p;

  try {
    p = new double[n];
  }
  catch (bad_alloc e) {
    cout << "Exception occurred: "
         << e.what() << endl;
  }

  return p;
}

void matrixFree(double** p, int n)
{
  for ( int i = 0; i < n; i++)
    delete[] p[i];
  delete[] p;
  p = 0;
}

void vectorFree(double* p)
{
  delete[] p;
  p = 0;
}

void matrixPrint(double** p, int n)
{
  for ( int i = 0; i < n; i++)
  {
    for ( int j = 0; j < n; j++)
      cout << setiosflags(ios::showpoint | ios::fixed | ios::right) 
           << setprecision(4) 
           << setw(12) << p[i][j];
    cout << endl;
  }
}

void vectorPrint(double* p, int n)
{
  for ( int i = 0; i < n; i++)
    cout << setiosflags(ios::showpoint | ios::fixed | ios::right) 
         << setprecision(4)
         << setw(12) << p[i];
  cout << endl << endl;
}

 

[Bill Simpson]

In the nested loops following
// convert to upper triangular form
ensure, if possible, that the k'th entry in the k'th row is nonzero.
One way this could be done, if possible, is by swapping a pair of rows.