---

matrix_double.h

Go to the documentation of this file.

// matrix.h

//  This matrix class is a C++ wrapper for the GNU Scientific Library
//  Copyright (C) 2001 Ramin Nakisa

//  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., 675 Mass Ave, Cambridge, MA 02139, USA.

#if !defined( _matrix_double_h )
#define _matrix_double_h

#include <iostream.h>
#include <fstream.h>
#include <iomanip.h>
#include <math.h>
#include <stdlib.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gslwrap/permutation.h>
#include <gslwrap/vector_double.h>

#define type_is
#ifdef  type_is
#define type_is_double
#endif

namespace gsl
{

class matrix
{
#ifdef type_is_double
        friend class matrix_float;
        friend class matrix_int;
#endif
public:
        matrix();
        matrix( size_t new_rows, size_t new_cols, bool clear = true );

        void copy(const matrix &other)
        {
                set_dimensions(other.size1(),other.size2());
                gsl_matrix_memcpy( m, other.m );
        }
        template<class oclass>
        void copy(const oclass &other)
        {
                set_dimensions( other.get_rows(), other.get_cols() );
                for ( size_t i = 0; i < get_rows(); i++ ) 
                {
                        for ( size_t j = 0; j < get_cols(); j++ ) 
                        {
                                gsl_matrix_set( m, i, j, (double)other(i,j));
                        }
                }
        }
        matrix( const matrix &other ):m(NULL) {copy(other);}
        template<class oclass>
        matrix( const oclass &other ):m(NULL) {copy(other);}

        ~matrix();
//      matrix( const char *Filename );
        size_t get_rows() const {return m->size1;}
        size_t get_cols() const {return m->size2;}
        size_t size1() const {return m->size1;}
        size_t size2() const {return m->size2;}
   

        void dimensions( size_t *num_rows, size_t *num_cols ) const;
double        get_element ( size_t row, size_t col ) const {return  gsl_matrix_get( m, row, col ) ;}
        const double &operator()( size_t row, size_t col ) const {return *gsl_matrix_ptr( m, row, col ) ;}
double       &operator()( size_t row, size_t col )       {return *gsl_matrix_ptr( m, row, col ) ;}
        void set_element( size_t row, size_t col, const double &v ){ gsl_matrix_set( m, row, col, v );}
        void set_elements( const double & new_value );
        void set_all ( const double & new_value ) {gsl_matrix_set_all ( m, new_value );}
        void set_zero() {gsl_matrix_set_zero( m );}
        void set_dimensions( size_t new_rows, size_t new_cols );
        void load( const char *filename );
        void save( const char *filename ) const;
        friend ostream& operator<< ( ostream& os, const matrix& m );
        //This function writes the elements of the matrix m to the stream stream in binary format. The return value is 0 for success and GSL_EFAILED if there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.
        int fwrite (FILE * stream) const {return gsl_matrix_fwrite (stream, m);}

//This function reads into the matrix m from the open stream stream in binary format. The matrix m must be preallocated with the correct dimensions since the function uses the size of m to determine how many bytes to read. The return value is 0 for success and GSL_EFAILED if there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture. 
        int fread (FILE * stream) {return gsl_matrix_fread (stream, m);}

        void load_binary( const char *filename );
        void save_binary( const char *filename ) const;
        bool operator==( const matrix &other ) const;
        bool operator!=( const matrix &other ) const {return !((*this)==other);}

        matrix &operator=( const matrix &other ){copy(other);return *this;}
        template<class omatrix>
        matrix &operator=( const omatrix &other ){copy(other);return *this;}
        matrix operator+( const matrix &other ) const;
        matrix operator+( const double &f ) const;
        friend matrix operator+( const double &f, const matrix &other );
        matrix &operator+=( const double &f );
        matrix &operator+=( const matrix &other );
        matrix operator-( const matrix &other ) const;
        matrix operator-( const double &f ) const;
        friend matrix operator-( const double &f, const matrix &other );
        matrix &operator-=( const double &f );
        matrix &operator-=( const matrix &other );
        matrix operator*( const matrix &other ) const;
        matrix operator*( const double &f ) const;
        friend matrix operator*( const double &f, const matrix &other );
        matrix &operator*=( const double &f );
        matrix &operator*=( const matrix &other );
        matrix operator/( const double &) const;
        matrix &operator/=( const double &);
        matrix transpose() const;
        matrix LU_decomp(gsl::permutation *perm=NULL,int *psign=NULL) const;
        matrix LU_invert() const;
private:
        void LU_decomp( gsl_matrix **a,
                                        gsl_permutation **permutation,
                                        int *sign ) const;
public:
    double sum() const;
        double LU_lndet() const;


        vector_view       row( size_t rowindex );
        const vector_view row( size_t rowindex ) const ;
        vector_view       column( size_t colindex );
        const vector_view column( size_t colindex ) const;

        matrix get_row( size_t rowindex ) const;
        matrix get_col( size_t colindex ) const;
        matrix row_sum() const;
        matrix column_sum() const;
        double trace() const;
        int cholesky_decomp( matrix &a ) const;
//      /** returns index of nearest row in matrix to vector argument. */
//      int nearest_row_index( const matrix &v ) const;
        matrix covariance() const;
        bool is_square() const;
        void set_diagonal( double f );
        void identity( size_t k );
        double norm( double n ) const;

        gsl_matrix       *gslobj()       {if (!m){cout << "matrix::gslobj ERROR, data not initialized!! " << endl; exit(-1);}return m;}
        const gsl_matrix *gslobj() const {if (!m){cout << "matrix::gslobj ERROR, data not initialized!! " << endl; exit(-1);}return m;}
 private:
   gsl_matrix *m;

};
}
#undef type_is
#undef type_is_double

#endif // _matrix_double_h

---