// *********************************************************************
//
// a C++ class of a matrix 
//
// In the case of overloaded mathematical operators, no size-checks are
// performed to gain speed.
//
// Variables:
//   cv        array of Complex objects
//   size		length of the vector
//
// *********************************************************************

#include <assert.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include "DVector.h"
#include "DMatrix.h"

#define size 3

    // ================================================================
	DMatrix::DMatrix(int r,int c) {
          init(r,c);
        }
    //
    // ================================================================
	DMatrix::DMatrix(double **x) {
	  //cout<<"In DMatrix::DMatrix(double **x)\n";
	  init(size,size);
          getValue(x);
        }
    // ================================================================
	DMatrix::DMatrix(double x[size][size]) {
	  //cout<<"In DMatrix::DMatrix(double x[3][3])\n";
	  init(size,size);
          getValue(x);
        }
	
    // ================================================================
	DMatrix::DMatrix(DVector &col1,DVector &col2, 
	                 DVector &col3)
        {
	  init(size,size);
	  for (int i=0; i<size; i++){
	  /*
	    cout << "col1[i]=" << col1[i] <<"\n";
	    cout << "col2[i]=" << col2[i] <<"\n";
	    cout << "col3[i]=" << col3[i] <<"\n";
	   */
	    cv[i][0]=col1[i];
	    //cout << "cv[0][1]="<<cv[0][1]<<"\n";
	    cv[i][1]=col2[i];
	    cv[i][2]=col3[i];
  	  }  
	}

    // ================================================================
        DMatrix::DMatrix(const DMatrix &matrix) {
	  init(size,size);
          getValue(matrix.cv);
     }

     // ================================================================
	void DMatrix::init(int nr, int nc) {
		//cout<<"In DMatrix::init(int nr, int nc)\n";
		cv = (double**) malloc(nr*nc* sizeof(double));
		*cv =new double[nr];
	        for (int i=0; i<nr; i++) {
		    cv[i] = new double[nc];
		}	
	}
	
     // ================================================================
     void DMatrix::getValue(double (*a)[3]) {
	  for (int i=0; i<size; i++) {
	    for (int j=0; j<size; j++) {
		cv[i][j] = a[i][j];
	    }
	  }
     }
     
     void DMatrix::getValue(double **array) {
          if (array != 0) {
		for (int i=0; i<size; i++) {
		 for (int j=0; j<size; j++) {
               	  cv[i][j] = array[i][j];
          	 }
		}
	  }
	  else {
		  for (int i=0; i<size; i++) {
		    for (int j=0; j<size; j++) {
			cv[i][j] = double(0.0);
		  }
		}
	  }	
     }

     // ================================================================
     DMatrix::~DMatrix() { 
        for (int i=0; i<size; i++){
	   delete cv[i];
	}   
        delete [] cv; 
     }      // destructor

     // ================================================================
     DMatrix& DMatrix::operator=( const DMatrix &cV) {
          if (this == &cV) return *this;
        for (int i=0; i<size; i++){
	   delete cv[i];
	}   
        delete [] cv; 
	init(size,size);
	getValue(cV.cv);
     }

     // ================================================================
     double& DMatrix::operator( )(int row, int column) {
          return cv[row][column];
     }


     // === Multiplication by a vector================================
	DVector DMatrix::operator * (DVector &vector) {
	        DVector result(size);
		
           /*	for(int i=0; i<size; i++) {
		   cout << vector[i]<<"\n";
	        }
		
	    */	
		for(int i=0; i<size; i++) {
		  for (int j=0; j<size; j++){
			result[i] += cv[i][j] * vector[j];
		  }	
		}
           /*	for(int i=0; i<size; i++) {
		   cout << result[i]<<"\n";
	        }
		
	    */	
		return result;
	}

     // === Multiplication by another matrix================================
	DMatrix DMatrix::operator * (DMatrix &matrix) {
	        DMatrix result(size,size);
/*		
       cout << "Test results of cv[0][0]=" << cv[0][0]<<"\n";
       cout << "Test results of cv[0][1]=" << cv[0][1]<<"\n";
       cout << "Test results of cv[0][2]=" << cv[0][2]<<"\n";
       cout << "Test results of cv[1][0]=" << cv[1][0]<<"\n";
       cout << "Test results of cv[1][1]=" << cv[1][1]<<"\n";
       cout << "Test results of cv[1][2]=" << cv[1][2]<<"\n";
       cout << "Test results of cv[2][0]=" << cv[2][0]<<"\n";
       cout << "Test results of cv[2][1]=" << cv[2][1]<<"\n";
       cout << "Test results of cv[2][2]=" << cv[2][2]<<"\n";
		
       cout << "Test results of matrix(0,0)=" << matrix(0,0)<<"\n";
       cout << "Test results of matrix(0,1)=" << matrix(0,1)<<"\n";
       cout << "Test results of matrix(0,2)=" << matrix(0,2)<<"\n";
       cout << "Test results of matrix(1,0)=" << matrix(1,0)<<"\n";
       cout << "Test results of matrix(1,1)=" << matrix(1,1)<<"\n";
       cout << "Test results of matrix(1,2)=" << matrix(1,2)<<"\n";
       cout << "Test results of matrix(2,0)=" << matrix(2,0)<<"\n";
       cout << "Test results of matrix(2,1)=" << matrix(2,1)<<"\n";
       cout << "Test results of matrix(2,2)=" << matrix(2,2)<<"\n";
     */  
       double sum=0.;
		for(int i=0; i<size; i++) {
		  for (int j=0; j<size; j++){
		  sum =0.;
		    for (int k=0; k<size; k++){
	//	        cout<<" i,k,j,cv[i][k] * matrix(k,j)="<<i<<k<<j<<","<<cv[i][k]<<"*"<<matrix(k,j)<<"\n";
			sum += cv[i][k] * matrix(k,j);
		    }
		   result(i,j)=sum; 
		   //cout<<"i,j,result(i,j)="<<i<<j<<" "<<result(i,j)<<"\n"; 
		  }
		}
		return result;
	}