trf++/include/trfutil
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