/***********************************************
* Name:       matrix_product-bugs.cxx
* Author:     Leo Liberti
* Source:     GNU C++
* Purpose:    multiply two matrices: bug (while sbuf.eof -> while !sbuf.eof)
* Build:      c++ -o matrix_product-bugs matrix_product-bugs.cxx
* History:    120117 work started
***********************************************/

#include<iostream> 
#include<fstream>
#include<cstdlib>
#include<vector>
#include<string>
#include<sstream>

// read rows x columns matrix A from file stream f
void readMatrix(std::ifstream& f, std::vector<std::vector<double> >& A) {
  using namespace std;
  string buffer;
  while(!f.eof()) {
    getline(f, buffer);
    if (buffer.size() > 0) {
      stringstream sbuf;
      sbuf << buffer;
      vector<double> theRow;    
      while(sbuf.eof()) {
	double t;
	sbuf >> t;
	theRow.push_back(t);
      }
      A.push_back(theRow);
    }
  }
}

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

// print a matrix on a stream
void printMatrix(std::ostream& fout, 
		 std::vector<std::vector<double> >& A) {
  using namespace std;
  for(int i = 0; i < A.size(); i++) {
    for(int j = 0; j < A[i].size() ; 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;
  f1.open(argv[1]);
  vector<vector<double> > A;
  readMatrix(f1, A);
  cout << "A(" << A.size() << "," << A[0].size() << ") = " << endl;
  printMatrix(cout, A);

  // second matrix
  ifstream f2;
  f2.open(argv[2]);
  vector<vector<double> > B;
  readMatrix(f2, B);
  cout << "B(" << B.size() << "," << B[0].size() << ") = " << endl;
  printMatrix(cout, B);

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

  // product
  vector<vector<double> > C;
  for(int i = 0; i < A.size(); i++) {
    vector<double> theRow;
    for(int j = 0; j < B[0].size(); j++) {
      theRow.push_back(scalarProduct(A, B, i, j));
    }
    C.push_back(theRow);
  }
  cout << "A*B(" << C.size() << "," << C[0].size() << ") = " << endl;
  printMatrix(cout, C);

  return ret;
}