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_BINDINGS)
56# include <boost/numeric/bindings/lapack/driver/gees.hpp>
57# include <boost/numeric/bindings/ublas.hpp>
60#if defined(_GLUCAT_USE_BLAZE)
61#include <blaze/Math.h>
62#include <blaze/math/DynamicMatrix.h>
63#include <blaze/math/DynamicVector.h>
66# if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
67# pragma GCC diagnostic pop
73namespace glucat {
namespace matrix
81 template<
typename LHS_T,
typename RHS_T >
83 kron(
const LHS_T& lhs,
const RHS_T& rhs) ->
const
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);
92 lhs_it1 = lhs.begin1();
93 lhs_it1 != lhs.end1();
96 lhs_it2 = lhs_it1.begin();
97 lhs_it2 != lhs_it1.end();
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;
104 rhs_it1 = rhs.begin1();
105 rhs_it1 != rhs.end1();
108 rhs_it2 = rhs_it1.begin();
109 rhs_it2 != rhs_it1.end();
111 result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
117 template<
typename LHS_T,
typename RHS_T >
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);
129 lhs_it1 = lhs.begin1();
130 lhs_it1 != lhs.end1();
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;
138 rhs_it1 = rhs.begin1();
139 rhs_it1 != rhs.end1();
142 const auto rhs_it2 = rhs_it1.begin();
143 result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
150 template<
typename LHS_T,
typename RHS_T >
153 const typename LHS_T::const_iterator2 lhs_it2,
155 const typename RHS_T::size_type res_s1,
156 const typename RHS_T::size_type res_s2)
159 const auto start1 = res_s1 * lhs_it2.index1();
160 const auto start2 = res_s2 * lhs_it2.index2();
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;
169 rhs_it1 = rhs_range.begin1();
170 rhs_it1 != rhs_range.end1();
173 rhs_it2 = rhs_it1.begin();
174 rhs_it2 != rhs_it1.end();
176 result(rhs_it2.index1(), rhs_it2.index2()) += lhs_val * *rhs_it2;
180 template<
typename LHS_T,
typename RHS_T >
182 nork(
const LHS_T& lhs,
const RHS_T& rhs,
const bool mono) ->
const
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;
199 throw error_t(
"matrix",
"nork: number of rows must not be 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");
209 auto result = RHS_T(res_s1, res_s2);
212 lhs_it1 = lhs.begin1();
213 lhs_it1 != lhs.end1();
216 lhs_it2 = lhs_it1.begin();
217 lhs_it2 != lhs_it1.end();
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;
226 template<
typename LHS_T,
typename RHS_T >
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);
244 lhs_it1 = lhs.begin1();
245 lhs_it1 != lhs.end1();
248 const auto lhs_it2 = lhs_it1.begin();
249 nork_range<LHS_T, RHS_T>(result, lhs_it2, rhs, res_s1, res_s2);
251 result /= norm_frob2_lhs;
256 template<
typename Matrix_T >
258 nnz(
const Matrix_T& m) ->
typename Matrix_T::size_type
260 using size_t =
typename Matrix_T::size_type;
261 auto result = size_t(0);
266 for (
auto& entry : it1)
273 template<
typename Matrix_T >
277 using Scalar_T =
typename Matrix_T::value_type;
282 for (
auto& entry : it1)
290 template<
typename Matrix_T >
294 using Scalar_T =
typename Matrix_T::value_type;
299 for (
auto& entry : it1)
307 template<
typename Matrix_T >
310 unit(
const typename Matrix_T::size_type dim) ->
const
313 using Scalar_T =
typename Matrix_T::value_type;
314 return ublas::identity_matrix<Scalar_T>(dim);
318 template<
typename LHS_T,
typename RHS_T >
321 const ublas::matrix_expression<RHS_T>& rhs) ->
const typename RHS_T::expression_type
323 using rhs_expression_t =
const RHS_T;
324 using matrix_row_t =
typename ublas::matrix_row<rhs_expression_t>;
326 const auto dim = lhs().size1();
328 auto result = RHS_T(dim, dim, dim);
330 lhs_row = lhs().begin1();
331 lhs_row != lhs().end1();
334 const auto& lhs_it = lhs_row.begin();
335 if (lhs_it != lhs_row.end())
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);
347 template<
typename LHS_T,
typename RHS_T >
351 const ublas::matrix_expression<RHS_T>& rhs) ->
const typename RHS_T::expression_type
353 using expression_t =
typename RHS_T::expression_type;
354 return ublas::sparse_prod<expression_t>(lhs(), rhs());
358 template<
typename LHS_T,
typename RHS_T >
361 prod(
const ublas::matrix_expression<LHS_T>& lhs,
362 const ublas::matrix_expression<RHS_T>& rhs) ->
const typename RHS_T::expression_type
364 const auto dim = lhs().size1();
365 RHS_T result(dim, dim);
366 ublas::axpy_prod(lhs, rhs, result,
true);
371 template<
typename Scalar_T,
typename LHS_T,
typename RHS_T >
373 inner(
const LHS_T& lhs,
const RHS_T& rhs) -> Scalar_T
375 auto result = Scalar_T(0);
377 lhs_it1 = lhs.begin1();
378 lhs_it1 != lhs.end1();
381 lhs_it2 = lhs_it1.begin();
382 lhs_it2 != lhs_it1.end();
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;
389 return result / lhs.size1();
393 template<
typename Matrix_T >
395 norm_frob2(
const Matrix_T& val) ->
typename Matrix_T::value_type
397 using Scalar_T =
typename Matrix_T::value_type;
399 auto result = Scalar_T(0);
401 val_it1 = val.begin1();
402 val_it1 != val.end1();
404 for (
auto& val_entry : val_it1)
408 result += val_entry * val_entry;
414 template<
typename Matrix_T >
416 trace(
const Matrix_T& val) ->
typename Matrix_T::value_type
418 using Scalar_T =
typename Matrix_T::value_type;
420 auto result = Scalar_T(0);
421 auto dim = val.size1();
423 ndx =
decltype(dim)(0);
427 const Scalar_T crd = val(ndx, ndx);
435#if defined(_GLUCAT_USE_BINDINGS)
437 template<
typename Matrix_T >
440 to_lapack(
const Matrix_T& val) -> ublas::matrix<double, ublas::column_major>
442 const auto s1 = val.size1();
443 const auto s2 = val.size2();
445 using lapack_matrix_t =
typename ublas::matrix<double, ublas::column_major>;
446 auto result = lapack_matrix_t(s1, s2);
449 using Scalar_T =
typename Matrix_T::value_type;
453 val_it1 = val.begin1();
454 val_it1 != val.end1();
457 val_it2 = val_it1.begin();
458 val_it2 != val_it1.end();
460 result(val_it2.index1(), val_it2.index2()) = traits_t::to_double(*val_it2);
466#if defined(_GLUCAT_USE_BLAZE)
468 template<
typename Matrix_T >
471 to_blaze(
const Matrix_T& val) -> blaze::DynamicMatrix<double,blaze::rowMajor>
473 const auto s1 = val.size1();
474 const auto s2 = val.size2();
476 using blaze_matrix_t =
typename blaze::DynamicMatrix<double,blaze::rowMajor>;
477 auto result = blaze_matrix_t(s1, s2);
479 using Scalar_T =
typename Matrix_T::value_type;
483 val_it1 = val.begin1();
484 val_it1 != val.end1();
487 val_it2 = val_it1.begin();
488 val_it2 != val_it1.end();
490 result(val_it2.index1(), val_it2.index2()) = traits_t::to_double(*val_it2);
498 template<
typename Matrix_T >
500 eigenvalues(
const Matrix_T& val) -> std::vector< std::complex<double> >
502 using complex_t = std::complex<double>;
503 using complex_vector_t =
typename std::vector<complex_t>;
505 const auto dim = val.size1();
506 auto lambda = complex_vector_t(dim);
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>;
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;
519 lapack::gees(
'N',
'N',
nullptr, T, sdim, real_lambda, imag_lambda, V);
522 k =
decltype(dim)(0);
525 lambda[k] = complex_t(real_lambda[k], imag_lambda[k]);
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>;
532 auto blaze_val = to_blaze(val);
533 auto blaze_lambda = blaze_complex_vector_t(dim);
534 blaze::geev(blaze_val, blaze_lambda);
537 k =
decltype(dim)(0);
540 lambda[k] = blaze_lambda[k];
546 template<
typename Matrix_T >
550 using Scalar_T =
typename Matrix_T::value_type;
553 using complex_t = std::complex<double>;
554 using complex_vector_t =
typename std::vector<complex_t>;
557 std::set<double> arg_set;
559 using vector_index_t =
typename complex_vector_t::size_type;
560 const auto dim = lambda.size();
562 std::max(std::numeric_limits<double>::epsilon(),
564 static const auto zero_eig_tol = 4096.0*
epsilon;
566 bool neg_real_eig_found =
false;
567 bool imag_eig_found =
false;
568 bool zero_eig_found =
false;
571 k =
decltype(dim)(0);
575 const auto lambda_k = lambda[k];
576 arg_set.insert(std::arg(lambda_k));
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);
582 if (!neg_real_eig_found &&
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 &&
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;
601 if (neg_real_eig_found)
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;
620 arg_it != arg_set.end();
623 const auto arg_diff = *arg_it - previous_arg;
624 if (arg_diff > best_diff)
626 best_diff = arg_diff;
627 best_arg = previous_arg;
629 previous_arg = *arg_it;
631 const auto arg_diff = first_arg + 2.0*pi - previous_arg;
632 if (arg_diff > best_diff)
634 best_diff = arg_diff;
635 best_arg = previous_arg;
637 result.
m_safe_arg = Scalar_T(pi - (best_arg + best_diff / 2.0));
Specific exception class.
Extra traits which extend numeric limits.
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 pi() -> Scalar_T
Pi.
static auto to_lapack(const Matrix_T &val) -> ublas::matrix< double, ublas::column_major >
Convert matrix to LAPACK format.
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.
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?