Matrix inversion through EROs

[tcdnewb]

Hello, I would appreciate some help with the following, the questions refer to the code below;


(i) Change the read matrix() routine so that it reads the dimensions m,n rather than being
passed them by the calling program. The specimen program overload.cpp (`Reference parameters.
. . ') will show the implications of this change.

(ii) Insert a routine
void print_matrix ( int m, int n, double a[10][10] )
which prints the matrix in a decent way. That is, the eld width should be 10, the precision 4, and the
notation ios::fixed

(iii) Delete the statement from the main program in which m and n are read. Now they are read
by the read matrix() routine.

(iv) Delete the existing output statements from main() and replace them with calls to your
print matrix() routine

#include <iostream>
#include <iomanip>

using namespace std;

void augment ( int m, double a[10][10], double b[10][20] )
{
  int i, j;
  for ( i=0; i<m; ++i )
  {
    for (j=0; j<m; ++j )
    {
      b[i][j] = a[i][j];
      b[i][j+m] = i == j;
    }
  }
}

void peeloff ( int m, double a[10][20], double c[10][10] )
{
  int i, j;
  for ( i=0; i<m; ++i )
  {
    for (j=0; j<m; ++j )
    {
      c[i][j] = a[i][j+m];
    }
  }
}


void swap ( int r, int s, int m, int n, double a[10][20] )
{
  int j;
  double temp;

  if ( r == s )
    return;

  for ( j=0; j<n; ++j )
  {
    temp = a[r][j];
    a[r][j] = a[s][j];
    a[s][j] = temp;
  }
}

void scale ( int r, double s, int m, int n, double a[10][20] )
{
  int j;
  for ( j=0; j<n; ++j )
  {
    a[r][j] *= s;
  }
}

void subtract ( int r, double s, int t, int m, int n, double a[10][20] )
{
  int j;
  for (j=0; j<n; ++j)
  {
    a[r][j] -= s * a[t][j];
  }
}

void multiply (int ell, int m, int n,  double a[10][10], double b[10][10],
        double c[10][10])
{
  int i,j,k;
  double tot;

  for (i=0; i<ell; ++i)
  {
    for (k=0; k < n; ++k )
    {
      tot = 0;
      for ( j=0; j<m; ++j )
      { tot += a[i][j] * b[j][k]; }
      c[i][k] = tot;
    }
  }
}

void read_matrix ( int m, int n,  double a[10][10] )
{
  int i,j;
  for ( i=0; i<m; ++i )
  {
    for ( j=0; j<n; ++j )
    {
      cin >> a[i][j];
    }
  }
}

void gje ( int m, int n, double a[10][20] )
{
  int i, j, k, pivot;
  double s;

  k = -1;

  for ( j=0; j<n; ++j )
  {
    pivot = -1;
    for ( i=k+1; pivot < 0 && i<m; ++i )
    {
      if ( a[i][j] != 0 )
	pivot = i;
    }

    if ( pivot >= 0 )
    {
      ++k;
      swap ( pivot, k, m, n,  a );
      s = 1 / a[k][j];
      scale ( k, s, m, n, a );
      for ( i=0; i<m; ++i )
      {
        if ( i != k )
	{
	  s = a[i][j];
	  subtract ( i, s, k, m, n, a );
        }
      }
    }
  } 

 cout << "gje finished k " << k << endl;
}

main()
{
  double a[10][10], b[10][20], c[10][10], d[10][10];

  int m, n, i,j;

  cin >> m >> n;

  read_matrix ( m,n,a );

  augment ( m, a, b );
  gje ( m, 2*m, b );
  peeloff ( m, b, c );
  multiply (m,m,m,a,c,d);

  cout << "Matrix\n";

  for ( i=0; i<m; ++i )
  {
    for ( j=0; j<m; ++j )
    {
      cout << setprecision ( 10 );
      cout <<  "  " << a[i][j];
    }
    cout << endl;
  }

  cout << "\nInverse\n";

  for ( i=0; i<m; ++i )
  {
    for ( j=0; j<m; ++j )
    {
      cout << setw ( 10 );
      cout << setprecision ( 10 );
      cout << "  " <<  c[i][j];
    }
    cout << endl;
  }

  cout << "\nProduct\n";
  for ( i=0; i<m; ++i )
  {
    for ( j=0; j<m; ++j )
    {
      cout << setw ( 10 );
      cout << setprecision ( 10 );
      cout << "  " <<  d[i][j];
    }
    cout << endl;
  }

  return 0;
}