Ignore:
Timestamp:
Dec 2, 2005, 1:15:50 AM (18 years ago)
Author:
Jari Häkkinen
Message:

Merged better_matrix_class changes r402:419 into the trunk.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/lib/gslapi/matrix.h

    r390 r420  
    11// $Id$
    22
    3 #ifndef _theplu_gslapi_matrix_ 
    4 #define _theplu_gslapi_matrix_ 
     3#ifndef _theplu_gslapi_matrix_
     4#define _theplu_gslapi_matrix_
    55
    66#include <c++_tools/gslapi/vector.h>
     7#include <c++_tools/utility/Exception.h>
    78
    89#include <gsl/gsl_matrix.h>
    910#include <iostream>
    10 
    1111
    1212namespace theplu {
    … …  
    2020  /// well in the future.
    2121  ///
    22   class matrix
    23   {
    24   public:
     22  /// \par[File streams] Reading and writing vectors to file streams
     23  /// are of course supported. These are implemented without using GSL
     24  /// functionality, and thus binary read and write to streams are not
     25  /// supported.
     26  ///
     27  /// \par[Matrix views] GSL matrix views are supported and these are
     28  /// disguised as ordinary gslapi::matrix'es. A support function is
     29  /// added, gslapi::matrix::isview(), that can be used to check if a
     30  /// matrix object is a view. Note that view matricess do not own the
     31  /// undelying data, and a view is not valid if the matrix owning the
     32  /// data is deallocated.
     33  ///
     34  /// @note All GSL matrix related functions are not implement but
     35  /// most functionality defined for GSL matrices can be achieved with
     36  /// this interface class. Most notable GSL functionality not
     37  /// supported are no binary file support and views on arrays,
     38  /// gslapi::vectors, gsl_vectors, diagonals, subdiagonals, and
     39  /// superdiagonals. If there is an interest from the user community,
     40  /// the omitted functionality could be included.
     41  ///
     42  /// @todo Maybe it would be smart to create temporary objects need
     43  /// for BLAS in constructors. As it is now temporary objects are
     44  /// called before BLAS functionality iss used, cf. const matrix&
     45  /// matrix::operator*=(const matrix& other)
     46  ///
     47  class matrix
     48  {
     49  public:
    2550
    2651    ///
    2752    /// The default constructor.
    28     ///
    29     matrix(void);
     53    ///
     54    /// This contructor does not initialize underlying (essential)
     55    /// structures.
     56    ///
     57    inline matrix(void) : m_(NULL), view_(NULL) {}
    3058
    3159    ///
    3260    /// Constructor. Allocates memory space for \a r times \a c
    3361    /// elements, and sets all elements to \a init_value.
    34     ///
    35     matrix(const size_t& r, const size_t& c, const double init_value=0);
    36  
    37     ///
     62    ///
     63    inline matrix(const size_t& r, const size_t& c, double init_value=0)
     64      : view_(NULL) { m_ = gsl_matrix_alloc(r,c); set_all(init_value); }
     65
     66    ///
    3867    /// The copy constructor.
    39     ///
    40     matrix(const matrix&);
     68    ///
     69    /// @note If the object to be copied is a matrix view, the values
     70    /// of the view will be copied, i.e. the view is not copied.
     71    ///
     72    inline matrix(const matrix& o) : view_(NULL) {m_=o.create_gsl_matrix_copy();}
     73
     74    ///
     75    /// The matrix view constructor.
     76    ///
     77    /// Create a view of matrix \a m, with starting row \a offset_row,
     78    /// starting column \a offset_column, row size \a n_row, and
     79    /// column size \a n_column.
     80    ///
     81    /// A matrix view can be used as any matrix with the difference
     82    /// that changes made to the view will also change the object that
     83    /// is viewed. Also, using the copy constructor will create a new
     84    /// matrix object that is a copy of whatever is viewed. If a copy
     85    /// of the view is needed then you should use this constructor to
     86    /// obtain a copy.
     87    ///
     88    /// @note If the object viewed by the view goes out of scope or is
     89    /// deleted, the view becomes invalid and the result of further
     90    /// use is undefined.
     91    ///
     92    matrix(matrix& m, size_t offset_row, size_t offset_column, size_t n_row,
     93           size_t n_column);
    4194
    4295    ///
    … …  
    4699    /// column elements in a row must be separated with white space
    47100    /// characters except the new line (\\n) character. Rows, and
    48     /// columns, are read consequently and only rectangular matrices
     101    /// columns, are read consecutively and only rectangular matrices
    49102    /// are supported. Empty lines are ignored. End of input is at end
    50103    /// of file marker.
    51104    ///
    52     explicit matrix(std::istream &);
    53 
    54     ///
     105    explicit matrix(std::istream &) throw (utility::IO_error,std::exception);
     106
     107    ///
    55108    /// The destructor.
    56109    ///
    57     ~matrix(void);
    58 
    59     ///
     110    ~matrix(void);
     111
     112    ///
     113    /// Elementwise addition of the elements of matrix \a b to the
     114    /// elements of the calling matrix ,\f$a_{ij} = a_{ij} + b_{ij} \;
     115    /// \forall i,j\f$. The result is stored into the calling matrix.
     116    ///
     117    /// @return Whatever GSL returns.
     118    ///
     119    inline int add(const matrix& b) { return gsl_matrix_add(m_,b.m_); }
     120
     121    ///
     122    /// Add the scalar value \a d to the elements of the calling
     123    /// matrix, \f$a_{ij} = a_{ij} + d \; \forall i,j\f$. The result
     124    /// is stored into the calling matrix.
     125    ///
     126    /// @return Whatever GSL returns.
     127    ///
     128    inline int
     129    add_constant(const double d) { return gsl_matrix_add_constant(m_,d); }
     130
     131    ///
    60132    /// @return The number of columns in the matrix.
    61133    ///
    62     inline size_t columns(void) const { return m_->size2; }
    63 
    64     ///
    65     /// Swap column \a i with \a j.
    66     ///
    67     inline void column_swap(const size_t& i,const size_t& j)
    68       { gsl_matrix_swap_columns(m_,i,j); }
    69 
    70     ///
    71     /// @return true if each element deviates less or equal than \a
    72     /// d. If any matrix contain a nan false is always returned.
    73     ///
    74     bool equal(const matrix&, const double d=0) const;
     134    inline size_t columns(void) const { return m_->size2; }
     135
     136    ///
     137    /// Elementwise division of the elemnts of the calling matrix by
     138    /// the elements of matrix \a b, \f$a_{ij} = a_{ij} / b_{ij} \;
     139    /// \forall i,j\f$. The result is stored into the calling matrix.
     140    ///
     141    /// @return Whatever GSL returns.
     142    ///
     143    inline int
     144    div_elements(const matrix& b) { return gsl_matrix_div_elements(m_,b.m_); }
     145
     146    ///
     147    /// Check whether matrices are equal within a user defined
     148    /// precision, set by \a precision.
     149    ///
     150    /// @return True if each element deviates less or equal than \a
     151    /// d. If any matrix contain a NaN, false is always returned.
     152    ///
     153    bool equal(const matrix&, const double precision=0) const;
     154
     155    ///
     156    /// @return A const pointer to the internal GSL matrix.
     157    ///
     158    inline const gsl_matrix* gsl_matrix_p(void) const { return m_; }
     159
     160    ///
     161    /// @return A pointer to the internal GSL matrix.
     162    ///
     163    // Jari, is this needed?
     164    inline gsl_matrix* gsl_matrix_p(void) { return m_; };
     165
     166    ///
     167    /// @return True if all elements in the matrix is zero, false
     168    /// othwerwise;
     169    ///
     170    inline bool isnull(void) const { return gsl_matrix_isnull(m_); }
     171
     172    ///
     173    /// Check if the matrix object is a view (sub-matrix) to another
     174    /// matrix.
     175    ///
     176    /// @return True if the object is a view, false othwerwise.
     177    ///
     178    inline bool isview(void) const { return view_; }
     179
     180    ///
     181    /// @return The maximum value of the matrix.
     182    ///
     183    inline double max(void) const { return gsl_matrix_max(m_); }
     184
     185    ///
     186    /// @return The index to the element with the maximum value of the
     187    /// matrix. The lowest index has precedence (searching in
     188    /// row-major order).
     189    ///
     190    inline std::pair<size_t,size_t> max_index(void) const;
     191
     192    ///
     193    /// @return The minimum value of the matrix.
     194    ///
     195    inline double min(void) const { return gsl_matrix_min(m_); }
     196
     197    ///
     198    /// @return The index to the element with the minimum value of the
     199    /// matrix. The lowest index has precedence (searching in
     200    /// row-major order).
     201    ///
     202    inline std::pair<size_t,size_t> min_index(void) const;
     203
     204    ///
     205    /// @return The indecies to the element with the minimum and
     206    /// maximum values of the matrix, respectively. The lowest index
     207    /// has precedence (searching in row-major order).
     208    ///
     209    inline void minmax_index(std::pair<size_t,size_t>& min,
     210                             std::pair<size_t,size_t>& max) const
     211      { gsl_matrix_minmax_index(m_,&min.first,&min.second,
     212                                &max.first,&max.second); }
     213
     214    ///
     215    /// @return The minimum and maximum values of the matrix, as the
     216    /// \a first and \a second member of the returned \a pair,
     217    /// respectively.
     218    ///
     219    std::pair<double,double> minmax(void) const;
     220
     221    ///
     222    /// Multiply the elements of matrix \a b with the elements of the
     223    /// calling matrix ,\f$a_{ij} = a_{ij} * b_{ij} \; \forall
     224    /// i,j\f$. The result is stored into the calling matrix.
     225    ///
     226    /// @return Whatever GSL returns.
     227    ///
     228    inline int
     229    mul_elements(const matrix& b) { return gsl_matrix_mul_elements(m_,b.m_); }
     230
     231    ///
     232    /// @return The number of rows in the matrix.
     233    ///
     234    inline size_t rows(void) const { return m_->size1; }
     235
     236    ///
     237    /// Multiply the elements of the calling matrix with a scalar \a
     238    /// d, \f$a_{ij} = d * a_{ij} \; \forall i,j\f$. The result is
     239    /// stored into the calling matrix.
     240    ///
     241    /// @return Whatever GSL returns.
     242    ///
     243    inline int scale(const double d) { return gsl_matrix_scale(m_,d); }
     244
     245    ///
     246    /// Set element values to values in \a mat. This function is
     247    /// needed for assignment of viewed elements.
     248    ///
     249    /// @return Whatever GSL returns.
     250    ///
     251    /// @note No check on size is done.
     252    ///
     253    /// @see const matrix& operator=(const matrix&)
     254    ///
     255    inline int set(const matrix& mat) { return gsl_matrix_memcpy(m_,mat.m_); }
     256
     257    ///
     258    /// Set all elements to \a value.
     259    ///
     260    inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }
     261
     262    ///
     263    /// Set \a column values to values in \a vec.
     264    ///
     265    /// @return Whatever GSL returns.
     266    ///
     267    /// @note No check on size is done.
     268    ///
     269    inline int set_column(const size_t column, const vector& vec)
     270      { return gsl_matrix_set_col(m_, column, vec.gsl_vector_p()); }
     271
     272    ///
     273    /// Set \a row values to values in \a vec.
     274    ///
     275    /// @return Whatever GSL returns.
     276    ///
     277    /// @note No check on size is done.
     278    ///
     279    inline int set_row(const size_t row, const vector& vec)
     280      { return gsl_matrix_set_row(m_, row, vec.gsl_vector_p()); }
     281
     282    ///
     283    /// Subtract the elements of matrix \a b from the elements of the
     284    /// calling matrix ,\f$a_{ij} = a_{ij} - b_{ij} \; \forall
     285    /// i,j\f$. The result is stored into the calling matrix.
     286    ///
     287    /// @return Whatever GSL returns.
     288    ///
     289    inline int sub(const matrix& b) { return gsl_matrix_sub(m_,b.m_); }
     290
     291    ///
     292    /// Exchange the elements of the this and \a other by copying. The
     293    /// two matrices must have the same size.
     294    ///
     295    /// @return Whatever GSL returns.
     296    ///
     297    inline int swap(matrix& other) { return gsl_matrix_swap(m_,other.m_); }
     298
     299    ///
     300    /// Swap columns \a i and \a j.
     301    ///
     302    inline int swap_columns(const size_t i,const size_t j)
     303      { return gsl_matrix_swap_columns(m_,i,j); }
     304
     305    ///
     306    /// Swap row \a i and column \a j.
     307    ///
     308    inline int swap_rowcol(const size_t i,const size_t j)
     309      { return gsl_matrix_swap_rowcol(m_,i,j); }
     310
     311    ///
     312    /// Swap rows \a i and \a j.
     313    ///
     314    inline int swap_rows(const size_t i, const size_t j)
     315      { return gsl_matrix_swap_rows(m_,i,j); }
     316
     317    ///
     318    /// Transpose the matrix.
     319    ///
     320    void transpose(void);
     321
     322    ///
     323    /// @return Reference to the element position (\a row, \a column).
     324    ///
     325    inline double& operator()(size_t row,size_t column)
     326      { return (*gsl_matrix_ptr(m_,row,column)); }
     327
     328    ///
     329    /// @return Const reference to the element position (\a row, \a
     330    /// column).
     331    ///
     332    inline const double& operator()(size_t row,size_t column) const
     333      { return (*gsl_matrix_const_ptr(m_,row,column)); }
     334
     335    ///
     336    /// Matrix-vector multiplication.
     337    ///
     338    /// @return The resulting vector.
     339    ///
     340    // Jari, where should this go?
     341    const vector operator*(const vector&) const;
     342
     343    ///
     344    /// Comparison operator.
     345    ///
     346    /// @return True if all elements are equal otherwise False.
     347    ///
     348    /// @see equal
     349    ///
     350    inline bool operator==(const matrix& other) const { return equal(other); }
     351
     352    ///
     353    /// Comparison operator.
     354    ///
     355    /// @return False if all elements are equal otherwise True.
     356    ///
     357    /// @see equal
     358    ///
     359    inline bool operator!=(const matrix& other) const { return !equal(other); }
     360
     361    ///
     362    /// The assignment operator. There is no requirements on
     363    /// dimensions, i.e. the matrix is remapped in memory if
     364    /// necessary. This implies that in general views cannot be
     365    /// assigned using this operator. Views will be mutated into
     366    /// normal matrices. The only exception to this behaviour on views
     367    /// is when self-assignemnt is done, since self-assignment is
     368    /// ignored.
     369    ///
     370    /// @return A const reference to the resulting matrix.
     371    ///
     372    /// @see int set(const matrix&)
     373    ///
     374    const matrix& operator=(const matrix& other);
     375
     376    ///
     377    /// Add and assign operator.
     378    ///
     379    inline const matrix& operator+=(const matrix& m) { add(m); return *this; }
     380
     381    ///
     382    /// Add and assign operator.
     383    ///
     384    inline const matrix&
     385    operator+=(const double d) { add_constant(d); return *this; }
     386
     387    ///
     388    /// Subtract and assign operator.
     389    ///
     390    inline const matrix& operator-=(const matrix& m) { sub(m); return *this; }
     391
     392    ///
     393    /// Multiply and assigment operator.
     394    ///
     395    /// @return Const reference to the resulting matrix.
     396    ///
     397    const matrix& operator*=(const matrix&);
     398
     399    ///
     400    /// Multiply and assign operator.
     401    ///
     402    inline const matrix& operator*=(const double d) { scale(d); return *this; }
     403
     404
     405  private:
    75406
    76407    ///
    … …  
    82413    /// @return A pointer to a copy of the internal GSL matrix.
    83414    ///
    84     gsl_matrix* gsl_matrix_copy(void) const;
    85 
    86     ///
    87     /// @return A const pointer to the internal GSL matrix.
    88     ///
    89     inline const gsl_matrix* gsl_matrix_pointer(void) const { return m_; }
    90 
    91     ///
    92     /// @return A pointer to the internal GSL matrix.
    93     ///
    94     inline gsl_matrix* gsl_matrix_pointer(void) { return m_; };
    95 
    96     ///
    97     /// This function multiplies all elements between two matrices and
    98     /// the result is stored back into the calling object, \f$a_{ij} =
    99     /// a_{ij} * b_{ij} \; \forall i,j\f$.
    100     ///
    101     inline void mul_elements(const matrix& b)
    102       { gsl_matrix_mul_elements(m_,b.m_); }
    103 
    104     ///
    105     /// Swap row \a i with \a j.
    106     ///
    107     inline void row_swap( const size_t& i, const size_t& j)
    108       { gsl_matrix_swap_rows(m_,i,j); }
    109 
    110     ///
    111     /// @return The number of rows in the matrix.
    112     ///
    113     inline size_t rows(void) const { return m_->size1; }
    114 
    115     ///
    116     /// Set all elements to \a value.
    117     ///
    118     inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }
    119 
    120     ///
    121     /// Set column \a i to \a vec.
    122     ///
    123     inline void set_column(const size_t i, const vector vec)
    124     {gsl_matrix_set_col(m_, i, vec.gsl_vector_pointer());}
    125 
    126     ///
    127     /// Set row \a i to \a vec.
    128     ///
    129     inline void set_row(const size_t i, const vector vec)
    130     {gsl_matrix_set_row(m_, i, vec.gsl_vector_pointer());}
    131 
    132     ///
    133     /// Set column \a i to \a vec
    134     ///
    135     //void matrix set_column(const size_t i, const vector vec) const;
    136 
    137     ///
    138     /// Transpose the matrix.
    139     ///
    140     /// @return Transpose of the matrix.
    141     ///
    142     matrix transpose(void) const;
    143    
    144     ///
    145     /// @return Reference to the element position (\a row, \a column).
    146     ///
    147     inline double& operator()(size_t row,size_t column)
    148       { return (*gsl_matrix_ptr(m_,row,column)); }
    149 
    150     ///
    151     /// @return Const reference to the element position (\a row, \a
    152     /// column).
    153     ///
    154     inline const double& operator()(size_t row,size_t column) const
    155       { return (*gsl_matrix_const_ptr(m_,row,column)); }
    156 
    157     ///
    158     /// @todo to be properly implemented. Should return a reference to a vector
    159     /// object. Underlying GSL supports only GSL_vector views, thus
    160     /// some redesign of the vector class is needed.
    161     ///
    162     /// inline vector& operator[](size_t i);
    163     /// So for now, stupid copying is done, and actually a matrix
    164     /// row is returned using this function.
    165     vector operator[](size_t row) const;
    166 
    167     ///
    168     /// @todo to be properly implemented. Should return a reference to
    169     /// a vector object (matrix row). Underlying GSL supports only GSL_vector
    170     /// views, thus some redesign of the vector class is needed.
    171     ///
    172     /// inline const vector& operator[](size_t i) const;
    173     ///
    174     /// So for now, stupid copying is done, and actually a matrix
    175     /// column is returned using this function.
    176     vector TEMP_col_return(size_t column) const;
    177 
    178     ///
    179     /// Matrix multiplication.
    180     ///
    181     /// @return The resulting matrix.
    182     ///
    183     matrix operator*(const matrix&) const;
    184 
    185     ///
    186     /// Matrix-vector multiplication.
    187     ///
    188     /// @return The resulting vector.
    189     ///
    190     vector operator*(const vector&) const;
    191 
    192     ///
    193     /// Matrix addition.
    194     ///
    195     /// @return The resulting matrix.
    196     ///
    197     matrix operator+(const matrix&) const;
    198 
    199     ///
    200     /// @todo Matrix addition.
    201     ///
    202     /// @return The resulting matrix.
    203     ///
    204     matrix operator+=(const matrix&);
    205 
    206     ///
    207     /// Matrix subtraction.
    208     ///
    209     /// @return The resulting matrix.
    210     ///
    211     matrix operator-(const matrix&) const;
    212 
    213     ///
    214     /// @todo Matrix subtraction.
    215     ///
    216     /// @return The resulting matrix.
    217     ///
    218     matrix operator-=(const matrix&) const;
    219 
    220     ///
    221     /// Comparison operator.
    222     ///
    223     /// @return True if all elements are equal otherwise False.
    224     ///
    225     /// @see equal
    226     ///
    227     inline bool operator==( const matrix& other ) const { return equal(other);}
    228 
    229     ///
    230     /// Comparison operator.
    231     ///
    232     /// @return False if all elements are equal otherwise True.
    233     ///
    234     /// @see equal
    235     ///
    236     inline bool operator!=( const matrix& other ) const {return !equal(other);}
    237 
    238     ///
    239     /// The assignment operator. There is no requirements on
    240     /// dimensions.
    241     ///
    242     matrix& operator=(const matrix& other);
    243 
    244     ///
    245     /// Subtract and assign operator.
    246     ///
    247     inline matrix& operator-=(const matrix& m)
    248       { gsl_matrix_sub(m_,m.m_); return *this; }
    249 
    250     ///
    251     /// Multiply and assign operator.
    252     ///
    253     inline matrix& operator*=(const double d)
    254       { gsl_matrix_scale(m_,d); return *this; }
    255 
    256 
    257   private:
    258 
    259     gsl_matrix* m_;
    260 
     415    gsl_matrix* create_gsl_matrix_copy(void) const;
     416
     417    gsl_matrix* m_;
     418    gsl_matrix_view* view_;
    261419  };
    262420
    263421  ///
    264422  /// The output operator for the matrix class.
    265   ///
    266   std::ostream& operator<< (std::ostream& s, const matrix&);
    267 
     423  ///
     424  std::ostream& operator<< (std::ostream& s, const matrix&);
    268425
    269426
    270427}} // of namespace gslapi and namespace theplu
    271428
    272 #endif 
     429#endif
Note: See TracChangeset for help on using the changeset viewer.