#include "matrix.h"

void error(const char* s)
{// remark and quit
  cerr << s;
  exit(0);
}



Matrix::Matrix(int r, int c){ // set up a matrix of zeros
   if (r<=0 || c<=0)error("bad matrix indeces \n");
   row=r;
   col=c;
   m=new double[r*c];
   for(int i=0; i<r*c; i++) m[i]=0.0;
 };

Matrix::Matrix(const Matrix& A){ // build a matrix from another one
   row = A.row;
   col = A.col;
   m = new double[row*col];
   for(int i=0; i< row*col; i++) m[i] =A.m[i];
 };

Matrix operator+(const Matrix& A, const Matrix& B){ // add two Matrices
 if (A.row!=B.row || A.col!=B.col)error("Matrix+: incomatable sizes \n");
 Matrix C(A.row, A.col);
 int i,rc=A.row*A.col;
 for(i=0; i<rc; i++)
    C.m[i] = A.m[i] + B.m[i];
 return C;
}

Matrix operator*(const Matrix& A, const Matrix& B){
 // Matrix multiplication
 if (A.col!=B.row)error("Matrix*: incomatable sizes \n");
 Matrix C(A.row, B.col);
 int i,j,k;
 for(i=0; i<A.row; i++)
    for(j=0; j<B.col; j++)
       for(k=0; k<A.col; k++)
          C.m[i*C.col+j]+=A.m[i*A.col + k] * B.m[k*B.col+j];
 return C;
}

Matrix operator-(const Matrix& A, const Matrix& B){
 // Matrix subtraction
 if (A.row!=B.row || A.col!=B.col)error("Matrix-: incomatable sizes \n");
 Matrix C(A.row, A.col);
 int i,rc=A.row*A.col;
 for(i=0; i<rc; i++)
    C.m[i]=A.m[i] - B.m[i];
 return C;
}

Matrix operator*(double a, const Matrix& A){
 // scalar times a matrix
 Matrix C(A.row, A.col);
 int i,rc=A.row*A.col;
 for(i=0; i<rc; i++)
    C.m[i]=a*A.m[i];
 return C;
}

Matrix operator*( const Matrix& A, double a){
 // scalar times a matrix
 return a*A;
}

Matrix operator&(const Matrix& A, const Matrix& B){
 if (A.row!=B.row || A.col!=B.col)error("Matrix&: incomatable sizes \n");
 Matrix C(A.row, A.col);
 int i,rc=A.row*A.col;
 for(i=0; i<rc; i++)
    C.m[i]=A.m[i]*B.m[i];
 return C;
}

ostream& operator<<(ostream& os,Matrix A){
 // output for small matrices
 int i, j;
 for(i=0; i<A.row; i++){
    for(j=0; j<A.col; j++)
       os << A.m[i*A.col+j] << "\t";
    os << "\n";
   }
 return os;
}

Matrix transpose (const Matrix& A){
 Matrix C(A.col, A.row);
 int i,j;
 for(i=0; i<A.row; i++)
    for(j=0; j<A.col; j++)
       C.m[j*C.col + i]=A.m[i*A.col + j];
    return C;   
}

Matrix solve(const Matrix& A, const Matrix& B){
 /*
   solve A x = B and return the answer.
   It is assumed that A can be factored without pivoting.
   The factors are not saved
 */
 if( A.row!=A.col || A.row!=B.row )
    error("Matrix solve: incompatible sizes\n");

 Matrix C(A),X(B); // local copies of A and B
 int i,j,k,cr=A.row,xc=B.col;
 double a,b;

 // First preprocess C
 for(i=0; i<cr; i++){ // for each row of C
    if( C.m[i*cr+i] == 0.0 ) error(" zero divide in solve ");
    a = (C.m[i*cr+i] = 1./C.m[i*cr+i]); 
    for( j=i+1; j<cr; j++ )  // normalize the current row
       C.m[i*cr+j] *=a;
    for( k=i+1; k<cr; k++ ){// for each row after i modify the row
       if((b = C.m[k*cr+i])!=0.0){ // but only if it is needed.
          for( j=i+1; j<cr; j++ ) 
             C.m[k*cr+j] -= b*C.m[i*cr+j];
       }
    }
 }

 // now solve the equations
 for(k=0; k<xc; k++){ // for each column of X
    for(i=0; i<cr; i++ ){ // forward solution
       for( j=0; j<i; j++ )
          X.m[i*xc+k] -= C.m[i*cr+j]*X.m[j*xc+k];
       X.m[i*xc+k] *= C.m[i*cr+i];
    }
    for(i=cr-2; i>=0; i--) // back solution
       for(j=i+1; j<cr;j++ )
          X.m[i*xc+k] -= C.m[i*cr+j]*X.m[j*xc+k];
 }

 return X;

}