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"
41 # if defined(_GLUCAT_HAVE_BOOST_SERIALIZATION_ARRAY_WRAPPER_H)
42 # include <boost/serialization/array_wrapper.hpp>
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>
53 #if defined(_GLUCAT_USE_BINDINGS)
54 # include <boost/numeric/bindings/lapack/driver/gees.hpp>
55 # include <boost/numeric/bindings/ublas.hpp>
58 # if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS)
59 # pragma GCC diagnostic pop
64 namespace glucat {
namespace matrix
72 template<
typename LHS_T,
typename RHS_T >
75 kron(
const LHS_T& lhs,
const RHS_T& rhs)
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);
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();
93 for (lhs_const_iterator2
94 lhs_it2 = lhs_it1.begin();
95 lhs_it2 != lhs_it1.end();
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();
105 for (rhs_const_iterator2
106 rhs_it2 = rhs_it1.begin();
107 rhs_it2 != rhs_it1.end();
109 result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
115 template<
typename LHS_T,
typename RHS_T >
118 mono_kron(
const LHS_T& lhs,
const RHS_T& rhs)
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);
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();
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();
147 const rhs_const_iterator2 rhs_it2 = rhs_it1.begin();
148 result(start1 + rhs_it2.index1(), start2 + rhs_it2.index2()) = lhs_val * *rhs_it2;
155 template<
typename LHS_T,
typename RHS_T >
158 const typename LHS_T::const_iterator2 lhs_it2,
160 const typename RHS_T::size_type res_s1,
161 const typename RHS_T::size_type res_s2)
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();
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;
175 for (
typename matrix_range_t::const_iterator1
176 rhs_it1 = rhs_range.begin1();
177 rhs_it1 != rhs_range.end1();
179 for (
typename matrix_range_t::const_iterator2
180 rhs_it2 = rhs_it1.begin();
181 rhs_it2 != rhs_it1.end();
183 result(rhs_it2.index1(), rhs_it2.index2()) += lhs_val * *rhs_it2;
187 template<
typename LHS_T,
typename RHS_T >
190 nork(
const LHS_T& lhs,
const RHS_T& rhs,
const bool mono)
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);
209 throw error_t(
"matrix",
"nork: number of rows must not be 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");
219 matrix_t result(res_s1, res_s2);
221 for (
typename LHS_T::const_iterator1
222 lhs_it1 = lhs.begin1();
223 lhs_it1 != lhs.end1();
225 for (
typename LHS_T::const_iterator2
226 lhs_it2 = lhs_it1.begin();
227 lhs_it2 != lhs_it1.end();
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;
236 template<
typename LHS_T,
typename RHS_T >
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);
256 for (
typename LHS_T::const_iterator1
257 lhs_it1 = lhs.begin1();
258 lhs_it1 != lhs.end1();
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);
264 result /= norm_frob2_lhs;
269 template<
typename Matrix_T >
270 typename Matrix_T::size_type
271 nnz(
const Matrix_T& m)
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;
292 template<
typename Matrix_T >
294 isnan(
const Matrix_T& m)
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;
315 template<
typename Matrix_T >
319 unit(
const typename Matrix_T::size_type dim)
321 typedef typename Matrix_T::value_type Scalar_T;
322 return ublas::identity_matrix<Scalar_T>(dim);
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)
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;
340 const matrix_index_t dim = lhs().size1();
342 matrix_t result = matrix_t(dim, dim, dim);
343 for (lhs_const_iterator1
344 lhs_row = lhs().begin1();
345 lhs_row != lhs().end1();
348 const lhs_const_iterator2& lhs_it = lhs_row.begin();
349 if (lhs_it != lhs_row.end())
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);
361 template<
typename LHS_T,
typename RHS_T >
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)
367 typedef typename RHS_T::expression_type expression_t;
368 return ublas::sparse_prod<expression_t>(lhs(), rhs());
372 template<
typename LHS_T,
typename RHS_T >
374 const typename RHS_T::expression_type
375 prod(
const ublas::matrix_expression<LHS_T>& lhs,
376 const ublas::matrix_expression<RHS_T>& rhs)
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);
385 typedef typename RHS_T::expression_type expression_t;
386 return ublas::sparse_prod<expression_t>(lhs(), rhs());
391 template<
typename Scalar_T,
typename LHS_T,
typename RHS_T >
393 inner(
const LHS_T& lhs,
const RHS_T& rhs)
395 Scalar_T result = Scalar_T(0);
396 for (
typename LHS_T::const_iterator1
397 lhs_it1 = lhs.begin1();
398 lhs_it1 != lhs.end1();
400 for (
typename LHS_T::const_iterator2
401 lhs_it2 = lhs_it1.begin();
402 lhs_it2 != lhs_it1.end();
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;
409 return result / lhs.size1();
413 template<
typename Matrix_T >
414 typename Matrix_T::value_type
417 typedef typename Matrix_T::value_type Scalar_T;
419 Scalar_T result = Scalar_T(0);
420 for (
typename Matrix_T::const_iterator1
421 val_it1 = val.begin1();
422 val_it1 != val.end1();
424 for (
typename Matrix_T::const_iterator2
425 val_it2 = val_it1.begin();
426 val_it2 != val_it1.end();
431 result += (*val_it2) * (*val_it2);
437 template<
typename Matrix_T >
438 typename Matrix_T::value_type
439 trace(
const Matrix_T& val)
441 typedef typename Matrix_T::value_type Scalar_T;
442 typedef typename Matrix_T::size_type matrix_index_t;
444 Scalar_T result = Scalar_T(0);
445 matrix_index_t dim = val.size1();
446 for (matrix_index_t ndx=0;
450 const Scalar_T crd = val(ndx, ndx);
458 #if defined(_GLUCAT_USE_BINDINGS)
459 template<
typename Matrix_T >
462 ublas::matrix<double, ublas::column_major>
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();
469 typedef typename ublas::matrix<double, ublas::column_major> lapack_matrix_t;
470 lapack_matrix_t result(s1, s2);
473 typedef typename Matrix_T::value_type Scalar_T;
474 typedef numeric_traits<Scalar_T> traits_t;
476 typedef typename Matrix_T::const_iterator1 const_iterator1;
477 typedef typename Matrix_T::const_iterator2 const_iterator2;
479 val_it1 = val.begin1();
480 val_it1 != val.end1();
483 val_it2 = val_it1.begin();
484 val_it2 != val_it1.end();
486 result(val_it2.index1(), val_it2.index2()) = traits_t::to_double(*val_it2);
493 template<
typename Matrix_T >
494 ublas::vector< std::complex<double> >
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;
501 const matrix_index_t dim = val.size1();
502 complex_vector_t lambda = complex_vector_t(dim);
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;
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;
516 lapack::gees (
'N',
'N', (external_fp)0, T, sdim, real_lambda, imag_lambda, V );
519 for (vector_t::size_type k=0; k!= dim; ++k)
520 lambda[k] = complex_t(real_lambda[k], imag_lambda[k]);
526 template<
typename Matrix_T >
530 typedef typename Matrix_T::value_type Scalar_T;
531 eig_genus<Matrix_T> result;
533 result.m_safe_arg = Scalar_T(0);
535 typedef std::complex<double> complex_t;
536 typedef typename ublas::vector<complex_t> complex_vector_t;
539 std::set<double> arg_set;
541 typedef typename complex_vector_t::size_type vector_index_t;
542 const vector_index_t dim = lambda.size();
547 bool negative_eig_found =
false;
548 bool imaginary_eig_found =
false;
549 for (vector_index_t k=0; k != dim; ++k)
551 const complex_t lambda_k = lambda[k];
552 arg_set.insert(std::arg(lambda_k));
554 const double real_lambda_k =
std::real(lambda_k);
558 if (!negative_eig_found &&
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 &&
565 (real_lambda_k == 0.0 ||
566 real_lambda_k * real_lambda_k < norm_tol))
567 imaginary_eig_found =
true;
571 if (negative_eig_found)
573 if (imaginary_eig_found)
578 result.m_safe_arg = -
pi/2.0;
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;
590 arg_it != arg_set.end();
593 const double arg_diff = *arg_it - previous_arg;
594 if (arg_diff > best_diff)
596 best_diff = arg_diff;
597 best_arg = previous_arg;
599 previous_arg = *arg_it;
601 const double arg_diff = first_arg + 2.0*
pi - previous_arg;
602 if (arg_diff > best_diff)
604 best_diff = arg_diff;
605 best_arg = previous_arg;
607 result.m_safe_arg =
pi - (best_arg + best_diff / 2.0);
613 #endif // _GLUCAT_MATRIX_IMP_H