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 
73 namespace 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)
436  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)
467  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 =
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
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 trace(const Matrix_T &val) -> typename Matrix_T::value_type
Matrix trace.
Definition: matrix_imp.h:416
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
static auto NaN() -> Scalar_T
Smart NaN.
Definition: scalar.h:115
float pi
Definition: PyClical.pyx:1907
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 norm(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Scalar_T norm == sum of norm of coordinates.
bool m_is_singular
Is the matrix singular?
Definition: matrix.h:144
Extra traits which extend numeric limits.
Definition: scalar.h:47
auto classify_eigenvalues(const Matrix_T &val) -> eig_genus< Matrix_T >
Classify the eigenvalues of a matrix.
Definition: matrix_imp.h:548
auto isnan(const Matrix_T &m) -> bool
Not a Number.
Definition: matrix_imp.h:292
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
const scalar_t epsilon
Definition: PyClical.h:150
Structure containing classification of eigenvalues.
Definition: matrix.h:140
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 eigenvalues(const Matrix_T &val) -> std::vector< std::complex< double > >
Eigenvalues of a matrix.
Definition: matrix_imp.h:500
eig_case_t m_eig_case
What kind of eigenvalues does the matrix contain?
Definition: matrix.h:146
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
Specific exception class.
Definition: errors.h:56
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
static auto pi() -> Scalar_T
Pi.
Definition: scalar.h:189
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 nnz(const Matrix_T &m) -> typename Matrix_T::size_type
Number of non-zeros.
Definition: matrix_imp.h:258
auto norm_frob2(const Matrix_T &val) -> typename Matrix_T::value_type
Square of Frobenius norm.
Definition: matrix_imp.h:395
static auto to_scalar_t(const Other_Scalar_T &val) -> Scalar_T
Cast to Scalar_T.
Definition: scalar.h:141
auto real(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Real part: synonym for scalar part.
Scalar_T m_safe_arg
Argument such that exp(pi-m_safe_arg) lies between arguments of eigenvalues.
Definition: matrix.h:148
auto imag(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Imaginary part: deprecated (always 0)
auto unit(const typename Matrix_T::size_type n) -> const Matrix_T
Unit matrix - as per Matlab eye.
Definition: matrix_imp.h:310
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
auto isinf(const Matrix_T &m) -> bool
Infinite.
Definition: matrix_imp.h:275