lib/maracasVisuLib/src/kernel/matrix.cxx

   1 /*=========================================================================
   2
   3   Program:   wxMaracas
   4   Module:    $RCSfile: matrix.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2008/10/31 16:32:56 $
   7   Version:   $Revision: 1.1 $
   8
   9   Copyright: (c) 2002, 2003
  10   License:
  11
  12      This software is distributed WITHOUT ANY WARRANTY; without even
  13      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  14      PURPOSE.  See the above copyright notice for more information.
  15
  16 =========================================================================*/
  17
  18 #include "gslobj.hxx"
  19 #include <string>
  20 #include <gsl/gsl_blas.h>
  21 #include <gsl/gsl_linalg.h>
  22
  23 // ---------------------------------------------------------------------
  24 std::ostream& operator<<( std::ostream& os, const gslobj_matrix& m )
  25 {
  26     for( int i = 0; i < m.size1( ); i++ ) {
  27
  28         for( int j = 0; j < m.size2( ); j++ )
  29             os << m( i, j ) << " ";
  30         os << std::endl;
  31
  32     } // rof
  33     return( os );
  34 }
  35
  36 // ---------------------------------------------------------------------
  37 gslobj_matrix::gslobj_matrix( size_t r, size_t c )
  38     : _copy( true )
  39 {
  40     _gsl = gsl_matrix_alloc( r, c );
  41 }
  42
  43 // ---------------------------------------------------------------------
  44 gslobj_matrix::gslobj_matrix( const gslobj_matrix& m )
  45     : _copy( true )
  46 {
  47     _gsl = gsl_matrix_alloc( m._gsl->size1, m._gsl->size2 );
  48     gsl_matrix_memcpy( _gsl, m._gsl );
  49 }
  50
  51 // ---------------------------------------------------------------------
  52 gslobj_matrix::gslobj_matrix( double* m, size_t r, size_t c )
  53     : _copy( false )
  54 {
  55     _view = gsl_matrix_view_array( m, r, c );
  56     _gsl = &( _view.matrix );
  57 }
  58
  59 // ---------------------------------------------------------------------
  60 gslobj_matrix::gslobj_matrix( double* m[], size_t r, size_t c )
  61     : _copy( false )
  62 {
  63     _view = gsl_matrix_view_array( m[ 0 ], r, c );
  64     _gsl = &( _view.matrix );
  65 }
  66
  67 // ---------------------------------------------------------------------
  68 gslobj_matrix::gslobj_matrix( gsl_matrix* m )
  69     : _copy( false )
  70 {
  71     _gsl = m;
  72 }
  73
  74 // ---------------------------------------------------------------------
  75 /*==PS========================================
  76 void gslobj_matrix::resize( size_t r, size_t c )
  77 {
  78     if( _copy ) {
  79
  80         if( r != rows( ) || c != columns( ) )
  81             gsl_matrix_free( _gsl );
  82         _gsl = gsl_matrix_alloc( r, c );
  83
  84     } // fi
  85 }
  86 ==PS========================================*/
  87 // ---------------------------------------------------------------------
  88 gslobj_matrix& gslobj_matrix::operator=( const gslobj_matrix& o )
  89 {
  90     gsl_matrix_memcpy( _gsl, o._gsl );
  91     return( *this );
  92 }
  93
  94 // ---------------------------------------------------------------------
  95 gslobj_matrix& gslobj_matrix::operator=( double o )
  96 {
  97     gsl_matrix_set_all( _gsl, o );
  98     return( *this );
  99 }
 100
 101 // ---------------------------------------------------------------------
 102 gslobj_matrix& gslobj_matrix::operator=( double* o )
 103 {
 104     memcpy( _gsl->data, o,
 105             _gsl->size1 * _gsl->size2 * sizeof( double ) );
 106     return( *this );
 107 }
 108
 109 // ---------------------------------------------------------------------
 110 gslobj_matrix& gslobj_matrix::operator=( double* o[] )
 111 {
 112     memcpy( _gsl->data, o[ 0 ],
 113             _gsl->size1 * _gsl->size2 * sizeof( double ) );
 114     return( *this );
 115 }
 116
 117 // ---------------------------------------------------------------------
 118 gslobj_matrix& gslobj_matrix::operator=( gsl_matrix* o )
 119 {
 120     gsl_matrix_memcpy( _gsl, o );
 121     return( *this );
 122 }
 123
 124 // ---------------------------------------------------------------------
 125 double& gslobj_matrix::operator()( size_t i, size_t j )
 126 {
 127     return( _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] );
 128 }
 129
 130 // ---------------------------------------------------------------------
 131 const double& gslobj_matrix::operator()( size_t i, size_t j ) const
 132 {
 133     return( _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] );
 134 }
 135
 136 // ---------------------------------------------------------------------
 137 /*==PS========================================
 138 void gslobj_matrix::setValue( size_t i, size_t j, double val )
 139 {
 140     _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] = val;
 141 }
 142 ==PS========================================*/
 143
 144 // ---------------------------------------------------------------------
 145 bool gslobj_matrix::operator==( const gslobj_matrix& o ) const
 146 {
 147     bool ret = true;
 148
 149     if( _gsl->size1 != o._gsl->size1 || _gsl->size2 != o._gsl->size2 )
 150         return( false );
 151
 152     for( int i = 0; i < _gsl->size1 && ret; i++ )
 153         for( int j = 0; j < _gsl->size2 && ret; j++ )
 154             ret = ret &&
 155                 ( _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] !=
 156                   o._gsl->data[ ( ( i * _gsl->size2 ) + j ) ] );
 157     return( ret );
 158 }
 159
 160 // ---------------------------------------------------------------------
 161 gslobj_matrix gslobj_matrix::operator+( const gslobj_matrix& o )
 162 {
 163     gslobj_matrix r( *this );
 164     gsl_matrix_add( r._gsl, o._gsl );
 165     return( r );
 166 }
 167
 168 // ---------------------------------------------------------------------
 169 gslobj_matrix gslobj_matrix::operator+( double o )
 170 {
 171     gslobj_matrix r( *this );
 172     gsl_matrix_add_constant( r._gsl, o );
 173     return( r );
 174 }
 175
 176 // ---------------------------------------------------------------------
 177 gslobj_matrix gslobj_matrix::operator+( double* o )
 178 {
 179     gslobj_matrix r( *this );
 180     gsl_matrix_view vw = gsl_matrix_view_array( o,
 181                                                 _gsl->size1,
 182                                                 _gsl->size2 );
 183     gsl_matrix_add( r._gsl, &( vw.matrix ) );
 184     return( r );
 185 }
 186
 187 // ---------------------------------------------------------------------
 188 gslobj_matrix gslobj_matrix::operator+( double* o[] )
 189 {
 190     gslobj_matrix r( *this );
 191     gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
 192                                                 _gsl->size1,
 193                                                 _gsl->size2 );
 194     gsl_matrix_add( r._gsl, &( vw.matrix ) );
 195     return( r );
 196 }
 197
 198 // ---------------------------------------------------------------------
 199 gslobj_matrix gslobj_matrix::operator+( gsl_matrix* o )
 200 {
 201     gslobj_matrix r( *this );
 202     gsl_matrix_add( r._gsl, o );
 203     return( r );
 204 }
 205
 206 // ---------------------------------------------------------------------
 207 gslobj_matrix& gslobj_matrix::operator+=( const gslobj_matrix& o )
 208 {
 209     gsl_matrix_add( _gsl, o._gsl );
 210     return( *this );
 211 }
 212
 213 // ---------------------------------------------------------------------
 214 gslobj_matrix& gslobj_matrix::operator+=( double o )
 215 {
 216     gsl_matrix_add_constant( _gsl, o );
 217     return( *this );
 218 }
 219
 220 // ---------------------------------------------------------------------
 221 gslobj_matrix& gslobj_matrix::operator+=( double* o )
 222 {
 223     gsl_matrix_view vw = gsl_matrix_view_array( o,
 224                                                 _gsl->size1,
 225                                                 _gsl->size2 );
 226     gsl_matrix_add( _gsl, &( vw.matrix ) );
 227     return( *this );
 228 }
 229
 230 // ---------------------------------------------------------------------
 231 gslobj_matrix& gslobj_matrix::operator+=( double* o[] )
 232 {
 233     gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
 234                                                 _gsl->size1,
 235                                                 _gsl->size2 );
 236     gsl_matrix_add( _gsl, &( vw.matrix ) );
 237     return( *this );
 238 }
 239
 240 // ---------------------------------------------------------------------
 241 gslobj_matrix& gslobj_matrix::operator+=( gsl_matrix* o )
 242 {
 243     gsl_matrix_add( _gsl, o );
 244     return( *this );
 245 }
 246
 247 // ---------------------------------------------------------------------
 248 gslobj_matrix gslobj_matrix::operator-( const gslobj_matrix& o )
 249 {
 250     gslobj_matrix r( *this );
 251     gsl_matrix_sub( r._gsl, o._gsl );
 252     return( r );
 253 }
 254
 255 // ---------------------------------------------------------------------
 256 gslobj_matrix gslobj_matrix::operator-( double o )
 257 {
 258     gslobj_matrix r( *this );
 259     gsl_matrix_add_constant( r._gsl, -o );
 260     return( r );
 261 }
 262
 263 // ---------------------------------------------------------------------
 264 gslobj_matrix gslobj_matrix::operator-( double* o )
 265 {
 266     gslobj_matrix r( *this );
 267     gsl_matrix_view vw = gsl_matrix_view_array( o,
 268                                                 _gsl->size1,
 269                                                 _gsl->size2 );
 270     gsl_matrix_sub( r._gsl, &( vw.matrix ) );
 271     return( r );
 272 }
 273
 274 // ---------------------------------------------------------------------
 275 gslobj_matrix gslobj_matrix::operator-( double* o[] )
 276 {
 277     gslobj_matrix r( *this );
 278     gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
 279                                                 _gsl->size1,
 280                                                 _gsl->size2 );
 281     gsl_matrix_sub( r._gsl, &( vw.matrix ) );
 282     return( r );
 283 }
 284
 285 // ---------------------------------------------------------------------
 286 gslobj_matrix gslobj_matrix::operator-( gsl_matrix* o )
 287 {
 288     gslobj_matrix r( *this );
 289     gsl_matrix_sub( r._gsl, o );
 290     return( r );
 291 }
 292
 293 // ---------------------------------------------------------------------
 294 gslobj_matrix& gslobj_matrix::operator-=( const gslobj_matrix& o )
 295 {
 296     gsl_matrix_sub( _gsl, o._gsl );
 297     return( *this );
 298 }
 299
 300 // ---------------------------------------------------------------------
 301 gslobj_matrix& gslobj_matrix::operator-=( double o )
 302 {
 303     gsl_matrix_add_constant( _gsl, -o );
 304     return( *this );
 305 }
 306
 307 // ---------------------------------------------------------------------
 308 gslobj_matrix& gslobj_matrix::operator-=( double* o )
 309 {
 310     gsl_matrix_view vw = gsl_matrix_view_array( o,
 311                                                 _gsl->size1,
 312                                                 _gsl->size2 );
 313     gsl_matrix_sub( _gsl, &( vw.matrix ) );
 314     return( *this );
 315 }
 316
 317 // ---------------------------------------------------------------------
 318 gslobj_matrix& gslobj_matrix::operator-=( double* o[] )
 319 {
 320     gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
 321                                                 _gsl->size1,
 322                                                 _gsl->size2 );
 323     gsl_matrix_sub( _gsl, &( vw.matrix ) );
 324     return( *this );
 325 }
 326
 327 // ---------------------------------------------------------------------
 328 gslobj_matrix& gslobj_matrix::operator-=( gsl_matrix* o )
 329 {
 330     gsl_matrix_sub( _gsl, o );
 331     return( *this );
 332 }
 333
 334 // ---------------------------------------------------------------------
 335 gslobj_matrix gslobj_matrix::operator*( const gslobj_matrix& o )
 336 {
 337     gslobj_matrix r( *this );
 338     gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
 339                     1.0, r._gsl, o._gsl,
 340                     0.0, _gsl );
 341     return( r );
 342 }
 343
 344 // ---------------------------------------------------------------------
 345 gslobj_matrix gslobj_matrix::operator*( double o )
 346 {
 347     gslobj_matrix r( *this );
 348     gsl_matrix_scale( r._gsl, o );
 349     return( r );
 350 }
 351
 352 // ---------------------------------------------------------------------
 353 gslobj_matrix gslobj_matrix::operator*( double* o )
 354 {
 355     gslobj_matrix r( *this );
 356     gsl_matrix_view vw = gsl_matrix_view_array( o,
 357                                                 _gsl->size1,
 358                                                 _gsl->size2 );
 359     gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
 360                     1.0, r._gsl, &( vw.matrix ),
 361                     0.0, _gsl );
 362     return( r );
 363 }
 364
 365 // ---------------------------------------------------------------------
 366 gslobj_matrix gslobj_matrix::operator*( double* o[] )
 367 {
 368     gslobj_matrix r( *this );
 369     gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
 370                                                 _gsl->size1,
 371                                                 _gsl->size2 );
 372     gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
 373                     1.0, r._gsl, &( vw.matrix ),
 374                     0.0, _gsl );
 375     return( r );
 376 }
 377
 378 // ---------------------------------------------------------------------
 379 gslobj_matrix gslobj_matrix::operator*( gsl_matrix* o )
 380 {
 381     gslobj_matrix r( *this );
 382     gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
 383                     1.0, r._gsl, o,
 384                     0.0, _gsl );
 385     return( r );
 386 }
 387
 388 // ---------------------------------------------------------------------
 389 gslobj_vector gslobj_matrix::operator*( const gslobj_vector& o )
 390 {
 391     gslobj_vector ret( rows( ) );
 392     gsl_blas_dgemv( CblasNoTrans, 1.0, _gsl,
 393                     ( gsl_vector* )o, 0.0,
 394                     ( gsl_vector* )ret );
 395     return( ret );
 396 }
 397
 398 // ---------------------------------------------------------------------
 399 gslobj_matrix& gslobj_matrix::operator*=( const gslobj_matrix& o )
 400 {
 401     gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
 402                     1.0, _gsl, o._gsl,
 403                     0.0, _gsl );
 404     return( *this );
 405 }
 406
 407 // ---------------------------------------------------------------------
 408 gslobj_matrix& gslobj_matrix::operator*=( double o )
 409 {
 410     gsl_matrix_scale( _gsl, o );
 411     return( *this );
 412 }
 413
 414 // ---------------------------------------------------------------------
 415 gslobj_matrix& gslobj_matrix::operator*=( double* o )
 416 {
 417     gsl_matrix_view vw = gsl_matrix_view_array( o,
 418                                                 _gsl->size1,
 419                                                 _gsl->size2 );
 420     gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
 421                     1.0, _gsl, &( vw.matrix ),
 422                     0.0, _gsl );
 423     return( *this );
 424 }
 425
 426 // ---------------------------------------------------------------------
 427 gslobj_matrix& gslobj_matrix::operator*=( double* o[] )
 428 {
 429     gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
 430                                                 _gsl->size1,
 431                                                 _gsl->size2 );
 432     gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
 433                     1.0, _gsl, &( vw.matrix ),
 434                     0.0, _gsl );
 435     return( *this );
 436 }
 437
 438 // ---------------------------------------------------------------------
 439 gslobj_matrix& gslobj_matrix::operator*=( gsl_matrix* o )
 440 {
 441     gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
 442                     1.0, _gsl, o,
 443                     0.0, _gsl );
 444     return( *this );
 445 }
 446
 447 // ---------------------------------------------------------------------
 448 /*==PS========================================
 449 gslobj_matrix gslobj_matrix::transpose( )
 450 {
 451     gslobj_matrix r( _gsl->size2, _gsl->size1 );
 452     gsl_matrix_transpose_memcpy( r._gsl, _gsl );
 453     return( r );
 454 }
 455 ==PS========================================*/
 456 // ---------------------------------------------------------------------
 457 /*==PS========================================
 458 gslobj_matrix gslobj_matrix::invertLU( )
 459 {
 460     gslobj_matrix r( *this ), lu( *this );
 461     gsl_permutation* perm = gsl_permutation_calloc( rows( ) );
 462     int sign;
 463
 464     gsl_linalg_LU_decomp( lu._gsl, perm, &sign );
 465     gsl_linalg_LU_invert( lu._gsl, perm, r._gsl );
 466     return( r );
 467 }
 468 ==PS========================================*/
 469 // ---------------------------------------------------------------------
 470 /*==PS========================================
 471 double gslobj_matrix::determinant( )
 472 {
 473     gslobj_matrix lu( *this );
 474     gsl_permutation* perm = gsl_permutation_calloc( rows( ) );
 475     int sign;
 476
 477     gsl_linalg_LU_decomp( lu._gsl, perm, &sign );
 478     return( gsl_linalg_LU_det( lu._gsl, sign ) );
 479 }
 480 ==PS========================================*/
 481 // ---------------------------------------------------------------------
 482 /*==PS========================================
 483 void gslobj_matrix::identity( )
 484 {
 485     gsl_matrix_set_identity( _gsl );
 486 }
 487 ==PS========================================*/
 488 // eof - gslobj_matrix.cxx