Gauss Elimination with Partial Pivoting


C++ SOURCE CODE FOR GAUSS ELIMINATION WITH PARTIAL PIVOTING



/**
    C++ code to solve system of linear equation of provided order
    by Gauss elimination with partial pivoting
*/


#include <iostream>
#include <algorithm>  //for swapping equation which is a horizontal vector
#include <cmath> //used for absolute value
#include <vector> //used as container for augmented matrix
#include <iomanip> //for formatting output

using namespace std;


/**
    This function read equation from input and assign it to 2-D vector

    @param augMatrix --> reference of 2-D vector with allocated size
    @return nothing
*/
void readSysEq(vector<vector<double> >& augMatrix) {
    cout << "Enter equation coeff in the form 'ax + by + ..... = c': " << endl;
    for(int i = 0; i < augMatrix.size(); i++) {
        cout << "Enter equation " << i + 1 << " :";
        for(int j = 0; j <= augMatrix.size(); j++)
            cin >> augMatrix[i][j];
    }
}

/**
    This function print the augmented matrix
    @param augMatrix --> copy of 2-D vector
*/
void printSysEq(vector<vector<double> > augMatrix) {
    cout << endl;
    for (int i=0; i < augMatrix.size(); i++) {
        for (int j=0; j <= augMatrix.size(); j++) {
            cout << setw(5) << augMatrix[i][j];
            if (j == augMatrix.size()-1) {
                cout << "| "; //separating rhs of equation from lhs
            }
        }
        cout << "\n";
    }
    cout << endl;
}

/**
    This function print the solution vector

    @param sol --> a copy of solution vector
*/
void printSolVector(vector<double> sol) {
    cout << "["; //printing in the format "[.....]"
    for(int i = 0; i < sol.size(); i++)
        cout << sol[i] << " ";
    cout << "]" << endl;
}


/**
    Our main function to solve given augmented matrix by gauss
    elimination with partial pivoting

    @param augMatrix --> 2-D container(in our case vector) as augmented matrix
    @return sol --> solution vector, a 1-D vector
*/
vector<double> gauss(vector<vector<double> > augMatrix) {
    int numOfEq = augMatrix.size(); //calculating number of provided equation

    //i is the iterator pointing a column
    for(int i = 0; i < numOfEq; i++) {

        //find max absolute coeff in the ith column
        int maxRowIndex = i;
        for(int j = i + 1; j < numOfEq; j++)
            if(abs(augMatrix[j][i]) > abs(augMatrix[maxRowIndex][i]))
                maxRowIndex = j;

        /**
            swap current row with max absolute pivot column coeff row
            One can use their own swapping technology instead of using
            swap function in algorithm class.
        */
        swap(augMatrix[i], augMatrix[maxRowIndex]);

        /**
            Reducing all coefficient below pivot element augMatrix[i][i] to zero
            Here, we don't need to worry about lower triangular portion of coefficient matrix
            which ultimately gonna be reduced to zero.
            So only compute in upper triangular matrix

            Now when we are in ith column we need jth row element in pivot column
            to be zero so multiply pivot row by factor (a[j][i] / a[i][i])
            subtract it with j row
        */
        for(int j = i + 1; j < numOfEq; j++)
            for(int k = i + 1; k <= numOfEq; k++)
                augMatrix[j][k] -= augMatrix[j][i] * augMatrix[i][k] / augMatrix[i][i];
    }

    //performing backward substitution
    vector<double> sol(numOfEq);
    for (int i = numOfEq - 1; i >= 0; i--) {
        sol[i] = augMatrix[i][numOfEq] / augMatrix[i][i];
        //substitute immediate solution backward and update RHS of each equation
        for (int j = i - 1; j >= 0; j--)
            augMatrix[j][numOfEq] -= augMatrix[j][i] * sol[i];
    }

    return sol;
}


int main() {
    int numOfEq;
    cout << "Enter number of equation(s): ";
    cin >> numOfEq;
    vector<vector<double> > augMatrix(numOfEq, vector<double>(numOfEq + 1));
    readSysEq(augMatrix);
    cout << "Given system of equation: ";
    printSysEq(augMatrix);
    cout << "The solution is: ";
    printSolVector(gauss(augMatrix));
    return 0;
}

Comments

Post a Comment

Need Help? Make a comment below & we'll help you out... :)

Popular posts from this blog

How to write a cover letter to register free .com.np in Nepal

Add a Comment Form Message on Blogger