1 /*=========================================================================
4 Module: $RCSfile: matrix.cxx,v $
6 Date: $Date: 2008/10/31 16:32:56 $
7 Version: $Revision: 1.1 $
9 Copyright: (c) 2002, 2003
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.
16 =========================================================================*/
20 #include <gsl/gsl_blas.h>
21 #include <gsl/gsl_linalg.h>
23 // ---------------------------------------------------------------------
24 std::ostream& operator<<( std::ostream& os, const gslobj_matrix& m )
26 for( int i = 0; i < m.size1( ); i++ ) {
28 for( int j = 0; j < m.size2( ); j++ )
29 os << m( i, j ) << " ";
36 // ---------------------------------------------------------------------
37 gslobj_matrix::gslobj_matrix( size_t r, size_t c )
40 _gsl = gsl_matrix_alloc( r, c );
43 // ---------------------------------------------------------------------
44 gslobj_matrix::gslobj_matrix( const gslobj_matrix& m )
47 _gsl = gsl_matrix_alloc( m._gsl->size1, m._gsl->size2 );
48 gsl_matrix_memcpy( _gsl, m._gsl );
51 // ---------------------------------------------------------------------
52 gslobj_matrix::gslobj_matrix( double* m, size_t r, size_t c )
55 _view = gsl_matrix_view_array( m, r, c );
56 _gsl = &( _view.matrix );
59 // ---------------------------------------------------------------------
60 gslobj_matrix::gslobj_matrix( double* m[], size_t r, size_t c )
63 _view = gsl_matrix_view_array( m[ 0 ], r, c );
64 _gsl = &( _view.matrix );
67 // ---------------------------------------------------------------------
68 gslobj_matrix::gslobj_matrix( gsl_matrix* m )
74 // ---------------------------------------------------------------------
75 /*==PS========================================
76 void gslobj_matrix::resize( size_t r, size_t c )
80 if( r != rows( ) || c != columns( ) )
81 gsl_matrix_free( _gsl );
82 _gsl = gsl_matrix_alloc( r, c );
86 ==PS========================================*/
87 // ---------------------------------------------------------------------
88 gslobj_matrix& gslobj_matrix::operator=( const gslobj_matrix& o )
90 gsl_matrix_memcpy( _gsl, o._gsl );
94 // ---------------------------------------------------------------------
95 gslobj_matrix& gslobj_matrix::operator=( double o )
97 gsl_matrix_set_all( _gsl, o );
101 // ---------------------------------------------------------------------
102 gslobj_matrix& gslobj_matrix::operator=( double* o )
104 memcpy( _gsl->data, o,
105 _gsl->size1 * _gsl->size2 * sizeof( double ) );
109 // ---------------------------------------------------------------------
110 gslobj_matrix& gslobj_matrix::operator=( double* o[] )
112 memcpy( _gsl->data, o[ 0 ],
113 _gsl->size1 * _gsl->size2 * sizeof( double ) );
117 // ---------------------------------------------------------------------
118 gslobj_matrix& gslobj_matrix::operator=( gsl_matrix* o )
120 gsl_matrix_memcpy( _gsl, o );
124 // ---------------------------------------------------------------------
125 double& gslobj_matrix::operator()( size_t i, size_t j )
127 return( _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] );
130 // ---------------------------------------------------------------------
131 const double& gslobj_matrix::operator()( size_t i, size_t j ) const
133 return( _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] );
136 // ---------------------------------------------------------------------
137 /*==PS========================================
138 void gslobj_matrix::setValue( size_t i, size_t j, double val )
140 _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] = val;
142 ==PS========================================*/
144 // ---------------------------------------------------------------------
145 bool gslobj_matrix::operator==( const gslobj_matrix& o ) const
149 if( _gsl->size1 != o._gsl->size1 || _gsl->size2 != o._gsl->size2 )
152 for( int i = 0; i < _gsl->size1 && ret; i++ )
153 for( int j = 0; j < _gsl->size2 && ret; j++ )
155 ( _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] !=
156 o._gsl->data[ ( ( i * _gsl->size2 ) + j ) ] );
160 // ---------------------------------------------------------------------
161 gslobj_matrix gslobj_matrix::operator+( const gslobj_matrix& o )
163 gslobj_matrix r( *this );
164 gsl_matrix_add( r._gsl, o._gsl );
168 // ---------------------------------------------------------------------
169 gslobj_matrix gslobj_matrix::operator+( double o )
171 gslobj_matrix r( *this );
172 gsl_matrix_add_constant( r._gsl, o );
176 // ---------------------------------------------------------------------
177 gslobj_matrix gslobj_matrix::operator+( double* o )
179 gslobj_matrix r( *this );
180 gsl_matrix_view vw = gsl_matrix_view_array( o,
183 gsl_matrix_add( r._gsl, &( vw.matrix ) );
187 // ---------------------------------------------------------------------
188 gslobj_matrix gslobj_matrix::operator+( double* o[] )
190 gslobj_matrix r( *this );
191 gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
194 gsl_matrix_add( r._gsl, &( vw.matrix ) );
198 // ---------------------------------------------------------------------
199 gslobj_matrix gslobj_matrix::operator+( gsl_matrix* o )
201 gslobj_matrix r( *this );
202 gsl_matrix_add( r._gsl, o );
206 // ---------------------------------------------------------------------
207 gslobj_matrix& gslobj_matrix::operator+=( const gslobj_matrix& o )
209 gsl_matrix_add( _gsl, o._gsl );
213 // ---------------------------------------------------------------------
214 gslobj_matrix& gslobj_matrix::operator+=( double o )
216 gsl_matrix_add_constant( _gsl, o );
220 // ---------------------------------------------------------------------
221 gslobj_matrix& gslobj_matrix::operator+=( double* o )
223 gsl_matrix_view vw = gsl_matrix_view_array( o,
226 gsl_matrix_add( _gsl, &( vw.matrix ) );
230 // ---------------------------------------------------------------------
231 gslobj_matrix& gslobj_matrix::operator+=( double* o[] )
233 gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
236 gsl_matrix_add( _gsl, &( vw.matrix ) );
240 // ---------------------------------------------------------------------
241 gslobj_matrix& gslobj_matrix::operator+=( gsl_matrix* o )
243 gsl_matrix_add( _gsl, o );
247 // ---------------------------------------------------------------------
248 gslobj_matrix gslobj_matrix::operator-( const gslobj_matrix& o )
250 gslobj_matrix r( *this );
251 gsl_matrix_sub( r._gsl, o._gsl );
255 // ---------------------------------------------------------------------
256 gslobj_matrix gslobj_matrix::operator-( double o )
258 gslobj_matrix r( *this );
259 gsl_matrix_add_constant( r._gsl, -o );
263 // ---------------------------------------------------------------------
264 gslobj_matrix gslobj_matrix::operator-( double* o )
266 gslobj_matrix r( *this );
267 gsl_matrix_view vw = gsl_matrix_view_array( o,
270 gsl_matrix_sub( r._gsl, &( vw.matrix ) );
274 // ---------------------------------------------------------------------
275 gslobj_matrix gslobj_matrix::operator-( double* o[] )
277 gslobj_matrix r( *this );
278 gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
281 gsl_matrix_sub( r._gsl, &( vw.matrix ) );
285 // ---------------------------------------------------------------------
286 gslobj_matrix gslobj_matrix::operator-( gsl_matrix* o )
288 gslobj_matrix r( *this );
289 gsl_matrix_sub( r._gsl, o );
293 // ---------------------------------------------------------------------
294 gslobj_matrix& gslobj_matrix::operator-=( const gslobj_matrix& o )
296 gsl_matrix_sub( _gsl, o._gsl );
300 // ---------------------------------------------------------------------
301 gslobj_matrix& gslobj_matrix::operator-=( double o )
303 gsl_matrix_add_constant( _gsl, -o );
307 // ---------------------------------------------------------------------
308 gslobj_matrix& gslobj_matrix::operator-=( double* o )
310 gsl_matrix_view vw = gsl_matrix_view_array( o,
313 gsl_matrix_sub( _gsl, &( vw.matrix ) );
317 // ---------------------------------------------------------------------
318 gslobj_matrix& gslobj_matrix::operator-=( double* o[] )
320 gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
323 gsl_matrix_sub( _gsl, &( vw.matrix ) );
327 // ---------------------------------------------------------------------
328 gslobj_matrix& gslobj_matrix::operator-=( gsl_matrix* o )
330 gsl_matrix_sub( _gsl, o );
334 // ---------------------------------------------------------------------
335 gslobj_matrix gslobj_matrix::operator*( const gslobj_matrix& o )
337 gslobj_matrix r( *this );
338 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
344 // ---------------------------------------------------------------------
345 gslobj_matrix gslobj_matrix::operator*( double o )
347 gslobj_matrix r( *this );
348 gsl_matrix_scale( r._gsl, o );
352 // ---------------------------------------------------------------------
353 gslobj_matrix gslobj_matrix::operator*( double* o )
355 gslobj_matrix r( *this );
356 gsl_matrix_view vw = gsl_matrix_view_array( o,
359 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
360 1.0, r._gsl, &( vw.matrix ),
365 // ---------------------------------------------------------------------
366 gslobj_matrix gslobj_matrix::operator*( double* o[] )
368 gslobj_matrix r( *this );
369 gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
372 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
373 1.0, r._gsl, &( vw.matrix ),
378 // ---------------------------------------------------------------------
379 gslobj_matrix gslobj_matrix::operator*( gsl_matrix* o )
381 gslobj_matrix r( *this );
382 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
388 // ---------------------------------------------------------------------
389 gslobj_vector gslobj_matrix::operator*( const gslobj_vector& o )
391 gslobj_vector ret( rows( ) );
392 gsl_blas_dgemv( CblasNoTrans, 1.0, _gsl,
393 ( gsl_vector* )o, 0.0,
394 ( gsl_vector* )ret );
398 // ---------------------------------------------------------------------
399 gslobj_matrix& gslobj_matrix::operator*=( const gslobj_matrix& o )
401 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
407 // ---------------------------------------------------------------------
408 gslobj_matrix& gslobj_matrix::operator*=( double o )
410 gsl_matrix_scale( _gsl, o );
414 // ---------------------------------------------------------------------
415 gslobj_matrix& gslobj_matrix::operator*=( double* o )
417 gsl_matrix_view vw = gsl_matrix_view_array( o,
420 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
421 1.0, _gsl, &( vw.matrix ),
426 // ---------------------------------------------------------------------
427 gslobj_matrix& gslobj_matrix::operator*=( double* o[] )
429 gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
432 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
433 1.0, _gsl, &( vw.matrix ),
438 // ---------------------------------------------------------------------
439 gslobj_matrix& gslobj_matrix::operator*=( gsl_matrix* o )
441 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
447 // ---------------------------------------------------------------------
448 /*==PS========================================
449 gslobj_matrix gslobj_matrix::transpose( )
451 gslobj_matrix r( _gsl->size2, _gsl->size1 );
452 gsl_matrix_transpose_memcpy( r._gsl, _gsl );
455 ==PS========================================*/
456 // ---------------------------------------------------------------------
457 /*==PS========================================
458 gslobj_matrix gslobj_matrix::invertLU( )
460 gslobj_matrix r( *this ), lu( *this );
461 gsl_permutation* perm = gsl_permutation_calloc( rows( ) );
464 gsl_linalg_LU_decomp( lu._gsl, perm, &sign );
465 gsl_linalg_LU_invert( lu._gsl, perm, r._gsl );
468 ==PS========================================*/
469 // ---------------------------------------------------------------------
470 /*==PS========================================
471 double gslobj_matrix::determinant( )
473 gslobj_matrix lu( *this );
474 gsl_permutation* perm = gsl_permutation_calloc( rows( ) );
477 gsl_linalg_LU_decomp( lu._gsl, perm, &sign );
478 return( gsl_linalg_LU_det( lu._gsl, sign ) );
480 ==PS========================================*/
481 // ---------------------------------------------------------------------
482 /*==PS========================================
483 void gslobj_matrix::identity( )
485 gsl_matrix_set_identity( _gsl );
487 ==PS========================================*/
488 // eof - gslobj_matrix.cxx