1#ifndef _GLUCAT_MATRIX_IMP_H
2#define _GLUCAT_MATRIX_IMP_H
39# if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
40# pragma GCC diagnostic push
41# pragma GCC diagnostic ignored "-Wunused-local-typedefs"
43# if defined(_GLUCAT_HAVE_BOOST_SERIALIZATION_ARRAY_WRAPPER_H)
44# include <boost/serialization/array_wrapper.hpp>
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>
55#if defined(_GLUCAT_USE_BLAZE)
56#include <blaze/Math.h>
57#include <blaze/math/DynamicMatrix.h>
58#include <blaze/math/DynamicVector.h>
61# if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
62# pragma GCC diagnostic pop
76 template<
typename LHS_T,
typename RHS_T >
78 kron(
const LHS_T& lhs,
const RHS_T& rhs) ->
const
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);
87 lhs_it1 = lhs.begin1();
88 lhs_it1 != lhs.end1();
91 lhs_it2 = lhs_it1.begin();
92 lhs_it2 != lhs_it1.end();
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;
99 rhs_it1 = rhs.begin1();
100 rhs_it1 != rhs.end1();
103 rhs_it2 = rhs_it1.begin();
104 rhs_it2 != rhs_it1.end();
106 result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
112 template<
typename LHS_T,
typename RHS_T >
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);
124 lhs_it1 = lhs.begin1();
125 lhs_it1 != lhs.end1();
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;
133 rhs_it1 = rhs.begin1();
134 rhs_it1 != rhs.end1();
137 const auto rhs_it2 = rhs_it1.begin();
138 result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
145 template<
typename LHS_T,
typename RHS_T >
148 const typename LHS_T::const_iterator2 lhs_it2,
150 const typename RHS_T::size_type res_s1,
151 const typename RHS_T::size_type res_s2)
154 const auto start1 = res_s1 * lhs_it2.index1();
155 const auto start2 = res_s2 * lhs_it2.index2();
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;
164 rhs_it1 = rhs_range.begin1();
165 rhs_it1 != rhs_range.end1();
168 rhs_it2 = rhs_it1.begin();
169 rhs_it2 != rhs_it1.end();
171 result(rhs_it2.index1(), rhs_it2.index2()) += lhs_val * *rhs_it2;
175 template<
typename LHS_T,
typename RHS_T >
177 nork(
const LHS_T& lhs,
const RHS_T& rhs,
const bool mono) ->
const
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;
194 throw error_t(
"matrix",
"nork: number of rows must not be 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");
204 auto result = RHS_T(res_s1, res_s2);
207 lhs_it1 = lhs.begin1();
208 lhs_it1 != lhs.end1();
211 lhs_it2 = lhs_it1.begin();
212 lhs_it2 != lhs_it1.end();
214 if (*lhs_it2 != Scalar_T(0))
216 result /= norm_frob2_lhs;
221 template<
typename LHS_T,
typename RHS_T >
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);
239 lhs_it1 = lhs.begin1();
240 lhs_it1 != lhs.end1();
243 const auto lhs_it2 = lhs_it1.begin();
246 result /= norm_frob2_lhs;
251 template<
typename Matrix_T >
253 nnz(
const Matrix_T& m) ->
typename Matrix_T::size_type
255 using size_t =
typename Matrix_T::size_type;
256 auto result = size_t(0);
261 for (
auto& entry : it1)
268 template<
typename Matrix_T >
272 using Scalar_T =
typename Matrix_T::value_type;
277 for (
auto& entry : it1)
285 template<
typename Matrix_T >
289 using Scalar_T =
typename Matrix_T::value_type;
294 for (
auto& entry : it1)
302 template<
typename Matrix_T >
305 unit(
const typename Matrix_T::size_type dim) ->
const
308 using Scalar_T =
typename Matrix_T::value_type;
309 return ublas::identity_matrix<Scalar_T>(dim);
313 template<
typename LHS_T,
typename RHS_T >
316 const ublas::matrix_expression<RHS_T>& rhs) ->
const typename RHS_T::expression_type
318 using rhs_expression_t =
const RHS_T;
319 using matrix_row_t =
typename ublas::matrix_row<rhs_expression_t>;
321 const auto dim = lhs().size1();
323 auto result = RHS_T(dim, dim, dim);
325 lhs_row = lhs().begin1();
326 lhs_row != lhs().end1();
329 const auto& lhs_it = lhs_row.begin();
330 if (lhs_it != lhs_row.end())
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);
342 template<
typename LHS_T,
typename RHS_T >
346 const ublas::matrix_expression<RHS_T>& rhs) ->
const typename RHS_T::expression_type
348 using expression_t =
typename RHS_T::expression_type;
349 return ublas::sparse_prod<expression_t>(lhs(), rhs());
353 template<
typename LHS_T,
typename RHS_T >
356 prod(
const ublas::matrix_expression<LHS_T>& lhs,
357 const ublas::matrix_expression<RHS_T>& rhs) ->
const typename RHS_T::expression_type
359 const auto dim = lhs().size1();
360 RHS_T result(dim, dim);
361 ublas::axpy_prod(lhs, rhs, result,
true);
366 template<
typename Scalar_T,
typename LHS_T,
typename RHS_T >
368 inner(
const LHS_T& lhs,
const RHS_T& rhs) -> Scalar_T
370 auto result = Scalar_T(0);
372 lhs_it1 = lhs.begin1();
373 lhs_it1 != lhs.end1();
376 lhs_it2 = lhs_it1.begin();
377 lhs_it2 != lhs_it1.end();
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;
384 return result / lhs.size1();
388 template<
typename Matrix_T >
390 norm_frob2(
const Matrix_T& val) ->
typename Matrix_T::value_type
392 using Scalar_T =
typename Matrix_T::value_type;
394 auto result = Scalar_T(0);
396 val_it1 = val.begin1();
397 val_it1 != val.end1();
399 for (
auto& val_entry : val_it1)
403 result += val_entry * val_entry;
409 template<
typename Matrix_T >
411 trace(
const Matrix_T& val) ->
typename Matrix_T::value_type
413 using Scalar_T =
typename Matrix_T::value_type;
415 auto result = Scalar_T(0);
416 auto dim = val.size1();
418 ndx =
decltype(dim)(0);
422 const Scalar_T crd = val(ndx, ndx);
430#if defined(_GLUCAT_USE_BLAZE)
432 template<
typename Matrix_T >
435 to_blaze(
const Matrix_T& val) -> blaze::DynamicMatrix<double,blaze::rowMajor>
437 const auto s1 = val.size1();
438 const auto s2 = val.size2();
440 using blaze_matrix_t =
typename blaze::DynamicMatrix<double,blaze::rowMajor>;
441 auto result = blaze_matrix_t(s1, s2);
443 using Scalar_T =
typename Matrix_T::value_type;
447 val_it1 = val.begin1();
448 val_it1 != val.end1();
451 val_it2 = val_it1.begin();
452 val_it2 != val_it1.end();
454 result(val_it2.index1(), val_it2.index2()) = traits_t::to_double(*val_it2);
462 template<
typename Matrix_T >
464 eigenvalues(
const Matrix_T& val) -> std::vector< std::complex<double> >
466 using complex_t = std::complex<double>;
467 using complex_vector_t =
typename std::vector<complex_t>;
469 const auto dim = val.size1();
470 auto lambda = complex_vector_t(dim);
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>;
477 auto blaze_lambda = blaze_complex_vector_t(dim);
478 blaze::geev(blaze_val, blaze_lambda);
481 k =
decltype(dim)(0);
484 lambda[k] = blaze_lambda[k];
490 template<
typename Matrix_T >
494 using Scalar_T =
typename Matrix_T::value_type;
499 std::set<double> arg_set;
501 const auto dim = lambda.size();
503 std::max(std::numeric_limits<double>::epsilon(),
505 static const auto zero_eig_tol = 4096.0*
epsilon;
507 bool neg_real_eig_found =
false;
508 bool imag_eig_found =
false;
509 bool zero_eig_found =
false;
512 k =
decltype(dim)(0);
516 const auto lambda_k = lambda[k];
517 arg_set.insert(std::arg(lambda_k));
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);
523 if (!neg_real_eig_found &&
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 &&
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;
542 if (neg_real_eig_found)
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;
561 arg_it != arg_set.end();
564 const auto arg_diff = *arg_it - previous_arg;
565 if (arg_diff > best_diff)
567 best_diff = arg_diff;
568 best_arg = previous_arg;
570 previous_arg = *arg_it;
572 const auto arg_diff = first_arg + 2.0*pi - previous_arg;
573 if (arg_diff > best_diff)
575 best_diff = arg_diff;
576 best_arg = previous_arg;
578 result.
m_safe_arg = Scalar_T(pi - (best_arg + best_diff / 2.0));
Specific exception class.
Extra traits which extend numeric limits.
static auto isNaN(const Scalar_T &val, bool_to_type< false >) -> bool
Smart isnan specialised for Scalar_T without quiet NaN.
static auto NaN() -> Scalar_T
Smart NaN.
static auto to_scalar_t(const Other_Scalar_T &val) -> Scalar_T
Cast to Scalar_T.
static auto to_double(const Scalar_T &val) -> double
Cast to double.
static auto isInf(const Scalar_T &val, bool_to_type< false >) -> bool
Smart isinf specialised for Scalar_T without infinity.
static auto pi() -> Scalar_T
Pi.
auto nork(const LHS_T &lhs, const RHS_T &rhs, const bool mono=true) -> const RHS_T
Left inverse of Kronecker product.
auto inner(const LHS_T &lhs, const RHS_T &rhs) -> Scalar_T
Inner product: sum(x(i,j)*y(i,j))/x.nrows()
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.
auto trace(const Matrix_T &val) -> typename Matrix_T::value_type
Matrix trace.
auto isinf(const Matrix_T &m) -> bool
Infinite.
auto nnz(const Matrix_T &m) -> typename Matrix_T::size_type
Number of non-zeros.
static auto to_blaze(const Matrix_T &val) -> blaze::DynamicMatrix< double, blaze::rowMajor >
Convert matrix to Blaze format.
auto mono_kron(const LHS_T &lhs, const RHS_T &rhs) -> const RHS_T
Sparse Kronecker tensor product of monomial matrices.
auto eigenvalues(const Matrix_T &val) -> std::vector< std::complex< double > >
Eigenvalues of a matrix.
auto norm_frob2(const Matrix_T &val) -> typename Matrix_T::value_type
Square of Frobenius norm.
auto classify_eigenvalues(const Matrix_T &val) -> eig_genus< Matrix_T >
Classify the eigenvalues of a matrix.
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.
auto kron(const LHS_T &lhs, const RHS_T &rhs) -> const RHS_T
Kronecker tensor product of matrices - as per Matlab kron.
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.
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.
auto unit(const typename Matrix_T::size_type n) -> const Matrix_T
Unit matrix - as per Matlab eye.
auto isnan(const Matrix_T &m) -> bool
Not a Number.
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.
Structure containing classification of eigenvalues.
bool m_is_singular
Is the matrix singular?
Scalar_T m_safe_arg
Argument such that exp(pi-m_safe_arg) lies between arguments of eigenvalues.
eig_case_t m_eig_case
What kind of eigenvalues does the matrix contain?