#include <iostream.h>
#include <stdlib.h>


// This program calculates the inverse of a square matrix.

int main()
{

 cout << "This program calculates the inverse of a square matrix."<<endl<< endl;

//////////////////////////////////////////////////////
// Get the dimensions of the matrix:
//////////////////////////////////////////////////////

  cout << endl << "Please enter the dimension of your matrix: ";
  int N;
  cin >> N;

  float A[N+1][N+1];

//////////////////////////////////////////////////////
// Get the matrix elements:
//////////////////////////////////////////////////////

   cout << endl << "Plese enter the matrix elements:" << endl;

   int k; int l;

   for(k=1;k<=N;k++)
   {
    for(l=1;l<=N;l++)
    {
    cout << "A[" << k << "," << l << "]=";
    cin >> A[k][l];
    }
   }
   cout << endl;

//////////////////////////////////////////////////////
// Print the matrix:
//////////////////////////////////////////////////////

  cout << "Your matrix is" << endl << endl;

    for(k=1;k<=N;k++)
   {
    for(l=1;l<=N;l++)
    {
    cout << A[k][l] << "   ";
    }
    cout << endl << endl;
   }

   system("PAUSE");
   cout << endl;

//////////////////////////////////////////////////////
// Form the augmented matrix
//////////////////////////////////////////////////////

   float Aug[N+1][2*N+1];

   for(k=1;k<=N;k++)
   {
     for(l=1;l<=N;l++)
     {
       Aug[k][l]=A[k][l];
       if(k==l){Aug[k][N+l]=1;}
       else{Aug[k][N+l]=0;}
     }
   }

//////////////////////////////////////////////////////
// Print the Augmented matrix:
//////////////////////////////////////////////////////

   cout << "The corresponding Augmented matrix is" << endl << endl;

    for(k=1;k<=N;k++)
   {
    for(l=1;l<=2*N;l++)
    {
    cout << Aug[k][l] << "   ";
    }
    cout << endl << endl;
   }

   system("PAUSE");
   cout << endl;

//////////////////////////////////////////////////////
// Now start the Gaussian algorithm:
//////////////////////////////////////////////////////

//a) First row echelon form:

   int d;

   float M[2*N+1];

   for(d=1;d<=N;d++)               // Start at the d th row.
   {                               //
     while(Aug[d][d]==0)           // If A[d][d]==0 then put that row
     {                             // at the bottom. This way one eventually
       for(l=1;l<=2*N;l++)         // gets a nonzero pivot element if the d th
       {                           // column is not zero altogether.
         M[l]=Aug[d][l];
       }
        for(k=d;k<N;k++)
        {
          for(l=1;l<=2*N;l++)
          {
            Aug[k][l]=Aug[k+1][l];
          }
        }
        for(l=1;l<=2*N;l++)
        {
          Aug[N][l]=M[l];
        }
     }                              // End of resorting.

     for(l=2*N; l>=d; l--)             // Divide the d th row by the pivot element.
     {
        Aug[d][l]=Aug[d][l]/Aug[d][d];
     }

     for(k=d+1;k<=N;k++)               // Use the pivot at the d th
     {                                 // row to cancel the terms
        for(l=2*N; l>=d; l--)          // under it.
        {
         Aug[k][l]=Aug[k][l]-Aug[k][d]*Aug[d][l];
        }
     }                                 // End of business with d th row.
                                       // Now go back and do it for d+1 st row.
   }


// b) Now the reduced row echelon form:

 for(d=N; d>1; d--)
 {
   for(k=d-1;k>=1;k--)           // Use the pivot at row d to cancel the terms
   {                             // above it.
     for(l=2*N; l>=1; l--)
     {
       Aug[k][l]=Aug[k][l]-Aug[k][d]*Aug[d][l];
     }
   }
 }

//////////////////////////////////////////////////////
// Write down the reduced matrix:
//////////////////////////////////////////////////////

   cout << "By applying Gaussian algorithm, we find" << endl << endl;

   for(k=1;k<=N;k++)
   {
     for(l=1;l<=2*N;l++)
     {
       cout << Aug[k][l] << "   ";
     }
     cout << endl << endl;
   }

   system("PAUSE");
   cout << endl;

//////////////////////////////////////////////////////
// Register the inverse of the matrix:
//////////////////////////////////////////////////////

   float Inv[N+1][N+1];

   for(k=1;k<=N;k++)
   {
     for(l=1;l<=N;l++)
     {
       Inv[k][l]=Aug[k][N+l];
     }
   }

//////////////////////////////////////////////////////
// Write down the inverse of the matrix:
//////////////////////////////////////////////////////

   cout << "The inverse of the matrix is" << endl << endl;

   for(k=1;k<=N;k++)
   {
     for(l=1;l<=N;l++)
     {
       cout << Inv[k][l] << "   ";
     }
     cout << endl << endl;
   }
   system("PAUSE");
   cout << endl;

//////////////////////////////////////////////////////
// Control the result:
//////////////////////////////////////////////////////


  cout << endl
  << "Multiplying your matrix with its inverse that I have found, we obtain"
  << endl << endl;

  float P[N+1][N+1];
  int m;

   for(k=1; k<=N; k++)
   {
     for(l=1; l<=N; l++)
     {
      P[k][l]=0;
      for(m=1; m<=N; m++)
      {
        P[k][l]+=A[k][m]*Inv[m][l];
      }
     }
   }

  for(k=1;k<=N;k++)
   {
     for(l=1;l<=N;l++)
     {
       cout << P[k][l] << "   ";
     }
     cout << endl << endl;
   }

//////////////////////////////////////////////////////
// End of Program:
//////////////////////////////////////////////////////

      cout << endl;
      system("PAUSE");
      return 0;
}