Changeset 403 for branches


Ignore:
Timestamp:
Nov 27, 2005, 1:26:25 AM (18 years ago)
Author:
Jari Häkkinen
Message:

Updated matrix class, but matrix view not considered yet.

Location:
branches/better_matrix_class/lib/gslapi
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/better_matrix_class/lib/gslapi/matrix.cc

    r357 r403  
    5353    std::vector<double> v;
    5454    for (nof_rows = 0; utility::read_to_double(is, v); nof_rows++) {
    55    
     55
    5656      // Ignoring empty lines
    5757      if(!v.size()){
    … …  
    5959        continue;
    6060      }
    61      
     61
    6262      if(nof_columns==0)
    6363        nof_columns=v.size();
    6464      else if(v.size()!=nof_columns) {
    6565        std::cerr << "matrix::matrix(std::istream&) data file error: "
    66                   << "line" << nof_rows+1 << " has " << v.size() 
     66                  << "line" << nof_rows+1 << " has " << v.size()
    6767                  << " columns; expected " << nof_columns
    6868                  << " columns"
    69                   << std::endl; 
     69                  << std::endl;
    7070        exit(1);
    7171      }
    … …  
    7777    // convert the data to a gsl matrix
    7878    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
    79     for(u_int i=0;i<nof_rows;i++) 
    80       for(u_int j=0;j<nof_columns;j++) 
     79    for(u_int i=0;i<nof_rows;i++)
     80      for(u_int j=0;j<nof_columns;j++)
    8181        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
    8282  }
    … …  
    9494
    9595
    96   gsl_matrix* matrix::gsl_matrix_copy(void) const
    97   {
    98     gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
    99     gsl_matrix_memcpy(m,m_);
    100     return m;
    101   }
    102 
    10396  bool matrix::equal(const matrix& other, const double d) const
    10497  {
    10598    if (columns()!=other.columns() || rows()!=other.rows())
    10699      return false;
    107     for (size_t i=0; i<rows(); i++) 
    108       for (size_t j=0; j<columns(); j++) 
     100    for (size_t i=0; i<rows(); i++)
     101      for (size_t j=0; j<columns(); j++)
    109102        // two last cond. are to check for nans
    110         if (fabs( (*this)(i,j)-other(i,j) ) > d || 
     103        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
    111104            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
    112105          return false;
    … …  
    114107    return true;
    115108  }
     109
     110
     111
     112  gsl_matrix* matrix::gsl_matrix_copy(void) const
     113  {
     114    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
     115    gsl_matrix_memcpy(m,m_);  // Jari, a GSL return value is ignored here
     116    return m;
     117  }
     118
     119
    116120
    117121  vector matrix::TEMP_col_return(size_t column) const
    … …  
    172176
    173177
    174  
     178
    175179  matrix& matrix::operator=( const matrix& other )
    176180  {
    … …  
    179183        gsl_matrix_free( m_ );
    180184      m_ = other.gsl_matrix_copy();
    181     }
     185    }
    182186    return *this;
    183187  }
    … …  
    185189
    186190
     191  // Jari, checkout GSL transpose support in GSL manual 8.4.9
    187192  matrix matrix::transpose(void) const
    188193  {
    … …  
    199204    s.setf( ios::dec );
    200205    s.precision(12);
    201     for(size_t i=0, j=0; i<m.rows(); i++) 
     206    for(size_t i=0, j=0; i<m.rows(); i++)
    202207      for (j=0; j<m.columns(); j++) {
    203208        s << m(i,j);
    … …  
    207212          s << "\n";
    208213      }
    209     return s;
     214    return s;
    210215  }
    211216

  • branches/better_matrix_class/lib/gslapi/matrix.h

    r390 r403  
    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>
    … …  
    2020  /// well in the future.
    2121  ///
    22   class matrix
    23   {
    24   public:
     22  /// @note All GSL matrix related functions are not implement but all
     23  /// functionality, except binary file support, defined for GSL
     24  /// matrices can be achieved with this interface class.
     25  ///
     26  class matrix
     27  {
     28  public:
    2529
    2630    ///
    2731    /// The default constructor.
    28     ///
    29     matrix(void);
     32    ///
     33    matrix(void);
    3034
    3135    ///
    3236    /// Constructor. Allocates memory space for \a r times \a c
    3337    /// 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     ///
     38    ///
     39    matrix(const size_t& r, const size_t& c, const double init_value=0);
     40
     41    ///
    3842    /// The copy constructor.
    39     ///
    40     matrix(const matrix&);
     43    ///
     44    /// @note If the object to be copied is a matrix view, the values
     45    /// of the view will be copied, i.e. the view is not copied.
     46    ///
     47    matrix(const matrix&);
    4148
    4249    ///
    … …  
    4653    /// column elements in a row must be separated with white space
    4754    /// characters except the new line (\\n) character. Rows, and
    48     /// columns, are read consequently and only rectangular matrices
     55    /// columns, are read consecutively and only rectangular matrices
    4956    /// are supported. Empty lines are ignored. End of input is at end
    5057    /// of file marker.
    … …  
    5259    explicit matrix(std::istream &);
    5360
    54     ///
     61    ///
    5562    /// The destructor.
    5663    ///
    57     ~matrix(void);
    58 
    59     ///
     64    ~matrix(void);
     65
     66    ///
     67    /// Elementwise addition of the elements of matrix \a b to the
     68    /// elements of the calling matrix ,\f$a_{ij} = a_{ij} + b_{ij} \;
     69    /// \forall i,j\f$. The result is stored into the calling matrix.
     70    ///
     71    /// @return Whatever GSL returns.
     72    ///
     73    inline int add(const matrix& b) { return gsl_matrix_add(m_,b.m_); }
     74
     75    ///
     76    /// Add the scalar value \a d to the elements of the calling
     77    /// matrix, \f$a_{ij} = a_{ij} + d \; \forall i,j\f$. The result
     78    /// is stored into the calling matrix.
     79    ///
     80    /// @return Whatever GSL returns.
     81    ///
     82    inline int add_constant(const double d)
     83      { return gsl_matrix_add_constant(m_,d); }
     84
     85    ///
    6086    /// @return The number of columns in the matrix.
    6187    ///
    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); }
     88    inline size_t columns(void) const { return m_->size2; }
     89
     90    ///
     91    /// Elementwise division of the elemnts of the calling matrix by
     92    /// the elements of matrix \a b, \f$a_{ij} = a_{ij} / b_{ij} \;
     93    /// \forall i,j\f$. The result is stored into the calling matrix.
     94    ///
     95    /// @return Whatever GSL returns.
     96    ///
     97    inline int div_elements(const matrix& b)
     98      { return gsl_matrix_div_elements(m_,b.m_); }
    6999
    70100    ///
    … …  
    75105
    76106    ///
    77     /// Create a new copy of the internal GSL matrix.
    78     ///
    79     /// Necessary memory for the new GSL matrix is allocated and the
    80     /// caller is responsible for freeing the allocated memory.
    81     ///
    82     /// @return A pointer to a copy of the internal GSL matrix.
    83     ///
    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     ///
     107    /// @return True if all elements in the matrix is zero, false
     108    /// othwerwise;
     109    ///
     110    inline bool isnull(void) const { return gsl_matrix_isnull(m_); }
     111
     112    ///
     113    /// @return The maximum value of the matrix.
     114    ///
     115    inline double max(void) const { return gsl_matrix_max(m_); }
     116
     117    ///
     118    /// @return The index to the element with the maximum value of the
     119    /// matrix. The lowest index has precedence (searching in
     120    /// row-major order).
     121    ///
     122    inline std::pair<size_t,size_t> max_index(void) const;
     123
     124    ///
     125    /// @return The minimum value of the matrix.
     126    ///
     127    inline double min(void) const { return gsl_matrix_min(m_); }
     128
     129    ///
     130    /// @return The index to the element with the minimum value of the
     131    /// matrix. The lowest index has precedence (searching in
     132    /// row-major order).
     133    ///
     134    inline std::pair<size_t,size_t> min_index(void) const;
     135
     136    ///
     137    /// @return The indecies to the element with the minimum and
     138    /// maximum values of the matrix, respectively. The lowest index
     139    /// has precedence (searching in row-major order).
     140    ///
     141    inline void minmax_index(std::pair<size_t,size_t>& min,
     142                             std::pair<size_t,size_t>& max) const
     143      { gsl_matrix_minmax_index(m_,&min.first,&min.second,
     144                                &max.first,&max.second); }
     145
     146    ///
     147    /// @return The minimum and maximum values of the matrix, as the
     148    /// \a first and \a second member of the returned \a pair,
     149    /// respectively.
     150    ///
     151    std::pair<double,double> minmax(void) const;
     152
     153    ///
     154    /// Multiply the elements of matrix \a b with the elements of the
     155    /// calling matrix ,\f$a_{ij} = a_{ij} * b_{ij} \; \forall
     156    /// i,j\f$. The result is stored into the calling matrix.
     157    ///
     158    /// @return Whatever GSL returns.
     159    ///
     160    inline int mul_elements(const matrix& b)
     161      { return gsl_matrix_mul_elements(m_,b.m_); }
     162
     163    ///
    111164    /// @return The number of rows in the matrix.
    112     ///
    113     inline size_t rows(void) const { return m_->size1; }
     165    ///
     166    inline size_t rows(void) const { return m_->size1; }
     167
     168    ///
     169    /// Multiply the elements of the calling matrix with a scalar \a
     170    /// d, \f$a_{ij} = d * a_{ij} \; \forall i,j\f$. The result is
     171    /// stored into the calling matrix.
     172    ///
     173    /// @return Whatever GSL returns.
     174    ///
     175    inline int scale(const double d) { return gsl_matrix_scale(m_,d); }
    114176
    115177    ///
    116178    /// 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     ///
     179    ///
     180    inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }
     181
     182    ///
     183    /// Set \a column to \a vec.
     184    ///
     185    /// @return Whatever GSL returns.
     186    ///
     187    inline int set_column(const size_t column, const vector& vec)
     188      { return gsl_matrix_set_col(m_, column, vec.gsl_vector_pointer()); }
     189
     190    ///
     191    /// Set \a row to \a vec.
     192    ///
     193    /// @return Whatever GSL returns.
     194    ///
     195    inline int set_row(const size_t row, const vector& vec)
     196      { return gsl_matrix_set_row(m_, row, vec.gsl_vector_pointer()); }
     197
     198    ///
     199    /// Subtract the elements of matrix \a b from the elements of the
     200    /// calling matrix ,\f$a_{ij} = a_{ij} - b_{ij} \; \forall
     201    /// i,j\f$. The result is stored into the calling matrix.
     202    ///
     203    /// @return Whatever GSL returns.
     204    ///
     205    inline int sub(const matrix& b) { return gsl_matrix_sub(m_,b.m_); }
     206
     207    ///
     208    /// Exchange the elements of the this and \a other by copying. The
     209    /// two matrices must have the same size.
     210    ///
     211    /// @return Whatever GSL returns.
     212    ///
     213    inline int swap(matrix& other) { return gsl_matrix_swap(m_,other.m_); }
     214
     215    ///
     216    /// Swap columns \a i and \a j.
     217    ///
     218    inline int swap_columns(const size_t i,const size_t j)
     219      { return gsl_matrix_swap_columns(m_,i,j); }
     220
     221    ///
     222    /// Swap row \a i and column \a j.
     223    ///
     224    inline int swap_rowcol(const size_t i,const size_t j)
     225      { return gsl_matrix_swap_rowcol(m_,i,j); }
     226
     227    ///
     228    /// Swap rows \a i and \a j.
     229    ///
     230    inline int swap_rows(const size_t i, const size_t j)
     231      { return gsl_matrix_swap_rows(m_,i,j); }
     232
     233    ///
    138234    /// Transpose the matrix.
    139235    ///
    140     /// @return Transpose of the matrix.
    141     ///
    142     matrix transpose(void) const;
    143    
     236    void transpose(void);
     237
    144238    ///
    145239    /// @return Reference to the element position (\a row, \a column).
    146     ///
     240    ///
    147241    inline double& operator()(size_t row,size_t column)
    148242      { return (*gsl_matrix_ptr(m_,row,column)); }
    … …  
    151245    /// @return Const reference to the element position (\a row, \a
    152246    /// column).
    153     ///
     247    ///
    154248    inline const double& operator()(size_t row,size_t column) const
    155249      { return (*gsl_matrix_const_ptr(m_,row,column)); }
    156250
    157251    ///
    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;
     252    /// Write me, should return a vector view.
     253    ///
     254    vector operator[](size_t row) const;
    166255
    167256    ///
    … …  
    169258    /// a vector object (matrix row). Underlying GSL supports only GSL_vector
    170259    /// views, thus some redesign of the vector class is needed.
    171     ///
    172     /// inline const vector& operator[](size_t i) const;
     260    ///
     261    /// inline const vector& operator[](size_t i) const;
    173262    ///
    174263    /// So for now, stupid copying is done, and actually a matrix
    175264    /// column is returned using this function.
     265    ///
     266    /// Jari, cf. the first two function in GSL manual section 8.4.8
    176267    vector TEMP_col_return(size_t column) const;
    177268
    178     ///
    179     /// Matrix multiplication.
    180     ///
    181     /// @return The resulting matrix.
    182     ///
    183     matrix operator*(const matrix&) const;
    184 
    185     ///
     269    ///
    186270    /// Matrix-vector multiplication.
    187271    ///
    188272    /// @return The resulting vector.
    189273    ///
    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;
     274    // Jari, where should this go?
     275    const vector operator*(const vector&) const;
    219276
    220277    ///
    … …  
    225282    /// @see equal
    226283    ///
    227     inline bool operator==( const matrix& other ) const { return equal(other);}
     284    inline bool operator==(const matrix& other) const { return equal(other); }
    228285
    229286    ///
    … …  
    234291    /// @see equal
    235292    ///
    236     inline bool operator!=( const matrix& other ) const {return !equal(other);}
     293    inline bool operator!=(const matrix& other) const { return !equal(other); }
    237294
    238295    ///
    … …  
    240297    /// dimensions.
    241298    ///
    242     matrix& operator=(const matrix& other);
     299    const matrix& operator=(const matrix& other);
     300
     301    ///
     302    /// Add and assign operator.
     303    ///
     304    inline const matrix& operator+=(const matrix& m) { add(m); return *this; }
     305
     306    ///
     307    /// Add and assign operator.
     308    ///
     309    inline const matrix& operator+=(const double d)
     310      { add_constant(d); return *this; }
    243311
    244312    ///
    245313    /// Subtract and assign operator.
    246314    ///
    247     inline matrix& operator-=(const matrix& m)
    248       { gsl_matrix_sub(m_,m.m_); return *this; }
     315    inline const matrix& operator-=(const matrix& m) { sub(m); return *this; }
     316
     317    ///
     318    /// Multiply and assigment operator.
     319    ///
     320    /// @return Const reference to the resulting matrix.
     321    ///
     322    const matrix& operator*=(const matrix&);
    249323
    250324    ///
    251325    /// Multiply and assign operator.
    252326    ///
    253     inline matrix& operator*=(const double d)
    254       { gsl_matrix_scale(m_,d); return *this; }
    255 
    256 
    257   private:
    258 
    259     gsl_matrix* m_;
     327    inline const matrix& operator*=(const double d) { scale(d); return *this; }
     328
     329
     330  private:
     331    ///
     332    /// Create a new copy of the internal GSL matrix.
     333    ///
     334    /// Necessary memory for the new GSL matrix is allocated and the
     335    /// caller is responsible for freeing the allocated memory.
     336    ///
     337    /// @return A pointer to a copy of the internal GSL matrix.
     338    ///
     339    gsl_matrix* gsl_matrix_copy(void) const;
     340
     341    ///
     342    /// @return A const pointer to the internal GSL matrix.
     343    ///
     344    // Jari, is this needed?
     345    inline const gsl_matrix* gsl_matrix_pointer(void) const { return m_; }
     346
     347    ///
     348    /// @return A pointer to the internal GSL matrix.
     349    ///
     350    // Jari, is this needed?
     351    inline gsl_matrix* gsl_matrix_pointer(void) { return m_; };
     352
     353    gsl_matrix* m_;
    260354
    261355  };
    … …  
    263357  ///
    264358  /// The output operator for the matrix class.
    265   ///
    266   std::ostream& operator<< (std::ostream& s, const matrix&);
    267 
    268 
     359  ///
     360  std::ostream& operator<< (std::ostream& s, const matrix&);
    269361
    270362}} // of namespace gslapi and namespace theplu
    271363
    272 #endif 
     364#endif
Note: See TracChangeset for help on using the changeset viewer.