Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Alterring a program

  1. Jun 6, 2010 #1
    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 (Text):
    #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;
    }
     
     
  2. jcsd
  3. Jun 7, 2010 #2
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Alterring a program
  1. Matlab programming (Replies: 0)

  2. Mathematical Programs (Replies: 1)

  3. Programming in Matlab (Replies: 4)

Loading...