// matrix.c

#ifndef matrix_C
#define matrix_C

#include "matrix.h"
#include "trfutil/trfstream.h"

// output stream
template <class T>
void matrix<T>::ostr(std::ostream& stream) const {
  stream << begin_object;
  int i = 0;
  for ( int irow=0; irow<nrow(); irow++ ) {
    if ( irow ) stream << new_line;
    for ( int icol=0; icol<ncol(); icol++ ) {
      if ( icol > 0 ) stream << ' ';
      stream << data[i++];
    }
  }
  stream << end_object;
}

//
// conversion from symmetric matrix to matrix
//
template <class T>
matrix<T>::matrix(const smatrix<T>& sma)
: array<T>(sma.nrow()*sma.nrow()) {
 _nrow = sma.nrow();
 _ncol = sma.nrow();
 for ( int irow=0; irow<_nrow; irow++ ) {
  for ( int icol=0; icol<_ncol; icol++ ) {
   data[ij(irow,icol)] = sma(irow,icol);
  }
 }
}

//
// mutiplication: matrix*matrix
//
template <class T>
matrix<T> operator*(const matrix<T>& mat1,const matrix<T>& mat2)
{
 if ( mat1.ncol() != mat2.nrow() ) {
  matrix<T> prod(0,0);
  return prod;
 }
 int nrow = mat1.nrow();
 int ncol = mat2.ncol();
 matrix<T> prod(nrow,ncol);
 for ( int row=0; row<nrow; row++ ) {
  for ( int col=0; col<ncol; col++ ) {
   T sum = 0.0;
   for ( int i=0; i<mat1.ncol(); i++ ) {
    sum += mat1(row,i) * mat2(i,col);
   }
   prod(row,col) = sum;
  }
 }
 return prod;
}

//
// mutiplication: matrix = smatrix*smatrix
//
template <class T>
matrix<T> operator*(const smatrix<T>& sma1,const smatrix<T>& sma2)
{
 if ( sma1.nrow() != sma2.nrow() ) {
  std::cerr << "operator*(sma,sma): nrow != ncol\n" << flush;
  matrix<T> prod(0,0);
  return prod;
 }
 int nrow = sma1.nrow();
 matrix<T> prod(nrow,nrow);
 for ( int row=0; row<nrow; row++ ) {
  for ( int col=0; col<nrow; col++ ) {
   T sum = 0.0;
   for ( int i=0; i<nrow; i++ ) {
    sum += sma1(row,i) * sma2(i,col);
   }
   prod(row,col) = sum;
  }
 }
 return prod;
}

//
// mutiplication: matrix*vector
//
template <class T>
nvector<T> operator*(const matrix<T>& mat1,const nvector<T>& vec2)
{
 if ( mat1.ncol() != vec2.length() ) {
  nvector<T> prod(0);
  return prod;
 }
 nvector<T> prod( mat1.nrow() );
 for ( int row=0; row<mat1.nrow(); row++ ) {
  T sum = 0.0;
  for ( int i=0; i<vec2.length(); i++ ) {
   sum += mat1(row,i) * vec2(i);
  }
  prod(row) = sum;
 }
 return prod;
}

//
// multiplication: M % S == M * S * MT
//
template <class T>
smatrix<T> operator%(const matrix<T>& mat, const smatrix<T>& sma)
{
  if ( mat.ncol() != sma.nrow() ) {
    smatrix<T> prod(0);
    return prod;
  }
  int dim1 = mat.ncol();
  int dim2 = mat.nrow();
  smatrix<T> prod( dim2 );
  if (dim2 == 1) {
    // N=1 case.
    T* smdata = sma.data;
    T* matdata = mat.data;
    int ii = 0;
    T sum = 0;
    for (int j=0; j<dim1; j++) {
      T sum1 = 0;
      for (int k=0; k<j; k++)
        sum1 += 2 * smdata[ii++] * matdata[k];
      sum += (sum1 + smdata[ii++] * matdata[j]) * matdata[j];
    }
    prod.data[0] = sum;
  }
  else {
    // General case.
    for ( int i=0; i<dim2; i++ ) {
      for ( int l=0; l<=i; l++ ) {
        T sum = 0.0;
        int smai = 0;
        T* mati = mat.data + dim1 * i;
        T* matl = mat.data + dim1 * l;
        for ( int j=0; j<dim1; j++ ) {
          for ( int k=0; k<j; k++ )
            sum += sma.data[smai++] * (mati[j]*matl[k] + mati[k]*matl[j]);
          sum += sma.data[smai++] * (mati[j]*matl[j]);
        }
        prod(i,l) = sum;
      }
    }
  }
  return prod;
}

#endif