/*=========================================================================

  Program:   wxMaracas
  Module:    $RCSfile: matrix.cxx,v $
  Language:  C++
  Date:      $Date: 2008/10/31 16:32:56 $
  Version:   $Revision: 1.1 $

  Copyright: (c) 2002, 2003
  License:
  
     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

#include "gslobj.hxx"
#include <string>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

// ---------------------------------------------------------------------
std::ostream& operator<<( std::ostream& os, const gslobj_matrix& m )
{
    for( int i = 0; i < m.size1( ); i++ ) {

	for( int j = 0; j < m.size2( ); j++ )
	    os << m( i, j ) << " ";
	os << std::endl;

    } // rof
    return( os );
}

// ---------------------------------------------------------------------
gslobj_matrix::gslobj_matrix( size_t r, size_t c )
    : _copy( true )
{
    _gsl = gsl_matrix_alloc( r, c );
}

// ---------------------------------------------------------------------
gslobj_matrix::gslobj_matrix( const gslobj_matrix& m )
    : _copy( true )
{
    _gsl = gsl_matrix_alloc( m._gsl->size1, m._gsl->size2 );
    gsl_matrix_memcpy( _gsl, m._gsl );
}

// ---------------------------------------------------------------------
gslobj_matrix::gslobj_matrix( double* m, size_t r, size_t c )
    : _copy( false )
{
    _view = gsl_matrix_view_array( m, r, c );
    _gsl = &( _view.matrix );
}

// ---------------------------------------------------------------------
gslobj_matrix::gslobj_matrix( double* m[], size_t r, size_t c )
    : _copy( false )
{
    _view = gsl_matrix_view_array( m[ 0 ], r, c );
    _gsl = &( _view.matrix );
}

// ---------------------------------------------------------------------
gslobj_matrix::gslobj_matrix( gsl_matrix* m )
    : _copy( false )
{
    _gsl = m;
}

// ---------------------------------------------------------------------
/*==PS========================================
void gslobj_matrix::resize( size_t r, size_t c )
{
    if( _copy ) {

	if( r != rows( ) || c != columns( ) )
	    gsl_matrix_free( _gsl );
	_gsl = gsl_matrix_alloc( r, c );

    } // fi
}
==PS========================================*/
// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator=( const gslobj_matrix& o )
{
    gsl_matrix_memcpy( _gsl, o._gsl );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator=( double o )
{
    gsl_matrix_set_all( _gsl, o );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator=( double* o )
{
    memcpy( _gsl->data, o,
	    _gsl->size1 * _gsl->size2 * sizeof( double ) );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator=( double* o[] )
{
    memcpy( _gsl->data, o[ 0 ],
	    _gsl->size1 * _gsl->size2 * sizeof( double ) );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator=( gsl_matrix* o )
{
    gsl_matrix_memcpy( _gsl, o );
    return( *this );
}

// ---------------------------------------------------------------------
double& gslobj_matrix::operator()( size_t i, size_t j )
{
    return( _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] );
}

// ---------------------------------------------------------------------
const double& gslobj_matrix::operator()( size_t i, size_t j ) const
{
    return( _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] );
}

// ---------------------------------------------------------------------
/*==PS========================================
void gslobj_matrix::setValue( size_t i, size_t j, double val )
{
    _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] = val;
}
==PS========================================*/

// ---------------------------------------------------------------------
bool gslobj_matrix::operator==( const gslobj_matrix& o ) const
{
    bool ret = true;

    if( _gsl->size1 != o._gsl->size1 || _gsl->size2 != o._gsl->size2 )
	return( false );

    for( int i = 0; i < _gsl->size1 && ret; i++ )
	for( int j = 0; j < _gsl->size2 && ret; j++ )
	    ret = ret &&
		( _gsl->data[ ( ( i * _gsl->size2 ) + j ) ] !=
		  o._gsl->data[ ( ( i * _gsl->size2 ) + j ) ] );
    return( ret );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator+( const gslobj_matrix& o )
{
    gslobj_matrix r( *this );
    gsl_matrix_add( r._gsl, o._gsl );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator+( double o )
{
    gslobj_matrix r( *this );
    gsl_matrix_add_constant( r._gsl, o );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator+( double* o )
{
    gslobj_matrix r( *this );
    gsl_matrix_view vw = gsl_matrix_view_array( o,
						_gsl->size1,
						_gsl->size2 );
    gsl_matrix_add( r._gsl, &( vw.matrix ) );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator+( double* o[] )
{
    gslobj_matrix r( *this );
    gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
						_gsl->size1,
						_gsl->size2 );
    gsl_matrix_add( r._gsl, &( vw.matrix ) );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator+( gsl_matrix* o )
{
    gslobj_matrix r( *this );
    gsl_matrix_add( r._gsl, o );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator+=( const gslobj_matrix& o )
{
    gsl_matrix_add( _gsl, o._gsl );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator+=( double o )
{
    gsl_matrix_add_constant( _gsl, o );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator+=( double* o )
{
    gsl_matrix_view vw = gsl_matrix_view_array( o,
						_gsl->size1,
						_gsl->size2 );
    gsl_matrix_add( _gsl, &( vw.matrix ) );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator+=( double* o[] )
{
    gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
						_gsl->size1,
						_gsl->size2 );
    gsl_matrix_add( _gsl, &( vw.matrix ) );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator+=( gsl_matrix* o )
{
    gsl_matrix_add( _gsl, o );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator-( const gslobj_matrix& o )
{
    gslobj_matrix r( *this );
    gsl_matrix_sub( r._gsl, o._gsl );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator-( double o )
{
    gslobj_matrix r( *this );
    gsl_matrix_add_constant( r._gsl, -o );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator-( double* o )
{
    gslobj_matrix r( *this );
    gsl_matrix_view vw = gsl_matrix_view_array( o,
						_gsl->size1,
						_gsl->size2 );
    gsl_matrix_sub( r._gsl, &( vw.matrix ) );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator-( double* o[] )
{
    gslobj_matrix r( *this );
    gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
						_gsl->size1,
						_gsl->size2 );
    gsl_matrix_sub( r._gsl, &( vw.matrix ) );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator-( gsl_matrix* o )
{
    gslobj_matrix r( *this );
    gsl_matrix_sub( r._gsl, o );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator-=( const gslobj_matrix& o )
{
    gsl_matrix_sub( _gsl, o._gsl );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator-=( double o )
{
    gsl_matrix_add_constant( _gsl, -o );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator-=( double* o )
{
    gsl_matrix_view vw = gsl_matrix_view_array( o,
						_gsl->size1,
						_gsl->size2 );
    gsl_matrix_sub( _gsl, &( vw.matrix ) );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator-=( double* o[] )
{
    gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
						_gsl->size1,
						_gsl->size2 );
    gsl_matrix_sub( _gsl, &( vw.matrix ) );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator-=( gsl_matrix* o )
{
    gsl_matrix_sub( _gsl, o );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator*( const gslobj_matrix& o )
{
    gslobj_matrix r( *this );
    gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
		    1.0, r._gsl, o._gsl,
		    0.0, _gsl );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator*( double o )
{
    gslobj_matrix r( *this );
    gsl_matrix_scale( r._gsl, o );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator*( double* o )
{
    gslobj_matrix r( *this );
    gsl_matrix_view vw = gsl_matrix_view_array( o,
						_gsl->size1,
						_gsl->size2 );
    gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
		    1.0, r._gsl, &( vw.matrix ),
		    0.0, _gsl );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator*( double* o[] )
{
    gslobj_matrix r( *this );
    gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
						_gsl->size1,
						_gsl->size2 );
    gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
		    1.0, r._gsl, &( vw.matrix ),
		    0.0, _gsl );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_matrix gslobj_matrix::operator*( gsl_matrix* o )
{
    gslobj_matrix r( *this );
    gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
		    1.0, r._gsl, o,
		    0.0, _gsl );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_vector gslobj_matrix::operator*( const gslobj_vector& o )
{
    gslobj_vector ret( rows( ) );
    gsl_blas_dgemv( CblasNoTrans, 1.0, _gsl,
		    ( gsl_vector* )o, 0.0,
		    ( gsl_vector* )ret );
    return( ret );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator*=( const gslobj_matrix& o )
{
    gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
		    1.0, _gsl, o._gsl,
		    0.0, _gsl );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator*=( double o )
{
    gsl_matrix_scale( _gsl, o );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator*=( double* o )
{
    gsl_matrix_view vw = gsl_matrix_view_array( o,
						_gsl->size1,
						_gsl->size2 );
    gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
		    1.0, _gsl, &( vw.matrix ),
		    0.0, _gsl );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator*=( double* o[] )
{
    gsl_matrix_view vw = gsl_matrix_view_array( o[ 0 ],
						_gsl->size1,
						_gsl->size2 );
    gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
		    1.0, _gsl, &( vw.matrix ),
		    0.0, _gsl );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_matrix& gslobj_matrix::operator*=( gsl_matrix* o )
{
    gsl_blas_dgemm( CblasNoTrans, CblasNoTrans,
		    1.0, _gsl, o,
		    0.0, _gsl );
    return( *this );
}

// ---------------------------------------------------------------------
/*==PS========================================
gslobj_matrix gslobj_matrix::transpose( )
{
    gslobj_matrix r( _gsl->size2, _gsl->size1 ); 
    gsl_matrix_transpose_memcpy( r._gsl, _gsl );
    return( r );
}
==PS========================================*/
// ---------------------------------------------------------------------
/*==PS========================================
gslobj_matrix gslobj_matrix::invertLU( )
{
    gslobj_matrix r( *this ), lu( *this );
    gsl_permutation* perm = gsl_permutation_calloc( rows( ) );
    int sign;

    gsl_linalg_LU_decomp( lu._gsl, perm, &sign );
    gsl_linalg_LU_invert( lu._gsl, perm, r._gsl );
    return( r );
}
==PS========================================*/
// ---------------------------------------------------------------------
/*==PS========================================
double gslobj_matrix::determinant( )
{
    gslobj_matrix lu( *this );
    gsl_permutation* perm = gsl_permutation_calloc( rows( ) );
    int sign;

    gsl_linalg_LU_decomp( lu._gsl, perm, &sign );
    return( gsl_linalg_LU_det( lu._gsl, sign ) );
}
==PS========================================*/
// ---------------------------------------------------------------------
/*==PS========================================
void gslobj_matrix::identity( )
{
    gsl_matrix_set_identity( _gsl );
}
==PS========================================*/
// eof - gslobj_matrix.cxx