Dogcows Code - chaz/yoink/blob - src/moof/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         typename cml::remove_const<Element>::type,
  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     /** Set the matrix row from the given vector. */
 225     void set_row(size_t i, const row_vector_type& row) {
 226       for(size_t j = 0; j < this->cols(); ++ j) (*this)(i,j) = row[j];
 227     }
 228
 229     /** Set the matrix column from the given vector. */
 230     void set_col(size_t j, const col_vector_type& col) {
 231       for(size_t i = 0; i < this->rows(); ++ i) (*this)(i,j) = col[i];
 232     }
 233
 234
 235   public:
 236
 237     CML_ASSIGN_MAT_22
 238     CML_ASSIGN_MAT_33
 239     CML_ASSIGN_MAT_44
 240
 241     /* Define class operators for external matrices. Note: external matrices
 242      * cannot be copy-constructed, but they can be assigned to:
 243      */
 244     CML_MAT_ASSIGN_FROM_MATTYPE
 245
 246     CML_MAT_ASSIGN_FROM_MAT(=, et::OpAssign)
 247     CML_MAT_ASSIGN_FROM_MAT(+=, et::OpAddAssign)
 248     CML_MAT_ASSIGN_FROM_MAT(-=, et::OpSubAssign)
 249
 250     CML_MAT_ASSIGN_FROM_MATXPR(=, et::OpAssign)
 251     CML_MAT_ASSIGN_FROM_MATXPR(+=, et::OpAddAssign)
 252     CML_MAT_ASSIGN_FROM_MATXPR(-=, et::OpSubAssign)
 253
 254     CML_MAT_ASSIGN_FROM_SCALAR(*=, et::OpMulAssign)
 255     CML_MAT_ASSIGN_FROM_SCALAR(/=, et::OpDivAssign)
 256
 257     CML_ACCUMULATED_MATRIX_MULT(const matrix_type&)
 258
 259     template<typename E, class AT, typename BO, typename L>
 260         CML_ACCUMULATED_MATRIX_MULT(const TEMPLATED_MATRIX_MACRO&)
 261
 262     template<class XprT>
 263         CML_ACCUMULATED_MATRIX_MULT(MATXPR_ARG_TYPE)
 264
 265
 266   protected:
 267
 268     value_type basis_element(size_t i, size_t j, row_basis) const {
 269         return (*this)(i,j);
 270     }
 271
 272     value_type basis_element(size_t i, size_t j, col_basis) const {
 273         return (*this)(j,i);
 274     }
 275
 276     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s, row_basis) {
 277         (*this)(i,j) = s;
 278     }
 279
 280     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s, col_basis) {
 281         (*this)(j,i) = s;
 282     }
 283
 284
 285   public:
 286
 287     /* Braces should only be used for testing: */
 288 #if defined(CML_ENABLE_MATRIX_BRACES)
 289     CML_MATRIX_BRACE_OPERATORS
 290 #endif
 291 };
 292
 293 /** Dynamic-size, external-memory matrix. */
 294 template<typename Element, typename BasisOrient, typename Layout>
 295 class matrix<Element,external<-1,-1>,BasisOrient,Layout>
 296 : public external_2D<Element,-1,-1,Layout>
 297 {
 298   public:
 299
 300     /* Shorthand for the generator: */
 301     typedef external<> generator_type;
 302
 303     /* Shorthand for the array type: */
 304     typedef external_2D<Element,-1,-1,Layout> array_type;
 305
 306     /* Shorthand for the type of this matrix: */
 307     typedef matrix<Element,generator_type,BasisOrient,Layout> matrix_type;
 308
 309     /* For integration into the expression template code: */
 310     typedef matrix_type expr_type;
 311
 312     /* For integration into the expression template code: */
 313     typedef matrix<Element,dynamic<>,BasisOrient,Layout> temporary_type;
 314     /* Note: this ensures that an external matrix is copied into the proper
 315      * temporary; external<> temporaries are not allowed.
 316      */
 317
 318     /* Standard: */
 319     typedef typename array_type::value_type value_type;
 320     typedef typename array_type::reference reference;
 321     typedef typename array_type::const_reference const_reference;
 322
 323     typedef matrix_type& expr_reference;
 324     typedef const matrix_type& expr_const_reference;
 325
 326     /* For matching by basis: */
 327     typedef BasisOrient basis_orient;
 328
 329     /* For matching by memory layout: */
 330     typedef typename array_type::layout layout;
 331
 332     /* For matching by storage type if necessary: */
 333     typedef typename array_type::memory_tag memory_tag;
 334
 335     /* For matching by size type if necessary: */
 336     typedef typename array_type::size_tag size_tag;
 337
 338     /* For matching by resizability: */
 339     typedef typename array_type::resizing_tag resizing_tag;
 340
 341     /* For matching by result-type: */
 342     typedef cml::et::matrix_result_tag result_tag;
 343
 344     /* For matching by assignability: */
 345     typedef cml::et::assignable_tag assignable_tag;
 346
 347     /* To simplify the matrix transpose operator: */
 348     typedef matrix<
 349         Element,
 350         typename array_type::transposed_type::generator_type,
 351         BasisOrient,
 352         Layout
 353     > transposed_type;
 354
 355     /* To simplify the matrix row and column operators: */
 356     typedef vector<
 357         Element,
 358         typename array_type::row_array_type::generator_type
 359     > row_vector_type;
 360
 361     typedef vector<
 362         Element,
 363         typename array_type::col_array_type::generator_type
 364     > col_vector_type;
 365
 366
 367   public:
 368
 369     /** Set this matrix to zero. */
 370     matrix_type& zero() {
 371         typedef cml::et::OpAssign<Element,Element> OpT;
 372         cml::et::UnrollAssignment<OpT>(*this,Element(0));
 373         return *this;
 374     }
 375
 376     /** Set this matrix to the identity.
 377      *
 378      * This only makes sense for a square matrix, but no error will be
 379      * signaled if the matrix is not square.
 380      */
 381     matrix_type& identity() {
 382         for(size_t i = 0; i < this->rows(); ++ i) {
 383             for(size_t j = 0; j < this->cols(); ++ j) {
 384                 (*this)(i,j) = value_type((i == j)?1:0);
 385             }
 386         }
 387         return *this;
 388     }
 389
 390     /** Set this matrix to its transpose.
 391      *
 392      * This only makes sense for a square matrix, but no error will be
 393      * signaled if the matrix is not square.
 394      */
 395     matrix_type& transpose() {
 396         /* transpose() returns a temporary: */
 397         *this = cml::transpose(*this);
 398         return *this;
 399     }
 400
 401     /** Set this matrix to its inverse.
 402      *
 403      * This only makes sense for a square matrix, but no error will be
 404      * signaled if the matrix is not square.
 405      */
 406     matrix_type& inverse() {
 407         /* inverse() returns a temporary: */
 408         *this = inverse(*this);
 409         return *this;
 410     }
 411
 412     /** Pairwise minimum of this matrix with another. */
 413     template<typename E, class AT, typename L>
 414     void minimize(const matrix<E,AT,basis_orient,L>& v) {
 415       /* XXX This should probably use ScalarPromote: */
 416       for (size_t i = 0; i < this->rows(); ++i) {
 417         for (size_t j = 0; j < this->cols(); ++j) {
 418           (*this)[i] = std::min((*this)(i,j),v(i,j));
 419         }
 420       }
 421     }
 422
 423     /** Pairwise maximum of this matrix with another. */
 424     template<typename E, class AT, class BO, typename L>
 425     void maximize(const matrix<E,AT,basis_orient,L>& v) {
 426       /* XXX This should probably use ScalarPromote: */
 427       for (size_t i = 0; i < this->rows(); ++i) {
 428         for (size_t j = 0; j < this->cols(); ++j) {
 429           (*this)[i] = std::max((*this)(i,j),v(i,j));
 430         }
 431       }
 432     }
 433
 434     /* Set each element to a random number in the range [min,max] */
 435     void random(ELEMENT_ARG_TYPE min, ELEMENT_ARG_TYPE max) {
 436       for(size_t i = 0; i < this->rows(); ++i) {
 437         for(size_t j = 0; j < this->cols(); ++j) {
 438           (*this)(i,j) = cml::random_real(min,max);
 439         }
 440       }
 441     }
 442
 443
 444   public:
 445
 446     /** Constructor for fixed-size external matrices.
 447      *
 448      * The array must be given as a pointer to Element*, not a
 449      * multi-dimensional array.  The caller owns the pointer, and is
 450      * responsible for doing any necessary memory management.
 451      *
 452      * @param ptr specify the external pointer.
 453      * @param rows the number of rows in the C array.
 454      * @param cols the number of columns in the C array.
 455      *
 456      * @throws same as the ArrayType constructor.
 457      */
 458     explicit matrix(value_type* const ptr, size_t rows, size_t cols)
 459         : array_type(ptr,rows,cols) {}
 460
 461
 462   public:
 463
 464     /** Return the matrix size as a pair. */
 465     matrix_size size() const {
 466         return matrix_size(this->rows(),this->cols());
 467     }
 468
 469     /** Return element j of basis vector i. */
 470     value_type basis_element(size_t i, size_t j) const {
 471         return basis_element(i,j,basis_orient());
 472     }
 473
 474     /** Set the given basis element. */
 475     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s) {
 476         set_basis_element(i,j,s,basis_orient());
 477     }
 478
 479
 480   public:
 481
 482     CML_ASSIGN_MAT_22
 483     CML_ASSIGN_MAT_33
 484     CML_ASSIGN_MAT_44
 485
 486     /* Define class operators for external matrices. Note: external matrices
 487      * cannot be copy-constructed, but they can be assigned to:
 488      */
 489     CML_MAT_ASSIGN_FROM_MATTYPE
 490
 491     CML_MAT_ASSIGN_FROM_MAT(=, et::OpAssign)
 492     CML_MAT_ASSIGN_FROM_MAT(+=, et::OpAddAssign)
 493     CML_MAT_ASSIGN_FROM_MAT(-=, et::OpSubAssign)
 494
 495     CML_MAT_ASSIGN_FROM_MATXPR(=, et::OpAssign)
 496     CML_MAT_ASSIGN_FROM_MATXPR(+=, et::OpAddAssign)
 497     CML_MAT_ASSIGN_FROM_MATXPR(-=, et::OpSubAssign)
 498
 499     CML_MAT_ASSIGN_FROM_SCALAR(*=, et::OpMulAssign)
 500     CML_MAT_ASSIGN_FROM_SCALAR(/=, et::OpDivAssign)
 501
 502     CML_ACCUMULATED_MATRIX_MULT(const matrix_type&)
 503
 504     template<typename E, class AT, typename BO, typename L>
 505         CML_ACCUMULATED_MATRIX_MULT(const TEMPLATED_MATRIX_MACRO&)
 506
 507     template<class XprT>
 508         CML_ACCUMULATED_MATRIX_MULT(MATXPR_ARG_TYPE)
 509
 510
 511   protected:
 512
 513     value_type basis_element(size_t i, size_t j, row_basis) const {
 514         return (*this)(i,j);
 515     }
 516
 517     value_type basis_element(size_t i, size_t j, col_basis) const {
 518         return (*this)(j,i);
 519     }
 520
 521     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s, row_basis) {
 522         (*this)(i,j) = s;
 523     }
 524
 525     void set_basis_element(size_t i, size_t j, ELEMENT_ARG_TYPE s, col_basis) {
 526         (*this)(j,i) = s;
 527     }
 528
 529
 530   public:
 531
 532     /* Braces should only be used for testing: */
 533 #if defined(CML_ENABLE_MATRIX_BRACES)
 534     CML_MATRIX_BRACE_OPERATORS
 535 #endif
 536 };
 537
 538 } // namespace cml
 539
 540 #endif
 541
 542 // -------------------------------------------------------------------------
 543 // vim:ft=cpp