Matrix Inversion (using partial pivoting)

//

//  main.cpp

//  2015 New Start

//

//  Created by zr9558 on 13/2/15.

//  Copyright (c) 2015 zr9558. All rights reserved.

//

#include <iostream>

using namespace std;

#include<math.h>

// 5.3 Algorithm Matrix_Inversion (using partial pivoting)

const double eps=1e-10;

int max1( double b[], int k, int n) // find the pivot from b[k] to b[n-1], absolute value

{

    double max=fabs(b[k]); int index=k;

    for( int i=k; i!=n; ++i)

        if( fabs(b[i])>max) { max=fabs(b[i]); index=i;}

    

    return index;

}

int main()

{

    double a[100][200], b[100];

    int n;

    cout<<"Enter the dimension of the matrix"<<endl;

    cin>>n;

    cout<<"Enter the matrix"<<endl;

    for( int i=0; i!=n; ++i)

        for( int j=0; j!=n; ++j)

            cin>>a[i][j];

    

    for( int i=0; i!=n; ++i) //Augment the matrix A

    {

        for( int j=n; j!=2*n; ++j)

        {

            if( j-i==n) a[i][j]=1;

            else a[i][j]=0;

        }

    }

    

    // Make upper triangular form;

    

    for( int k=0; k!=n; ++k)

    {

        for( int i=k; i!=n; ++i)

            b[i]=a[i][k];

        

        int j=max1(b,k,n); // find the pivot from the elements a[k][k],...,a[n-1][k] in the kth column of A and let a[j][k] be the pivot;

        

        if( fabs(a[j][k])<=eps) // if a[j][k]==0, then not invertible;

        {

            cout<<"not invertible"<<endl;

            exit(0);

        }

        

        if( j!=k)

        {

            for( int i=0; i!=2*n; ++i)

            {

                double temp=a[j][i]; a[j][i]=a[k][i]; a[k][i]=temp;

            }

        }

        

        if( a[k][k]!=1) // if a[k][k]!=1, then divide all the elements of kth row by a[k][k],

        {

            double temp=a[k][k];

            for( int i=0; i!=2*n; ++i)

                a[k][i]/=temp;

        }

        

        for( int j=k+1; j!=n; ++j) // substract a[j][k] times the kth row to the jth row for j=k+1,...,n-1;

        {

            double temp=a[j][k];

            for( int i=0; i!=2*n; ++i)

                a[j][i]-=temp*a[k][i];

        }

    }

    

    

    // Make the left half of A a unit matrix;

    

    for( int k=1; k!=n; ++k)

    {

        for( int j=k-1; j>=0; --j) // substract a[j][k] times the kth row to the jth row for j=k-1,...,0;

        {

            double temp=a[j][k];

            for( int i=0; i!=2*n; ++i)

                a[j][i]-=temp*a[k][i];

        }

    }

    

    

    //output the inverse matrix

    cout<<"Inverse Matrix is"<<endl;

    for( int j=0; j!=n; ++j)

    {

        for( int i=n; i!=2*n; ++i)

        {

            cout<<a[j][i]<<" ";

        }

        cout<<endl;

    }

    

    

    

    return 0;

}