/***********************************************
* Name:       matrix_product-ptr-bugs.cxx
* Author:     Leo Liberti
* Source:     GNU C++
* Purpose:    multiply two matrices, use pointers rather than std template lib
*             show bugs: in scalarProduct (A[i][k]*B[j][k] -> A[i][k]*B[k][j])
*                        in deallocateMemory (loop on cols instead of rows)
* Build:      c++ -o matrix_product-ptr matrix_product-ptr.cxx
* History:    120117 work started
***********************************************/

#include<iostream> 
#include<fstream>
#include<cstdlib>

// allocate rows x columns matrix memory to A
void allocateMemory(double**& A, int rows, int columns) {
  A = new double* [rows];
  for(int i = 0; i < rows; i++) {
    A[i] = new double [columns];
  }
}

// free rows x columns matrix memory from A
void deallocateMemory(double**& A, int rows, int columns) {
  for(int i = columns; i > 0; i--) {
    delete[] A[i-1];
  }
  delete [] A;
}

// read rows x columns matrix A from file stream f
void readMatrix(std::ifstream& f, int rows, int columns, double** A) {
  for(int i = 0; i < rows; i++) {
    for(int j = 0; j < columns; j++) {
      f >> A[i][j];
    }
  }
}

// scalar product of vectors row(A,i) and col(B,j) of size n
double scalarProduct(double** A, double** B, int i, int j, int n) {
  double ret = 0;
  for(int k = 0; k < n; k++) {
    ret += A[i][k] * B[j][k];
  }
  return ret;
}

// print a matrix on a stream
void printMatrix(std::ostream& fout, double** A, int rows, int columns) {
  using namespace std;
  for(int i = 0; i < rows; i++) {
    for(int j = 0; j < columns ; j++) {
      fout << A[i][j] << " ";
    }
    fout << endl;
  }
}

// main procedure
int main(int argc, char** argv) {
  using namespace std;
  int ret = 0;

  // test command line
  if (argc < 3) {
    cerr << argv[0] << ": multiply two matrices found in files given on "
	 << "cmd line" << endl;
    cerr << "  syntax is " << argv[0] << " file1 file2" << endl;
    exit(1);
  }

  // first matrix
  ifstream f1;
  int rows1 = 0;
  int columns1 = 0;
  f1.open(argv[1]);
  f1 >> rows1 >> columns1;
  double** A;
  allocateMemory(A, rows1, columns1);
  readMatrix(f1, rows1, columns1, A);
  cout << "A(" << rows1 << "," << columns1 << ") = " << endl;
  printMatrix(cout, A, rows1, columns1);


  // second matrix
  ifstream f2;
  int rows2 = 0;
  int columns2 = 0;
  f2.open(argv[2]);
  f2 >> rows2 >> columns2;
  double** B;
  allocateMemory(B, rows2, columns2);
  readMatrix(f2, rows2, columns2, B);
  cout << "B(" << rows2 << "," << columns2 << ") = " << endl;
  printMatrix(cout, B, rows2, columns2);

  // verify size consistency
  if (columns1 != rows2) {
    cerr << argv[0] << ": A with " << columns1 << " cols and B with " 
	 << rows2 << " rows can't be multiplied" << endl;
    exit(4);
  }

  // product
  int n = columns1;
  double** C;
  allocateMemory(C, rows1, columns2);
  for(int i = 0; i < rows1; i++) {
    for(int j = 0; j < columns2; j++) {
      C[i][j] = scalarProduct(A, B, i, j, n);
    }
  }
  cout << "A*B(" << rows1 << "," << columns2 << ") = " << endl;
  printMatrix(cout, C, rows1, columns2);

  // deallocate space
  deallocateMemory(A, rows1, columns1);
  deallocateMemory(B, rows2, columns2);
  deallocateMemory(C, rows1, columns2);  

  return ret;
}