glucat 0.13.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_BLAZE)
56#include <blaze/Math.h>
57#include <blaze/math/DynamicMatrix.h>
58#include <blaze/math/DynamicVector.h>
59#endif
60
61# if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
62# pragma GCC diagnostic pop
63# endif
64
65#include <set>
66#include <vector>
67
68namespace glucat { namespace matrix
69{
70 // References for algorithms:
71 // [v]: C. F. van Loan and N. Pitsianis, "Approximation with Kronecker products",
72 // in Linear Algebra for Large Scale and Real-Time Applications, Marc S. Moonen,
73 // Gene H. Golub, and Bart L. R. Moor (eds.), 1993, pp. 293--314.
74
76 template< typename LHS_T, typename RHS_T >
77 auto
78 kron(const LHS_T& lhs, const RHS_T& rhs) -> const
79 RHS_T
80 {
81 const auto rhs_s1 = rhs.size1();
82 const auto rhs_s2 = rhs.size2();
83 auto result = RHS_T(lhs.size1()*rhs_s1, lhs.size2()*rhs_s2);
84 result.clear();
85
86 for (auto
87 lhs_it1 = lhs.begin1();
88 lhs_it1 != lhs.end1();
89 ++lhs_it1)
90 for (auto
91 lhs_it2 = lhs_it1.begin();
92 lhs_it2 != lhs_it1.end();
93 ++lhs_it2)
94 {
95 const auto start1 = rhs_s1 * lhs_it2.index1();
96 const auto start2 = rhs_s2 * lhs_it2.index2();
97 const auto& lhs_val = *lhs_it2;
98 for (auto
99 rhs_it1 = rhs.begin1();
100 rhs_it1 != rhs.end1();
101 ++rhs_it1)
102 for (auto
103 rhs_it2 = rhs_it1.begin();
104 rhs_it2 != rhs_it1.end();
105 ++rhs_it2)
106 result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
107 }
108 return result;
109 }
110
112 template< typename LHS_T, typename RHS_T >
113 auto
114 mono_kron(const LHS_T& lhs, const RHS_T& rhs) -> const
115 RHS_T
116 {
117 const auto rhs_s1 = rhs.size1();
118 const auto rhs_s2 = rhs.size2();
119 const auto dim = lhs.size1()*rhs_s1;
120 auto result = RHS_T(dim, dim, dim);
121 result.clear();
122
123 for (auto
124 lhs_it1 = lhs.begin1();
125 lhs_it1 != lhs.end1();
126 ++lhs_it1)
127 {
128 const auto lhs_it2 = lhs_it1.begin();
129 const auto start1 = rhs_s1 * lhs_it2.index1();
130 const auto start2 = rhs_s2 * lhs_it2.index2();
131 const auto& lhs_val = *lhs_it2;
132 for (auto
133 rhs_it1 = rhs.begin1();
134 rhs_it1 != rhs.end1();
135 ++rhs_it1)
136 {
137 const auto rhs_it2 = rhs_it1.begin();
138 result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
139 }
140 }
141 return result;
142 }
143
145 template< typename LHS_T, typename RHS_T >
146 void
147 nork_range(RHS_T& result,
148 const typename LHS_T::const_iterator2 lhs_it2,
149 const RHS_T& rhs,
150 const typename RHS_T::size_type res_s1,
151 const typename RHS_T::size_type res_s2)
152 {
153 // Definition matches [v] Section 4, Theorem 4.1.
154 const auto start1 = res_s1 * lhs_it2.index1();
155 const auto start2 = res_s2 * lhs_it2.index2();
156 using ublas::range;
157 const auto& range1 = range(start1, start1 + res_s1);
158 const auto& range2 = range(start2, start2 + res_s2);
159 using matrix_range_t = ublas::matrix_range<const RHS_T>;
160 const auto& rhs_range = matrix_range_t(rhs, range1, range2);
161 using Scalar_T = typename RHS_T::value_type;
162 const auto lhs_val = numeric_traits<Scalar_T>::to_scalar_t(*lhs_it2);
163 for (auto
164 rhs_it1 = rhs_range.begin1();
165 rhs_it1 != rhs_range.end1();
166 ++rhs_it1)
167 for (auto
168 rhs_it2 = rhs_it1.begin();
169 rhs_it2 != rhs_it1.end();
170 ++rhs_it2)
171 result(rhs_it2.index1(), rhs_it2.index2()) += lhs_val * *rhs_it2;
172 }
173
175 template< typename LHS_T, typename RHS_T >
176 auto
177 nork(const LHS_T& lhs, const RHS_T& rhs, const bool mono) -> const
178 RHS_T
179 {
180 // nork(A, kron(A, B)) is close to B
181 // Definition matches [v] Section 4, Theorem 4.1.
182 const auto lhs_s1 = lhs.size1();
183 const auto lhs_s2 = lhs.size2();
184 const auto rhs_s1 = rhs.size1();
185 const auto rhs_s2 = rhs.size2();
186 const auto res_s1 = rhs_s1 / lhs_s1;
187 const auto res_s2 = rhs_s2 / lhs_s2;
188 using Scalar_T = typename RHS_T::value_type;
189 const auto norm_frob2_lhs = norm_frob2(lhs);
190 if (!mono)
191 {
192 using error_t = error<RHS_T>;
193 if (rhs_s1 == 0)
194 throw error_t("matrix", "nork: number of rows must not be 0");
195 if (rhs_s2 == 0)
196 throw error_t("matrix", "nork: number of cols must not be 0");
197 if (res_s1 * lhs_s1 != rhs_s1)
198 throw error_t("matrix", "nork: incompatible numbers of rows");
199 if (res_s2 * lhs_s2 != rhs_s2)
200 throw error_t("matrix", "nork: incompatible numbers of cols");
201 if (norm_frob2_lhs == Scalar_T(0))
202 throw error_t("matrix", "nork: LHS must not be 0");
203 }
204 auto result = RHS_T(res_s1, res_s2);
205 result.clear();
206 for (auto
207 lhs_it1 = lhs.begin1();
208 lhs_it1 != lhs.end1();
209 ++lhs_it1)
210 for (auto
211 lhs_it2 = lhs_it1.begin();
212 lhs_it2 != lhs_it1.end();
213 ++lhs_it2)
214 if (*lhs_it2 != Scalar_T(0))
215 nork_range<LHS_T, RHS_T>(result, lhs_it2, rhs, res_s1, res_s2);
216 result /= norm_frob2_lhs;
217 return result;
218 }
219
221 template< typename LHS_T, typename RHS_T >
222 auto
223 signed_perm_nork(const LHS_T& lhs, const RHS_T& rhs) -> const
224 RHS_T
225 {
226 // signed_perm_nork(A, kron(A, B)) is close to B
227 // Definition matches [v] Section 4, Theorem 4.1.
228 const auto lhs_s1 = lhs.size1();
229 const auto lhs_s2 = lhs.size2();
230 const auto rhs_s1 = rhs.size1();
231 const auto rhs_s2 = rhs.size2();
232 const auto res_s1 = rhs_s1 / lhs_s1;
233 const auto res_s2 = rhs_s2 / lhs_s2;
234 using Scalar_T = typename RHS_T::value_type;
235 const auto norm_frob2_lhs = Scalar_T( double(lhs_s1) );
236 auto result = RHS_T(res_s1, res_s2);
237 result.clear();
238 for (auto
239 lhs_it1 = lhs.begin1();
240 lhs_it1 != lhs.end1();
241 ++lhs_it1)
242 {
243 const auto lhs_it2 = lhs_it1.begin();
244 nork_range<LHS_T, RHS_T>(result, lhs_it2, rhs, res_s1, res_s2);
245 }
246 result /= norm_frob2_lhs;
247 return result;
248 }
249
251 template< typename Matrix_T >
252 auto
253 nnz(const Matrix_T& m) -> typename Matrix_T::size_type
254 {
255 using size_t = typename Matrix_T::size_type;
256 auto result = size_t(0);
257 for (auto
258 it1 = m.begin1();
259 it1 != m.end1();
260 ++it1)
261 for (auto& entry : it1)
262 if (entry != 0)
263 ++result;
264 return result;
265 }
266
268 template< typename Matrix_T >
269 auto
270 isinf(const Matrix_T& m) -> bool
271 {
272 using Scalar_T = typename Matrix_T::value_type;
273 for (auto
274 it1 = m.begin1();
275 it1 != m.end1();
276 ++it1)
277 for (auto& entry : it1)
279 return true;
280
281 return false;
282 }
283
285 template< typename Matrix_T >
286 auto
287 isnan(const Matrix_T& m) -> bool
288 {
289 using Scalar_T = typename Matrix_T::value_type;
290 for (auto
291 it1 = m.begin1();
292 it1 != m.end1();
293 ++it1)
294 for (auto& entry : it1)
296 return true;
297
298 return false;
299 }
300
302 template< typename Matrix_T >
303 inline
304 auto
305 unit(const typename Matrix_T::size_type dim) -> const
306 Matrix_T
307 {
308 using Scalar_T = typename Matrix_T::value_type;
309 return ublas::identity_matrix<Scalar_T>(dim);
310 }
311
313 template< typename LHS_T, typename RHS_T >
314 auto
315 mono_prod(const ublas::matrix_expression<LHS_T>& lhs,
316 const ublas::matrix_expression<RHS_T>& rhs) -> const typename RHS_T::expression_type
317 {
318 using rhs_expression_t = const RHS_T;
319 using matrix_row_t = typename ublas::matrix_row<rhs_expression_t>;
320
321 const auto dim = lhs().size1();
322 // The following assumes that RHS_T is a sparse matrix type.
323 auto result = RHS_T(dim, dim, dim);
324 for (auto
325 lhs_row = lhs().begin1();
326 lhs_row != lhs().end1();
327 ++lhs_row)
328 {
329 const auto& lhs_it = lhs_row.begin();
330 if (lhs_it != lhs_row.end())
331 {
332 const auto& rhs_row = matrix_row_t(rhs(), lhs_it.index2());
333 const auto& rhs_it = rhs_row.begin();
334 if (rhs_it != rhs_row.end())
335 result(lhs_it.index1(), rhs_it.index()) = (*lhs_it) * (*rhs_it);
336 }
337 }
338 return result;
339 }
340
342 template< typename LHS_T, typename RHS_T >
343 inline
344 auto
345 sparse_prod(const ublas::matrix_expression<LHS_T>& lhs,
346 const ublas::matrix_expression<RHS_T>& rhs) -> const typename RHS_T::expression_type
347 {
348 using expression_t = typename RHS_T::expression_type;
349 return ublas::sparse_prod<expression_t>(lhs(), rhs());
350 }
351
353 template< typename LHS_T, typename RHS_T >
354 inline
355 auto
356 prod(const ublas::matrix_expression<LHS_T>& lhs,
357 const ublas::matrix_expression<RHS_T>& rhs) -> const typename RHS_T::expression_type
358 {
359 const auto dim = lhs().size1();
360 RHS_T result(dim, dim);
361 ublas::axpy_prod(lhs, rhs, result, true);
362 return result;
363 }
364
366 template< typename Scalar_T, typename LHS_T, typename RHS_T >
367 auto
368 inner(const LHS_T& lhs, const RHS_T& rhs) -> Scalar_T
369 {
370 auto result = Scalar_T(0);
371 for (auto
372 lhs_it1 = lhs.begin1();
373 lhs_it1 != lhs.end1();
374 ++lhs_it1)
375 for (auto
376 lhs_it2 = lhs_it1.begin();
377 lhs_it2 != lhs_it1.end();
378 ++lhs_it2)
379 {
380 const auto& rhs_val = rhs(lhs_it2.index1(),lhs_it2.index2());
381 if (rhs_val != Scalar_T(0))
382 result += (*lhs_it2) * rhs_val;
383 }
384 return result / lhs.size1();
385 }
386
388 template< typename Matrix_T >
389 auto
390 norm_frob2(const Matrix_T& val) -> typename Matrix_T::value_type
391 {
392 using Scalar_T = typename Matrix_T::value_type;
393
394 auto result = Scalar_T(0);
395 for (auto
396 val_it1 = val.begin1();
397 val_it1 != val.end1();
398 ++val_it1)
399 for (auto& val_entry : val_it1)
400 {
401 if (numeric_traits<Scalar_T>::isNaN(val_entry))
403 result += val_entry * val_entry;
404 }
405 return result;
406 }
407
409 template< typename Matrix_T >
410 auto
411 trace(const Matrix_T& val) -> typename Matrix_T::value_type
412 {
413 using Scalar_T = typename Matrix_T::value_type;
414
415 auto result = Scalar_T(0);
416 auto dim = val.size1();
417 for (auto
418 ndx = decltype(dim)(0);
419 ndx != dim;
420 ++ndx)
421 {
422 const Scalar_T crd = val(ndx, ndx);
425 result += crd;
426 }
427 return result;
428 }
429
430#if defined(_GLUCAT_USE_BLAZE)
432 template< typename Matrix_T >
433 static
434 auto
435 to_blaze(const Matrix_T& val) -> blaze::DynamicMatrix<double,blaze::rowMajor>
436 {
437 const auto s1 = val.size1();
438 const auto s2 = val.size2();
439
440 using blaze_matrix_t = typename blaze::DynamicMatrix<double,blaze::rowMajor>;
441 auto result = blaze_matrix_t(s1, s2);
442
443 using Scalar_T = typename Matrix_T::value_type;
444 using traits_t = numeric_traits<Scalar_T>;
445
446 for (auto
447 val_it1 = val.begin1();
448 val_it1 != val.end1();
449 ++val_it1)
450 for (auto
451 val_it2 = val_it1.begin();
452 val_it2 != val_it1.end();
453 ++val_it2)
454 result(val_it2.index1(), val_it2.index2()) = traits_t::to_double(*val_it2);
455
456 return result;
457 }
458
459#endif
460
462 template< typename Matrix_T >
463 auto
464 eigenvalues(const Matrix_T& val) -> std::vector< std::complex<double> >
465 {
466 using complex_t = std::complex<double>;
467 using complex_vector_t = typename std::vector<complex_t>;
468
469 const auto dim = val.size1();
470 auto lambda = complex_vector_t(dim);
471
472#if defined(_GLUCAT_USE_BLAZE)
473 using complex_t = std::complex<double>;
474 using blaze_complex_vector_t = blaze::DynamicVector<complex_t, blaze::columnVector>;
475
476 auto blaze_val = to_blaze(val);
477 auto blaze_lambda = blaze_complex_vector_t(dim);
478 blaze::geev(blaze_val, blaze_lambda);
479
480 for (auto
481 k = decltype(dim)(0);
482 k != dim;
483 ++k)
484 lambda[k] = blaze_lambda[k];
485#endif
486 return lambda;
487 }
488
490 template< typename Matrix_T >
491 auto
493 {
494 using Scalar_T = typename Matrix_T::value_type;
495 eig_genus<Matrix_T> result;
496
497 auto lambda = eigenvalues(val);
498
499 std::set<double> arg_set;
500
501 const auto dim = lambda.size();
502 static const auto epsilon =
503 std::max(std::numeric_limits<double>::epsilon(),
504 numeric_traits<Scalar_T>::to_double(std::numeric_limits<Scalar_T>::epsilon()));
505 static const auto zero_eig_tol = 4096.0*epsilon;
506
507 bool neg_real_eig_found = false;
508 bool imag_eig_found = false;
509 bool zero_eig_found = false;
510
511 for (auto
512 k = decltype(dim)(0);
513 k != dim;
514 ++k)
515 {
516 const auto lambda_k = lambda[k];
517 arg_set.insert(std::arg(lambda_k));
518
519 const auto real_lambda_k = std::real(lambda_k);
520 const auto imag_lambda_k = std::imag(lambda_k);
521 const auto norm_tol = 4096.0*epsilon*std::norm(lambda_k);
522
523 if (!neg_real_eig_found &&
524 real_lambda_k < -epsilon &&
525 (imag_lambda_k == 0.0 ||
526 imag_lambda_k * imag_lambda_k < norm_tol))
527 neg_real_eig_found = true;
528 if (!imag_eig_found &&
529 imag_lambda_k > epsilon &&
530 (real_lambda_k == 0.0 ||
531 real_lambda_k * real_lambda_k < norm_tol))
532 imag_eig_found = true;
533 if (!zero_eig_found &&
534 std::norm(lambda_k) < zero_eig_tol)
535 zero_eig_found = true;
536 }
537
538 if (zero_eig_found)
539 result.m_is_singular = true;
540
541 static const auto pi = numeric_traits<double>::pi();
542 if (neg_real_eig_found)
543 {
544 if (imag_eig_found)
545 result.m_eig_case = both_eigs;
546 else
547 {
548 result.m_eig_case = neg_real_eigs;
549 result.m_safe_arg = Scalar_T(-pi / 2.0);
550 }
551 }
552
553 if (result.m_eig_case == both_eigs)
554 {
555 auto arg_it = arg_set.begin();
556 auto first_arg = *arg_it;
557 auto best_arg = first_arg;
558 auto best_diff = 0.0;
559 auto previous_arg = first_arg;
560 for (++arg_it;
561 arg_it != arg_set.end();
562 ++arg_it)
563 {
564 const auto arg_diff = *arg_it - previous_arg;
565 if (arg_diff > best_diff)
566 {
567 best_diff = arg_diff;
568 best_arg = previous_arg;
569 }
570 previous_arg = *arg_it;
571 }
572 const auto arg_diff = first_arg + 2.0*pi - previous_arg;
573 if (arg_diff > best_diff)
574 {
575 best_diff = arg_diff;
576 best_arg = previous_arg;
577 }
578 result.m_safe_arg = Scalar_T(pi - (best_arg + best_diff / 2.0));
579 }
580 return result;
581 }
582} }
583
584#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
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:177
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:368
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:223
auto trace(const Matrix_T &val) -> typename Matrix_T::value_type
Matrix trace.
Definition matrix_imp.h:411
auto isinf(const Matrix_T &m) -> bool
Infinite.
Definition matrix_imp.h:270
auto nnz(const Matrix_T &m) -> typename Matrix_T::size_type
Number of non-zeros.
Definition matrix_imp.h:253
static auto to_blaze(const Matrix_T &val) -> blaze::DynamicMatrix< double, blaze::rowMajor >
Convert matrix to Blaze format.
Definition matrix_imp.h:435
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:114
auto eigenvalues(const Matrix_T &val) -> std::vector< std::complex< double > >
Eigenvalues of a matrix.
Definition matrix_imp.h:464
auto norm_frob2(const Matrix_T &val) -> typename Matrix_T::value_type
Square of Frobenius norm.
Definition matrix_imp.h:390
auto classify_eigenvalues(const Matrix_T &val) -> eig_genus< Matrix_T >
Classify the eigenvalues of a matrix.
Definition matrix_imp.h:492
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:356
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:78
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:315
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:345
auto unit(const typename Matrix_T::size_type n) -> const Matrix_T
Unit matrix - as per Matlab eye.
Definition matrix_imp.h:305
auto isnan(const Matrix_T &m) -> bool
Not a Number.
Definition matrix_imp.h:287
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:147
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