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

  Program:   wxMaracas
  Module:    $RCSfile: vector.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>

// ---------------------------------------------------------------------
std::ostream& operator<<( std::ostream& os, const gslobj_vector& v )
{
    for( int i = 0; i < v.size( ); i++ ) os << " " << v( i );
    return( os );
}

// ---------------------------------------------------------------------
gslobj_vector::gslobj_vector( size_t s )
    : _copy( true )
{
    _gsl = gsl_vector_alloc( s );
}

// ---------------------------------------------------------------------
gslobj_vector::gslobj_vector( const gslobj_vector& v )
    : _copy( true )
{
    _gsl = gsl_vector_alloc( v._gsl->size );
    gsl_vector_memcpy( _gsl, v._gsl );
}

// ---------------------------------------------------------------------
gslobj_vector::gslobj_vector( double* v, size_t s )
    : _copy( false )
{
    _view = gsl_vector_view_array( v, s );
    _gsl = &( _view.vector );
}

// ---------------------------------------------------------------------
gslobj_vector::gslobj_vector( gsl_vector* v )
    : _copy( false )
{
    _gsl = v;
}

// ---------------------------------------------------------------------
/*==PS========================================
void gslobj_vector::resize( size_t s )
{
    if( _copy ) {

	if( s != size( ) )
	    gsl_vector_free( _gsl );
	_gsl = gsl_vector_alloc( s );

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

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

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

// ---------------------------------------------------------------------
gslobj_vector& gslobj_vector::operator=( gsl_vector* o )
{
    gsl_vector_memcpy( _gsl, o );
    return( *this );
}

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

    if( _gsl->size != o._gsl->size )
	return( false );

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

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

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

// ---------------------------------------------------------------------
gslobj_vector gslobj_vector::operator+( double* o )
{
    gslobj_vector r( *this );
    gsl_vector_view vw = gsl_vector_view_array( o, _gsl->size );
    gsl_vector_add( r._gsl, &( vw.vector ) );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_vector gslobj_vector::operator+( gsl_vector* o )
{
    gslobj_vector r( *this );
    gsl_vector_add( r._gsl, o );
    return( r );
}

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

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

// ---------------------------------------------------------------------
gslobj_vector& gslobj_vector::operator+=( double* o )
{
    gsl_vector_view vw = gsl_vector_view_array( o, _gsl->size );
    gsl_vector_add( _gsl, &( vw.vector ) );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_vector& gslobj_vector::operator+=( gsl_vector* o )
{
    gsl_vector_add( _gsl, o );
    return( *this );
}

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

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

// ---------------------------------------------------------------------
gslobj_vector gslobj_vector::operator-( double* o )
{
    gslobj_vector r( *this );
    gsl_vector_view vw = gsl_vector_view_array( o, _gsl->size );
    gsl_vector_sub( r._gsl, &( vw.vector ) );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_vector gslobj_vector::operator-( gsl_vector* o )
{
    gslobj_vector r( *this );
    gsl_vector_sub( r._gsl, o );
    return( r );
}

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

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

// ---------------------------------------------------------------------
gslobj_vector& gslobj_vector::operator-=( double* o )
{
    gsl_vector_view vw = gsl_vector_view_array( o, _gsl->size );
    gsl_vector_sub( _gsl, &( vw.vector ) );
    return( *this );
}

// ---------------------------------------------------------------------
gslobj_vector& gslobj_vector::operator-=( gsl_vector* o )
{
    gsl_vector_sub( _gsl, o );
    return( *this );
}

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

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

// ---------------------------------------------------------------------
double gslobj_vector::dot( const gslobj_vector& o )
{
    double r;
    gsl_blas_ddot( _gsl, o._gsl, &r );
    return( r );
}

// ---------------------------------------------------------------------
double gslobj_vector::dot( double* o )
{
    double r;
    gsl_vector_view vw = gsl_vector_view_array( o, _gsl->size );
    gsl_blas_ddot( _gsl, &( vw.vector ), &r );
    return( r );
}

// ---------------------------------------------------------------------
double gslobj_vector::dot( gsl_vector* o )
{
    double r;
    gsl_blas_ddot( _gsl, o, &r );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_vector gslobj_vector::cross( const gslobj_vector& o )
{
    return( cross( o._gsl->data ) );
}

// ---------------------------------------------------------------------
gslobj_vector gslobj_vector::cross( double* o )
{
    gslobj_vector r( *this );
    double S, s;
	
    // 1st element
    S = _gsl->data[ 1 ] * o[ 2 ];
    s = _gsl->data[ 2 ] * o[ 1 ];
    r._gsl->data[ 0 ] = S - s;
		
    // 2nd element
    S = _gsl->data[ 2 ] * o[ 0 ];
    s = _gsl->data[ 0 ] * o[ 2 ];
    r._gsl->data[ 1 ] = S - s;
		
    // 3rd element
    S = _gsl->data[ 0 ] * o[ 1 ];
    s = _gsl->data[ 1 ] * o[ 0 ];
    r._gsl->data[ 2 ] = S - s;
    return( r );
}

// ---------------------------------------------------------------------
gslobj_vector gslobj_vector::cross( gsl_vector* o )
{
    return( cross( o->data ) );
}

// ---------------------------------------------------------------------
int gslobj_vector::scross( const gslobj_vector& o )
{
    *this = cross( o );
    return( 1 );
}

// ---------------------------------------------------------------------
int gslobj_vector::scross( double* o )
{
    *this = cross( o );
    return( 1 );
}

// ---------------------------------------------------------------------
int gslobj_vector::scross( gsl_vector* o )
{
    *this = cross( o );
    return( 1 );
}

// ---------------------------------------------------------------------
gslobj_vector gslobj_vector::operator/( double o )
{
    gslobj_vector r( *this );
    gsl_vector_scale( r._gsl, 1.0 / o );
    return( r );
}

// ---------------------------------------------------------------------
gslobj_vector& gslobj_vector::operator/=( double o )
{
    gsl_vector_scale( _gsl, 1.0 / o );
    return( *this );
}

// ---------------------------------------------------------------------
/*==PS========================================
double gslobj_vector::norm1( )
{
    return( gsl_blas_dasum( _gsl ) );
}
==PS========================================*/
// ---------------------------------------------------------------------
double gslobj_vector::norm2( )
{
    return( gsl_blas_dnrm2( _gsl ) );
}

// ---------------------------------------------------------------------
gslobj_vector gslobj_vector::normalize( )
{
    gslobj_vector r( *this );
    gsl_vector_scale( r._gsl, 1.0 / gsl_blas_dnrm2( _gsl ) );
    return( r );
}

// ---------------------------------------------------------------------
int gslobj_vector::snormalize( )
{
    return( gsl_vector_scale( _gsl, 1.0 / gsl_blas_dnrm2( _gsl ) ) );
}

// eof - gslobj_vector.cxx