//
// 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;
}