1 /* -*- C++ -*- ------------------------------------------------------------
3 Copyright (c) 2007 Jesse Anders and Demian Nave http://cmldev.net/
5 The Configurable Math Library (CML) is distributed under the terms of the
6 Boost Software License, v1.0 (see cml/LICENSE for details).
8 *-----------------------------------------------------------------------*/
13 #ifndef external_matrix_h
14 #define external_matrix_h
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>
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
>
32 /* Shorthand for the generator: */
33 typedef external
<Rows
,Cols
> generator_type
;
35 /* Shorthand for the array type: */
36 typedef external_2D
<Element
,Rows
,Cols
,Layout
> array_type
;
38 /* Shorthand for the type of this matrix: */
39 typedef matrix
<Element
,generator_type
,BasisOrient
,Layout
> matrix_type
;
41 /* For integration into the expression template code: */
42 typedef matrix_type expr_type
;
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.
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
;
55 typedef matrix_type
& expr_reference
;
56 typedef const matrix_type
& expr_const_reference
;
58 /* For matching by basis: */
59 typedef BasisOrient basis_orient
;
61 /* For matching by memory layout: */
62 typedef typename
array_type::layout layout
;
64 /* For matching by storage type if necessary: */
65 typedef typename
array_type::memory_tag memory_tag
;
67 /* For matching by size type if necessary: */
68 typedef typename
array_type::size_tag size_tag
;
70 /* For matching by resizability: */
71 typedef typename
array_type::resizing_tag resizing_tag
;
73 /* For matching by result-type: */
74 typedef cml::et::matrix_result_tag result_tag
;
76 /* For matching by assignability: */
77 typedef cml::et::assignable_tag assignable_tag
;
79 /* To simplify the matrix transpose operator: */
81 typename
cml::remove_const
<Element
>::type
,
82 typename
array_type::transposed_type::generator_type
,
87 /* To simplify the matrix row and column operators: */
90 typename
array_type::row_array_type::generator_type
95 typename
array_type::col_array_type::generator_type
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));
108 /** Set this matrix to the identity.
110 * This only makes sense for a square matrix, but no error will be
111 * signaled if the matrix is not square.
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);
122 /** Set this matrix to its transpose.
124 * This only makes sense for a square matrix, but no error will be
125 * signaled if the matrix is not square.
127 matrix_type
& transpose() {
128 /* transpose() returns a temporary: */
129 *this = transpose(*this);
133 /** Set this matrix to its inverse.
135 * This only makes sense for a square matrix, but no error will be
136 * signaled if the matrix is not square.
138 matrix_type
& inverse() {
139 /* inverse() returns a temporary: */
140 *this = cml::inverse(*this);
144 /* NOTE: minimize() and maximize() no longer supported (Jesse) */
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
));
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
));
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
);
182 /** Constructor for fixed-size external matrices.
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.
188 * @param ptr specify the external pointer.
190 * @throws same as the ArrayType constructor.
192 explicit matrix(value_type ptr
[Rows
][Cols
]) : array_type(ptr
) {}
194 /** Constructor for fixed-size external matrices.
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.
200 * @param ptr specify the external pointer.
202 * @throws same as the ArrayType constructor.
204 explicit matrix(value_type
* ptr
) : array_type(ptr
) {}
209 /** Return the matrix size as a pair. */
210 matrix_size
size() const {
211 return matrix_size(this->rows(),this->cols());
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());
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());
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
];
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
];
241 /* Define class operators for external matrices. Note: external matrices
242 * cannot be copy-constructed, but they can be assigned to:
244 CML_MAT_ASSIGN_FROM_MATTYPE
246 CML_MAT_ASSIGN_FROM_MAT(=, et::OpAssign
)
247 CML_MAT_ASSIGN_FROM_MAT(+=, et::OpAddAssign
)
248 CML_MAT_ASSIGN_FROM_MAT(-=, et::OpSubAssign
)
250 CML_MAT_ASSIGN_FROM_MATXPR(=, et::OpAssign
)
251 CML_MAT_ASSIGN_FROM_MATXPR(+=, et::OpAddAssign
)
252 CML_MAT_ASSIGN_FROM_MATXPR(-=, et::OpSubAssign
)
254 CML_MAT_ASSIGN_FROM_SCALAR(*=, et::OpMulAssign
)
255 CML_MAT_ASSIGN_FROM_SCALAR(/=, et::OpDivAssign
)
257 CML_ACCUMULATED_MATRIX_MULT(const matrix_type
&)
259 template<typename E
, class AT
, typename BO
, typename L
>
260 CML_ACCUMULATED_MATRIX_MULT(const TEMPLATED_MATRIX_MACRO
&)
263 CML_ACCUMULATED_MATRIX_MULT(MATXPR_ARG_TYPE
)
268 value_type
basis_element(size_t i
, size_t j
, row_basis
) const {
272 value_type
basis_element(size_t i
, size_t j
, col_basis
) const {
276 void set_basis_element(size_t i
, size_t j
, ELEMENT_ARG_TYPE s
, row_basis
) {
280 void set_basis_element(size_t i
, size_t j
, ELEMENT_ARG_TYPE s
, col_basis
) {
287 /* Braces should only be used for testing: */
288 #if defined(CML_ENABLE_MATRIX_BRACES)
289 CML_MATRIX_BRACE_OPERATORS
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
>
300 /* Shorthand for the generator: */
301 typedef external
<> generator_type
;
303 /* Shorthand for the array type: */
304 typedef external_2D
<Element
,-1,-1,Layout
> array_type
;
306 /* Shorthand for the type of this matrix: */
307 typedef matrix
<Element
,generator_type
,BasisOrient
,Layout
> matrix_type
;
309 /* For integration into the expression template code: */
310 typedef matrix_type expr_type
;
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.
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
;
323 typedef matrix_type
& expr_reference
;
324 typedef const matrix_type
& expr_const_reference
;
326 /* For matching by basis: */
327 typedef BasisOrient basis_orient
;
329 /* For matching by memory layout: */
330 typedef typename
array_type::layout layout
;
332 /* For matching by storage type if necessary: */
333 typedef typename
array_type::memory_tag memory_tag
;
335 /* For matching by size type if necessary: */
336 typedef typename
array_type::size_tag size_tag
;
338 /* For matching by resizability: */
339 typedef typename
array_type::resizing_tag resizing_tag
;
341 /* For matching by result-type: */
342 typedef cml::et::matrix_result_tag result_tag
;
344 /* For matching by assignability: */
345 typedef cml::et::assignable_tag assignable_tag
;
347 /* To simplify the matrix transpose operator: */
350 typename
array_type::transposed_type::generator_type
,
355 /* To simplify the matrix row and column operators: */
358 typename
array_type::row_array_type::generator_type
363 typename
array_type::col_array_type::generator_type
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));
376 /** Set this matrix to the identity.
378 * This only makes sense for a square matrix, but no error will be
379 * signaled if the matrix is not square.
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);
390 /** Set this matrix to its transpose.
392 * This only makes sense for a square matrix, but no error will be
393 * signaled if the matrix is not square.
395 matrix_type
& transpose() {
396 /* transpose() returns a temporary: */
397 *this = cml::transpose(*this);
401 /** Set this matrix to its inverse.
403 * This only makes sense for a square matrix, but no error will be
404 * signaled if the matrix is not square.
406 matrix_type
& inverse() {
407 /* inverse() returns a temporary: */
408 *this = inverse(*this);
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
));
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
));
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
);
446 /** Constructor for fixed-size external matrices.
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.
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.
456 * @throws same as the ArrayType constructor.
458 explicit matrix(value_type
* const ptr
, size_t rows
, size_t cols
)
459 : array_type(ptr
,rows
,cols
) {}
464 /** Return the matrix size as a pair. */
465 matrix_size
size() const {
466 return matrix_size(this->rows(),this->cols());
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());
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());
486 /* Define class operators for external matrices. Note: external matrices
487 * cannot be copy-constructed, but they can be assigned to:
489 CML_MAT_ASSIGN_FROM_MATTYPE
491 CML_MAT_ASSIGN_FROM_MAT(=, et::OpAssign
)
492 CML_MAT_ASSIGN_FROM_MAT(+=, et::OpAddAssign
)
493 CML_MAT_ASSIGN_FROM_MAT(-=, et::OpSubAssign
)
495 CML_MAT_ASSIGN_FROM_MATXPR(=, et::OpAssign
)
496 CML_MAT_ASSIGN_FROM_MATXPR(+=, et::OpAddAssign
)
497 CML_MAT_ASSIGN_FROM_MATXPR(-=, et::OpSubAssign
)
499 CML_MAT_ASSIGN_FROM_SCALAR(*=, et::OpMulAssign
)
500 CML_MAT_ASSIGN_FROM_SCALAR(/=, et::OpDivAssign
)
502 CML_ACCUMULATED_MATRIX_MULT(const matrix_type
&)
504 template<typename E
, class AT
, typename BO
, typename L
>
505 CML_ACCUMULATED_MATRIX_MULT(const TEMPLATED_MATRIX_MACRO
&)
508 CML_ACCUMULATED_MATRIX_MULT(MATXPR_ARG_TYPE
)
513 value_type
basis_element(size_t i
, size_t j
, row_basis
) const {
517 value_type
basis_element(size_t i
, size_t j
, col_basis
) const {
521 void set_basis_element(size_t i
, size_t j
, ELEMENT_ARG_TYPE s
, row_basis
) {
525 void set_basis_element(size_t i
, size_t j
, ELEMENT_ARG_TYPE s
, col_basis
) {
532 /* Braces should only be used for testing: */
533 #if defined(CML_ENABLE_MATRIX_BRACES)
534 CML_MATRIX_BRACE_OPERATORS
542 // -------------------------------------------------------------------------