1 /*=========================================================================
3 /*# ---------------------------------------------------------------------
5 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
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
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.
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
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 # ------------------------------------------------------------------------ */
29 Module: $RCSfile: matrix.cxx,v $
31 Date: $Date: 2012/11/15 14:15:31 $
32 Version: $Revision: 1.2 $
34 Copyright: (c) 2002, 2003
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.
41 =========================================================================*/
45 #include <gsl/gsl_blas.h>
46 #include <gsl/gsl_linalg.h>
48 // ---------------------------------------------------------------------
49 std::ostream& operator<<( std::ostream& os, const gslobj_matrix& m )
51 for( int i = 0; i < m.size1( ); i++ ) {
53 for( int j = 0; j < m.size2( ); j++ )
54 os << m( i, j ) << " ";
61 // ---------------------------------------------------------------------
62 gslobj_matrix::gslobj_matrix( size_t r, size_t c )
65 _gsl = gsl_matrix_alloc( r, c );
68 // ---------------------------------------------------------------------
69 gslobj_matrix::gslobj_matrix( const gslobj_matrix& m )
72 _gsl = gsl_matrix_alloc( m._gsl->size1, m._gsl->size2 );
73 gsl_matrix_memcpy( _gsl, m._gsl );
76 // ---------------------------------------------------------------------
77 gslobj_matrix::gslobj_matrix( double* m, size_t r, size_t c )
80 _view = gsl_matrix_view_array( m, r, c );
81 _gsl = &( _view.matrix );
84 // ---------------------------------------------------------------------
85 gslobj_matrix::gslobj_matrix( double* m[], size_t r, size_t c )
88 _view = gsl_matrix_view_array( m[ 0 ], r, c );
89 _gsl = &( _view.matrix );
92 // ---------------------------------------------------------------------
93 gslobj_matrix::gslobj_matrix( gsl_matrix* m )
99 // ---------------------------------------------------------------------
100 /*==PS========================================
101 void gslobj_matrix::resize( size_t r, size_t c )
105 if( r != rows( ) || c != columns( ) )
106 gsl_matrix_free( _gsl );
107 _gsl = gsl_matrix_alloc( r, c );
111 ==PS========================================*/
112 // ---------------------------------------------------------------------
113 gslobj_matrix& gslobj_matrix::operator=( const gslobj_matrix& o )
115 gsl_matrix_memcpy( _gsl, o._gsl );
119 // ---------------------------------------------------------------------
120 gslobj_matrix& gslobj_matrix::operator=( double o )
122 gsl_matrix_set_all( _gsl, o );
126 // ---------------------------------------------------------------------
127 gslobj_matrix& gslobj_matrix::operator=( double* o )
129 memcpy( _gsl->data, o,
130 _gsl->size1 * _gsl->size2 * sizeof( double ) );
134 // ---------------------------------------------------------------------
135 gslobj_matrix& gslobj_matrix::operator=( double* o[] )
137 memcpy( _gsl->data, o[ 0 ],
138 _gsl->size1 * _gsl->size2 * sizeof( double ) );
142 // ---------------------------------------------------------------------
143 gslobj_matrix& gslobj_matrix::operator=( gsl_matrix* o )
145 gsl_matrix_memcpy( _gsl, o );
149 // ---------------------------------------------------------------------
150 double& gslobj_matrix::operator()( size_t i, size_t j )
152 return( _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] );
155 // ---------------------------------------------------------------------
156 const double& gslobj_matrix::operator()( size_t i, size_t j ) const
158 return( _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] );
161 // ---------------------------------------------------------------------
162 /*==PS========================================
163 void gslobj_matrix::setValue( size_t i, size_t j, double val )
165 _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] = val;
167 ==PS========================================*/
169 // ---------------------------------------------------------------------
170 bool gslobj_matrix::operator==( const gslobj_matrix& o ) const
174 if( _gsl->size1 != o._gsl->size1 || _gsl->size2 != o._gsl->size2 )
177 for( int i = 0; i < _gsl->size1 && ret; i++ )
178 for( int j = 0; j < _gsl->size2 && ret; j++ )
180 ( _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] !=
181 o._gsl->data[ ( ( i * _gsl->size2 ) + j ) ] );
185 // ---------------------------------------------------------------------
186 gslobj_matrix gslobj_matrix::operator+( const gslobj_matrix& o )
188 gslobj_matrix r( *this );
189 gsl_matrix_add( r._gsl, o._gsl );
193 // ---------------------------------------------------------------------
194 gslobj_matrix gslobj_matrix::operator+( double o )
196 gslobj_matrix r( *this );
197 gsl_matrix_add_constant( r._gsl, o );
201 // ---------------------------------------------------------------------
202 gslobj_matrix gslobj_matrix::operator+( double* o )
204 gslobj_matrix r( *this );
205 gsl_matrix_view vw = gsl_matrix_view_array( o,
208 gsl_matrix_add( r._gsl, &( vw.matrix ) );
212 // ---------------------------------------------------------------------
213 gslobj_matrix gslobj_matrix::operator+( double* o[] )
215 gslobj_matrix r( *this );
216 gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
219 gsl_matrix_add( r._gsl, &( vw.matrix ) );
223 // ---------------------------------------------------------------------
224 gslobj_matrix gslobj_matrix::operator+( gsl_matrix* o )
226 gslobj_matrix r( *this );
227 gsl_matrix_add( r._gsl, o );
231 // ---------------------------------------------------------------------
232 gslobj_matrix& gslobj_matrix::operator+=( const gslobj_matrix& o )
234 gsl_matrix_add( _gsl, o._gsl );
238 // ---------------------------------------------------------------------
239 gslobj_matrix& gslobj_matrix::operator+=( double o )
241 gsl_matrix_add_constant( _gsl, o );
245 // ---------------------------------------------------------------------
246 gslobj_matrix& gslobj_matrix::operator+=( double* o )
248 gsl_matrix_view vw = gsl_matrix_view_array( o,
251 gsl_matrix_add( _gsl, &( vw.matrix ) );
255 // ---------------------------------------------------------------------
256 gslobj_matrix& gslobj_matrix::operator+=( double* o[] )
258 gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
261 gsl_matrix_add( _gsl, &( vw.matrix ) );
265 // ---------------------------------------------------------------------
266 gslobj_matrix& gslobj_matrix::operator+=( gsl_matrix* o )
268 gsl_matrix_add( _gsl, o );
272 // ---------------------------------------------------------------------
273 gslobj_matrix gslobj_matrix::operator-( const gslobj_matrix& o )
275 gslobj_matrix r( *this );
276 gsl_matrix_sub( r._gsl, o._gsl );
280 // ---------------------------------------------------------------------
281 gslobj_matrix gslobj_matrix::operator-( double o )
283 gslobj_matrix r( *this );
284 gsl_matrix_add_constant( r._gsl, -o );
288 // ---------------------------------------------------------------------
289 gslobj_matrix gslobj_matrix::operator-( double* o )
291 gslobj_matrix r( *this );
292 gsl_matrix_view vw = gsl_matrix_view_array( o,
295 gsl_matrix_sub( r._gsl, &( vw.matrix ) );
299 // ---------------------------------------------------------------------
300 gslobj_matrix gslobj_matrix::operator-( double* o[] )
302 gslobj_matrix r( *this );
303 gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
306 gsl_matrix_sub( r._gsl, &( vw.matrix ) );
310 // ---------------------------------------------------------------------
311 gslobj_matrix gslobj_matrix::operator-( gsl_matrix* o )
313 gslobj_matrix r( *this );
314 gsl_matrix_sub( r._gsl, o );
318 // ---------------------------------------------------------------------
319 gslobj_matrix& gslobj_matrix::operator-=( const gslobj_matrix& o )
321 gsl_matrix_sub( _gsl, o._gsl );
325 // ---------------------------------------------------------------------
326 gslobj_matrix& gslobj_matrix::operator-=( double o )
328 gsl_matrix_add_constant( _gsl, -o );
332 // ---------------------------------------------------------------------
333 gslobj_matrix& gslobj_matrix::operator-=( double* o )
335 gsl_matrix_view vw = gsl_matrix_view_array( o,
338 gsl_matrix_sub( _gsl, &( vw.matrix ) );
342 // ---------------------------------------------------------------------
343 gslobj_matrix& gslobj_matrix::operator-=( double* o[] )
345 gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
348 gsl_matrix_sub( _gsl, &( vw.matrix ) );
352 // ---------------------------------------------------------------------
353 gslobj_matrix& gslobj_matrix::operator-=( gsl_matrix* o )
355 gsl_matrix_sub( _gsl, o );
359 // ---------------------------------------------------------------------
360 gslobj_matrix gslobj_matrix::operator*( const gslobj_matrix& o )
362 gslobj_matrix r( *this );
363 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
369 // ---------------------------------------------------------------------
370 gslobj_matrix gslobj_matrix::operator*( double o )
372 gslobj_matrix r( *this );
373 gsl_matrix_scale( r._gsl, o );
377 // ---------------------------------------------------------------------
378 gslobj_matrix gslobj_matrix::operator*( double* o )
380 gslobj_matrix r( *this );
381 gsl_matrix_view vw = gsl_matrix_view_array( o,
384 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
385 1.0, r._gsl, &( vw.matrix ),
390 // ---------------------------------------------------------------------
391 gslobj_matrix gslobj_matrix::operator*( double* o[] )
393 gslobj_matrix r( *this );
394 gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
397 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
398 1.0, r._gsl, &( vw.matrix ),
403 // ---------------------------------------------------------------------
404 gslobj_matrix gslobj_matrix::operator*( gsl_matrix* o )
406 gslobj_matrix r( *this );
407 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
413 // ---------------------------------------------------------------------
414 gslobj_vector gslobj_matrix::operator*( const gslobj_vector& o )
416 gslobj_vector ret( rows( ) );
417 gsl_blas_dgemv( CblasNoTrans, 1.0, _gsl,
418 ( gsl_vector* )o, 0.0,
419 ( gsl_vector* )ret );
423 // ---------------------------------------------------------------------
424 gslobj_matrix& gslobj_matrix::operator*=( const gslobj_matrix& o )
426 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
432 // ---------------------------------------------------------------------
433 gslobj_matrix& gslobj_matrix::operator*=( double o )
435 gsl_matrix_scale( _gsl, o );
439 // ---------------------------------------------------------------------
440 gslobj_matrix& gslobj_matrix::operator*=( double* o )
442 gsl_matrix_view vw = gsl_matrix_view_array( o,
445 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
446 1.0, _gsl, &( vw.matrix ),
451 // ---------------------------------------------------------------------
452 gslobj_matrix& gslobj_matrix::operator*=( double* o[] )
454 gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
457 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
458 1.0, _gsl, &( vw.matrix ),
463 // ---------------------------------------------------------------------
464 gslobj_matrix& gslobj_matrix::operator*=( gsl_matrix* o )
466 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
472 // ---------------------------------------------------------------------
473 /*==PS========================================
474 gslobj_matrix gslobj_matrix::transpose( )
476 gslobj_matrix r( _gsl->size2, _gsl->size1 );
477 gsl_matrix_transpose_memcpy( r._gsl, _gsl );
480 ==PS========================================*/
481 // ---------------------------------------------------------------------
482 /*==PS========================================
483 gslobj_matrix gslobj_matrix::invertLU( )
485 gslobj_matrix r( *this ), lu( *this );
486 gsl_permutation* perm = gsl_permutation_calloc( rows( ) );
489 gsl_linalg_LU_decomp( lu._gsl, perm, &sign );
490 gsl_linalg_LU_invert( lu._gsl, perm, r._gsl );
493 ==PS========================================*/
494 // ---------------------------------------------------------------------
495 /*==PS========================================
496 double gslobj_matrix::determinant( )
498 gslobj_matrix lu( *this );
499 gsl_permutation* perm = gsl_permutation_calloc( rows( ) );
502 gsl_linalg_LU_decomp( lu._gsl, perm, &sign );
503 return( gsl_linalg_LU_det( lu._gsl, sign ) );
505 ==PS========================================*/
506 // ---------------------------------------------------------------------
507 /*==PS========================================
508 void gslobj_matrix::identity( )
510 gsl_matrix_set_identity( _gsl );
512 ==PS========================================*/
513 // eof - gslobj_matrix.cxx