glucat  0.8.2
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/matrix.h"
36 
37 # if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
38 # pragma GCC diagnostic push
39 # pragma GCC diagnostic ignored "-Wunused-local-typedefs"
40 # endif
41 
42 #include <boost/numeric/ublas/vector.hpp>
43 #include <boost/numeric/ublas/vector_proxy.hpp>
44 #include <boost/numeric/ublas/matrix.hpp>
45 #include <boost/numeric/ublas/matrix_expression.hpp>
46 #include <boost/numeric/ublas/matrix_proxy.hpp>
47 #include <boost/numeric/ublas/matrix_sparse.hpp>
48 #include <boost/numeric/ublas/operation.hpp>
49 #include <boost/numeric/ublas/operation_sparse.hpp>
50 
51 #if defined(_GLUCAT_USE_BINDINGS)
52 # include <boost/numeric/bindings/lapack/driver/gees.hpp>
53 # include <boost/numeric/bindings/ublas.hpp>
54 #endif
55 
56 # if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
57 # pragma GCC diagnostic pop
58 # endif
59 
60 #include <set>
61 
62 namespace glucat { namespace matrix
63 {
64  // References for algorithms:
65  // [v]: C. F. van Loan and N. Pitsianis, "Approximation with Kronecker products",
66  // in Linear Algebra for Large Scale and Real-Time Applications, Marc S. Moonen,
67  // Gene H. Golub, and Bart L. R. Moor (eds.), 1993, pp. 293--314.
68 
70  template< typename LHS_T, typename RHS_T >
71  const
72  RHS_T
73  kron(const LHS_T& lhs, const RHS_T& rhs)
74  {
75  typedef RHS_T matrix_t;
76  typedef typename matrix_t::size_type matrix_index_t;
77  const matrix_index_t rhs_s1 = rhs.size1();
78  const matrix_index_t rhs_s2 = rhs.size2();
79  matrix_t result(lhs.size1()*rhs_s1, lhs.size2()*rhs_s2);
80  result.clear();
81 
82  typedef typename matrix_t::value_type Scalar_T;
83  typedef typename LHS_T::const_iterator1 lhs_const_iterator1;
84  typedef typename LHS_T::const_iterator2 lhs_const_iterator2;
85  typedef typename RHS_T::const_iterator1 rhs_const_iterator1;
86  typedef typename RHS_T::const_iterator2 rhs_const_iterator2;
87  for (lhs_const_iterator1
88  lhs_it1 = lhs.begin1();
89  lhs_it1 != lhs.end1();
90  ++lhs_it1)
91  for (lhs_const_iterator2
92  lhs_it2 = lhs_it1.begin();
93  lhs_it2 != lhs_it1.end();
94  ++lhs_it2)
95  {
96  const matrix_index_t start1 = rhs_s1 * lhs_it2.index1();
97  const matrix_index_t start2 = rhs_s2 * lhs_it2.index2();
98  const Scalar_T& lhs_val = *lhs_it2;
99  for (rhs_const_iterator1
100  rhs_it1 = rhs.begin1();
101  rhs_it1 != rhs.end1();
102  ++rhs_it1)
103  for (rhs_const_iterator2
104  rhs_it2 = rhs_it1.begin();
105  rhs_it2 != rhs_it1.end();
106  ++rhs_it2)
107  result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
108  }
109  return result;
110  }
111 
113  template< typename LHS_T, typename RHS_T >
114  const
115  RHS_T
116  mono_kron(const LHS_T& lhs, const RHS_T& rhs)
117  {
118  typedef RHS_T matrix_t;
119  typedef typename matrix_t::size_type matrix_index_t;
120  const matrix_index_t rhs_s1 = rhs.size1();
121  const matrix_index_t rhs_s2 = rhs.size2();
122  const matrix_index_t dim = lhs.size1()*rhs_s1;
123  matrix_t result(dim, dim, dim);
124  result.clear();
125 
126  typedef typename matrix_t::value_type Scalar_T;
127  typedef typename LHS_T::const_iterator1 lhs_const_iterator1;
128  typedef typename LHS_T::const_iterator2 lhs_const_iterator2;
129  typedef typename RHS_T::const_iterator1 rhs_const_iterator1;
130  typedef typename RHS_T::const_iterator2 rhs_const_iterator2;
131  for (lhs_const_iterator1
132  lhs_it1 = lhs.begin1();
133  lhs_it1 != lhs.end1();
134  ++lhs_it1)
135  {
136  const lhs_const_iterator2 lhs_it2 = lhs_it1.begin();
137  const matrix_index_t start1 = rhs_s1 * lhs_it2.index1();
138  const matrix_index_t start2 = rhs_s2 * lhs_it2.index2();
139  const Scalar_T& lhs_val = *lhs_it2;
140  for (rhs_const_iterator1
141  rhs_it1 = rhs.begin1();
142  rhs_it1 != rhs.end1();
143  ++rhs_it1)
144  {
145  const rhs_const_iterator2 rhs_it2 = rhs_it1.begin();
146  result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
147  }
148  }
149  return result;
150  }
151 
153  template< typename LHS_T, typename RHS_T >
154  void
155  nork_range(RHS_T& result,
156  const typename LHS_T::const_iterator2 lhs_it2,
157  const RHS_T& rhs,
158  const typename RHS_T::size_type res_s1,
159  const typename RHS_T::size_type res_s2)
160  {
161  // Definition matches [v] Section 4, Theorem 4.1.
162  typedef RHS_T matrix_t;
163  typedef typename matrix_t::size_type matrix_index_t;
164  const matrix_index_t start1 = res_s1 * lhs_it2.index1();
165  const matrix_index_t start2 = res_s2 * lhs_it2.index2();
166  using ublas::range;
167  const range& range1 = range(start1, start1 + res_s1);
168  const range& range2 = range(start2, start2 + res_s2);
169  typedef ublas::matrix_range<const matrix_t> matrix_range_t;
170  const matrix_range_t& rhs_range = matrix_range_t(rhs, range1, range2);
171  typedef typename matrix_t::value_type Scalar_T;
172  const Scalar_T lhs_val = numeric_traits<Scalar_T>::to_scalar_t(*lhs_it2);
173  for (typename matrix_range_t::const_iterator1
174  rhs_it1 = rhs_range.begin1();
175  rhs_it1 != rhs_range.end1();
176  ++rhs_it1)
177  for (typename matrix_range_t::const_iterator2
178  rhs_it2 = rhs_it1.begin();
179  rhs_it2 != rhs_it1.end();
180  ++rhs_it2)
181  result(rhs_it2.index1(), rhs_it2.index2()) += lhs_val * *rhs_it2;
182  }
183 
185  template< typename LHS_T, typename RHS_T >
186  const
187  RHS_T
188  nork(const LHS_T& lhs, const RHS_T& rhs, const bool mono)
189  {
190  // nork(A, kron(A, B)) is close to B
191  // Definition matches [v] Section 4, Theorem 4.1.
192  typedef RHS_T matrix_t;
193  typedef typename LHS_T::size_type lhs_index_t;
194  typedef typename matrix_t::size_type matrix_index_t;
195  const lhs_index_t lhs_s1 = lhs.size1();
196  const lhs_index_t lhs_s2 = lhs.size2();
197  const matrix_index_t rhs_s1 = rhs.size1();
198  const matrix_index_t rhs_s2 = rhs.size2();
199  const matrix_index_t res_s1 = rhs_s1 / lhs_s1;
200  const matrix_index_t res_s2 = rhs_s2 / lhs_s2;
201  typedef typename matrix_t::value_type Scalar_T;
202  const Scalar_T norm_frob2_lhs = norm_frob2(lhs);
203  if (!mono)
204  {
205  typedef error<matrix_t> error_t;
206  if (rhs_s1 == 0)
207  throw error_t("matrix", "nork: number of rows must not be 0");
208  if (rhs_s2 == 0)
209  throw error_t("matrix", "nork: number of cols must not be 0");
210  if (res_s1 * lhs_s1 != rhs_s1)
211  throw error_t("matrix", "nork: incompatible numbers of rows");
212  if (res_s2 * lhs_s2 != rhs_s2)
213  throw error_t("matrix", "nork: incompatible numbers of cols");
214  if (norm_frob2_lhs == Scalar_T(0))
215  throw error_t("matrix", "nork: LHS must not be 0");
216  }
217  matrix_t result(res_s1, res_s2);
218  result.clear();
219  for (typename LHS_T::const_iterator1
220  lhs_it1 = lhs.begin1();
221  lhs_it1 != lhs.end1();
222  ++lhs_it1)
223  for (typename LHS_T::const_iterator2
224  lhs_it2 = lhs_it1.begin();
225  lhs_it2 != lhs_it1.end();
226  ++lhs_it2)
227  if (*lhs_it2 != Scalar_T(0))
228  nork_range<LHS_T, RHS_T>(result, lhs_it2, rhs, res_s1, res_s2);
229  result /= norm_frob2_lhs;
230  return result;
231  }
232 
234  template< typename LHS_T, typename RHS_T >
235  const
236  RHS_T
237  signed_perm_nork(const LHS_T& lhs, const RHS_T& rhs)
238  {
239  // signed_perm_nork(A, kron(A, B)) is close to B
240  // Definition matches [v] Section 4, Theorem 4.1.
241  typedef RHS_T matrix_t;
242  typedef typename LHS_T::size_type lhs_index_t;
243  typedef typename matrix_t::size_type matrix_index_t;
244  const lhs_index_t lhs_s1 = lhs.size1();
245  const lhs_index_t lhs_s2 = lhs.size2();
246  const matrix_index_t rhs_s1 = rhs.size1();
247  const matrix_index_t rhs_s2 = rhs.size2();
248  const matrix_index_t res_s1 = rhs_s1 / lhs_s1;
249  const matrix_index_t res_s2 = rhs_s2 / lhs_s2;
250  typedef typename matrix_t::value_type Scalar_T;
251  const Scalar_T norm_frob2_lhs = Scalar_T( double(lhs_s1) );
252  matrix_t result(res_s1, res_s2);
253  result.clear();
254  for (typename LHS_T::const_iterator1
255  lhs_it1 = lhs.begin1();
256  lhs_it1 != lhs.end1();
257  ++lhs_it1)
258  {
259  const typename LHS_T::const_iterator2 lhs_it2 = lhs_it1.begin();
260  nork_range<LHS_T, RHS_T>(result, lhs_it2, rhs, res_s1, res_s2);
261  }
262  result /= norm_frob2_lhs;
263  return result;
264  }
265 
267  template< typename Matrix_T >
268  typename Matrix_T::size_type
269  nnz(const Matrix_T& m)
270  {
271  typedef Matrix_T matrix_t;
272  typedef typename matrix_t::size_type matrix_index_t;
273  typedef typename matrix_t::const_iterator1 const_iterator1;
274  typedef typename matrix_t::const_iterator2 const_iterator2;
275  matrix_index_t result = 0;
276  for (const_iterator1
277  it1 = m.begin1();
278  it1 != m.end1();
279  ++it1)
280  for (const_iterator2
281  it2 = it1.begin();
282  it2 != it1.end();
283  ++it2)
284  if (*it2 != 0)
285  ++result;
286  return result;
287  }
288 
290  template< typename Matrix_T >
291  bool
292  isnan(const Matrix_T& m)
293  {
294  typedef Matrix_T matrix_t;
295  typedef typename matrix_t::value_type Scalar_T;
296  typedef typename matrix_t::const_iterator1 const_iterator1;
297  typedef typename matrix_t::const_iterator2 const_iterator2;
298  for (const_iterator1
299  it1 = m.begin1();
300  it1 != m.end1();
301  ++it1)
302  for (const_iterator2
303  it2 = it1.begin();
304  it2 != it1.end();
305  ++it2)
307  return true;
308 
309  return false;
310  }
311 
313  template< typename Matrix_T >
314  inline
315  const
316  Matrix_T
317  unit(const typename Matrix_T::size_type dim)
318  {
319  typedef typename Matrix_T::value_type Scalar_T;
320  return ublas::identity_matrix<Scalar_T>(dim);
321  }
322 
324  template< typename LHS_T, typename RHS_T >
325  const typename RHS_T::expression_type
326  mono_prod(const ublas::matrix_expression<LHS_T>& lhs,
327  const ublas::matrix_expression<RHS_T>& rhs)
328  {
329  typedef const LHS_T lhs_expression_t;
330  typedef const RHS_T rhs_expression_t;
331  typedef typename RHS_T::expression_type matrix_t;
332  typedef typename matrix_t::size_type matrix_index_t;
333  typedef typename lhs_expression_t::const_iterator1 lhs_const_iterator1;
334  typedef typename lhs_expression_t::const_iterator2 lhs_const_iterator2;
335  typedef typename ublas::matrix_row<rhs_expression_t> matrix_row_t;
336  typedef typename matrix_row_t::const_iterator row_const_iterator;
337 
338  const matrix_index_t dim = lhs().size1();
339  // The following assumes that matrix_t is a sparse matrix type.
340  matrix_t result = matrix_t(dim, dim, dim);
341  for (lhs_const_iterator1
342  lhs_row = lhs().begin1();
343  lhs_row != lhs().end1();
344  ++lhs_row)
345  {
346  const lhs_const_iterator2& lhs_it = lhs_row.begin();
347  if (lhs_it != lhs_row.end())
348  {
349  const matrix_row_t& rhs_row = matrix_row_t(rhs(), lhs_it.index2());
350  const row_const_iterator& rhs_it = rhs_row.begin();
351  if (rhs_it != rhs_row.end())
352  result(lhs_it.index1(), rhs_it.index()) = (*lhs_it) * (*rhs_it);
353  }
354  }
355  return result;
356  }
357 
359  template< typename LHS_T, typename RHS_T >
360  inline
361  const typename RHS_T::expression_type
362  sparse_prod(const ublas::matrix_expression<LHS_T>& lhs,
363  const ublas::matrix_expression<RHS_T>& rhs)
364  {
365  typedef typename RHS_T::expression_type expression_t;
366  return ublas::sparse_prod<expression_t>(lhs(), rhs());
367  }
368 
370  template< typename LHS_T, typename RHS_T >
371  inline
372  const typename RHS_T::expression_type
373  prod(const ublas::matrix_expression<LHS_T>& lhs,
374  const ublas::matrix_expression<RHS_T>& rhs)
375  {
376 #if defined(_GLUCAT_USE_DENSE_MATRICES)
377  typedef typename RHS_T::size_type matrix_index_t;
378  const matrix_index_t dim = lhs().size1();
379  RHS_T result(dim, dim);
380  ublas::axpy_prod(lhs, rhs, result, true);
381  return result;
382 #else
383  typedef typename RHS_T::expression_type expression_t;
384  return ublas::sparse_prod<expression_t>(lhs(), rhs());
385 #endif
386  }
387 
389  template< typename Scalar_T, typename LHS_T, typename RHS_T >
390  Scalar_T
391  inner(const LHS_T& lhs, const RHS_T& rhs)
392  {
393  Scalar_T result = Scalar_T(0);
394  for (typename LHS_T::const_iterator1
395  lhs_it1 = lhs.begin1();
396  lhs_it1 != lhs.end1();
397  ++lhs_it1)
398  for (typename LHS_T::const_iterator2
399  lhs_it2 = lhs_it1.begin();
400  lhs_it2 != lhs_it1.end();
401  ++lhs_it2)
402  {
403  const Scalar_T& rhs_val = rhs(lhs_it2.index1(),lhs_it2.index2());
404  if (rhs_val != Scalar_T(0))
405  result += (*lhs_it2) * rhs_val;
406  }
407  return result / lhs.size1();
408  }
409 
411  template< typename Matrix_T >
412  typename Matrix_T::value_type
413  norm_frob2(const Matrix_T& val)
414  {
415  typedef typename Matrix_T::value_type Scalar_T;
416 
417  Scalar_T result = Scalar_T(0);
418  for (typename Matrix_T::const_iterator1
419  val_it1 = val.begin1();
420  val_it1 != val.end1();
421  ++val_it1)
422  for (typename Matrix_T::const_iterator2
423  val_it2 = val_it1.begin();
424  val_it2 != val_it1.end();
425  ++val_it2)
426  {
427  if (numeric_traits<Scalar_T>::isNaN(*val_it2))
429  result += (*val_it2) * (*val_it2);
430  }
431  return result;
432  }
433 
435  template< typename Matrix_T >
436  typename Matrix_T::value_type
437  trace(const Matrix_T& val)
438  {
439  typedef typename Matrix_T::value_type Scalar_T;
440  typedef typename Matrix_T::size_type matrix_index_t;
441 
442  Scalar_T result = Scalar_T(0);
443  matrix_index_t dim = val.size1();
444  for (matrix_index_t ndx=0;
445  ndx != dim;
446  ++ndx)
447  {
448  const Scalar_T crd = val(ndx, ndx);
451  result += crd;
452  }
453  return result;
454  }
455 
456 #if defined(_GLUCAT_USE_BINDINGS)
457  template< typename Matrix_T >
459  static
460  ublas::matrix<double, ublas::column_major>
461  to_lapack(const Matrix_T& val)
462  {
463  typedef typename Matrix_T::size_type matrix_index_t;
464  const matrix_index_t s1 = val.size1();
465  const matrix_index_t s2 = val.size2();
466 
467  typedef typename ublas::matrix<double, ublas::column_major> lapack_matrix_t;
468  lapack_matrix_t result(s1, s2);
469  result.clear();
470 
471  typedef typename Matrix_T::value_type Scalar_T;
472  typedef numeric_traits<Scalar_T> traits_t;
473 
474  typedef typename Matrix_T::const_iterator1 const_iterator1;
475  typedef typename Matrix_T::const_iterator2 const_iterator2;
476  for (const_iterator1
477  val_it1 = val.begin1();
478  val_it1 != val.end1();
479  ++val_it1)
480  for (const_iterator2
481  val_it2 = val_it1.begin();
482  val_it2 != val_it1.end();
483  ++val_it2)
484  result(val_it2.index1(), val_it2.index2()) = traits_t::to_double(*val_it2);
485 
486  return result;
487  }
488 #endif
489 
491  template< typename Matrix_T >
492  ublas::vector< std::complex<double> >
493  eigenvalues(const Matrix_T& val)
494  {
495  typedef std::complex<double> complex_t;
496  typedef typename ublas::vector<complex_t> complex_vector_t;
497  typedef typename Matrix_T::size_type matrix_index_t;
498 
499  const matrix_index_t dim = val.size1();
500  complex_vector_t lambda = complex_vector_t(dim);
501  lambda.clear();
502 
503 #if defined(_GLUCAT_USE_BINDINGS)
504  namespace lapack = boost::numeric::bindings::lapack;
505  typedef typename ublas::matrix<double, ublas::column_major> lapack_matrix_t;
506 
507  lapack_matrix_t T = to_lapack(val);
508  lapack_matrix_t V = T;
509  typedef typename ublas::vector<double> vector_t;
510  vector_t real_lambda = vector_t(dim);
511  vector_t imag_lambda = vector_t(dim);
512  fortran_int_t sdim = 0;
513 
514  lapack::gees ('N', 'N', (external_fp)0, T, sdim, real_lambda, imag_lambda, V );
515 
516  lambda.clear();
517  for (vector_t::size_type k=0; k!= dim; ++k)
518  lambda[k] = complex_t(real_lambda[k], imag_lambda[k]);
519 #endif
520  return lambda;
521  }
522 
524  template< typename Matrix_T >
525  eig_genus<Matrix_T>
526  classify_eigenvalues(const Matrix_T& val)
527  {
528  typedef typename Matrix_T::value_type Scalar_T;
529  eig_genus<Matrix_T> result;
530  result.m_eig_case = safe_eig_case;
531  result.m_safe_arg = Scalar_T(0);
532 
533  typedef std::complex<double> complex_t;
534  typedef typename ublas::vector<complex_t> complex_vector_t;
535  complex_vector_t lambda = eigenvalues(val);
536 
537  std::set<double> arg_set;
538 
539  typedef typename complex_vector_t::size_type vector_index_t;
540  const vector_index_t dim = lambda.size();
541  static const double epsilon =
542  std::max(std::numeric_limits<double>::epsilon(),
543  numeric_traits<Scalar_T>::to_double(std::numeric_limits<Scalar_T>::epsilon()));
544 
545  bool negative_eig_found = false;
546  bool imaginary_eig_found = false;
547  for (vector_index_t k=0; k != dim; ++k)
548  {
549  const complex_t lambda_k = lambda[k];
550  arg_set.insert(std::arg(lambda_k));
551 
552  const double real_lambda_k = std::real(lambda_k);
553  const double imag_lambda_k = std::imag(lambda_k);
554  const double norm_tol = 4096.0*epsilon*std::norm(lambda_k);
555 
556  if (!negative_eig_found &&
557  real_lambda_k < -epsilon &&
558  (imag_lambda_k == 0.0 ||
559  imag_lambda_k * imag_lambda_k < norm_tol))
560  negative_eig_found = true;
561  if (!imaginary_eig_found &&
562  std::abs(imag_lambda_k) > epsilon &&
563  (real_lambda_k == 0.0 ||
564  real_lambda_k * real_lambda_k < norm_tol))
565  imaginary_eig_found = true;
566  }
567 
568  static const double pi = numeric_traits<double>::pi();
569  if (negative_eig_found)
570  {
571  if (imaginary_eig_found)
572  result.m_eig_case = both_eig_case;
573  else
574  {
575  result.m_eig_case = negative_eig_case;
576  result.m_safe_arg = -pi/2.0;
577  }
578  }
579 
580  if (result.m_eig_case == both_eig_case)
581  {
582  typename std::set<double>::const_iterator arg_it = arg_set.begin();
583  double first_arg = *arg_it;
584  double best_arg = first_arg;
585  double best_diff = 0.0;
586  double previous_arg = first_arg;
587  for (++arg_it;
588  arg_it != arg_set.end();
589  ++arg_it)
590  {
591  const double arg_diff = *arg_it - previous_arg;
592  if (arg_diff > best_diff)
593  {
594  best_diff = arg_diff;
595  best_arg = previous_arg;
596  }
597  previous_arg = *arg_it;
598  }
599  const double arg_diff = first_arg + 2.0*pi - previous_arg;
600  if (arg_diff > best_diff)
601  {
602  best_diff = arg_diff;
603  best_arg = previous_arg;
604  }
605  result.m_safe_arg = pi - (best_arg + best_diff / 2.0);
606  }
607  return result;
608  }
609 } }
610 
611 #endif // _GLUCAT_MATRIX_IMP_H
const RHS_T nork(const LHS_T &lhs, const RHS_T &rhs, const bool mono=true)
Left inverse of Kronecker product.
Definition: matrix_imp.h:188
const RHS_T::expression_type prod(const ublas::matrix_expression< LHS_T > &lhs, const ublas::matrix_expression< RHS_T > &rhs)
Product of matrices.
Definition: matrix_imp.h:373
static Scalar_T to_scalar_t(const Other_Scalar_T &val)
Cast to Scalar_T.
Definition: scalar.h:140
Scalar_T abs(const Multivector< Scalar_T, LO, HI > &val)
Absolute value == sqrt(norm)
const RHS_T kron(const LHS_T &lhs, const RHS_T &rhs)
Kronecker tensor product of matrices - as per Matlab kron.
Definition: matrix_imp.h:73
float pi
Definition: PyClical.pyx:1857
const RHS_T mono_kron(const LHS_T &lhs, const RHS_T &rhs)
Sparse Kronecker tensor product of monomial matrices.
Definition: matrix_imp.h:116
Scalar_T norm(const Multivector< Scalar_T, LO, HI > &val)
Scalar_T norm == sum of norm of coordinates.
static Scalar_T pi()
Pi.
Definition: scalar.h:188
Extra traits which extend numeric limits.
Definition: scalar.h:46
Matrix_T::value_type norm_frob2(const Matrix_T &val)
Square of Frobenius norm.
Definition: matrix_imp.h:413
bool isnan(const Matrix_T &m)
Not a Number.
Definition: matrix_imp.h:292
static Scalar_T NaN()
Smart NaN.
Definition: scalar.h:114
ublas::vector< std::complex< double > > eigenvalues(const Matrix_T &val)
Eigenvalues of a matrix.
Definition: matrix_imp.h:493
Scalar_T imag(const Multivector< Scalar_T, LO, HI > &val)
Imaginary part: deprecated (always 0)
Structure containing classification of eigenvalues.
Definition: matrix.h:131
const Matrix_T unit(const typename Matrix_T::size_type n)
Unit matrix - as per Matlab eye.
Definition: matrix_imp.h:317
eig_case_t m_eig_case
What kind of eigenvalues does the matrix contain?
Definition: matrix.h:135
static ublas::matrix< double, ublas::column_major > to_lapack(const Matrix_T &val)
Convert matrix to LAPACK format.
Definition: matrix_imp.h:461
Specific exception class.
Definition: errors.h:57
Scalar_T real(const Multivector< Scalar_T, LO, HI > &val)
Real part: synonym for scalar part.
const RHS_T::expression_type mono_prod(const ublas::matrix_expression< LHS_T > &lhs, const ublas::matrix_expression< RHS_T > &rhs)
Product of monomial matrices.
Definition: matrix_imp.h:326
const RHS_T::expression_type sparse_prod(const ublas::matrix_expression< LHS_T > &lhs, const ublas::matrix_expression< RHS_T > &rhs)
Product of sparse matrices.
Definition: matrix_imp.h:362
Matrix_T::value_type trace(const Matrix_T &val)
Matrix trace.
Definition: matrix_imp.h:437
Scalar_T m_safe_arg
Argument such that exp(pi-m_safe_arg) lies between arguments of eigenvalues.
Definition: matrix.h:137
const RHS_T signed_perm_nork(const LHS_T &lhs, const RHS_T &rhs)
Left inverse of Kronecker product where lhs is a signed permutation matrix.
Definition: matrix_imp.h:237
eig_genus< Matrix_T > classify_eigenvalues(const Matrix_T &val)
Classify the eigenvalues of a matrix.
Definition: matrix_imp.h:526
Scalar_T inner(const LHS_T &lhs, const RHS_T &rhs)
Inner product: sum(x(i,j)*y(i,j))/x.nrows()
Definition: matrix_imp.h:391
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:155
Matrix_T::size_type nnz(const Matrix_T &m)
Number of non-zeros.
Definition: matrix_imp.h:269