matrix_source.h Source File

00001 // matrix.h
00002 
00003 //  This matrix class is a C++ wrapper for the GNU Scientific Library
00004 //  Copyright (C) 2001 Ramin Nakisa
00005 
00006 //  This program is free software; you can redistribute it and/or modify
00007 //  it under the terms of the GNU General Public License as published by
00008 //  the Free Software Foundation; either version 2 of the License, or
00009 //  (at your option) any later version.
00010 
00011 //  This program is distributed in the hope that it will be useful,
00012 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014 //  GNU General Public License for more details.
00015 
00016 //  You should have received a copy of the GNU General Public License
00017 //  along with this program; if not, write to the Free Software
00018 //  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00019 
00020 #if !defined( _matrix_#type#_h )
00021 #define _matrix_#type#_h
00022 
00023 #include <iostream.h>
00024 #include <fstream.h>
00025 #include <iomanip.h>
00026 #include <math.h>
00027 #include <stdlib.h>
00028 
00030 #include <gsl/gsl_math.h>
00031 #include <gsl/gsl_matrix#typeext#.h>
00032 #include <gsl/gsl_linalg.h>
00033 #include <gslwrap/permutation.h>
00034 #include <gslwrap/vector_#type#.h>
00035 
00036 #define type_is#typeext#
00037 #ifdef  type_is
00038 #define type_is_double
00039 #endif
00040 
00041 namespace gsl
00042 {
00043 
00045 class matrix#typeext#
00046 {
00047 #ifdef type_is_double
00048         friend class matrix_float;
00049         friend class matrix_int;
00050 #endif
00051 public:
00053         matrix#typeext#();
00055         matrix#typeext#( size_t new_rows, size_t new_cols, bool clear = true );
00056 
00057         void copy(const matrix#typeext# &other)
00058         {
00059                 set_dimensions(other.size1(),other.size2());
00060                 gsl_matrix#typeext#_memcpy( m, other.m );
00061         }
00062         template<class oclass>
00063         void copy(const oclass &other)
00064         {
00065                 set_dimensions( other.get_rows(), other.get_cols() );
00066                 for ( size_t i = 0; i < get_rows(); i++ ) 
00067                 {
00068                         for ( size_t j = 0; j < get_cols(); j++ ) 
00069                         {
00070                                 gsl_matrix#typeext#_set( m, i, j, (#type#)other(i,j));
00071                         }
00072                 }
00073         }
00074         matrix#typeext#( const matrix#typeext# &other ):m(NULL) {copy(other);}
00076         template<class oclass>
00077         matrix#typeext#( const oclass &other ):m(NULL) {copy(other);}
00078 
00080         ~matrix#typeext#();
00082 //      matrix#typeext#( const char *Filename );
00084         size_t get_rows() const {return m->size1;}
00086         size_t get_cols() const {return m->size2;}
00088         size_t size1() const {return m->size1;}
00090         size_t size2() const {return m->size2;}
00091    
00092 
00094         void dimensions( size_t *num_rows, size_t *num_cols ) const;
00096 #type#        get_element ( size_t row, size_t col ) const {return  gsl_matrix#typeext#_get( m, row, col ) ;}
00097         const #type# &operator()( size_t row, size_t col ) const {return *gsl_matrix#typeext#_ptr( m, row, col ) ;}
00098 #type#       &operator()( size_t row, size_t col )       {return *gsl_matrix#typeext#_ptr( m, row, col ) ;}
00099 
00100         void set_element( size_t row, size_t col, const #type# &v ){ gsl_matrix#typeext#_set( m, row, col, v );}
00102         void set_elements( const #type# & new_value );
00103         void set_all ( const #type# & new_value ) {gsl_matrix#typeext#_set_all ( m, new_value );}
00104         void set_zero() {gsl_matrix#typeext#_set_zero( m );}
00106         void set_dimensions( size_t new_rows, size_t new_cols );
00108         void load( const char *filename );
00110         void save( const char *filename ) const;
00112         friend ostream& operator<< ( ostream& os, const matrix#typeext#& m );
00113         //This function writes the elements of the matrix m to the stream stream in binary format. The return value is 0 for success and GSL_EFAILED if there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.
00114         int fwrite (FILE * stream) const {return gsl_matrix#typeext#_fwrite (stream, m);}
00115 
00116 //This function reads into the matrix m from the open stream stream in binary format. The matrix m must be preallocated with the correct dimensions since the function uses the size of m to determine how many bytes to read. The return value is 0 for success and GSL_EFAILED if there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture. 
00117         int fread (FILE * stream) {return gsl_matrix#typeext#_fread (stream, m);}
00118 
00120         void load_binary( const char *filename );
00122         void save_binary( const char *filename ) const;
00124         bool operator==( const matrix#typeext# &other ) const;
00125         bool operator!=( const matrix#typeext# &other ) const {return !((*this)==other);}
00126 
00127         matrix#typeext# &operator=( const matrix#typeext# &other ){copy(other);return *this;}
00129         template<class omatrix>
00130         matrix#typeext# &operator=( const omatrix &other ){copy(other);return *this;}
00132         matrix#typeext# operator+( const matrix#typeext# &other ) const;
00134         matrix#typeext# operator+( const #type# &f ) const;
00136         friend matrix#typeext# operator+( const #type# &f, const matrix#typeext# &other );
00138         matrix#typeext# &operator+=( const #type# &f );
00140         matrix#typeext# &operator+=( const matrix#typeext# &other );
00142         matrix#typeext# operator-( const matrix#typeext# &other ) const;
00144         matrix#typeext# operator-( const #type# &f ) const;
00146         friend matrix#typeext# operator-( const #type# &f, const matrix#typeext# &other );
00148         matrix#typeext# &operator-=( const #type# &f );
00150         matrix#typeext# &operator-=( const matrix#typeext# &other );
00152         matrix#typeext# operator*( const matrix#typeext# &other ) const;
00154         matrix#typeext# operator*( const #type# &f ) const;
00156         friend matrix#typeext# operator*( const #type# &f, const matrix#typeext# &other );
00158         matrix#typeext# &operator*=( const #type# &f );
00160         matrix#typeext# &operator*=( const matrix#typeext# &other );
00162         matrix#typeext# operator/( const #type# &) const;
00164         matrix#typeext# &operator/=( const #type# &);
00166         matrix#typeext# transpose() const;
00168         matrix#typeext# LU_decomp(gsl::permutation *perm=NULL,int *psign=NULL) const;
00170         matrix#typeext# LU_invert() const;
00171 private:
00173         void LU_decomp( gsl_matrix#typeext# **a,
00174                                         gsl_permutation **permutation,
00175                                         int *sign ) const;
00176 public:
00178     #type# sum() const;
00179 
00180         double LU_lndet() const;
00181 
00182 
00184         vector#typeext#_view       row( size_t rowindex );
00185         const vector#typeext#_view row( size_t rowindex ) const ;
00187         vector#typeext#_view       column( size_t colindex );
00188         const vector#typeext#_view column( size_t colindex ) const;
00189 
00191         matrix#typeext# get_row( size_t rowindex ) const;
00193         matrix#typeext# get_col( size_t colindex ) const;
00195         matrix#typeext# row_sum() const;
00197         matrix#typeext# column_sum() const;
00199         double trace() const;
00201         int cholesky_decomp( matrix#typeext# &a ) const;
00202 //      /** returns index of nearest row in matrix to vector argument. */
00203 //      int nearest_row_index( const matrix#typeext# &v ) const;
00205         matrix#typeext# covariance() const;
00207         bool is_square() const;
00209         void set_diagonal( #type# f );
00211         void identity( size_t k );
00213         double norm( double n ) const;
00214 
00216         gsl_matrix#typeext#       *gslobj()       {if (!m){cout << "matrix#typeext#::gslobj ERROR, data not initialized!! " << endl; exit(-1);}return m;}
00217         const gsl_matrix#typeext# *gslobj() const {if (!m){cout << "matrix#typeext#::gslobj ERROR, data not initialized!! " << endl; exit(-1);}return m;}
00218  private:
00220    gsl_matrix#typeext# *m;
00221 
00222 };
00223 }
00224 #undef type_is#typeext#
00225 #undef type_is_double
00226 
00227 #endif // _matrix_#type#_h