glucat 0.12.0
matrix_imp.h
Go to the documentation of this file.
1#ifndef _GLUCAT_MATRIX_IMP_H
2#define _GLUCAT_MATRIX_IMP_H
3/***************************************************************************
4 GluCat : Generic library of universal Clifford algebra templates
5 matrix_imp.h : Implement common matrix functions
6 -------------------
7 begin : Sun 2001-12-09
8 copyright : (C) 2001-2012 by Paul C. Leopardi
9 : uBLAS interface contributed by Joerg Walter
10 ***************************************************************************
11
12 This library is free software: you can redistribute it and/or modify
13 it under the terms of the GNU Lesser General Public License as published
14 by the Free Software Foundation, either version 3 of the License, or
15 (at your option) any later version.
16
17 This library is distributed in the hope that it will be useful,
18 but WITHOUT ANY WARRANTY; without even the implied warranty of
19 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 GNU Lesser General Public License for more details.
21
22 You should have received a copy of the GNU Lesser General Public License
23 along with this library. If not, see <http://www.gnu.org/licenses/>.
24
25 ***************************************************************************
26 This library is based on a prototype written by Arvind Raja and was
27 licensed under the LGPL with permission of the author. See Arvind Raja,
28 "Object-oriented implementations of Clifford algebras in C++: a prototype",
29 in Ablamowicz, Lounesto and Parra (eds.)
30 "Clifford algebras with numeric and symbolic computations", Birkhauser, 1996.
31 ***************************************************************************
32 See also Arvind Raja's original header comments in glucat.h
33 ***************************************************************************/
34
35#include "glucat/errors.h"
36#include "glucat/scalar.h"
37#include "glucat/matrix.h"
38
39# if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
40# pragma GCC diagnostic push
41# pragma GCC diagnostic ignored "-Wunused-local-typedefs"
42# endif
43# if defined(_GLUCAT_HAVE_BOOST_SERIALIZATION_ARRAY_WRAPPER_H)
44# include <boost/serialization/array_wrapper.hpp>
45# endif
46#include <boost/numeric/ublas/vector.hpp>
47#include <boost/numeric/ublas/vector_proxy.hpp>
48#include <boost/numeric/ublas/matrix.hpp>
49#include <boost/numeric/ublas/matrix_expression.hpp>
50#include <boost/numeric/ublas/matrix_proxy.hpp>
51#include <boost/numeric/ublas/matrix_sparse.hpp>
52#include <boost/numeric/ublas/operation.hpp>
53#include <boost/numeric/ublas/operation_sparse.hpp>
54
55#if defined(_GLUCAT_USE_BINDINGS)
56# include <boost/numeric/bindings/lapack/driver/gees.hpp>
57# include <boost/numeric/bindings/ublas.hpp>
58#endif
59
60#if defined(_GLUCAT_USE_BLAZE)
61#include <blaze/Math.h>
62#include <blaze/math/DynamicMatrix.h>
63#include <blaze/math/DynamicVector.h>
64#endif
65
66# if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
67# pragma GCC diagnostic pop
68# endif
69
70#include <set>
71#include <vector>
72
73namespace glucat { namespace matrix
74{
75 // References for algorithms:
76 // [v]: C. F. van Loan and N. Pitsianis, "Approximation with Kronecker products",
77 // in Linear Algebra for Large Scale and Real-Time Applications, Marc S. Moonen,
78 // Gene H. Golub, and Bart L. R. Moor (eds.), 1993, pp. 293--314.
79
81 template< typename LHS_T, typename RHS_T >
82 auto
83 kron(const LHS_T& lhs, const RHS_T& rhs) -> const
84 RHS_T
85 {
86 const auto rhs_s1 = rhs.size1();
87 const auto rhs_s2 = rhs.size2();
88 auto result = RHS_T(lhs.size1()*rhs_s1, lhs.size2()*rhs_s2);
89 result.clear();
90
91 for (auto
92 lhs_it1 = lhs.begin1();
93 lhs_it1 != lhs.end1();
94 ++lhs_it1)
95 for (auto
96 lhs_it2 = lhs_it1.begin();
97 lhs_it2 != lhs_it1.end();
98 ++lhs_it2)
99 {
100 const auto start1 = rhs_s1 * lhs_it2.index1();
101 const auto start2 = rhs_s2 * lhs_it2.index2();
102 const auto& lhs_val = *lhs_it2;
103 for (auto
104 rhs_it1 = rhs.begin1();
105 rhs_it1 != rhs.end1();
106 ++rhs_it1)
107 for (auto
108 rhs_it2 = rhs_it1.begin();
109 rhs_it2 != rhs_it1.end();
110 ++rhs_it2)
111 result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
112 }
113 return result;
114 }
115
117 template< typename LHS_T, typename RHS_T >
118 auto
119 mono_kron(const LHS_T& lhs, const RHS_T& rhs) -> const
120 RHS_T
121 {
122 const auto rhs_s1 = rhs.size1();
123 const auto rhs_s2 = rhs.size2();
124 const auto dim = lhs.size1()*rhs_s1;
125 auto result = RHS_T(dim, dim, dim);
126 result.clear();
127
128 for (auto
129 lhs_it1 = lhs.begin1();
130 lhs_it1 != lhs.end1();
131 ++lhs_it1)
132 {
133 const auto lhs_it2 = lhs_it1.begin();
134 const auto start1 = rhs_s1 * lhs_it2.index1();
135 const auto start2 = rhs_s2 * lhs_it2.index2();
136 const auto& lhs_val = *lhs_it2;
137 for (auto
138 rhs_it1 = rhs.begin1();
139 rhs_it1 != rhs.end1();
140 ++rhs_it1)
141 {
142 const auto rhs_it2 = rhs_it1.begin();
143 result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
144 }
145 }
146 return result;
147 }
148
150 template< typename LHS_T, typename RHS_T >
151 void
152 nork_range(RHS_T& result,
153 const typename LHS_T::const_iterator2 lhs_it2,
154 const RHS_T& rhs,
155 const typename RHS_T::size_type res_s1,
156 const typename RHS_T::size_type res_s2)
157 {
158 // Definition matches [v] Section 4, Theorem 4.1.
159 const auto start1 = res_s1 * lhs_it2.index1();
160 const auto start2 = res_s2 * lhs_it2.index2();
161 using ublas::range;
162 const auto& range1 = range(start1, start1 + res_s1);
163 const auto& range2 = range(start2, start2 + res_s2);
164 using matrix_range_t = ublas::matrix_range<const RHS_T>;
165 const auto& rhs_range = matrix_range_t(rhs, range1, range2);
166 using Scalar_T = typename RHS_T::value_type;
167 const auto lhs_val = numeric_traits<Scalar_T>::to_scalar_t(*lhs_it2);
168 for (auto
169 rhs_it1 = rhs_range.begin1();
170 rhs_it1 != rhs_range.end1();
171 ++rhs_it1)
172 for (auto
173 rhs_it2 = rhs_it1.begin();
174 rhs_it2 != rhs_it1.end();
175 ++rhs_it2)
176 result(rhs_it2.index1(), rhs_it2.index2()) += lhs_val * *rhs_it2;
177 }
178
180 template< typename LHS_T, typename RHS_T >
181 auto
182 nork(const LHS_T& lhs, const RHS_T& rhs, const bool mono) -> const
183 RHS_T
184 {
185 // nork(A, kron(A, B)) is close to B
186 // Definition matches [v] Section 4, Theorem 4.1.
187 const auto lhs_s1 = lhs.size1();
188 const auto lhs_s2 = lhs.size2();
189 const auto rhs_s1 = rhs.size1();
190 const auto rhs_s2 = rhs.size2();
191 const auto res_s1 = rhs_s1 / lhs_s1;
192 const auto res_s2 = rhs_s2 / lhs_s2;
193 using Scalar_T = typename RHS_T::value_type;
194 const auto norm_frob2_lhs = norm_frob2(lhs);
195 if (!mono)
196 {
197 using error_t = error<RHS_T>;
198 if (rhs_s1 == 0)
199 throw error_t("matrix", "nork: number of rows must not be 0");
200 if (rhs_s2 == 0)
201 throw error_t("matrix", "nork: number of cols must not be 0");
202 if (res_s1 * lhs_s1 != rhs_s1)
203 throw error_t("matrix", "nork: incompatible numbers of rows");
204 if (res_s2 * lhs_s2 != rhs_s2)
205 throw error_t("matrix", "nork: incompatible numbers of cols");
206 if (norm_frob2_lhs == Scalar_T(0))
207 throw error_t("matrix", "nork: LHS must not be 0");
208 }
209 auto result = RHS_T(res_s1, res_s2);
210 result.clear();
211 for (auto
212 lhs_it1 = lhs.begin1();
213 lhs_it1 != lhs.end1();
214 ++lhs_it1)
215 for (auto
216 lhs_it2 = lhs_it1.begin();
217 lhs_it2 != lhs_it1.end();
218 ++lhs_it2)
219 if (*lhs_it2 != Scalar_T(0))
220 nork_range<LHS_T, RHS_T>(result, lhs_it2, rhs, res_s1, res_s2);
221 result /= norm_frob2_lhs;
222 return result;
223 }
224
226 template< typename LHS_T, typename RHS_T >
227 auto
228 signed_perm_nork(const LHS_T& lhs, const RHS_T& rhs) -> const
229 RHS_T
230 {
231 // signed_perm_nork(A, kron(A, B)) is close to B
232 // Definition matches [v] Section 4, Theorem 4.1.
233 const auto lhs_s1 = lhs.size1();
234 const auto lhs_s2 = lhs.size2();
235 const auto rhs_s1 = rhs.size1();
236 const auto rhs_s2 = rhs.size2();
237 const auto res_s1 = rhs_s1 / lhs_s1;
238 const auto res_s2 = rhs_s2 / lhs_s2;
239 using Scalar_T = typename RHS_T::value_type;
240 const auto norm_frob2_lhs = Scalar_T( double(lhs_s1) );
241 auto result = RHS_T(res_s1, res_s2);
242 result.clear();
243 for (auto
244 lhs_it1 = lhs.begin1();
245 lhs_it1 != lhs.end1();
246 ++lhs_it1)
247 {
248 const auto lhs_it2 = lhs_it1.begin();
249 nork_range<LHS_T, RHS_T>(result, lhs_it2, rhs, res_s1, res_s2);
250 }
251 result /= norm_frob2_lhs;
252 return result;
253 }
254
256 template< typename Matrix_T >
257 auto
258 nnz(const Matrix_T& m) -> typename Matrix_T::size_type
259 {
260 using size_t = typename Matrix_T::size_type;
261 auto result = size_t(0);
262 for (auto
263 it1 = m.begin1();
264 it1 != m.end1();
265 ++it1)
266 for (auto& entry : it1)
267 if (entry != 0)
268 ++result;
269 return result;
270 }
271
273 template< typename Matrix_T >
274 auto
275 isinf(const Matrix_T& m) -> bool
276 {
277 using Scalar_T = typename Matrix_T::value_type;
278 for (auto
279 it1 = m.begin1();
280 it1 != m.end1();
281 ++it1)
282 for (auto& entry : it1)
284 return true;
285
286 return false;
287 }
288
290 template< typename Matrix_T >
291 auto
292 isnan(const Matrix_T& m) -> bool
293 {
294 using Scalar_T = typename Matrix_T::value_type;
295 for (auto
296 it1 = m.begin1();
297 it1 != m.end1();
298 ++it1)
299 for (auto& entry : it1)
301 return true;
302
303 return false;
304 }
305
307 template< typename Matrix_T >
308 inline
309 auto
310 unit(const typename Matrix_T::size_type dim) -> const
311 Matrix_T
312 {
313 using Scalar_T = typename Matrix_T::value_type;
314 return ublas::identity_matrix<Scalar_T>(dim);
315 }
316
318 template< typename LHS_T, typename RHS_T >
319 auto
320 mono_prod(const ublas::matrix_expression<LHS_T>& lhs,
321 const ublas::matrix_expression<RHS_T>& rhs) -> const typename RHS_T::expression_type
322 {
323 using rhs_expression_t = const RHS_T;
324 using matrix_row_t = typename ublas::matrix_row<rhs_expression_t>;
325
326 const auto dim = lhs().size1();
327 // The following assumes that RHS_T is a sparse matrix type.
328 auto result = RHS_T(dim, dim, dim);
329 for (auto
330 lhs_row = lhs().begin1();
331 lhs_row != lhs().end1();
332 ++lhs_row)
333 {
334 const auto& lhs_it = lhs_row.begin();
335 if (lhs_it != lhs_row.end())
336 {
337 const auto& rhs_row = matrix_row_t(rhs(), lhs_it.index2());
338 const auto& rhs_it = rhs_row.begin();
339 if (rhs_it != rhs_row.end())
340 result(lhs_it.index1(), rhs_it.index()) = (*lhs_it) * (*rhs_it);
341 }
342 }
343 return result;
344 }
345
347 template< typename LHS_T, typename RHS_T >
348 inline
349 auto
350 sparse_prod(const ublas::matrix_expression<LHS_T>& lhs,
351 const ublas::matrix_expression<RHS_T>& rhs) -> const typename RHS_T::expression_type
352 {
353 using expression_t = typename RHS_T::expression_type;
354 return ublas::sparse_prod<expression_t>(lhs(), rhs());
355 }
356
358 template< typename LHS_T, typename RHS_T >
359 inline
360 auto
361 prod(const ublas::matrix_expression<LHS_T>& lhs,
362 const ublas::matrix_expression<RHS_T>& rhs) -> const typename RHS_T::expression_type
363 {
364 const auto dim = lhs().size1();
365 RHS_T result(dim, dim);
366 ublas::axpy_prod(lhs, rhs, result, true);
367 return result;
368 }
369
371 template< typename Scalar_T, typename LHS_T, typename RHS_T >
372 auto
373 inner(const LHS_T& lhs, const RHS_T& rhs) -> Scalar_T
374 {
375 auto result = Scalar_T(0);
376 for (auto
377 lhs_it1 = lhs.begin1();
378 lhs_it1 != lhs.end1();
379 ++lhs_it1)
380 for (auto
381 lhs_it2 = lhs_it1.begin();
382 lhs_it2 != lhs_it1.end();
383 ++lhs_it2)
384 {
385 const auto& rhs_val = rhs(lhs_it2.index1(),lhs_it2.index2());
386 if (rhs_val != Scalar_T(0))
387 result += (*lhs_it2) * rhs_val;
388 }
389 return result / lhs.size1();
390 }
391
393 template< typename Matrix_T >
394 auto
395 norm_frob2(const Matrix_T& val) -> typename Matrix_T::value_type
396 {
397 using Scalar_T = typename Matrix_T::value_type;
398
399 auto result = Scalar_T(0);
400 for (auto
401 val_it1 = val.begin1();
402 val_it1 != val.end1();
403 ++val_it1)
404 for (auto& val_entry : val_it1)
405 {
406 if (numeric_traits<Scalar_T>::isNaN(val_entry))
408 result += val_entry * val_entry;
409 }
410 return result;
411 }
412
414 template< typename Matrix_T >
415 auto
416 trace(const Matrix_T& val) -> typename Matrix_T::value_type
417 {
418 using Scalar_T = typename Matrix_T::value_type;
419
420 auto result = Scalar_T(0);
421 auto dim = val.size1();
422 for (auto
423 ndx = decltype(dim)(0);
424 ndx != dim;
425 ++ndx)
426 {
427 const Scalar_T crd = val(ndx, ndx);
430 result += crd;
431 }
432 return result;
433 }
434
435#if defined(_GLUCAT_USE_BINDINGS)
437 template< typename Matrix_T >
438 static
439 auto
440 to_lapack(const Matrix_T& val) -> ublas::matrix<double, ublas::column_major>
441 {
442 const auto s1 = val.size1();
443 const auto s2 = val.size2();
444
445 using lapack_matrix_t = typename ublas::matrix<double, ublas::column_major>;
446 auto result = lapack_matrix_t(s1, s2);
447 result.clear();
448
449 using Scalar_T = typename Matrix_T::value_type;
450 using traits_t = numeric_traits<Scalar_T>;
451
452 for (auto
453 val_it1 = val.begin1();
454 val_it1 != val.end1();
455 ++val_it1)
456 for (auto
457 val_it2 = val_it1.begin();
458 val_it2 != val_it1.end();
459 ++val_it2)
460 result(val_it2.index1(), val_it2.index2()) = traits_t::to_double(*val_it2);
461
462 return result;
463 }
464#endif
465
466#if defined(_GLUCAT_USE_BLAZE)
468 template< typename Matrix_T >
469 static
470 auto
471 to_blaze(const Matrix_T& val) -> blaze::DynamicMatrix<double,blaze::rowMajor>
472 {
473 const auto s1 = val.size1();
474 const auto s2 = val.size2();
475
476 using blaze_matrix_t = typename blaze::DynamicMatrix<double,blaze::rowMajor>;
477 auto result = blaze_matrix_t(s1, s2);
478
479 using Scalar_T = typename Matrix_T::value_type;
480 using traits_t = numeric_traits<Scalar_T>;
481
482 for (auto
483 val_it1 = val.begin1();
484 val_it1 != val.end1();
485 ++val_it1)
486 for (auto
487 val_it2 = val_it1.begin();
488 val_it2 != val_it1.end();
489 ++val_it2)
490 result(val_it2.index1(), val_it2.index2()) = traits_t::to_double(*val_it2);
491
492 return result;
493 }
494
495#endif
496
498 template< typename Matrix_T >
499 auto
500 eigenvalues(const Matrix_T& val) -> std::vector< std::complex<double> >
501 {
502 using complex_t = std::complex<double>;
503 using complex_vector_t = typename std::vector<complex_t>;
504
505 const auto dim = val.size1();
506 auto lambda = complex_vector_t(dim);
507
508#if defined(_GLUCAT_USE_BINDINGS)
509 namespace lapack = boost::numeric::bindings::lapack;
510 using lapack_matrix_t = typename ublas::matrix<double, ublas::column_major>;
511
512 auto T = to_lapack(val);
513 auto V = T;
514 using vector_t = typename ublas::vector<double>;
515 auto real_lambda = vector_t(dim);
516 auto imag_lambda = vector_t(dim);
517 fortran_int_t sdim = 0;
518
519 lapack::gees('N', 'N', nullptr, T, sdim, real_lambda, imag_lambda, V);
520
521 for (auto
522 k = decltype(dim)(0);
523 k != dim;
524 ++k)
525 lambda[k] = complex_t(real_lambda[k], imag_lambda[k]);
526#endif
527#if defined(_GLUCAT_USE_BLAZE)
528 using blaze_matrix_t = typename blaze::DynamicMatrix<double, blaze::rowMajor>;
529 using complex_t = std::complex<double>;
530 using blaze_complex_vector_t = blaze::DynamicVector<complex_t, blaze::columnVector>;
531
532 auto blaze_val = to_blaze(val);
533 auto blaze_lambda = blaze_complex_vector_t(dim);
534 blaze::geev(blaze_val, blaze_lambda);
535
536 for (auto
537 k = decltype(dim)(0);
538 k != dim;
539 ++k)
540 lambda[k] = blaze_lambda[k];
541#endif
542 return lambda;
543 }
544
546 template< typename Matrix_T >
547 auto
549 {
550 using Scalar_T = typename Matrix_T::value_type;
551 eig_genus<Matrix_T> result;
552
553 using complex_t = std::complex<double>;
554 using complex_vector_t = typename std::vector<complex_t>;
555 auto lambda = eigenvalues(val);
556
557 std::set<double> arg_set;
558
559 using vector_index_t = typename complex_vector_t::size_type;
560 const auto dim = lambda.size();
561 static const auto epsilon =
562 std::max(std::numeric_limits<double>::epsilon(),
563 numeric_traits<Scalar_T>::to_double(std::numeric_limits<Scalar_T>::epsilon()));
564 static const auto zero_eig_tol = 4096.0*epsilon;
565
566 bool neg_real_eig_found = false;
567 bool imag_eig_found = false;
568 bool zero_eig_found = false;
569
570 for (auto
571 k = decltype(dim)(0);
572 k != dim;
573 ++k)
574 {
575 const auto lambda_k = lambda[k];
576 arg_set.insert(std::arg(lambda_k));
577
578 const auto real_lambda_k = std::real(lambda_k);
579 const auto imag_lambda_k = std::imag(lambda_k);
580 const auto norm_tol = 4096.0*epsilon*std::norm(lambda_k);
581
582 if (!neg_real_eig_found &&
583 real_lambda_k < -epsilon &&
584 (imag_lambda_k == 0.0 ||
585 imag_lambda_k * imag_lambda_k < norm_tol))
586 neg_real_eig_found = true;
587 if (!imag_eig_found &&
588 imag_lambda_k > epsilon &&
589 (real_lambda_k == 0.0 ||
590 real_lambda_k * real_lambda_k < norm_tol))
591 imag_eig_found = true;
592 if (!zero_eig_found &&
593 std::norm(lambda_k) < zero_eig_tol)
594 zero_eig_found = true;
595 }
596
597 if (zero_eig_found)
598 result.m_is_singular = true;
599
600 static const auto pi = numeric_traits<double>::pi();
601 if (neg_real_eig_found)
602 {
603 if (imag_eig_found)
604 result.m_eig_case = both_eigs;
605 else
606 {
607 result.m_eig_case = neg_real_eigs;
608 result.m_safe_arg = Scalar_T(-pi / 2.0);
609 }
610 }
611
612 if (result.m_eig_case == both_eigs)
613 {
614 auto arg_it = arg_set.begin();
615 auto first_arg = *arg_it;
616 auto best_arg = first_arg;
617 auto best_diff = 0.0;
618 auto previous_arg = first_arg;
619 for (++arg_it;
620 arg_it != arg_set.end();
621 ++arg_it)
622 {
623 const auto arg_diff = *arg_it - previous_arg;
624 if (arg_diff > best_diff)
625 {
626 best_diff = arg_diff;
627 best_arg = previous_arg;
628 }
629 previous_arg = *arg_it;
630 }
631 const auto arg_diff = first_arg + 2.0*pi - previous_arg;
632 if (arg_diff > best_diff)
633 {
634 best_diff = arg_diff;
635 best_arg = previous_arg;
636 }
637 result.m_safe_arg = Scalar_T(pi - (best_arg + best_diff / 2.0));
638 }
639 return result;
640 }
641} }
642
643#endif // _GLUCAT_MATRIX_IMP_H
const scalar_t epsilon
Definition PyClical.h:150
Specific exception class.
Definition errors.h:57
Extra traits which extend numeric limits.
Definition scalar.h:48
static auto isNaN(const Scalar_T &val, bool_to_type< false >) -> bool
Smart isnan specialised for Scalar_T without quiet NaN.
Definition scalar.h:68
static auto NaN() -> Scalar_T
Smart NaN.
Definition scalar.h:115
static auto to_scalar_t(const Other_Scalar_T &val) -> Scalar_T
Cast to Scalar_T.
Definition scalar.h:141
static auto to_double(const Scalar_T &val) -> double
Cast to double.
Definition scalar.h:133
static auto isInf(const Scalar_T &val, bool_to_type< false >) -> bool
Smart isinf specialised for Scalar_T without infinity.
Definition scalar.h:54
static auto pi() -> Scalar_T
Pi.
Definition scalar.h:189
static auto to_lapack(const Matrix_T &val) -> ublas::matrix< double, ublas::column_major >
Convert matrix to LAPACK format.
Definition matrix_imp.h:440
auto nork(const LHS_T &lhs, const RHS_T &rhs, const bool mono=true) -> const RHS_T
Left inverse of Kronecker product.
Definition matrix_imp.h:182
auto inner(const LHS_T &lhs, const RHS_T &rhs) -> Scalar_T
Inner product: sum(x(i,j)*y(i,j))/x.nrows()
Definition matrix_imp.h:373
auto signed_perm_nork(const LHS_T &lhs, const RHS_T &rhs) -> const RHS_T
Left inverse of Kronecker product where lhs is a signed permutation matrix.
Definition matrix_imp.h:228
auto trace(const Matrix_T &val) -> typename Matrix_T::value_type
Matrix trace.
Definition matrix_imp.h:416
auto isinf(const Matrix_T &m) -> bool
Infinite.
Definition matrix_imp.h:275
auto nnz(const Matrix_T &m) -> typename Matrix_T::size_type
Number of non-zeros.
Definition matrix_imp.h:258
auto mono_kron(const LHS_T &lhs, const RHS_T &rhs) -> const RHS_T
Sparse Kronecker tensor product of monomial matrices.
Definition matrix_imp.h:119
auto eigenvalues(const Matrix_T &val) -> std::vector< std::complex< double > >
Eigenvalues of a matrix.
Definition matrix_imp.h:500
auto norm_frob2(const Matrix_T &val) -> typename Matrix_T::value_type
Square of Frobenius norm.
Definition matrix_imp.h:395
auto classify_eigenvalues(const Matrix_T &val) -> eig_genus< Matrix_T >
Classify the eigenvalues of a matrix.
Definition matrix_imp.h:548
auto prod(const ublas::matrix_expression< LHS_T > &lhs, const ublas::matrix_expression< RHS_T > &rhs) -> const typename RHS_T::expression_type
Product of matrices.
Definition matrix_imp.h:361
auto kron(const LHS_T &lhs, const RHS_T &rhs) -> const RHS_T
Kronecker tensor product of matrices - as per Matlab kron.
Definition matrix_imp.h:83
auto mono_prod(const ublas::matrix_expression< LHS_T > &lhs, const ublas::matrix_expression< RHS_T > &rhs) -> const typename RHS_T::expression_type
Product of monomial matrices.
Definition matrix_imp.h:320
auto sparse_prod(const ublas::matrix_expression< LHS_T > &lhs, const ublas::matrix_expression< RHS_T > &rhs) -> const typename RHS_T::expression_type
Product of sparse matrices.
Definition matrix_imp.h:350
auto unit(const typename Matrix_T::size_type n) -> const Matrix_T
Unit matrix - as per Matlab eye.
Definition matrix_imp.h:310
auto isnan(const Matrix_T &m) -> bool
Not a Number.
Definition matrix_imp.h:292
void nork_range(RHS_T &result, const typename LHS_T::const_iterator2 lhs_it2, const RHS_T &rhs, const typename RHS_T::size_type res_s1, const typename RHS_T::size_type res_s2)
Utility routine for nork: calculate result for a range of indices.
Definition matrix_imp.h:152
Structure containing classification of eigenvalues.
Definition matrix.h:141
bool m_is_singular
Is the matrix singular?
Definition matrix.h:144
Scalar_T m_safe_arg
Argument such that exp(pi-m_safe_arg) lies between arguments of eigenvalues.
Definition matrix.h:148
eig_case_t m_eig_case
What kind of eigenvalues does the matrix contain?
Definition matrix.h:146