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