SawMatrix.cpp

/***************************************************************************
*   Copyright (C) 2005 by Nestor Aguirre                                  *
*   nfaguirrec@unal.edu.co                                                *
*                                                                         *
*   This program is free software; you can redistribute it and/or modify  *
*   it under the terms of the GNU General Public License as published by  *
*   the Free Software Foundation; either version 2 of the License, or     *
*   (at your option) any later version.                                   *
*                                                                         *
*   This program is distributed in the hope that it will be useful,       *
*   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
*   GNU General Public License for more details.                          *
*                                                                         *
*   You should have received a copy of the GNU General Public License     *
*   along with this program; if not, write to the                         *
*   Free Software Foundation, Inc.,                                       *
*   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
***************************************************************************/
#include "SawMatrix.h"

/***********************************
* Project
*/
#include <SawVector.h> 
#include <SawVectorList.h>
//**********************************

SawMatrix::SawMatrix( uint rows, uint cols, double a )
                : SawObject()
{
        this->rows = rows ;
        this->cols = cols ;
        matrix = gsl_matrix_alloc( rows, cols ) ;
        gsl_matrix_set_all( matrix, a ) ;
}

SawMatrix::SawMatrix( double* array, uint rows, uint cols )
{
        this->rows = rows ;
        this->cols = cols ;
        matrix = gsl_matrix_alloc( rows, cols ) ;

        for ( uint i = 0; i < rows; i++ ) {
                for ( uint j = 0; j < cols; j++ ) {
                        gsl_matrix_set( matrix, i, j, array[ cols * i + j ] ) ;
                }
        }
}

SawMatrix::SawMatrix( const SawMatrix& m ):
SawObject(m) 
{
        this->cols = m.getColumns() ;
        this->rows = m.getRows() ;

        this->matrix = gsl_matrix_alloc( rows, cols ) ;
        gsl_matrix_memcpy( this->matrix, m.getGSLMatrix() ) ;
        
}

void SawMatrix::operator = ( const SawMatrix& m )
{
        this->cols = m.getColumns() ;
        this->rows = m.getRows() ;
        
        this->matrix = gsl_matrix_alloc( cols, rows ) ;
        gsl_matrix_memcpy( this->matrix, m.getGSLMatrix() ) ;
        
}

SawMatrix::~SawMatrix()
{
        gsl_matrix_free( matrix ) ;
}

void SawMatrix::swapRows( uint i, uint j )
{
        gsl_matrix_swap_rows( matrix, i, j ) ;
}

void SawMatrix::swapColumns( uint i, uint j )
{
        gsl_matrix_swap_columns( matrix, i, j ) ;
}

gsl_matrix* SawMatrix::getGSLMatrix() const
{
        return matrix ;
}

void SawMatrix::setGSLMatrix( gsl_matrix* matrix )
{
        this->matrix = matrix ;
}

uint SawMatrix::getRows() const
{
        return rows ;
}

uint SawMatrix::getColumns() const
{
        return cols ;
}

double SawMatrix::get( uint i, uint j ) const
        {
                return gsl_matrix_get( matrix, i, j ) ;
        }

void SawMatrix::set( uint i, uint j, double value )
{
        return gsl_matrix_set( matrix, i, j, value ) ;
}

void SawMatrix::setIdentity()
{
        gsl_matrix_set_identity( matrix ) ;
}

void SawMatrix::setNull()
{
        gsl_matrix_set_zero( matrix ) ;
}

bool SawMatrix::operator == ( SawMatrix m )
{
        for ( uint i = 0; i < this->rows; i++ ) {
                for ( uint j = 0; j < this->cols; j++ ) {
                        if ( gsl_fcmp( get( i, j ), m.get( i, j ), thresholdDoubleComparation ) != 0 )
                                        return false ;
                }
        }


        return true ;
}

bool SawMatrix::operator != ( SawMatrix m )
{
        for ( uint i = 0; i < this->rows; i++ ) {
                for ( uint j = 0; j < this->cols; j++ ) {
                        if ( gsl_fcmp( get( i, j ), m.get( i, j ), thresholdDoubleComparation ) == 0 )
                                        return false ;
                }
        }


        return true ;
}

double SawMatrix::getMax( uint* i, uint* j ) const
{
        if ( i && j ) {
                gsl_matrix_max_index( matrix, i, j ) ;
        }

        return gsl_matrix_max( matrix ) ;
}

double SawMatrix::getMin( uint* i, uint* j ) const
{
        if ( i && j ) {
                gsl_matrix_min_index( matrix, i, j ) ;
        }

        return gsl_matrix_min( matrix ) ;
}

double SawMatrix::getDeterminant()
{
        gsl_matrix* tmp = gsl_matrix_alloc( this->rows, this->cols ) ;
        gsl_matrix_memcpy( tmp, matrix ) ;
        
        gsl_permutation* p = gsl_permutation_alloc( this->rows ) ;
        int s ;
        
        gsl_linalg_LU_decomp( tmp, p, &s ) ;
        
        return gsl_linalg_LU_det( tmp, s ) ;
}

SawMatrix SawMatrix::getInverse()
{
        SawMatrix output( rows, cols ) ;
        
        gsl_matrix* tmp = gsl_matrix_alloc( this->rows, this->cols ) ;
        gsl_matrix_memcpy( tmp, matrix ) ;
        
        gsl_matrix* invert = gsl_matrix_alloc( this->rows, this->cols ) ;
        gsl_permutation* p = gsl_permutation_alloc( this->rows ) ;
        
        int s ;
        
        gsl_linalg_LU_decomp( tmp, p, &s ) ;
        gsl_linalg_LU_invert( tmp, p, invert ) ;

        output.setGSLMatrix( invert ) ;
        
        gsl_permutation_free (p);
        
        return output ;
}

SawVector SawMatrix::getEigenValues( SawVectorList* eigenVectors )
{
        gsl_matrix* tmp = gsl_matrix_alloc( this->rows, this->cols ) ;
        gsl_matrix_memcpy( tmp, matrix ) ;

        // Se necesita que la matriz sea cuadrada
        gsl_vector * eval = gsl_vector_alloc( rows );
        gsl_matrix *evec = gsl_matrix_alloc( rows, rows );
        gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc( rows );

        gsl_eigen_symmv( tmp, eval, evec, w );
        gsl_eigen_symmv_free( w );
        gsl_eigen_symmv_sort( eval, evec, GSL_EIGEN_SORT_ABS_ASC );

        SawVector eigenValues( rows ) ;

        if( eigenVectors ){
                for ( int i=0; i<rows ; i++ ) {
                        /***************************************************************
                        * Se calcula el valor propio y se agrega a eigenValues
                        */
                        double eval_i = gsl_vector_get( eval, i );
                        eigenValues.set( i, eval_i ) ;
                        
                        /***************************************************************
                        * Se calcula el vector propio y se agrega a eigenVectors
                        */
                        gsl_vector * v = gsl_vector_alloc( rows ) ;
                        gsl_vector_memcpy( v, &gsl_matrix_column( evec, i ).vector );
                        SawVector tmp( rows ) ;
                        tmp.setGSLVector( v ) ;
                        
                        (*eigenVectors)[i] = tmp ;
                }
        }else{
                for ( int i=0; i<rows; i++ ) {
                        double eval_i = gsl_vector_get( eval, i );
                        eigenValues.set( i, eval_i ) ;
                }
        }


        return eigenValues ;
}

bool SawMatrix::isNull() const
{
        if ( gsl_matrix_isnull( matrix ) == 1 )
                return true ;
        else
                return false ;
}

SawMatrix SawMatrix::getTransPose()
{
        gsl_matrix * tmp = gsl_matrix_alloc( this->rows, this->cols ) ;
        gsl_matrix_transpose_memcpy( tmp, matrix ) ;

        SawMatrix output( this->rows, this->cols ) ;
        output.setGSLMatrix( tmp ) ;

        return output ;
}

SawMatrix operator + ( SawMatrix a, SawMatrix b )
{
        gsl_matrix * tmp = gsl_matrix_alloc( a.rows, a.cols ) ;
        gsl_matrix_memcpy( tmp, a.matrix ) ;
        gsl_matrix_add( tmp, b.matrix ) ;

        SawMatrix output( a.rows, a.cols ) ;
        output.setGSLMatrix( tmp ) ;

        return output ;
}

SawMatrix operator - ( SawMatrix a, SawMatrix b )
{
        gsl_matrix * tmp = gsl_matrix_alloc( a.rows, a.cols ) ;
        gsl_matrix_memcpy( tmp, a.matrix ) ;
        gsl_matrix_sub( tmp, b.matrix ) ;

        SawMatrix output( a.rows, a.cols ) ;
        output.setGSLMatrix( tmp ) ;

        return output ;
}

SawMatrix operator * ( SawMatrix a, SawMatrix b )
{
        gsl_matrix * tmp = gsl_matrix_alloc( a.rows, b.cols ) ;

        gsl_blas_dgemm ( CblasNoTrans, CblasNoTrans, 1.0, a.matrix, b.matrix, 0.0, tmp );

        SawMatrix output( a.rows, b.cols ) ;
        output.setGSLMatrix( tmp ) ;

        return output ;
}

SawMatrix operator * ( SawMatrix m, double a )
{
        gsl_matrix * tmp = gsl_matrix_alloc( m.rows, m.cols ) ;
        gsl_matrix_memcpy( tmp, m.getGSLMatrix() ) ;
        gsl_matrix_scale( tmp, a ) ;

        SawMatrix output( m.rows, m.cols ) ;
        output.setGSLMatrix( tmp ) ;

        return output ;
}

SawMatrix operator * ( double a, SawMatrix m )
{
        return m * a ;
}

ostream& operator << ( ostream& os, const SawMatrix& m )
{
        switch ( m.getFormatNumberOutputStream() ) {
                        case SawObject::FIXED:
                        os.setf( ios::fixed , ios::floatfield ) ;
                        break ;

                        case SawObject::SCIENTIFIC:
                        os.setf( ios::scientific, ios::floatfield ) ;
                        break ;
        }

        os.precision( m.precisionOutputStream ) ;
        os << endl ;
        //      os << "[ " ;
        for ( uint i = 0; i < m.rows; i++ ) {

                switch ( m.getFormatNumberOutputStream() ) {
                                case SawObject::FIXED:

                                os << setw( 4 + m.widthOutputStream - m.precisionOutputStream ) << "|" ;

                                for ( uint j = 0; j < m.cols; j++ )
                                        os << setw( m.widthOutputStream ) << gsl_matrix_get( m.matrix, i, j ) ;

                                os << setw( m.widthOutputStream - m.precisionOutputStream - 1 ) << "|" ;
                                os << endl ;

                                break ;

                                case SawObject::SCIENTIFIC:

                                os << setw( m.widthOutputStream - m.precisionOutputStream ) << "|" ;

                                for ( uint j = 0; j < m.cols; j++ )
                                        os << setw( m.widthOutputStream + 3 ) << gsl_matrix_get( m.matrix, i, j ) ;

                                os << setw( 4 ) << "|" ;
                                os << endl ;

                                break ;
                }
                //              if( i< (m.rows-1) ){
                //                      os << endl ;
                //                      os << "  |" ;
                //              }
        }
        //      os << "]" ;
        return os ;
}

istream& operator >> ( istream& is, const SawMatrix& m )
{
        double tmp ;
        for ( uint i = 0; i < m.rows; i++ ) {
                for ( uint j = 0; j < m.cols; j++ ) {
                        is >> tmp ;
                        gsl_matrix_set( m.matrix, i, j, tmp ) ;
                }
        }
        return is ;
}

double matrixSin( double x ){ return sin(x) ; }
double matrixCos( double x ){ return cos(x) ; }
double matrixTan( double x ){ return tan(x) ; }
double matrixAsin( double x ){ return asin(x) ; }
double matrixAcos( double x ){ return acos(x) ; }
double matrixAtan( double x ){ return atan(x) ; }
double matrixSinh( double x ){ return sinh(x) ; }
double matrixCosh( double x ){ return cosh(x) ; }
double matrixExp( double x ){ return exp(x) ; }
double matrixLn( double x ){ return log(x) ; }
double matrixLog( double x ){ return log10(x) ; }
double matrixSqrt( double x ){ return sqrt(x) ; }
double matrixPow( double x, double y ){ return pow(x, y) ; }

SawMatrix SawMatrix::genericFunction( SawMatrix a, double(* func)(double) )
{
        // Se necesita que la matriz sea cuadrada
        gsl_vector *eval = gsl_vector_alloc( a.getRows() );
        gsl_matrix *evec = gsl_matrix_alloc( a.getRows(), a.getColumns() );
        gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc( a.getRows() );
        
        // Se calculan los valores y vectores propios
        gsl_eigen_symmv( a.getGSLMatrix(), eval, evec, w );
        gsl_eigen_symmv_free( w );
        gsl_eigen_symmv_sort( eval, evec, GSL_EIGEN_SORT_ABS_ASC );
        
        // Se construye la matriz diagonal de los valores propios evalueando la funci�
        gsl_matrix* diagonal = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
        gsl_matrix_set_zero( diagonal )  ;
        for( int i=0; i<a.getRows(); i++ )
                gsl_matrix_set( diagonal, i, i, (*func)( gsl_vector_get( eval, i ) ) ) ;
                
        // Se calcula la matriz de vectores propios transpuesta
        gsl_matrix* evec_t = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
        gsl_matrix_transpose_memcpy( evec_t, evec ) ;
        
        // Se hacen las respectivas multiplicaciones
        gsl_matrix* work = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
        gsl_matrix_set_zero( work )  ;
        
        // A*B
        gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, evec, diagonal, 0.0, work );
        
        gsl_matrix* work2 = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
        gsl_matrix_set_zero( work2 )  ;
        
        // (A*B)*A^t
        gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, work, evec_t, 0.0, work2 );
        
        // Se almacena en la matriz de salida
        SawMatrix output( a.getRows(), a.getColumns() ) ;
        output.setGSLMatrix( work2 ) ;

//      SawMatrix d( a.getRows(), a.getColumns() ) ;
//      SawVector eigenValues( a.getRows() ) ;
//      SawMatrix eigenVectors( a.getRows(), a.getColumns() ) ;
//      eigenValues = a.getEigenValues( &eigenVectors ) ;
//      
//      for( int i=0; i<a.getRows(); i++ ){
//              d.set( i, i, (*func)( eigenValues.get(i) ) ) ;
//      }
//      
//      SawMatrix output( a.getRows(), a.getColumns() ) ;
//      output = eigenVectors*d*eigenVectors.getTransPose() ;
        
        return output ;
}

SawMatrix SawMatrix::genericFunction( SawMatrix a, double(* func)(double, double), double param )
{
        // Se necesita que la matriz sea cuadrada
        gsl_vector *eval = gsl_vector_alloc( a.getRows() );
        gsl_matrix *evec = gsl_matrix_alloc( a.getRows(), a.getColumns() );
        gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc( a.getRows() );
        
        // Se calculan los valores y vectores propios
        gsl_eigen_symmv( a.getGSLMatrix(), eval, evec, w );
        gsl_eigen_symmv_free( w );
        gsl_eigen_symmv_sort( eval, evec, GSL_EIGEN_SORT_ABS_ASC );
        
        // Se construye la matriz diagonal de los valores propios evalueando la funci�
        gsl_matrix* diagonal = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
        gsl_matrix_set_zero( diagonal )  ;
        for( int i=0; i<a.getRows(); i++ )
                gsl_matrix_set( diagonal, i, i, (*func)( gsl_vector_get( eval, i ), param ) ) ;
                
        // Se calcula la matriz de vectores propios transpuesta
        gsl_matrix* evec_t = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
        gsl_matrix_transpose_memcpy( evec_t, evec ) ;
        
        // Se hacen las respectivas multiplicaciones
        gsl_matrix* work = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
        gsl_matrix_set_zero( work )  ;
        
        // A*B
        gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, evec, diagonal, 0.0, work );
        
        gsl_matrix* work2 = gsl_matrix_alloc( a.getRows(), a.getColumns() ) ;
        gsl_matrix_set_zero( work2 )  ;
        
        // (A*B)*A^t
        gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, work, evec_t, 0.0, work2 );
        
        // Se almacena en la matriz de salida
        SawMatrix output( a.getRows(), a.getColumns() ) ;
        output.setGSLMatrix( work2 ) ;
        
        return output ;
}

SawMatrix SawMatrix::sin( SawMatrix a )
{
        return SawMatrix::genericFunction( a, *matrixSin ) ;
}

SawMatrix SawMatrix::cos( SawMatrix a )
{
        return SawMatrix::genericFunction( a, *matrixCos ) ;
}

SawMatrix SawMatrix::tan( SawMatrix a )
{
        return SawMatrix::genericFunction( a, *matrixTan ) ;
}

SawMatrix SawMatrix::asin( SawMatrix a )
{
        return SawMatrix::genericFunction( a, *matrixAsin ) ;
}

SawMatrix SawMatrix::acos( SawMatrix a )
{
        return SawMatrix::genericFunction( a, *matrixAcos ) ;
}

SawMatrix SawMatrix::atan( SawMatrix a )
{
        return SawMatrix::genericFunction( a, *matrixAtan ) ;
}

SawMatrix SawMatrix::sinh( SawMatrix a )
{
        return SawMatrix::genericFunction( a, *matrixSinh ) ;
}

SawMatrix SawMatrix::cosh( SawMatrix a )
{
        return SawMatrix::genericFunction( a, *matrixCosh ) ;
}

SawMatrix SawMatrix::exp( SawMatrix a )
{
        return SawMatrix::genericFunction( a, *matrixExp ) ;
}

SawMatrix SawMatrix::ln( SawMatrix a )
{
        return SawMatrix::genericFunction( a, *matrixLn ) ;
}

SawMatrix SawMatrix::log( SawMatrix a )
{
        return SawMatrix::genericFunction( a, *matrixLog ) ;
}

SawMatrix SawMatrix::sqrt( SawMatrix a )
{
        return SawMatrix::genericFunction( a, *matrixSqrt ) ;
}

SawMatrix SawMatrix::pow( SawMatrix a, double e )
{
        return SawMatrix::genericFunction( a, *matrixPow, e ) ;
}

---