00001 00002 #include<gslwrap/matrix_double.h> 00003 #include<gslwrap/matrix_float.h> 00004 #include<gslwrap/vector_float.h> 00005 #include<iomanip.h> 00006 00007 #define type_is_float 00008 #ifdef type_is 00009 #define type_is_double 00010 #endif 00011 00012 namespace gsl 00013 { 00014 00077 matrix_float::matrix_float():m(NULL) 00078 { 00079 } 00080 00081 matrix_float::matrix_float( size_t rows, size_t cols , bool clear) 00082 { 00083 if(clear){m = gsl_matrix_float_calloc( rows, cols );} 00084 else {m = gsl_matrix_float_alloc ( rows, cols );} 00085 } 00086 00087 00088 matrix_float::~matrix_float() 00089 { 00090 if ( m ) {gsl_matrix_float_free( m );} 00091 } 00092 00093 // matrix_float::matrix_float( const char *filename ) 00094 // { 00095 // ; 00096 // } 00097 00098 void matrix_float::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_float::set_elements( const float & new_value ) 00106 { 00107 gsl_matrix_float_set_all( m, new_value ); 00108 } 00109 00110 void matrix_float::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_float_free( m ); 00115 } 00116 00117 // allocate 00118 m = gsl_matrix_float_calloc( new_rows, new_cols ); 00119 } 00120 00121 00122 void matrix_float::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_float::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_float & 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_float::load_binary( const char *filename ) 00166 { 00167 ; 00168 } 00169 00170 void matrix_float::save_binary( const char *filename ) const 00171 { 00172 ; 00173 } 00174 00175 bool matrix_float::operator==( const matrix_float &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_float matrix_float::operator+( const matrix_float &other ) const 00197 { 00198 matrix_float result(*this); 00199 gsl_matrix_float_add( result.m, other.m ); 00200 return result; 00201 } 00202 00203 matrix_float matrix_float::operator+( const float &f ) const 00204 { 00205 matrix_float result( *this ); 00206 gsl_matrix_float_add_constant( result.m, f ); 00207 00208 return( result ); 00209 } 00210 00211 matrix_float operator+( const float &f, const matrix_float &other ) 00212 { 00213 matrix_float result( other ); 00214 gsl_matrix_float_add_constant( result.m, f ); 00215 00216 return( result ); 00217 } 00218 00219 matrix_float &matrix_float::operator+=( const float &f ) 00220 { 00221 gsl_matrix_float_add_constant( m, f ); 00222 00223 return( *this ); 00224 } 00225 00226 matrix_float &matrix_float::operator+=( const matrix_float &other ) 00227 { 00228 gsl_matrix_float_add( m, other.m ); 00229 00230 return( *this ); 00231 } 00232 00233 matrix_float matrix_float::operator-( const matrix_float &other ) const 00234 { 00235 matrix_float result( *this ); 00236 gsl_matrix_float_sub( result.m, other.m ); 00237 00238 return( result ); 00239 } 00240 00241 matrix_float matrix_float::operator-( const float &f ) const 00242 { 00243 matrix_float result( *this ); 00244 gsl_matrix_float_add_constant( result.m, -f ); 00245 00246 return( result ); 00247 } 00248 00249 matrix_float operator-( const float &f, const matrix_float &other ) 00250 { 00251 matrix_float result( -1 * other ); 00252 gsl_matrix_float_add_constant( result.m, f ); 00253 00254 return( result ); 00255 } 00256 00257 matrix_float &matrix_float::operator-=( const float &f ) 00258 { 00259 gsl_matrix_float_add_constant( m, -f ); 00260 00261 return( *this ); 00262 } 00263 00264 matrix_float &matrix_float::operator-=( const matrix_float &other ) 00265 { 00266 gsl_matrix_float_sub( m, other.m ); 00267 00268 return( *this ); 00269 } 00270 00271 00272 matrix_float matrix_float::operator*( const matrix_float &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_float matrix_float::operator*( const float &f ) const 00286 { 00287 matrix_float result( *this ); 00288 gsl_matrix_float_scale( result.m, f ); 00289 return( result ); 00290 } 00291 00292 matrix_float operator*( const float &f, const matrix_float &other ) 00293 { 00294 matrix_float result( other ); 00295 gsl_matrix_float_scale( result.m, f ); 00296 00297 return( result ); 00298 } 00299 00300 matrix_float &matrix_float::operator*=( const float &f ) 00301 { 00302 gsl_matrix_float_scale( m, f ); 00303 00304 return( *this ); 00305 } 00306 00307 matrix_float &matrix_float::operator*=( const matrix_float &other ) 00308 { 00309 *this=(*this)*other; 00310 return( *this ); 00311 } 00312 00313 matrix_float matrix_float::operator/( const float &f ) const 00314 { 00315 matrix_float result( *this ); 00316 00317 if ( f != 0.0 ) { 00318 gsl_matrix_float_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_float &matrix_float::operator/=( const float &f ) 00328 { 00329 if ( f != 0.0 ) { 00330 gsl_matrix_float_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_float matrix_float::transpose() const 00340 { 00341 int i, j; 00342 matrix_float 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_float_set( result.m, j, i, gsl_matrix_float_get( m, i, j ) ); 00347 } 00348 } 00349 00350 return( result ); 00351 } 00352 00353 00354 matrix_float matrix_float::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_float matrix_float::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 float matrix_float::sum() const 00383 { 00384 int i, j; 00385 float sum = 0; 00386 00387 for ( i = 0; i < get_rows(); i++ ) { 00388 for ( j = 0; j < get_cols(); j++ ) { 00389 sum += gsl_matrix_float_get( m, i, j ); 00390 } 00391 } 00392 00393 return( sum ); 00394 } 00395 00396 // returns logarithm of the determinant of the matrix. 00397 double matrix_float::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_float_view 00407 matrix_float::row( size_t rowindex ) 00408 { 00409 gsl_vector_float_view view=gsl_matrix_float_row(m, rowindex); 00410 return vector_float_view::create_vector_view(view); 00411 } 00412 00413 const 00414 vector_float_view 00415 matrix_float::row( size_t rowindex ) const 00416 { 00417 gsl_vector_float_view view=gsl_matrix_float_row(m, rowindex); 00418 return vector_float_view::create_vector_view(view); 00419 } 00420 00421 vector_float_view 00422 matrix_float::column( size_t colindex ) 00423 { 00424 gsl_vector_float_view view=gsl_matrix_float_column(m, colindex); 00425 return vector_float_view::create_vector_view(view); 00426 } 00427 00428 const 00429 vector_float_view 00430 matrix_float::column( size_t colindex ) const 00431 { 00432 gsl_vector_float_view view=gsl_matrix_float_column(m, colindex); 00433 return vector_float_view::create_vector_view(view); 00434 } 00435 00437 matrix_float matrix_float::get_row( size_t rowindex ) const 00438 { 00439 matrix_float rowmatrix( 1, get_cols() ); 00440 gsl_vector_float *tempvector = gsl_vector_float_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_float_get_row( tempvector, m, rowindex ); 00449 gsl_matrix_float_set_row( rowmatrix.m, 0, tempvector ); 00450 00451 // tidy up 00452 gsl_vector_float_free( tempvector ); 00453 00454 return( rowmatrix ); 00455 } 00456 00458 matrix_float matrix_float::get_col( size_t colindex ) const 00459 { 00460 matrix_float columnmatrix( get_rows(), 1 ); 00461 gsl_vector_float *tempvector = gsl_vector_float_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_float_get_col( tempvector, m, colindex ); 00470 gsl_matrix_float_set_col( columnmatrix.m, 0, tempvector ); 00471 for ( int i = 0; i < get_rows(); i++ ) 00472 cout << gsl_vector_float_get( tempvector, i ) << endl; 00473 00474 // tidy up 00475 gsl_vector_float_free( tempvector ); 00476 00477 return( columnmatrix ); 00478 } 00479 00481 matrix_float matrix_float::row_sum() const 00482 { 00483 int i; 00484 matrix_float 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_float matrix_float::column_sum() const 00496 { 00497 int i; 00498 matrix_float 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_float::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_float::cholesky_decomp( matrix_float &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_float::identity( size_t k ) 00543 { 00544 set_dimensions( k, k ); 00545 set_zero(); 00546 set_diagonal( 1 ); 00547 } 00548 00550 void matrix_float::set_diagonal( float 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_float::is_square() const 00560 { 00561 if ( get_rows() == get_cols() ) 00562 return true; 00563 return false; 00564 } 00565 00567 double matrix_float::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 }