gsl_matrix_double.h

Go to the documentation of this file.

//

//  This matrix class has been adapted from the gslwrap project.
//  It 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.

#ifndef _GSL_MATRIX_DOUBLE_H
#define _GSL_MATRIX_DOUBLE_H

#include <iostream>
#include <fstream>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h> 
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include "gsl_permutation.h"
#include "gsl_vector_double.h"

#define type_is_int
#ifdef  type_is
#define type_is_double
#endif

//

class matrix 
{

 friend class matrix_int;

 private: 
 
//

 gsl_matrix *m;
 
 public:
 
 typedef double value_type;
 typedef vector vector_type;

//---------------------------------------------------------------------------

//
//

  matrix();
  
//
//
  
  matrix( size_t new_rows, size_t new_cols, bool clear = true);
  
//
//

  ~matrix();

//
//
        
  void set_element( size_t row, size_t col, const double &v );
  
//
//
        
  double get_element( size_t row, size_t col ) const;
  
//
//
        
  void dimensions( size_t *num_rows, size_t *num_cols ) const;
  
//
//
        
  void set_dimensions( size_t new_rows, size_t new_cols );

//
//
        
  const double &operator()( size_t row, size_t col ) const;
  
//
//
        
  double &operator()( size_t row, size_t col );
  
//
//

  template<class oclass>
  void copy(const oclass &other)
  {
                if ( static_cast<const void *>( this ) == static_cast<const void *>( &other ) )
                        return;

                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 submatrix(size_t row_min, size_t row_max, size_t col_min, size_t col_max) const;         

// 

    void set_zero();

//              

        size_t get_rows() const; 

//

        size_t get_cols() const; 

//      
        size_t size1() const; 
        
//
        size_t size2() const; 

//      

        int fwrite (FILE * stream) const;

// 

        int fread (FILE * stream);

// 
//

        int fprintf(FILE * stream, const char * format) const;

// 
//

        int fscanf (FILE * stream);
        
//

        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;

//

        vector_view       diagonal();

//

        const vector_view diagonal() const;

//
//

        matrix &operator=( const matrix& other );

//
//
        
        template<class omatrix>
        matrix &operator=( const omatrix& other )
        {
                        copy(other);
                        return *this;
        }
        
//
//      
        
        bool operator==( const matrix &other ) const;

//
//
        
        bool operator!=( const matrix &other ) const {return !((*this)==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 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 &);
        
//
//      
        
        void diag(const vector& v);
        
//
//

        void set_tridiag(const vector& v1, const vector& v2, const vector& v3);

//
        
        matrix transpose() const;
        
//

        matrix LU_decomp(permutation *perm=NULL,int *psign=NULL) const;

//
        
        matrix LU_invert() const;

//
        
        double sum() const;

//
//

        double row_sum( size_t rowindex ) const;

//
//
        
        double col_sum( size_t colindex ) const;

//
        
        double LU_lndet() const;

//
        
        int cholesky_decomp( matrix &a ) const;
         
//      

        gsl_matrix       *gslobj()       {assert(m);return m;} 
        
//

        const gsl_matrix *gslobj() const {assert(m);return m;}
        
};

#endif // GSL_MATRIX_DOUBLE_H