SawMatrix.cpp

00001 /***************************************************************************
00002 *   Copyright (C) 2005 by Nestor Aguirre                                  *
00003 *   nfaguirrec@unal.edu.co                                                *
00004 *                                                                         *
00005 *   This program is free software; you can redistribute it and/or modify  *
00006 *   it under the terms of the GNU General Public License as published by  *
00007 *   the Free Software Foundation; either version 2 of the License, or     *
00008 *   (at your option) any later version.                                   *
00009 *                                                                         *
00010 *   This program is distributed in the hope that it will be useful,       *
00011 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00012 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00013 *   GNU General Public License for more details.                          *
00014 *                                                                         *
00015 *   You should have received a copy of the GNU General Public License     *
00016 *   along with this program; if not, write to the                         *
00017 *   Free Software Foundation, Inc.,                                       *
00018 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
00019 ***************************************************************************/
00020 #include "SawMatrix.h"
00021 
00022 /***********************************
00023 * Project
00024 */
00025 #include <SawVector.h> 
00026 #include <SawVectorList.h>
00027 //**********************************
00028 
00032 SawMatrix::SawMatrix( uint rows, uint cols, double a )
00033                 : SawObject()
00034 {
00035         this->rows = rows ;
00036         this->cols = cols ;
00037         matrix = gsl_matrix_alloc( rows, cols ) ;
00038         gsl_matrix_set_all( matrix, a ) ;
00039 }
00040 
00047 SawMatrix::SawMatrix( double* array, uint rows, uint cols )
00048 {
00049         this->rows = rows ;
00050         this->cols = cols ;
00051         matrix = gsl_matrix_alloc( rows, cols ) ;
00052 
00053         for ( uint i = 0; i < rows; i++ ) {
00054                 for ( uint j = 0; j < cols; j++ ) {
00055                         gsl_matrix_set( matrix, i, j, array[ cols * i + j ] ) ;
00056                 }
00057         }
00059 }
00060 
00065 SawMatrix::SawMatrix( const SawMatrix& m ):
00066 SawObject(m) 
00067 {
00068         this->cols = m.getColumns() ;
00069         this->rows = m.getRows() ;
00070 
00071         this->matrix = gsl_matrix_alloc( rows, cols ) ;
00072         gsl_matrix_memcpy( this->matrix, m.getGSLMatrix() ) ;
00073         
00075 }
00076 
00081 void SawMatrix::operator = ( const SawMatrix& m )
00082 {
00083         this->cols = m.getColumns() ;
00084         this->rows = m.getRows() ;
00085         
00086         this->matrix = gsl_matrix_alloc( cols, rows ) ;
00087         gsl_matrix_memcpy( this->matrix, m.getGSLMatrix() ) ;
00088         
00090 }
00091 
00095 SawMatrix::~SawMatrix()
00096 {
00097         gsl_matrix_free( matrix ) ;
00098 }
00099 
00105 void SawMatrix::swapRows( uint i, uint j )
00106 {
00107         gsl_matrix_swap_rows( matrix, i, j ) ;
00109 }
00110 
00116 void SawMatrix::swapColumns( uint i, uint j )
00117 {
00118         gsl_matrix_swap_columns( matrix, i, j ) ;
00120 }
00121 
00125 gsl_matrix* SawMatrix::getGSLMatrix() const
00126 {
00127         return matrix ;
00128 }
00129 
00133 void SawMatrix::setGSLMatrix( gsl_matrix* matrix )
00134 {
00135         this->matrix = matrix ;
00136 }
00137 
00142 uint SawMatrix::getRows() const
00143 {
00144         return rows ;
00145 }
00146 
00151 uint SawMatrix::getColumns() const
00152 {
00153         return cols ;
00154 }
00155 
00162 double SawMatrix::get( uint i, uint j ) const
00163         {
00164                 return gsl_matrix_get( matrix, i, j ) ;
00166         }
00167 
00174 void SawMatrix::set( uint i, uint j, double value )
00175 {
00176         return gsl_matrix_set( matrix, i, j, value ) ;
00178 }
00179 
00184 void SawMatrix::setIdentity()
00185 {
00186         gsl_matrix_set_identity( matrix ) ;
00187 }
00188 
00193 void SawMatrix::setNull()
00194 {
00195         gsl_matrix_set_zero( matrix ) ;
00196 }
00197 
00203 bool SawMatrix::operator == ( SawMatrix m )
00204 {
00205         for ( uint i = 0; i < this->rows; i++ ) {
00206                 for ( uint j = 0; j < this->cols; j++ ) {
00207                         if ( gsl_fcmp( get( i, j ), m.get( i, j ), thresholdDoubleComparation ) != 0 )
00208                                         return false ;
00209                 }
00210         }
00211 
00213 
00214         return true ;
00215 }
00216 
00222 bool SawMatrix::operator != ( SawMatrix m )
00223 {
00224         for ( uint i = 0; i < this->rows; i++ ) {
00225                 for ( uint j = 0; j < this->cols; j++ ) {
00226                         if ( gsl_fcmp( get( i, j ), m.get( i, j ), thresholdDoubleComparation ) == 0 )
00227                                         return false ;
00228                 }
00229         }
00230 
00232 
00233         return true ;
00234 }
00235 
00242 double SawMatrix::getMax( uint* i, uint* j ) const
00243 {
00244         if ( i && j ) {
00245                 gsl_matrix_max_index( matrix, i, j ) ;
00246         }
00247 
00248         return gsl_matrix_max( matrix ) ;
00249 }
00250 
00257 double SawMatrix::getMin( uint* i, uint* j ) const
00258 {
00259         if ( i && j ) {
00260                 gsl_matrix_min_index( matrix, i, j ) ;
00261         }
00262 
00263         return gsl_matrix_min( matrix ) ;
00264 }
00265 
00270 double SawMatrix::getDeterminant()
00271 {
00272         gsl_matrix* tmp = gsl_matrix_alloc( this->rows, this->cols ) ;
00273         gsl_matrix_memcpy( tmp, matrix ) ;
00274         
00275         gsl_permutation* p = gsl_permutation_alloc( this->rows ) ;
00276         int s ;
00277         
00278         gsl_linalg_LU_decomp( tmp, p, &s ) ;
00279         
00280         return gsl_linalg_LU_det( tmp, s ) ;
00281 }
00282 
00287 SawMatrix SawMatrix::getInverse()
00288 {
00289         SawMatrix output( rows, cols ) ;
00290         
00291         gsl_matrix* tmp = gsl_matrix_alloc( this->rows, this->cols ) ;
00292         gsl_matrix_memcpy( tmp, matrix ) ;
00293         
00294         gsl_matrix* invert = gsl_matrix_alloc( this->rows, this->cols ) ;
00295         gsl_permutation* p = gsl_permutation_alloc( this->rows ) ;
00296         
00297         int s ;
00298         
00299         gsl_linalg_LU_decomp( tmp, p, &s ) ;
00300         gsl_linalg_LU_invert( tmp, p, invert ) ;
00301 
00302         output.setGSLMatrix( invert ) ;
00303         
00304         gsl_permutation_free (p);
00305         
00307         return output ;
00308 }
00309 
00314 SawVector SawMatrix::getEigenValues( SawVectorList* eigenVectors )
00315 {
00316         gsl_matrix* tmp = gsl_matrix_alloc( this->rows, this->cols ) ;
00317         gsl_matrix_memcpy( tmp, matrix ) ;
00318 
00319         // Se necesita que la matriz sea cuadrada
00320         gsl_vector * eval = gsl_vector_alloc( rows );
00321         gsl_matrix *evec = gsl_matrix_alloc( rows, rows );
00322         gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc( rows );
00323 
00324         gsl_eigen_symmv( tmp, eval, evec, w );
00325         gsl_eigen_symmv_free( w );
00326         gsl_eigen_symmv_sort( eval, evec, GSL_EIGEN_SORT_ABS_ASC );
00327 
00328         SawVector eigenValues( rows ) ;
00329 
00330         if( eigenVectors ){
00331                 for ( int i=0; i<rows ; i++ ) {
00332                         /***************************************************************
00333                         * Se calcula el valor propio y se agrega a eigenValues
00334                         */
00335                         double eval_i = gsl_vector_get( eval, i );
00336                         eigenValues.set( i, eval_i ) ;
00337                         
00338                         /***************************************************************
00339                         * Se calcula el vector propio y se agrega a eigenVectors
00340                         */
00341                         gsl_vector * v = gsl_vector_alloc( rows ) ;
00342                         gsl_vector_memcpy( v, &gsl_matrix_column( evec, i ).vector );
00343                         SawVector tmp( rows ) ;
00344                         tmp.setGSLVector( v ) ;
00345                         
00346                         (*eigenVectors)[i] = tmp ;
00347                 }
00348         }else{
00349                 for ( int i=0; i<rows; i++ ) {
00350                         double eval_i = gsl_vector_get( eval, i );
00351                         eigenValues.set( i, eval_i ) ;
00352                 }
00353         }
00354 
00357 
00358         return eigenValues ;
00359 }
00360 
00365 bool SawMatrix::isNull() const
00366 {
00367         if ( gsl_matrix_isnull( matrix ) == 1 )
00368                 return true ;
00369         else
00370                 return false ;
00371 }
00372 
00377 SawMatrix SawMatrix::getTransPose()
00378 {
00379         gsl_matrix * tmp = gsl_matrix_alloc( this->rows, this->cols ) ;
00380         gsl_matrix_transpose_memcpy( tmp, matrix ) ;
00381 
00382         SawMatrix output( this->rows, this->cols ) ;
00383         output.setGSLMatrix( tmp ) ;
00384 
00385         return output ;
00386 }
00387 
00393 SawMatrix operator + ( SawMatrix a, SawMatrix b )
00394 {
00395         gsl_matrix * tmp = gsl_matrix_alloc( a.rows, a.cols ) ;
00396         gsl_matrix_memcpy( tmp, a.matrix ) ;
00397         gsl_matrix_add( tmp, b.matrix ) ;
00398 
00399         SawMatrix output( a.rows, a.cols ) ;
00400         output.setGSLMatrix( tmp ) ;
00401 
00403         return output ;
00404 }
00405 
00411 SawMatrix operator - ( SawMatrix a, SawMatrix b )
00412 {
00413         gsl_matrix * tmp = gsl_matrix_alloc( a.rows, a.cols ) ;
00414         gsl_matrix_memcpy( tmp, a.matrix ) ;
00415         gsl_matrix_sub( tmp, b.matrix ) ;
00416 
00417         SawMatrix output( a.rows, a.cols ) ;
00418         output.setGSLMatrix( tmp ) ;
00419 
00421         return output ;
00422 }
00423 
00429 SawMatrix operator * ( SawMatrix a, SawMatrix b )
00430 {
00431         gsl_matrix * tmp = gsl_matrix_alloc( a.rows, b.cols ) ;
00432 
00433         gsl_blas_dgemm ( CblasNoTrans, CblasNoTrans, 1.0, a.matrix, b.matrix, 0.0, tmp );
00434 
00435         SawMatrix output( a.rows, b.cols ) ;
00436         output.setGSLMatrix( tmp ) ;
00437 
00439         return output ;
00440 }
00441 
00448 SawMatrix operator * ( SawMatrix m, double a )
00449 {
00450         gsl_matrix * tmp = gsl_matrix_alloc( m.rows, m.cols ) ;
00451         gsl_matrix_memcpy( tmp, m.getGSLMatrix() ) ;
00452         gsl_matrix_scale( tmp, a ) ;
00453 
00454         SawMatrix output( m.rows, m.cols ) ;
00455         output.setGSLMatrix( tmp ) ;
00456 
00457         return output ;
00458 }
00459 
00466 SawMatrix operator * ( double a, SawMatrix m )
00467 {
00468         return m * a ;
00469 }
00470 
00477 ostream& operator << ( ostream& os, const SawMatrix& m )
00478 {
00479         switch ( m.getFormatNumberOutputStream() ) {
00480                         case SawObject::FIXED:
00481                         os.setf( ios::fixed , ios::floatfield ) ;
00482                         break ;
00483 
00484                         case SawObject::SCIENTIFIC:
00485                         os.setf( ios::scientific, ios::floatfield ) ;
00486                         break ;
00487         }
00488 
00489         os.precision( m.precisionOutputStream ) ;
00490         os << endl ;
00491         //      os << "[ " ;
00492         for ( uint i = 0; i < m.rows; i++ ) {
00493 
00494                 switch ( m.getFormatNumberOutputStream() ) {
00495                                 case SawObject::FIXED:
00496 
00497                                 os << setw( 4 + m.widthOutputStream - m.precisionOutputStream ) << "|" ;
00498 
00499                                 for ( uint j = 0; j < m.cols; j++ )
00500                                         os << setw( m.widthOutputStream ) << gsl_matrix_get( m.matrix, i, j ) ;
00501 
00502                                 os << setw( m.widthOutputStream - m.precisionOutputStream - 1 ) << "|" ;
00503                                 os << endl ;
00504 
00505                                 break ;
00506 
00507                                 case SawObject::SCIENTIFIC:
00508 
00509                                 os << setw( m.widthOutputStream - m.precisionOutputStream ) << "|" ;
00510 
00511                                 for ( uint j = 0; j < m.cols; j++ )
00512                                         os << setw( m.widthOutputStream + 3 ) << gsl_matrix_get( m.matrix, i, j ) ;
00513 
00514                                 os << setw( 4 ) << "|" ;
00515                                 os << endl ;
00516 
00517                                 break ;
00518                 }
00519                 //              if( i< (m.rows-1) ){
00520                 //                      os << endl ;
00521                 //                      os << "  |" ;
00522                 //              }
00523         }
00524         //      os << "]" ;
00525         return os ;
00526 }
00527 
00534 istream& operator >> ( istream& is, const SawMatrix& m )
00535 {
00536         double tmp ;
00537         for ( uint i = 0; i < m.rows; i++ ) {
00538                 for ( uint j = 0; j < m.cols; j++ ) {
00539                         is >> tmp ;
00540                         gsl_matrix_set( m.matrix, i, j, tmp ) ;
00541                 }
00542         }
00543         return is ;
00544 }
00545 
00546 double matrixSin( double x ){ return sin(x) ; }
00547 double matrixCos( double x ){ return cos(x) ; }
00548 double matrixTan( double x ){ return tan(x) ; }
00549 double matrixAsin( double x ){ return asin(x) ; }
00550 double matrixAcos( double x ){ return acos(x) ; }
00551 double matrixAtan( double x ){ return atan(x) ; }
00552 double matrixSinh( double x ){ return sinh(x) ; }
00553 double matrixCosh( double x ){ return cosh(x) ; }
00554 double matrixExp( double x ){ return exp(x) ; }
00555 double matrixLn( double x ){ return log(x) ; }
00556 double matrixLog( double x ){ return log10(x) ; }
00557 double matrixSqrt( double x ){ return sqrt(x) ; }
00558 double matrixPow( double x, double y ){ return pow(x, y) ; }
00559 
00560 SawMatrix SawMatrix::genericFunction( SawMatrix a, double(* func)(double) )
00561 {
00562         // Se necesita que la matriz sea cuadrada
00563         gsl_vector *eval = gsl_vector_alloc( a.getRows() );
00564         gsl_matrix *evec = gsl_matrix_alloc( a.getRows(), a.getColumns() );
00565         gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc( a.getRows() );
00566         
00567         // Se calculan los valores y vectores propios
00568         gsl_eigen_symmv( a.getGSLMatrix(), eval, evec, w );
00569         gsl_eigen_symmv_free( w );
00570         gsl_eigen_symmv_sort( eval, evec, GSL_EIGEN_SORT_ABS_ASC );
00571         
00572         // Se construye la matriz diagonal de los valores propios evalueando la funci�
00573         gsl_matrix* diagonal = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
00574         gsl_matrix_set_zero( diagonal )  ;
00575         for( int i=0; i<a.getRows(); i++ )
00576                 gsl_matrix_set( diagonal, i, i, (*func)( gsl_vector_get( eval, i ) ) ) ;
00577                 
00578         // Se calcula la matriz de vectores propios transpuesta
00579         gsl_matrix* evec_t = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
00580         gsl_matrix_transpose_memcpy( evec_t, evec ) ;
00581         
00582         // Se hacen las respectivas multiplicaciones
00583         gsl_matrix* work = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
00584         gsl_matrix_set_zero( work )  ;
00585         
00586         // A*B
00587         gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, evec, diagonal, 0.0, work );
00588         
00589         gsl_matrix* work2 = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
00590         gsl_matrix_set_zero( work2 )  ;
00591         
00592         // (A*B)*A^t
00593         gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, work, evec_t, 0.0, work2 );
00594         
00595         // Se almacena en la matriz de salida
00596         SawMatrix output( a.getRows(), a.getColumns() ) ;
00597         output.setGSLMatrix( work2 ) ;
00598 
00599 //      SawMatrix d( a.getRows(), a.getColumns() ) ;
00600 //      SawVector eigenValues( a.getRows() ) ;
00601 //      SawMatrix eigenVectors( a.getRows(), a.getColumns() ) ;
00602 //      eigenValues = a.getEigenValues( &eigenVectors ) ;
00603 //      
00604 //      for( int i=0; i<a.getRows(); i++ ){
00605 //              d.set( i, i, (*func)( eigenValues.get(i) ) ) ;
00606 //      }
00607 //      
00608 //      SawMatrix output( a.getRows(), a.getColumns() ) ;
00609 //      output = eigenVectors*d*eigenVectors.getTransPose() ;
00610         
00613         return output ;
00614 }
00615 
00616 SawMatrix SawMatrix::genericFunction( SawMatrix a, double(* func)(double, double), double param )
00617 {
00618         // Se necesita que la matriz sea cuadrada
00619         gsl_vector *eval = gsl_vector_alloc( a.getRows() );
00620         gsl_matrix *evec = gsl_matrix_alloc( a.getRows(), a.getColumns() );
00621         gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc( a.getRows() );
00622         
00623         // Se calculan los valores y vectores propios
00624         gsl_eigen_symmv( a.getGSLMatrix(), eval, evec, w );
00625         gsl_eigen_symmv_free( w );
00626         gsl_eigen_symmv_sort( eval, evec, GSL_EIGEN_SORT_ABS_ASC );
00627         
00628         // Se construye la matriz diagonal de los valores propios evalueando la funci�
00629         gsl_matrix* diagonal = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
00630         gsl_matrix_set_zero( diagonal )  ;
00631         for( int i=0; i<a.getRows(); i++ )
00632                 gsl_matrix_set( diagonal, i, i, (*func)( gsl_vector_get( eval, i ), param ) ) ;
00633                 
00634         // Se calcula la matriz de vectores propios transpuesta
00635         gsl_matrix* evec_t = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
00636         gsl_matrix_transpose_memcpy( evec_t, evec ) ;
00637         
00638         // Se hacen las respectivas multiplicaciones
00639         gsl_matrix* work = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
00640         gsl_matrix_set_zero( work )  ;
00641         
00642         // A*B
00643         gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, evec, diagonal, 0.0, work );
00644         
00645         gsl_matrix* work2 = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
00646         gsl_matrix_set_zero( work2 )  ;
00647         
00648         // (A*B)*A^t
00649         gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, work, evec_t, 0.0, work2 );
00650         
00651         // Se almacena en la matriz de salida
00652         SawMatrix output( a.getRows(), a.getColumns() ) ;
00653         output.setGSLMatrix( work2 ) ;
00654         
00657         return output ;
00658 }
00659 
00665 SawMatrix SawMatrix::sin( SawMatrix a )
00666 {
00667         return SawMatrix::genericFunction( a, *matrixSin ) ;
00668 }
00669 
00675 SawMatrix SawMatrix::cos( SawMatrix a )
00676 {
00677         return SawMatrix::genericFunction( a, *matrixCos ) ;
00678 }
00679 
00685 SawMatrix SawMatrix::tan( SawMatrix a )
00686 {
00687         return SawMatrix::genericFunction( a, *matrixTan ) ;
00688 }
00689 
00695 SawMatrix SawMatrix::asin( SawMatrix a )
00696 {
00697         return SawMatrix::genericFunction( a, *matrixAsin ) ;
00698 }
00699 
00705 SawMatrix SawMatrix::acos( SawMatrix a )
00706 {
00707         return SawMatrix::genericFunction( a, *matrixAcos ) ;
00708 }
00709 
00715 SawMatrix SawMatrix::atan( SawMatrix a )
00716 {
00717         return SawMatrix::genericFunction( a, *matrixAtan ) ;
00718 }
00719 
00725 SawMatrix SawMatrix::sinh( SawMatrix a )
00726 {
00727         return SawMatrix::genericFunction( a, *matrixSinh ) ;
00728 }
00729 
00735 SawMatrix SawMatrix::cosh( SawMatrix a )
00736 {
00737         return SawMatrix::genericFunction( a, *matrixCosh ) ;
00738 }
00739 
00745 SawMatrix SawMatrix::exp( SawMatrix a )
00746 {
00747         return SawMatrix::genericFunction( a, *matrixExp ) ;
00748 }
00749 
00755 SawMatrix SawMatrix::ln( SawMatrix a )
00756 {
00757         return SawMatrix::genericFunction( a, *matrixLn ) ;
00758 }
00759 
00765 SawMatrix SawMatrix::log( SawMatrix a )
00766 {
00767         return SawMatrix::genericFunction( a, *matrixLog ) ;
00768 }
00769 
00775 SawMatrix SawMatrix::sqrt( SawMatrix a )
00776 {
00777         return SawMatrix::genericFunction( a, *matrixSqrt ) ;
00778 }
00779 
00785 SawMatrix SawMatrix::pow( SawMatrix a, double e )
00786 {
00787         return SawMatrix::genericFunction( a, *matrixPow, e ) ;
00788 }

Generado el Sun Jul 22 18:05:41 2007 para SAW por  doxygen 1.5.1