diff -N matrix.c
--- matrix.c	7 Jul 2010 16:10:46 -0000	1.1
+++ /dev/null	1 Jan 1970 00:00:00 -0000
@@ -1,157 +0,0 @@

-// 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 << this->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++ ) {
-   this->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" << std::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