Main Page   Namespace List   Class Hierarchy   Compound List   File List   Namespace Members   Compound Members   File Members  

matrix_double.cc

Go to the documentation of this file.
00001 
00002 #include<gslwrap/matrix_double.h>
00003 #include<gslwrap/matrix_double.h>
00004 #include<gslwrap/vector_double.h>
00005 #include<iomanip.h>
00006 
00007 #define type_is
00008 #ifdef type_is
00009 #define type_is_double
00010 #endif
00011 
00012 namespace gsl
00013 {
00014 
00077 matrix::matrix():m(NULL)
00078 {
00079 }
00080 
00081 matrix::matrix( size_t rows, size_t cols , bool clear)
00082 {
00083    if(clear){m = gsl_matrix_calloc( rows, cols );}
00084    else     {m = gsl_matrix_alloc ( rows, cols );}
00085 }
00086 
00087 
00088 matrix::~matrix()
00089 {
00090    if ( m ) {gsl_matrix_free( m );}
00091 }
00092 
00093 //  matrix::matrix( const char *filename )
00094 //  {
00095 //     ;
00096 //  }
00097 
00098 void matrix::dimensions( size_t *num_rows, size_t *num_cols ) const
00099 {
00100    *num_rows = m->size1;
00101    *num_cols = m->size2;
00102 }
00103 
00104 
00105 void matrix::set_elements( const double & new_value  )
00106 {
00107    gsl_matrix_set_all( m, new_value );
00108 }
00109 
00110 void matrix::set_dimensions( size_t new_rows, size_t new_cols )
00111 {
00112    // if dimensions have changed re-allocate matrix
00113    if ( m && (get_rows() != new_rows || get_cols() != new_cols )) {
00114               gsl_matrix_free( m );
00115    }
00116 
00117    // allocate
00118    m = gsl_matrix_calloc( new_rows, new_cols );
00119 }
00120 
00121 
00122 void matrix::load( const char *filename )
00123 {
00124    FILE * f = fopen( filename, "r" ) ;
00125    matrix temp;
00126    gsl_matrix_fread ( f, temp.m );
00127    fclose (f);
00128    *this = temp;
00129 }
00130 
00131 void matrix::save( const char *filename ) const
00132 {
00133    FILE * f = fopen( filename, "w" ) ;
00134    matrix temp = *this;
00135    gsl_matrix_fwrite ( f, temp.m );
00136    fclose ( f );
00137 }
00138 
00139 ostream& operator<< ( ostream& os, const matrix & m )
00140 {
00141    size_t i, j;
00142 
00143    os.setf( ios::fixed );
00144 
00145 //FIXME for aCC (doesn't find correct outstream function
00146 //     for ( i = 0; i < m.get_rows(); i++ ) {
00147 //         for ( j = 0; j < m.get_cols() - 1; j++ ) {
00148 //                 os << setprecision( 6 ) << setw( 11 ) ;//<< m.get_element( i, j ) << " ";
00149 //         }
00150 //         os << setprecision( 6 ) << setw( 11 ) ;//<< m.get_element( i, j ) << endl;
00151 //     }
00152 
00153    for ( i = 0; i < m.get_rows(); i++ ) {
00154            for ( j = 0; j < m.get_cols() - 1; j++ ) {
00155                    os << m.get_element( i, j ) << " ";
00156            }
00157            os << m.get_element( i, j ) << endl;
00158    }
00159 
00160    return os;
00161 }
00162 
00163 
00164 
00165 void matrix::load_binary( const char *filename )
00166 {
00167    ;
00168 }
00169 
00170 void matrix::save_binary( const char *filename ) const
00171 {
00172    ;
00173 }
00174 
00175 bool matrix::operator==( const matrix &other ) const
00176 {
00177    size_t i, j;
00178    
00179    // first check for same dimensions
00180    if ( size1() != other.size1() || size2() != other.size2() )
00181    {
00182       return false;
00183    }
00184 
00185    for ( i = 0; i < size1(); i++ ) {
00186       for ( j = 0; j < size2(); j++ ) {
00187          if ( this->get_element( i, j ) != other.get_element( i, j ) ) {
00188             return false;
00189          }
00190       }
00191    }
00192 
00193    return true;
00194 }
00195 
00196 matrix matrix::operator+( const matrix &other ) const
00197 {
00198         matrix result(*this);
00199         gsl_matrix_add( result.m, other.m );
00200         return result;
00201 }
00202 
00203 matrix matrix::operator+( const double &f ) const
00204 {
00205    matrix result( *this );
00206    gsl_matrix_add_constant( result.m, f );
00207 
00208    return( result );
00209 }
00210 
00211 matrix operator+( const double &f, const matrix &other )
00212 {
00213    matrix result( other );
00214    gsl_matrix_add_constant( result.m, f );
00215 
00216    return( result );
00217 }
00218 
00219 matrix &matrix::operator+=( const double &f )
00220 {
00221    gsl_matrix_add_constant( m, f );
00222 
00223    return( *this );
00224 }
00225 
00226 matrix &matrix::operator+=( const matrix &other )
00227 {
00228    gsl_matrix_add( m, other.m );
00229 
00230    return( *this );
00231 }
00232 
00233 matrix matrix::operator-( const matrix &other ) const
00234 {
00235    matrix result( *this );
00236    gsl_matrix_sub( result.m, other.m );
00237 
00238    return( result );
00239 }
00240 
00241 matrix matrix::operator-( const double &f ) const
00242 {
00243    matrix result( *this );
00244    gsl_matrix_add_constant( result.m, -f );
00245 
00246    return( result );
00247 }
00248 
00249 matrix operator-( const double &f, const matrix &other )
00250 {
00251    matrix result( -1 * other );
00252    gsl_matrix_add_constant( result.m, f );
00253 
00254    return( result );
00255 }
00256 
00257 matrix &matrix::operator-=( const double &f )
00258 {
00259    gsl_matrix_add_constant( m, -f );
00260 
00261    return( *this );
00262 }
00263 
00264 matrix &matrix::operator-=( const matrix &other )
00265 {
00266    gsl_matrix_sub( m, other.m );
00267 
00268    return( *this );
00269 }
00270 
00271 
00272 matrix matrix::operator*( const matrix &other ) const
00273 {
00274         matrix result( get_rows(), other.get_cols() );
00275 #ifdef type_is_double
00276         gsl_linalg_matmult(m,other.m,result.m);
00277 #else //type_is_double
00278         matrix a=*this;
00279         matrix b=other;
00280         gsl_linalg_matmult(a.m,b.m,result.m);
00281 #endif //type_is_double
00282         return result ;
00283 }
00284 
00285 matrix matrix::operator*( const double &f ) const
00286 {
00287    matrix result( *this );
00288    gsl_matrix_scale( result.m, f );
00289    return( result );
00290 }
00291 
00292 matrix operator*( const double &f, const matrix &other )
00293 {
00294    matrix result( other );
00295    gsl_matrix_scale( result.m, f );
00296 
00297    return( result );
00298 }
00299 
00300 matrix &matrix::operator*=( const double &f )
00301 {
00302    gsl_matrix_scale( m, f );
00303 
00304    return( *this );
00305 }
00306 
00307 matrix &matrix::operator*=( const matrix &other )
00308 {
00309    *this=(*this)*other;
00310    return( *this );
00311 }
00312 
00313 matrix matrix::operator/( const double &f ) const
00314 {
00315    matrix result( *this );
00316 
00317    if ( f != 0.0 ) {
00318       gsl_matrix_scale( result.m, 1.0 / f );
00319    } else {
00320       cout << "e_r_r_o_r: division by zero." << endl;
00321       return( result );
00322    }
00323 
00324    return( result );
00325 }
00326 
00327 matrix &matrix::operator/=( const double &f )
00328 {
00329    if ( f != 0.0 ) {
00330       gsl_matrix_scale( m, 1.0 / f );
00331    } else {
00332       cout << "e_rr_or: division by zero." << endl;
00333       return( *this );
00334    }
00335 
00336    return( *this );
00337 }
00338 
00339 matrix matrix::transpose() const 
00340 {
00341    int i, j;
00342    matrix result( get_cols(), get_rows() );
00343 
00344    for ( i = 0; i < get_rows(); i++ ) {
00345       for ( j = 0; j < get_cols(); j++ ) {
00346          gsl_matrix_set( result.m, j, i, gsl_matrix_get( m, i, j ) );
00347       }
00348    }
00349    
00350    return( result );
00351 }
00352 
00353 
00354 matrix matrix::LU_decomp(gsl::permutation *perm,int *psign) const
00355 {
00356         bool retPerm=perm!=NULL;
00357         if(!perm){perm = new permutation();}
00358         int sign;
00359         perm->resize( get_rows() );
00360         matrix result=*this;// this does conversion  if necessary
00361         gsl_linalg_LU_decomp(  result.m, perm->gsldata, &sign );
00362 
00363         if(!retPerm){delete perm;}
00364         if(psign){*psign=sign;}
00365 
00366    return result;// this does conversion  if necessary
00367 }
00368 
00369 matrix matrix::LU_invert() const
00370 {
00371    int i, j;
00372    permutation p;
00373    matrix a=*this;
00374    a=a.LU_decomp(&p);
00375    matrix inverse(size1(),size2());
00376    gsl_linalg_LU_invert( a.m, p.gsldata, inverse.m );
00377    return inverse;
00378 }
00379 
00380 
00381 // returns sum of all the elements.
00382 double matrix::sum() const 
00383 {
00384         int i, j;
00385         double sum = 0;
00386 
00387         for ( i = 0; i < get_rows(); i++ ) {
00388                 for ( j = 0; j < get_cols(); j++ ) {
00389                         sum += gsl_matrix_get( m, i, j );
00390                 }
00391         }
00392 
00393         return( sum );
00394 }
00395 
00396 // returns logarithm of the determinant of the matrix.
00397 double matrix::LU_lndet() const 
00398 {
00399         matrix a=*this;
00400         a=a.LU_decomp();
00401 
00402         // calculate log determinant from LU decomposed matrix "a"
00403         return gsl_linalg_LU_lndet( a.m );
00404 }
00405 
00406 vector_view 
00407 matrix::row( size_t rowindex )
00408 {
00409         gsl_vector_view view=gsl_matrix_row(m, rowindex);
00410         return vector_view::create_vector_view(view);
00411 }
00412 
00413 const 
00414 vector_view 
00415 matrix::row( size_t rowindex ) const
00416 {
00417         gsl_vector_view view=gsl_matrix_row(m, rowindex);
00418         return vector_view::create_vector_view(view);
00419 }
00420 
00421 vector_view 
00422 matrix::column( size_t colindex )
00423 {
00424         gsl_vector_view view=gsl_matrix_column(m, colindex);
00425         return vector_view::create_vector_view(view);
00426 }
00427 
00428 const 
00429 vector_view 
00430 matrix::column( size_t colindex ) const
00431 {
00432         gsl_vector_view view=gsl_matrix_column(m, colindex);
00433         return vector_view::create_vector_view(view);
00434 }
00435 
00437 matrix matrix::get_row( size_t rowindex ) const 
00438 {
00439         matrix rowmatrix( 1, get_cols() );
00440         gsl_vector *tempvector = gsl_vector_calloc( get_cols() );
00441         
00442         if ( rowindex < 0 || rowindex >= get_rows() )
00443         {
00444                 cerr << "row index must be in range 0 to " << get_rows() - 1 << endl;
00445                 exit( 1 );
00446         }
00447 
00448         gsl_matrix_get_row( tempvector, m, rowindex );
00449         gsl_matrix_set_row( rowmatrix.m, 0, tempvector );
00450 
00451         // tidy up
00452         gsl_vector_free( tempvector );
00453         
00454         return( rowmatrix );
00455 }
00456 
00458 matrix matrix::get_col( size_t colindex ) const 
00459 {
00460         matrix columnmatrix( get_rows(), 1 );
00461         gsl_vector *tempvector = gsl_vector_calloc( get_rows() );
00462         
00463         if ( colindex < 0 || colindex >= get_cols() )
00464         {
00465                 cerr << "column index must be in range 0 to " << get_cols() - 1 << endl;
00466                 exit( 1 );
00467         }
00468         
00469         gsl_matrix_get_col( tempvector, m, colindex );
00470         gsl_matrix_set_col( columnmatrix.m, 0, tempvector );
00471         for ( int i = 0; i < get_rows(); i++ )
00472                 cout << gsl_vector_get( tempvector, i ) << endl;
00473 
00474         // tidy up
00475         gsl_vector_free( tempvector );
00476         
00477         return( columnmatrix );
00478 }
00479 
00481 matrix matrix::row_sum() const 
00482 {
00483         int     i;
00484         matrix sum( get_rows(), 1 );
00485         
00486         sum.set_zero();
00487         for ( i = 0; i < get_rows(); i++ ) {
00488                 sum.set_element( i, 0, get_row( i ).sum() );
00489         }
00490         
00491         return( sum );
00492 }
00493 
00495 matrix matrix::column_sum() const 
00496 {
00497         int     i;
00498         matrix sum( 1, get_cols() );
00499         
00500         sum.set_zero( );
00501         for ( i = 0; i < get_cols(); i++ ) {
00502                 sum.set_element( 0, i, get_col( i ).sum() );
00503         }
00504                 
00505         return( sum );
00506 }
00507 
00509 double matrix::trace() const 
00510 {
00511         int     i;
00512         double sum = 0.0;
00513         
00514         if ( get_rows() != get_cols() ) {
00515                 cerr << "e_r_r_o_r: cannot calculate trace of non-square matrix.";
00516                 cerr << endl;
00517                 exit( 1 );
00518         }
00519 
00520         // calculate sum of diagonal elements
00521         for ( i = 0; i < get_rows(); i++ ) {
00522                 sum += get_element( i, i );
00523         }
00524 
00525         return( sum );
00526 }
00527 
00529 // don't forget to de-allocate a, which is allocated in this method.
00530 int matrix::cholesky_decomp( matrix &a ) const 
00531 {
00532         int i, j;
00533         int error;
00534         matrix result=*this;
00535         // do decomposition with call to g_s_l
00536         error = gsl_linalg_cholesky_decomp( result.m );
00537         a=result;
00538         return error;
00539 }
00540 
00542 void matrix::identity( size_t k ) 
00543 {
00544         set_dimensions( k, k );
00545         set_zero();
00546         set_diagonal( 1 );
00547 }
00548 
00550 void matrix::set_diagonal( double f ) 
00551 {
00552         size_t i;
00553         int mn=(get_rows()<get_cols() ? get_rows() : get_cols());
00554         for ( i = 0; i < mn; i++ )
00555                         set_element( i, i, f );
00556 }
00557 
00559 bool matrix::is_square() const 
00560 {
00561         if ( get_rows() == get_cols() )
00562                 return true;
00563         return false;
00564 }
00565 
00567 double matrix::norm( double n ) const 
00568 {
00569         int i, j;
00570         double sum = 0.0;
00571         
00572         for ( i = 0; i < get_rows(); i++ ) {
00573                 for ( j = 0; j < get_cols(); j++ ) {
00574                         sum += pow( get_element( i, j ), n );
00575                 }
00576         }
00577         
00578         return sum;
00579 }
00580 
00581 }

Generated at Sun Dec 16 23:44:43 2001 for gslwrap by doxygen1.2.8.1 written by Dimitri van Heesch, © 1997-2001