- #1
philosafari
- 1
- 0
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.
Code:
#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;
}