1 #ifndef _GLUCAT_MATRIX_IMP_H 2 #define _GLUCAT_MATRIX_IMP_H 37 # if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS) 38 # pragma GCC diagnostic push 39 # pragma GCC diagnostic ignored "-Wunused-local-typedefs" 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> 51 #if defined(_GLUCAT_USE_BINDINGS) 52 # include <boost/numeric/bindings/lapack/driver/gees.hpp> 53 # include <boost/numeric/bindings/ublas.hpp> 56 # if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS) 57 # pragma GCC diagnostic pop 62 namespace glucat {
namespace matrix
70 template<
typename LHS_T,
typename RHS_T >
73 kron(
const LHS_T& lhs,
const RHS_T& rhs)
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);
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();
91 for (lhs_const_iterator2
92 lhs_it2 = lhs_it1.begin();
93 lhs_it2 != lhs_it1.end();
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();
103 for (rhs_const_iterator2
104 rhs_it2 = rhs_it1.begin();
105 rhs_it2 != rhs_it1.end();
107 result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
113 template<
typename LHS_T,
typename RHS_T >
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);
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();
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();
145 const rhs_const_iterator2 rhs_it2 = rhs_it1.begin();
146 result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
153 template<
typename LHS_T,
typename RHS_T >
156 const typename LHS_T::const_iterator2 lhs_it2,
158 const typename RHS_T::size_type res_s1,
159 const typename RHS_T::size_type res_s2)
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();
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;
173 for (
typename matrix_range_t::const_iterator1
174 rhs_it1 = rhs_range.begin1();
175 rhs_it1 != rhs_range.end1();
177 for (
typename matrix_range_t::const_iterator2
178 rhs_it2 = rhs_it1.begin();
179 rhs_it2 != rhs_it1.end();
181 result(rhs_it2.index1(), rhs_it2.index2()) += lhs_val * *rhs_it2;
185 template<
typename LHS_T,
typename RHS_T >
188 nork(
const LHS_T& lhs,
const RHS_T& rhs,
const bool mono)
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);
207 throw error_t(
"matrix",
"nork: number of rows must not be 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");
217 matrix_t result(res_s1, res_s2);
219 for (
typename LHS_T::const_iterator1
220 lhs_it1 = lhs.begin1();
221 lhs_it1 != lhs.end1();
223 for (
typename LHS_T::const_iterator2
224 lhs_it2 = lhs_it1.begin();
225 lhs_it2 != lhs_it1.end();
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;
234 template<
typename LHS_T,
typename RHS_T >
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);
254 for (
typename LHS_T::const_iterator1
255 lhs_it1 = lhs.begin1();
256 lhs_it1 != lhs.end1();
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);
262 result /= norm_frob2_lhs;
267 template<
typename Matrix_T >
268 typename Matrix_T::size_type
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;
290 template<
typename Matrix_T >
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;
313 template<
typename Matrix_T >
317 unit(
const typename Matrix_T::size_type dim)
319 typedef typename Matrix_T::value_type Scalar_T;
320 return ublas::identity_matrix<Scalar_T>(dim);
324 template<
typename LHS_T,
typename RHS_T >
325 const typename RHS_T::expression_type
327 const ublas::matrix_expression<RHS_T>& rhs)
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;
338 const matrix_index_t dim = lhs().size1();
340 matrix_t result = matrix_t(dim, dim, dim);
341 for (lhs_const_iterator1
342 lhs_row = lhs().begin1();
343 lhs_row != lhs().end1();
346 const lhs_const_iterator2& lhs_it = lhs_row.begin();
347 if (lhs_it != lhs_row.end())
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);
359 template<
typename LHS_T,
typename RHS_T >
361 const typename RHS_T::expression_type
363 const ublas::matrix_expression<RHS_T>& rhs)
365 typedef typename RHS_T::expression_type expression_t;
366 return ublas::sparse_prod<expression_t>(lhs(), rhs());
370 template<
typename LHS_T,
typename RHS_T >
372 const typename RHS_T::expression_type
373 prod(
const ublas::matrix_expression<LHS_T>& lhs,
374 const ublas::matrix_expression<RHS_T>& rhs)
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);
383 typedef typename RHS_T::expression_type expression_t;
384 return ublas::sparse_prod<expression_t>(lhs(), rhs());
389 template<
typename Scalar_T,
typename LHS_T,
typename RHS_T >
391 inner(
const LHS_T& lhs,
const RHS_T& rhs)
393 Scalar_T result = Scalar_T(0);
394 for (
typename LHS_T::const_iterator1
395 lhs_it1 = lhs.begin1();
396 lhs_it1 != lhs.end1();
398 for (
typename LHS_T::const_iterator2
399 lhs_it2 = lhs_it1.begin();
400 lhs_it2 != lhs_it1.end();
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;
407 return result / lhs.size1();
411 template<
typename Matrix_T >
412 typename Matrix_T::value_type
415 typedef typename Matrix_T::value_type Scalar_T;
417 Scalar_T result = Scalar_T(0);
418 for (
typename Matrix_T::const_iterator1
419 val_it1 = val.begin1();
420 val_it1 != val.end1();
422 for (
typename Matrix_T::const_iterator2
423 val_it2 = val_it1.begin();
424 val_it2 != val_it1.end();
429 result += (*val_it2) * (*val_it2);
435 template<
typename Matrix_T >
436 typename Matrix_T::value_type
439 typedef typename Matrix_T::value_type Scalar_T;
440 typedef typename Matrix_T::size_type matrix_index_t;
442 Scalar_T result = Scalar_T(0);
443 matrix_index_t dim = val.size1();
444 for (matrix_index_t ndx=0;
448 const Scalar_T crd = val(ndx, ndx);
456 #if defined(_GLUCAT_USE_BINDINGS) 457 template<
typename Matrix_T >
460 ublas::matrix<double, ublas::column_major>
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();
467 typedef typename ublas::matrix<double, ublas::column_major> lapack_matrix_t;
468 lapack_matrix_t result(s1, s2);
471 typedef typename Matrix_T::value_type Scalar_T;
474 typedef typename Matrix_T::const_iterator1 const_iterator1;
475 typedef typename Matrix_T::const_iterator2 const_iterator2;
477 val_it1 = val.begin1();
478 val_it1 != val.end1();
481 val_it2 = val_it1.begin();
482 val_it2 != val_it1.end();
484 result(val_it2.index1(), val_it2.index2()) = traits_t::to_double(*val_it2);
491 template<
typename Matrix_T >
492 ublas::vector< std::complex<double> >
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;
499 const matrix_index_t dim = val.size1();
500 complex_vector_t lambda = complex_vector_t(dim);
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;
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;
514 lapack::gees (
'N',
'N', (external_fp)0, T, sdim, real_lambda, imag_lambda, V );
517 for (vector_t::size_type k=0; k!= dim; ++k)
518 lambda[k] = complex_t(real_lambda[k], imag_lambda[k]);
524 template<
typename Matrix_T >
528 typedef typename Matrix_T::value_type Scalar_T;
533 typedef std::complex<double> complex_t;
534 typedef typename ublas::vector<complex_t> complex_vector_t;
537 std::set<double> arg_set;
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(),
545 bool negative_eig_found =
false;
546 bool imaginary_eig_found =
false;
547 for (vector_index_t k=0; k != dim; ++k)
549 const complex_t lambda_k = lambda[k];
550 arg_set.insert(std::arg(lambda_k));
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);
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;
569 if (negative_eig_found)
571 if (imaginary_eig_found)
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;
588 arg_it != arg_set.end();
591 const double arg_diff = *arg_it - previous_arg;
592 if (arg_diff > best_diff)
594 best_diff = arg_diff;
595 best_arg = previous_arg;
597 previous_arg = *arg_it;
599 const double arg_diff = first_arg + 2.0*pi - previous_arg;
600 if (arg_diff > best_diff)
602 best_diff = arg_diff;
603 best_arg = previous_arg;
605 result.
m_safe_arg = pi - (best_arg + best_diff / 2.0);
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.
const RHS_T::expression_type prod(const ublas::matrix_expression< LHS_T > &lhs, const ublas::matrix_expression< RHS_T > &rhs)
Product of matrices.
static Scalar_T to_scalar_t(const Other_Scalar_T &val)
Cast to Scalar_T.
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.
const RHS_T mono_kron(const LHS_T &lhs, const RHS_T &rhs)
Sparse Kronecker tensor product of monomial matrices.
Scalar_T norm(const Multivector< Scalar_T, LO, HI > &val)
Scalar_T norm == sum of norm of coordinates.
Extra traits which extend numeric limits.
Matrix_T::value_type norm_frob2(const Matrix_T &val)
Square of Frobenius norm.
bool isnan(const Matrix_T &m)
Not a Number.
static Scalar_T NaN()
Smart NaN.
ublas::vector< std::complex< double > > eigenvalues(const Matrix_T &val)
Eigenvalues of a matrix.
Scalar_T imag(const Multivector< Scalar_T, LO, HI > &val)
Imaginary part: deprecated (always 0)
Structure containing classification of eigenvalues.
const Matrix_T unit(const typename Matrix_T::size_type n)
Unit matrix - as per Matlab eye.
eig_case_t m_eig_case
What kind of eigenvalues does the matrix contain?
static ublas::matrix< double, ublas::column_major > to_lapack(const Matrix_T &val)
Convert matrix to LAPACK format.
Specific exception class.
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.
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.
Matrix_T::value_type trace(const Matrix_T &val)
Matrix trace.
Scalar_T m_safe_arg
Argument such that exp(pi-m_safe_arg) lies between arguments of eigenvalues.
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.
eig_genus< Matrix_T > classify_eigenvalues(const Matrix_T &val)
Classify the eigenvalues of a matrix.
Scalar_T inner(const LHS_T &lhs, const RHS_T &rhs)
Inner product: sum(x(i,j)*y(i,j))/x.nrows()
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.
Matrix_T::size_type nnz(const Matrix_T &m)
Number of non-zeros.