Dogcows Code - chaz/yoink/blob - src/cml/matrix/external.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 external_matrix_h
  14 #define external_matrix_h
  15
  16 #include <cml/core/external_2D.h>
  17 #include <cml/matrix/matrix_expr.h>
  18 #include <cml/matrix/class_ops.h>
  19 #include <cml/matrix/matrix_unroller.h>
  20 #include <cml/matrix/dynamic.h>
  21
  22 namespace cml {
  23
  24 /** Fixed-size, external-memory matrix. */
  25 template<typename Element, int Rows, int Cols,
  26     typename BasisOrient, typename Layout>
  27 class matrix<Element,external<Rows,Cols>,BasisOrient,Layout>
  28 : public external_2D<Element,Rows,Cols,Layout>
  29 {
  30   public:
  31
  32     /* Shorthand for the generator: */
  33     typedef external<Rows,Cols> generator_type;
  34
  35     /* Shorthand for the array type: */
  36     typedef external_2D<Element,Rows,Cols,Layout> array_type;
  37
  38     /* Shorthand for the type of this matrix: */
  39     typedef matrix<Element,generator_type,BasisOrient,Layout> matrix_type;
  40
  41     /* For integration into the expression template code: */
  42     typedef matrix_type expr_type;
  43
  44     /* For integration into the expression template code: */
  45     typedef matrix<Element,fixed<Rows,Cols>,BasisOrient,Layout> temporary_type;
  46     /* Note: this ensures that an external matrix is copied into the proper
  47      * temporary; external<> temporaries are not allowed.
  48      */
  49
  50     /* Standard: */
  51     typedef typename array_type::value_type value_type;
  52     typedef typename array_type::reference reference;
  53     typedef typename array_type::const_reference const_reference;
  54
  55     typedef matrix_type& expr_reference;
  56     typedef const matrix_type& expr_const_reference;
  57
  58     /* For matching by basis: */
  59     typedef BasisOrient basis_orient;
  60
  61     /* For matching by memory layout: */
  62     typedef typename array_type::layout layout;
  63
  64     /* For matching by storage type if necessary: */
  65     typedef typename array_type::memory_tag memory_tag;
  66
  67     /* For matching by size type if necessary: */
  68     typedef typename array_type::size_tag size_tag;
  69
  70     /* For matching by resizability: */
  71     typedef typename array_type::resizing_tag resizing_tag;
  72
  73     /* For matching by result-type: */
  74     typedef cml::et::matrix_result_tag result_tag;
  75
  76     /* For matching by assignability: */
  77     typedef cml::et::assignable_tag assignable_tag;
  78
  79     /* To simplify the matrix transpose operator: */
  80     typedef matrix<
  81         Element,
  82         typename array_type::transposed_type::generator_type,
  83         BasisOrient,
  84         Layout
  85     > transposed_type;
  86
  87     /* To simplify the matrix row and column operators: */
  88     typedef vector<
  89         Element,
  90         typename array_type::row_array_type::generator_type
  91     > row_vector_type;
  92
  93     typedef vector<
  94         Element,
  95         typename array_type::col_array_type::generator_type
  96     > col_vector_type;
  97
  98
  99   public:
 100
 101     /** Set this matrix to zero. */
 102     matrix_type& zero() {
 103         typedef cml::et::OpAssign<Element,Element> OpT;
 104         cml::et::UnrollAssignment<OpT>(*this,Element(0));
 105         return *this;
 106     }
 107
 108     /** Set this matrix to the identity.
 109      *
 110      * This only makes sense for a square matrix, but no error will be
 111      * signaled if the matrix is not square.
 112      */
 113     matrix_type& identity() {
 114         for(size_t i = 0; i < this->rows(); ++ i) {
 115             for(size_t j = 0; j < this->cols(); ++ j) {
 116                 (*this)(i,j) = value_type((i == j)?1:0);
 117             }
 118         }
 119         return *this;
 120     }
 121
 122     /** Set this matrix to its transpose.
 123      *
 124      * This only makes sense for a square matrix, but no error will be
 125      * signaled if the matrix is not square.
 126      */
 127     matrix_type& transpose() {
 128         /* transpose() returns a temporary: */
 129         *this = transpose(*this);
 130         return *this;
 131     }
 132
 133     /** Set this matrix to its inverse.
 134      *
 135      * This only makes sense for a square matrix, but no error will be
 136      * signaled if the matrix is not square.
 137      */
 138     matrix_type& inverse() {
 139         /* inverse() returns a temporary: */
 140         *this = cml::inverse(*this);
 141         return *this;
 142     }
 143
 144     /* NOTE: minimize() and maximize() no longer supported (Jesse) */
 145
 146     #if 0
 147     /** Pairwise minimum of this matrix with another. */
 148     template<typename E, class AT, typename L>
 149     void minimize(const matrix<E,AT,basis_orient,L>& v) {
 150       /* XXX This should probably use ScalarPromote: */
 151       for (size_t i = 0; i < this->rows(); ++i) {
 152         for (size_t j = 0; j < this->cols(); ++j) {
 153           (*this)(i,j) = std::min((*this)(i,j),v(i,j));
 154         }
 155       }
 156     }
 157
 158     /** Pairwise maximum of this matrix with another. */
 159     template<typename E, class AT, typename L>
 160     void maximize(const matrix<E,AT,basis_orient,L>& v) {
 161       /* XXX This should probably use ScalarPromote: */
 162       for (size_t i = 0; i < this->rows(); ++i) {
 163         for (size_t j = 0; j < this->cols(); ++j) {
 164           (*this)(i,j) = std::max((*this)(i,j),v(i,j));
 165         }
 166       }
 167     }
 168     #endif
 169
 170     /* Set each element to a random number in the range [min,max] */
 171     void random(ELEMENT_ARG_TYPE min, ELEMENT_ARG_TYPE max) {
 172       for(size_t i = 0; i < this->rows(); ++i) {
 173         for(size_t j = 0; j < this->cols(); ++j) {
 174           (*this)(i,j) = random_real(min,max);
 175         }
 176       }
 177     }
 178
 179
 180   public:
 181
 182     /** Constructor for fixed-size external matrices.
 183      *
 184      * The array must be given as a pointer to Element*, not a
 185      * multi-dimensional array.  The caller owns the pointer, and is
 186      * responsible for doing any necessary memory management.
 187      *
 188      * @param ptr specify the external pointer.
 189      *
 190      * @throws same as the ArrayType constructor.
 191      */
 192     explicit matrix(value_type ptr[Rows][Cols]) : array_type(ptr) {}
 193
 194     /** Constructor for fixed-size external matrices.
 195      *
 196      * The array must be given as a pointer to Element*, not a
 197      * multi-dimensional array.  The caller owns the pointer, and is
 198      * responsible for doing any necessary memory management.
 199      *
 200      * @param ptr specify the external pointer.
 201      *
 202      * @throws same as the ArrayType constructor.
 203      */
 204     explicit matrix(value_type* ptr) : array_type(ptr) {}
 205
 206
 207   public:
 208
 209     /** Return the matrix size as a pair. */
 210     matrix_size size() const {
 211         return matrix_size(this->rows(),this->cols());
 212     }
 213
 214     /** Return element j of basis vector i. */
 215     value_type basis_element(size_t i, size_t j) const {
 216         return basis_element(i,j,basis_orient());
 217     }
 218
 219     /** Set the given basis element. */
 220     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s) {
 221         set_basis_element(i,j,s,basis_orient());
 222     }
 223
 224
 225   public:
 226
 227     CML_ASSIGN_MAT_22
 228     CML_ASSIGN_MAT_33
 229     CML_ASSIGN_MAT_44
 230
 231     /* Define class operators for external matrices. Note: external matrices
 232      * cannot be copy-constructed, but they can be assigned to:
 233      */
 234     CML_MAT_ASSIGN_FROM_MATTYPE
 235
 236     CML_MAT_ASSIGN_FROM_MAT(=, et::OpAssign)
 237     CML_MAT_ASSIGN_FROM_MAT(+=, et::OpAddAssign)
 238     CML_MAT_ASSIGN_FROM_MAT(-=, et::OpSubAssign)
 239
 240     CML_MAT_ASSIGN_FROM_MATXPR(=, et::OpAssign)
 241     CML_MAT_ASSIGN_FROM_MATXPR(+=, et::OpAddAssign)
 242     CML_MAT_ASSIGN_FROM_MATXPR(-=, et::OpSubAssign)
 243
 244     CML_MAT_ASSIGN_FROM_SCALAR(*=, et::OpMulAssign)
 245     CML_MAT_ASSIGN_FROM_SCALAR(/=, et::OpDivAssign)
 246
 247     CML_ACCUMULATED_MATRIX_MULT(const matrix_type&)
 248
 249     template<typename E, class AT, typename BO, typename L>
 250         CML_ACCUMULATED_MATRIX_MULT(const TEMPLATED_MATRIX_MACRO&)
 251
 252     template<class XprT>
 253         CML_ACCUMULATED_MATRIX_MULT(MATXPR_ARG_TYPE)
 254
 255
 256   protected:
 257
 258     value_type basis_element(size_t i, size_t j, row_basis) const {
 259         return (*this)(i,j);
 260     }
 261
 262     value_type basis_element(size_t i, size_t j, col_basis) const {
 263         return (*this)(j,i);
 264     }
 265
 266     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s, row_basis) {
 267         (*this)(i,j) = s;
 268     }
 269
 270     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s, col_basis) {
 271         (*this)(j,i) = s;
 272     }
 273
 274
 275   public:
 276
 277     /* Braces should only be used for testing: */
 278 #if defined(CML_ENABLE_MATRIX_BRACES)
 279     CML_MATRIX_BRACE_OPERATORS
 280 #endif
 281 };
 282
 283 /** Dynamic-size, external-memory matrix. */
 284 template<typename Element, typename BasisOrient, typename Layout>
 285 class matrix<Element,external<-1,-1>,BasisOrient,Layout>
 286 : public external_2D<Element,-1,-1,Layout>
 287 {
 288   public:
 289
 290     /* Shorthand for the generator: */
 291     typedef external<> generator_type;
 292
 293     /* Shorthand for the array type: */
 294     typedef external_2D<Element,-1,-1,Layout> array_type;
 295
 296     /* Shorthand for the type of this matrix: */
 297     typedef matrix<Element,generator_type,BasisOrient,Layout> matrix_type;
 298
 299     /* For integration into the expression template code: */
 300     typedef matrix_type expr_type;
 301
 302     /* For integration into the expression template code: */
 303     typedef matrix<Element,dynamic<>,BasisOrient,Layout> temporary_type;
 304     /* Note: this ensures that an external matrix is copied into the proper
 305      * temporary; external<> temporaries are not allowed.
 306      */
 307
 308     /* Standard: */
 309     typedef typename array_type::value_type value_type;
 310     typedef typename array_type::reference reference;
 311     typedef typename array_type::const_reference const_reference;
 312
 313     typedef matrix_type& expr_reference;
 314     typedef const matrix_type& expr_const_reference;
 315
 316     /* For matching by basis: */
 317     typedef BasisOrient basis_orient;
 318
 319     /* For matching by memory layout: */
 320     typedef typename array_type::layout layout;
 321
 322     /* For matching by storage type if necessary: */
 323     typedef typename array_type::memory_tag memory_tag;
 324
 325     /* For matching by size type if necessary: */
 326     typedef typename array_type::size_tag size_tag;
 327
 328     /* For matching by resizability: */
 329     typedef typename array_type::resizing_tag resizing_tag;
 330
 331     /* For matching by result-type: */
 332     typedef cml::et::matrix_result_tag result_tag;
 333
 334     /* For matching by assignability: */
 335     typedef cml::et::assignable_tag assignable_tag;
 336
 337     /* To simplify the matrix transpose operator: */
 338     typedef matrix<
 339         Element,
 340         typename array_type::transposed_type::generator_type,
 341         BasisOrient,
 342         Layout
 343     > transposed_type;
 344
 345     /* To simplify the matrix row and column operators: */
 346     typedef vector<
 347         Element,
 348         typename array_type::row_array_type::generator_type
 349     > row_vector_type;
 350
 351     typedef vector<
 352         Element,
 353         typename array_type::col_array_type::generator_type
 354     > col_vector_type;
 355
 356
 357   public:
 358
 359     /** Set this matrix to zero. */
 360     matrix_type& zero() {
 361         typedef cml::et::OpAssign<Element,Element> OpT;
 362         cml::et::UnrollAssignment<OpT>(*this,Element(0));
 363         return *this;
 364     }
 365
 366     /** Set this matrix to the identity.
 367      *
 368      * This only makes sense for a square matrix, but no error will be
 369      * signaled if the matrix is not square.
 370      */
 371     matrix_type& identity() {
 372         for(size_t i = 0; i < this->rows(); ++ i) {
 373             for(size_t j = 0; j < this->cols(); ++ j) {
 374                 (*this)(i,j) = value_type((i == j)?1:0);
 375             }
 376         }
 377         return *this;
 378     }
 379
 380     /** Set this matrix to its transpose.
 381      *
 382      * This only makes sense for a square matrix, but no error will be
 383      * signaled if the matrix is not square.
 384      */
 385     matrix_type& transpose() {
 386         /* transpose() returns a temporary: */
 387         *this = cml::transpose(*this);
 388         return *this;
 389     }
 390
 391     /** Set this matrix to its inverse.
 392      *
 393      * This only makes sense for a square matrix, but no error will be
 394      * signaled if the matrix is not square.
 395      */
 396     matrix_type& inverse() {
 397         /* inverse() returns a temporary: */
 398         *this = inverse(*this);
 399         return *this;
 400     }
 401
 402     /** Pairwise minimum of this matrix with another. */
 403     template<typename E, class AT, typename L>
 404     void minimize(const matrix<E,AT,basis_orient,L>& v) {
 405       /* XXX This should probably use ScalarPromote: */
 406       for (size_t i = 0; i < this->rows(); ++i) {
 407         for (size_t j = 0; j < this->cols(); ++j) {
 408           (*this)[i] = std::min((*this)(i,j),v(i,j));
 409         }
 410       }
 411     }
 412
 413     /** Pairwise maximum of this matrix with another. */
 414     template<typename E, class AT, class BO, typename L>
 415     void maximize(const matrix<E,AT,basis_orient,L>& v) {
 416       /* XXX This should probably use ScalarPromote: */
 417       for (size_t i = 0; i < this->rows(); ++i) {
 418         for (size_t j = 0; j < this->cols(); ++j) {
 419           (*this)[i] = std::max((*this)(i,j),v(i,j));
 420         }
 421       }
 422     }
 423
 424     /* Set each element to a random number in the range [min,max] */
 425     void random(ELEMENT_ARG_TYPE min, ELEMENT_ARG_TYPE max) {
 426       for(size_t i = 0; i < this->rows(); ++i) {
 427         for(size_t j = 0; j < this->cols(); ++j) {
 428           (*this)(i,j) = cml::random_real(min,max);
 429         }
 430       }
 431     }
 432
 433
 434   public:
 435
 436     /** Constructor for fixed-size external matrices.
 437      *
 438      * The array must be given as a pointer to Element*, not a
 439      * multi-dimensional array.  The caller owns the pointer, and is
 440      * responsible for doing any necessary memory management.
 441      *
 442      * @param ptr specify the external pointer.
 443      *
 444      * @throws same as the ArrayType constructor.
 445      */
 446     explicit matrix(value_type* const ptr, size_t rows, size_t cols)
 447         : array_type(ptr,rows,cols) {}
 448
 449
 450   public:
 451
 452     /** Return the matrix size as a pair. */
 453     matrix_size size() const {
 454         return matrix_size(this->rows(),this->cols());
 455     }
 456
 457     /** Return element j of basis vector i. */
 458     value_type basis_element(size_t i, size_t j) const {
 459         return basis_element(i,j,basis_orient());
 460     }
 461
 462     /** Set the given basis element. */
 463     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s) {
 464         set_basis_element(i,j,s,basis_orient());
 465     }
 466
 467
 468   public:
 469
 470     CML_ASSIGN_MAT_22
 471     CML_ASSIGN_MAT_33
 472     CML_ASSIGN_MAT_44
 473
 474     /* Define class operators for external matrices. Note: external matrices
 475      * cannot be copy-constructed, but they can be assigned to:
 476      */
 477     CML_MAT_ASSIGN_FROM_MATTYPE
 478
 479     CML_MAT_ASSIGN_FROM_MAT(=, et::OpAssign)
 480     CML_MAT_ASSIGN_FROM_MAT(+=, et::OpAddAssign)
 481     CML_MAT_ASSIGN_FROM_MAT(-=, et::OpSubAssign)
 482
 483     CML_MAT_ASSIGN_FROM_MATXPR(=, et::OpAssign)
 484     CML_MAT_ASSIGN_FROM_MATXPR(+=, et::OpAddAssign)
 485     CML_MAT_ASSIGN_FROM_MATXPR(-=, et::OpSubAssign)
 486
 487     CML_MAT_ASSIGN_FROM_SCALAR(*=, et::OpMulAssign)
 488     CML_MAT_ASSIGN_FROM_SCALAR(/=, et::OpDivAssign)
 489
 490     CML_ACCUMULATED_MATRIX_MULT(const matrix_type&)
 491
 492     template<typename E, class AT, typename BO, typename L>
 493         CML_ACCUMULATED_MATRIX_MULT(const TEMPLATED_MATRIX_MACRO&)
 494
 495     template<class XprT>
 496         CML_ACCUMULATED_MATRIX_MULT(MATXPR_ARG_TYPE)
 497
 498
 499   protected:
 500
 501     value_type basis_element(size_t i, size_t j, row_basis) const {
 502         return (*this)(i,j);
 503     }
 504
 505     value_type basis_element(size_t i, size_t j, col_basis) const {
 506         return (*this)(j,i);
 507     }
 508
 509     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s, row_basis) {
 510         (*this)(i,j) = s;
 511     }
 512
 513     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s, col_basis) {
 514         (*this)(j,i) = s;
 515     }
 516
 517
 518   public:
 519
 520     /* Braces should only be used for testing: */
 521 #if defined(CML_ENABLE_MATRIX_BRACES)
 522     CML_MATRIX_BRACE_OPERATORS
 523 #endif
 524 };
 525
 526 } // namespace cml
 527
 528 #endif
 529
 530 // -------------------------------------------------------------------------
 531 // vim:ft=cpp