lib/maracasVisuLib/src/kernel/include/matrix.cxx

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