Dogcows Code - chaz/yoink/blob - src/cml/matrix/dynamic.h

   1 /* -*- C++ -*- ------------------------------------------------------------
   2
   3 Copyright (c) 2007 Jesse Anders and Demian Nave http://cmldev.net/
   4
   5 The Configurable Math Library (CML) is distributed under the terms of the
   6 Boost Software License, v1.0 (see cml/LICENSE for details).
   7
   8  *-----------------------------------------------------------------------*/
   9 /** @file
  10  *  @brief
  11  */
  12
  13 #ifndef dynamic_matrix_h
  14 #define dynamic_matrix_h
  15
  16 #include <cml/core/dynamic_2D.h>
  17 #include <cml/matrix/matrix_expr.h>
  18 #include <cml/matrix/class_ops.h>
  19 #include <cml/matrix/matrix_unroller.h>
  20
  21 namespace cml {
  22
  23 /** Resizeable, dynamic-memory matrix. */
  24 template<typename Element, typename Alloc,
  25     typename BasisOrient, typename Layout>
  26 class matrix<Element,dynamic<Alloc>,BasisOrient,Layout>
  27 : public dynamic_2D<Element,Layout,Alloc>
  28 {
  29   public:
  30
  31     /* Shorthand for the generator: */
  32     typedef dynamic<Alloc> generator_type;
  33
  34     /* Shorthand for the array type: */
  35     typedef dynamic_2D<Element,Layout,Alloc> array_type;
  36
  37     /* Shorthand for the type of this matrix: */
  38     typedef matrix<Element,generator_type,BasisOrient,Layout> matrix_type;
  39
  40     /* For integration into the expression template code: */
  41     typedef matrix_type expr_type;
  42
  43     /* For integration into the expression template code: */
  44     typedef matrix_type temporary_type;
  45
  46     /* Standard: */
  47     typedef typename array_type::value_type value_type;
  48     typedef typename array_type::reference reference;
  49     typedef typename array_type::const_reference const_reference;
  50
  51     /* For integration into the expression templates code: */
  52     typedef matrix_type& expr_reference;
  53     typedef const matrix_type& expr_const_reference;
  54
  55     /* For matching by basis: */
  56     typedef BasisOrient basis_orient;
  57
  58     /* For matching by memory layout: */
  59     typedef typename array_type::layout layout;
  60
  61     /* For matching by storage type: */
  62     typedef typename array_type::memory_tag memory_tag;
  63
  64     /* For matching by size type if necessary: */
  65     typedef typename array_type::size_tag size_tag;
  66
  67     /* For matching by resizability: */
  68     typedef typename array_type::resizing_tag resizing_tag;
  69
  70     /* For matching by result type: */
  71     typedef cml::et::matrix_result_tag result_tag;
  72
  73     /* For matching by assignability: */
  74     typedef cml::et::assignable_tag assignable_tag;
  75
  76     /* To simplify the matrix transpose operator: */
  77     typedef matrix<
  78         Element,
  79         typename array_type::transposed_type::generator_type,
  80         BasisOrient,
  81         Layout
  82     > transposed_type;
  83
  84     /* To simplify the matrix row and column operators: */
  85     typedef vector<
  86         Element,
  87         typename array_type::row_array_type::generator_type
  88     > row_vector_type;
  89
  90     typedef vector<
  91         Element,
  92         typename array_type::col_array_type::generator_type
  93     > col_vector_type;
  94
  95
  96   public:
  97
  98     /** Set this matrix to zero. */
  99     matrix_type& zero() {
 100         typedef cml::et::OpAssign<Element,Element> OpT;
 101         cml::et::UnrollAssignment<OpT>(*this,Element(0));
 102         return *this;
 103     }
 104
 105     /** Set this matrix to the identity.
 106      *
 107      * This only makes sense for a square matrix, but no error will be
 108      * signaled if the matrix is not square.
 109      */
 110     matrix_type& identity() {
 111         for(size_t i = 0; i < this->rows(); ++ i) {
 112             for(size_t j = 0; j < this->cols(); ++ j) {
 113                 (*this)(i,j) = value_type((i == j)?1:0);
 114             }
 115         }
 116         return *this;
 117     }
 118
 119     /** Set this matrix to its transpose.
 120      *
 121      * This only makes sense for a square matrix, but no error will be
 122      * signaled if the matrix is not square.
 123      */
 124     matrix_type& transpose() {
 125         /* transpose() returns a temporary: */
 126         *this = cml::transpose(*this);
 127         return *this;
 128     }
 129
 130     /** Set this matrix to its inverse.
 131      *
 132      * This only makes sense for a square matrix, but no error will be
 133      * signaled if the matrix is not square.
 134      */
 135     matrix_type& inverse() {
 136         /* inverse() returns a temporary: */
 137         *this = cml::inverse(*this);
 138         return *this;
 139     }
 140
 141     /* NOTE: minimize() and maximize() no longer supported (Jesse) */
 142
 143     #if 0
 144     /** Pairwise minimum of this matrix with another. */
 145     template<typename E, class AT, typename L>
 146     void minimize(const matrix<E,AT,basis_orient,L>& v) {
 147       /* XXX This should probably use ScalarPromote: */
 148       for (size_t i = 0; i < this->rows(); ++i) {
 149         for (size_t j = 0; j < this->cols(); ++j) {
 150           (*this)(i,j) = std::min((*this)(i,j),v(i,j));
 151         }
 152       }
 153     }
 154
 155     /** Pairwise maximum of this matrix with another. */
 156     template<typename E, class AT, typename L>
 157     void maximize(const matrix<E,AT,basis_orient,L>& v) {
 158       /* XXX This should probably use ScalarPromote: */
 159       for (size_t i = 0; i < this->rows(); ++i) {
 160         for (size_t j = 0; j < this->cols(); ++j) {
 161           (*this)(i,j) = std::max((*this)(i,j),v(i,j));
 162         }
 163       }
 164     }
 165     #endif
 166
 167     /* Set each element to a random number in the range [min,max] */
 168     void random(ELEMENT_ARG_TYPE min, ELEMENT_ARG_TYPE max) {
 169       for(size_t i = 0; i < this->rows(); ++i) {
 170         for(size_t j = 0; j < this->cols(); ++j) {
 171           (*this)(i,j) = cml::random_real(min,max);
 172         }
 173       }
 174     }
 175
 176
 177   public:
 178
 179     /** Default constructor. */
 180     matrix() {}
 181
 182     /** Constructor for dynamically-sized arrays.
 183      *
 184      * @param rows specify the number of rows.
 185      * @param cols specify the number of cols.
 186      */
 187     explicit matrix(size_t rows, size_t cols)
 188         : array_type(rows,cols) {}
 189
 190
 191   public:
 192
 193     /** Return the matrix size as a pair. */
 194     matrix_size size() const {
 195         return matrix_size(this->rows(),this->cols());
 196     }
 197
 198     /** Return element j of basis vector i. */
 199     value_type basis_element(size_t i, size_t j) const {
 200         return basis_element(i,j,basis_orient());
 201     }
 202
 203     /** Set the given basis element. */
 204     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s) {
 205         set_basis_element(i,j,s,basis_orient());
 206     }
 207
 208     /** Set the matrix row from the given vector. */
 209     void set_row(size_t i, const row_vector_type& row) {
 210       for(size_t j = 0; j < this->cols(); ++ j) (*this)(i,j) = row[j];
 211     }
 212
 213     /** Set the matrix column from the given vector. */
 214     void set_col(size_t j, const col_vector_type& col) {
 215       for(size_t i = 0; i < this->rows(); ++ i) (*this)(i,j) = col[i];
 216     }
 217
 218
 219   public:
 220
 221     /* Define common class operators: */
 222
 223     CML_CONSTRUCT_MAT_22
 224     CML_CONSTRUCT_MAT_33
 225     CML_CONSTRUCT_MAT_44
 226
 227     CML_MAT_COPY_FROM_ARRAY(: array_type())
 228     CML_MAT_COPY_FROM_MATTYPE
 229     CML_MAT_COPY_FROM_MAT
 230     CML_MAT_COPY_FROM_MATXPR
 231
 232     CML_ASSIGN_MAT_22
 233     CML_ASSIGN_MAT_33
 234     CML_ASSIGN_MAT_44
 235
 236     CML_MAT_ASSIGN_FROM_MATTYPE
 237
 238     CML_MAT_ASSIGN_FROM_MAT(=, et::OpAssign)
 239     CML_MAT_ASSIGN_FROM_MAT(+=, et::OpAddAssign)
 240     CML_MAT_ASSIGN_FROM_MAT(-=, et::OpSubAssign)
 241
 242     CML_MAT_ASSIGN_FROM_MATXPR(=, et::OpAssign)
 243     CML_MAT_ASSIGN_FROM_MATXPR(+=, et::OpAddAssign)
 244     CML_MAT_ASSIGN_FROM_MATXPR(-=, et::OpSubAssign)
 245
 246     CML_MAT_ASSIGN_FROM_SCALAR(*=, et::OpMulAssign)
 247     CML_MAT_ASSIGN_FROM_SCALAR(/=, et::OpDivAssign)
 248
 249     /** Accumulated matrix multiplication.
 250      *
 251      * This only makes sense for a square matrix, but no error will be
 252      * signaled if the matrix is not square.
 253      */
 254     matrix_type& operator*=(const matrix_type& m) {
 255         /* Matrix multiplication returns a temporary: */
 256         *this = (*this)*m;
 257         return *this;
 258     }
 259
 260     /** Accumulated matrix multiplication.
 261      *
 262      * This only makes sense for a square matrix, but no error will be
 263      * signaled if the matrix is not square.
 264      */
 265     template<typename E, class AT, typename BO, typename L> matrix_type&
 266     operator*=(const matrix<E,AT,BO,L>& m) {
 267         /* Matrix multiplication returns a temporary: */
 268         *this = (*this)*m;
 269         return *this;
 270     }
 271
 272     /** Accumulated matrix multiplication.
 273      *
 274      * This only makes sense for a square matrix, but no error will be
 275      * signaled if the matrix is not square.
 276      */
 277     template<class XprT> matrix_type&
 278     operator*=(MATXPR_ARG_TYPE e) {
 279         /* Verify that a promotion exists at compile time: */
 280         typedef typename et::MatrixPromote<
 281             matrix_type, typename XprT::result_type>::type result_type;
 282         /* Matrix multiplication returns a temporary: */
 283         *this = (*this)*e;
 284         return *this;
 285     }
 286
 287
 288   protected:
 289
 290     value_type basis_element(size_t i, size_t j, row_basis) const {
 291         return (*this)(i,j);
 292     }
 293
 294     value_type basis_element(size_t i, size_t j, col_basis) const {
 295         return (*this)(j,i);
 296     }
 297
 298     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s, row_basis) {
 299         (*this)(i,j) = s;
 300     }
 301
 302     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s, col_basis) {
 303         (*this)(j,i) = s;
 304     }
 305
 306
 307   public:
 308
 309     /* Braces should only be used for testing: */
 310 #if defined(CML_ENABLE_MATRIX_BRACES)
 311     CML_MATRIX_BRACE_OPERATORS
 312 #endif
 313 };
 314
 315 } // namespace cml
 316
 317 #endif
 318
 319 // -------------------------------------------------------------------------
 320 // vim:ft=cpp