This is an automated email from the git hooks/post-receive script. logari81 pushed a commit to branch master in repository getfem.
The following commit(s) were added to refs/heads/master by this push: new 57d84b08 Whitespace only 57d84b08 is described below commit 57d84b087d1e2a0927f18016a00923ee1d1a2ab7 Author: Konstantinos Poulios <logar...@gmail.com> AuthorDate: Thu Mar 23 20:11:18 2023 +0100 Whitespace only --- src/gmm/gmm_blas.h | 409 +++++++++++++++++++++-------------------- src/gmm/gmm_dense_lu.h | 48 ++--- src/gmm/gmm_lapack_interface.h | 24 +-- 3 files changed, 248 insertions(+), 233 deletions(-) diff --git a/src/gmm/gmm_blas.h b/src/gmm/gmm_blas.h index e41b5da1..4426bb7d 100644 --- a/src/gmm/gmm_blas.h +++ b/src/gmm/gmm_blas.h @@ -45,14 +45,14 @@ namespace gmm { /* ******************************************************************** */ - /* */ - /* Generic algorithms */ - /* */ + /* */ + /* Generic algorithms */ + /* */ /* ******************************************************************** */ /* ******************************************************************** */ - /* Miscellaneous */ + /* Miscellaneous */ /* ******************************************************************** */ /** clear (fill with zeros) a vector or matrix. */ @@ -78,7 +78,7 @@ namespace gmm { template <typename L> inline size_type nnz(const L& l, abstract_matrix) { return nnz(l, typename principal_orientation_type<typename - linalg_traits<L>::sub_orientation>::potype()); + linalg_traits<L>::sub_orientation>::potype()); } template <typename L> inline size_type nnz(const L& l, row_major) { @@ -113,16 +113,16 @@ namespace gmm { template <typename L> inline // to be optimized for dense vectors ... void fill(L& l, typename gmm::linalg_traits<L>::value_type x, - abstract_vector) { + abstract_vector) { for (size_type i = 0; i < vect_size(l); ++i) l[i] = x; } template <typename L> inline // to be optimized for dense matrices ... void fill(L& l, typename gmm::linalg_traits<L>::value_type x, - abstract_matrix) { + abstract_matrix) { for (size_type i = 0; i < mat_nrows(l); ++i) for (size_type j = 0; j < mat_ncols(l); ++j) - l(i,j) = x; + l(i,j) = x; } /** fill a vector or matrix with random value (uniform [-1,1]). */ @@ -132,7 +132,7 @@ namespace gmm { ///@cond DOXY_SHOW_ALL_FUNCTIONS template <typename L> inline void fill_random(const L& l) { fill_random(linalg_const_cast(l), - typename linalg_traits<L>::linalg_type()); + typename linalg_traits<L>::linalg_type()); } template <typename L> inline void fill_random(L& l, abstract_vector) { @@ -143,7 +143,7 @@ namespace gmm { template <typename L> inline void fill_random(L& l, abstract_matrix) { for (size_type i = 0; i < mat_nrows(l); ++i) for (size_type j = 0; j < mat_ncols(l); ++j) - l(i,j) = gmm::random(typename linalg_traits<L>::value_type()); + l(i,j) = gmm::random(typename linalg_traits<L>::value_type()); } ///@endcond @@ -157,19 +157,19 @@ namespace gmm { template <typename L> inline void fill_random(const L& l, double cfill) { fill_random(linalg_const_cast(l), cfill, - typename linalg_traits<L>::linalg_type()); + typename linalg_traits<L>::linalg_type()); } template <typename L> inline void fill_random(L& l, double cfill, abstract_vector) { typedef typename linalg_traits<L>::value_type T; size_type ntot = std::min(vect_size(l), - size_type(double(vect_size(l))*cfill) + 1); + size_type(double(vect_size(l))*cfill) + 1); for (size_type nb = 0; nb < ntot;) { size_type i = gmm::irandom(vect_size(l)); if (l[i] == T(0)) { - l[i] = gmm::random(typename linalg_traits<L>::value_type()); - ++nb; + l[i] = gmm::random(typename linalg_traits<L>::value_type()); + ++nb; } } } @@ -177,7 +177,7 @@ namespace gmm { template <typename L> inline void fill_random(L& l, double cfill, abstract_matrix) { fill_random(l, cfill, typename principal_orientation_type<typename - linalg_traits<L>::sub_orientation>::potype()); + linalg_traits<L>::sub_orientation>::potype()); } template <typename L> inline @@ -253,7 +253,7 @@ namespace gmm { /* ******************************************************************** */ - /* Scalar product */ + /* Scalar product */ /* ******************************************************************** */ ///@endcond @@ -264,8 +264,8 @@ namespace gmm { GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch, " << vect_size(v1) << " !=" << vect_size(v2)); return vect_sp(v1, v2, - typename linalg_traits<V1>::storage_type(), - typename linalg_traits<V2>::storage_type()); + typename linalg_traits<V1>::storage_type(), + typename linalg_traits<V2>::storage_type()); } /** scalar product between two vectors, using a matrix. @@ -277,7 +277,7 @@ namespace gmm { typename strongest_value_type3<V1,V2,MATSP>::value_type vect_sp(const MATSP &ps, const V1 &v1, const V2 &v2) { return vect_sp_with_mat(ps, v1, v2, - typename linalg_traits<MATSP>::sub_orientation()); + typename linalg_traits<MATSP>::sub_orientation()); } ///@cond DOXY_SHOW_ALL_FUNCTIONS @@ -285,13 +285,13 @@ namespace gmm { typename strongest_value_type3<V1,V2,MATSP>::value_type vect_sp_with_mat(const MATSP &ps, const V1 &v1, const V2 &v2, row_major) { return vect_sp_with_matr(ps, v1, v2, - typename linalg_traits<V2>::storage_type()); + typename linalg_traits<V2>::storage_type()); } template <typename MATSP, typename V1, typename V2> inline typename strongest_value_type3<V1,V2,MATSP>::value_type vect_sp_with_matr(const MATSP &ps, const V1 &v1, const V2 &v2, - abstract_sparse) { + abstract_sparse) { GMM_ASSERT2(vect_size(v1) == mat_ncols(ps) && vect_size(v2) == mat_nrows(ps), "dimensions mismatch"); size_type nr = mat_nrows(ps); @@ -306,13 +306,13 @@ namespace gmm { template <typename MATSP, typename V1, typename V2> inline typename strongest_value_type3<V1,V2,MATSP>::value_type vect_sp_with_matr(const MATSP &ps, const V1 &v1, const V2 &v2, - abstract_skyline) + abstract_skyline) { return vect_sp_with_matr(ps, v1, v2, abstract_sparse()); } template <typename MATSP, typename V1, typename V2> inline typename strongest_value_type3<V1,V2,MATSP>::value_type vect_sp_with_matr(const MATSP &ps, const V1 &v1, const V2 &v2, - abstract_dense) { + abstract_dense) { GMM_ASSERT2(vect_size(v1) == mat_ncols(ps) && vect_size(v2) == mat_nrows(ps), "dimensions mismatch"); typename linalg_traits<V2>::const_iterator @@ -332,13 +332,13 @@ namespace gmm { typename strongest_value_type3<V1,V2,MATSP>::value_type vect_sp_with_mat(const MATSP &ps, const V1 &v1, const V2 &v2,col_major){ return vect_sp_with_matc(ps, v1, v2, - typename linalg_traits<V1>::storage_type()); + typename linalg_traits<V1>::storage_type()); } template <typename MATSP, typename V1, typename V2> inline typename strongest_value_type3<V1,V2,MATSP>::value_type vect_sp_with_matc(const MATSP &ps, const V1 &v1, const V2 &v2, - abstract_sparse) { + abstract_sparse) { GMM_ASSERT2(vect_size(v1) == mat_ncols(ps) && vect_size(v2) == mat_nrows(ps), "dimensions mismatch"); typename linalg_traits<V1>::const_iterator @@ -352,13 +352,13 @@ namespace gmm { template <typename MATSP, typename V1, typename V2> inline typename strongest_value_type3<V1,V2,MATSP>::value_type vect_sp_with_matc(const MATSP &ps, const V1 &v1, const V2 &v2, - abstract_skyline) + abstract_skyline) { return vect_sp_with_matc(ps, v1, v2, abstract_sparse()); } template <typename MATSP, typename V1, typename V2> inline typename strongest_value_type3<V1,V2,MATSP>::value_type vect_sp_with_matc(const MATSP &ps, const V1 &v1, const V2 &v2, - abstract_dense) { + abstract_dense) { GMM_ASSERT2(vect_size(v1) == mat_ncols(ps) && vect_size(v2) == mat_nrows(ps), "dimensions mismatch"); typename linalg_traits<V1>::const_iterator @@ -377,7 +377,7 @@ namespace gmm { template <typename MATSP, typename V1, typename V2> inline typename strongest_value_type3<V1,V2,MATSP>::value_type vect_sp_with_mat(const MATSP &ps, const V1 &v1, const V2 &v2, - abstract_null_type) { + abstract_null_type) { typename temporary_vector<V1>::vector_type w(mat_nrows(ps)); GMM_WARNING2("Warning, a temporary is used in scalar product\n"); mult(ps, v1, w); @@ -386,7 +386,7 @@ namespace gmm { template <typename IT1, typename IT2> inline typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type, - typename std::iterator_traits<IT2>::value_type>::T + typename std::iterator_traits<IT2>::value_type>::T vect_sp_dense_(IT1 it, IT1 ite, IT2 it2) { typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type, typename std::iterator_traits<IT2>::value_type>::T res(0); @@ -396,10 +396,10 @@ namespace gmm { template <typename IT1, typename V> inline typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type, - typename linalg_traits<V>::value_type>::T + typename linalg_traits<V>::value_type>::T vect_sp_sparse_(IT1 it, IT1 ite, const V &v) { typename strongest_numeric_type<typename std::iterator_traits<IT1>::value_type, - typename linalg_traits<V>::value_type>::T res(0); + typename linalg_traits<V>::value_type>::T res(0); for (; it != ite; ++it) res += (*it) * v[it.index()]; return res; } @@ -408,7 +408,7 @@ namespace gmm { typename strongest_value_type<V1,V2>::value_type vect_sp(const V1 &v1, const V2 &v2, abstract_dense, abstract_dense) { return vect_sp_dense_(vect_const_begin(v1), vect_const_end(v1), - vect_const_begin(v2)); + vect_const_begin(v2)); } template <typename V1, typename V2> inline @@ -481,7 +481,7 @@ namespace gmm { while (it1 != ite1 && it2 != ite2) { if (it1.index() == it2.index()) - { res += (*it1) * *it2; ++it1; ++it2; } + { res += (*it1) * *it2; ++it1; ++it2; } else if (it1.index() < it2.index()) ++it1; else ++it2; } return res; @@ -497,12 +497,12 @@ namespace gmm { typename strongest_value_type<V1,V2>::value_type vect_sp(const V1 &v1, const V2 &v2,abstract_sparse,abstract_sparse) { return vect_sp_sparse_sparse(v1, v2, - typename linalg_and<typename linalg_traits<V1>::index_sorted, - typename linalg_traits<V2>::index_sorted>::bool_type()); + typename linalg_and<typename linalg_traits<V1>::index_sorted, + typename linalg_traits<V2>::index_sorted>::bool_type()); } /* ******************************************************************** */ - /* Hermitian product */ + /* Hermitian product */ /* ******************************************************************** */ ///@endcond /** Hermitian product. */ @@ -519,7 +519,7 @@ namespace gmm { } /* ******************************************************************** */ - /* Trace of a matrix */ + /* Trace of a matrix */ /* ******************************************************************** */ /** Trace of a matrix */ @@ -534,7 +534,7 @@ namespace gmm { } /* ******************************************************************** */ - /* Euclidean norm */ + /* Euclidean norm */ /* ******************************************************************** */ /** squared Euclidean norm of a vector. */ @@ -571,18 +571,18 @@ namespace gmm { R res(0); while (it1 != ite1 && it2 != ite2) { size_type i1 = index_of_it(it1, k1, - typename linalg_traits<V1>::storage_type()); + typename linalg_traits<V1>::storage_type()); size_type i2 = index_of_it(it2, k2, - typename linalg_traits<V2>::storage_type()); + typename linalg_traits<V2>::storage_type()); if (i1 == i2) { - res += gmm::abs_sqr(*it2 - *it1); ++it1; ++k1; ++it2; ++k2; + res += gmm::abs_sqr(*it2 - *it1); ++it1; ++k1; ++it2; ++k2; } else if (i1 < i2) { - res += gmm::abs_sqr(*it1); ++it1; ++k1; + res += gmm::abs_sqr(*it1); ++it1; ++k1; } else { - res += gmm::abs_sqr(*it2); ++it2; ++k2; + res += gmm::abs_sqr(*it2); ++it2; ++k2; } } while (it1 != ite1) { res += gmm::abs_sqr(*it1); ++it1; } @@ -625,8 +625,8 @@ namespace gmm { ::magnitude_type mat_euclidean_norm_sqr(const M &m) { return mat_euclidean_norm_sqr(m, - typename principal_orientation_type<typename - linalg_traits<M>::sub_orientation>::potype()); + typename principal_orientation_type<typename + linalg_traits<M>::sub_orientation>::potype()); } /** Euclidean norm of a matrix. */ @@ -637,7 +637,7 @@ namespace gmm { { return gmm::sqrt(mat_euclidean_norm_sqr(m)); } /* ******************************************************************** */ - /* vector norm1 */ + /* vector norm1 */ /* ******************************************************************** */ /** 1-norm of a vector */ template <typename V> @@ -646,7 +646,7 @@ namespace gmm { vect_norm1(const V &v) { auto it = vect_const_begin(v), ite = vect_const_end(v); typename number_traits<typename linalg_traits<V>::value_type> - ::magnitude_type res(0); + ::magnitude_type res(0); for (; it != ite; ++it) res += gmm::abs(*it); return res; } @@ -664,18 +664,18 @@ namespace gmm { R res(0); while (it1 != ite1 && it2 != ite2) { size_type i1 = index_of_it(it1, k1, - typename linalg_traits<V1>::storage_type()); + typename linalg_traits<V1>::storage_type()); size_type i2 = index_of_it(it2, k2, - typename linalg_traits<V2>::storage_type()); + typename linalg_traits<V2>::storage_type()); if (i1 == i2) { - res += gmm::abs(*it2 - *it1); ++it1; ++k1; ++it2; ++k2; + res += gmm::abs(*it2 - *it1); ++it1; ++k1; ++it2; ++k2; } else if (i1 < i2) { - res += gmm::abs(*it1); ++it1; ++k1; + res += gmm::abs(*it1); ++it1; ++k1; } else { - res += gmm::abs(*it2); ++it2; ++k2; + res += gmm::abs(*it2); ++it2; ++k2; } } while (it1 != ite1) { res += gmm::abs(*it1); ++it1; } @@ -684,7 +684,7 @@ namespace gmm { } /* ******************************************************************** */ - /* vector Infinity norm */ + /* vector Infinity norm */ /* ******************************************************************** */ /** Infinity norm of a vector. */ template <typename V> @@ -711,18 +711,18 @@ namespace gmm { R res(0); while (it1 != ite1 && it2 != ite2) { size_type i1 = index_of_it(it1, k1, - typename linalg_traits<V1>::storage_type()); + typename linalg_traits<V1>::storage_type()); size_type i2 = index_of_it(it2, k2, - typename linalg_traits<V2>::storage_type()); + typename linalg_traits<V2>::storage_type()); if (i1 == i2) { - res = std::max(res, gmm::abs(*it2 - *it1)); ++it1; ++k1; ++it2; ++k2; + res = std::max(res, gmm::abs(*it2 - *it1)); ++it1; ++k1; ++it2; ++k2; } else if (i1 < i2) { - res = std::max(res, gmm::abs(*it1)); ++it1; ++k1; + res = std::max(res, gmm::abs(*it1)); ++it1; ++k1; } else { - res = std::max(res, gmm::abs(*it2)); ++it2; ++k2; + res = std::max(res, gmm::abs(*it2)); ++it2; ++k2; } } while (it1 != ite1) { res = std::max(res, gmm::abs(*it1)); ++it1; } @@ -731,7 +731,7 @@ namespace gmm { } /* ******************************************************************** */ - /* matrix norm1 */ + /* matrix norm1 */ /* ******************************************************************** */ ///@cond DOXY_SHOW_ALL_FUNCTIONS template <typename M> @@ -758,7 +758,7 @@ namespace gmm { typename linalg_traits<M>::const_sub_row_type row = mat_const_row(m, i); auto it = vect_const_begin(row), ite = vect_const_end(row); for (size_type k = 0; it != ite; ++it, ++k) - aux[index_of_it(it, k, store_type())] += gmm::abs(*it); + aux[index_of_it(it, k, store_type())] += gmm::abs(*it); } return vect_norminf(aux); } @@ -785,7 +785,7 @@ namespace gmm { /* ******************************************************************** */ - /* matrix Infinity norm */ + /* matrix Infinity norm */ /* ******************************************************************** */ ///@cond DOXY_SHOW_ALL_FUNCTIONS template <typename M> @@ -812,7 +812,7 @@ namespace gmm { typename linalg_traits<M>::const_sub_col_type col = mat_const_col(m, i); auto it = vect_const_begin(col), ite = vect_const_end(col); for (size_type k = 0; it != ite; ++it, ++k) - aux[index_of_it(it, k, store_type())] += gmm::abs(*it); + aux[index_of_it(it, k, store_type())] += gmm::abs(*it); } return vect_norminf(aux); } @@ -838,7 +838,7 @@ namespace gmm { } /* ******************************************************************** */ - /* Max norm for matrices */ + /* Max norm for matrices */ /* ******************************************************************** */ ///@cond DOXY_SHOW_ALL_FUNCTIONS template <typename M> @@ -868,13 +868,13 @@ namespace gmm { typename number_traits<typename linalg_traits<M>::value_type> ::magnitude_type mat_maxnorm(const M &m) { - return mat_maxnorm(m, - typename principal_orientation_type<typename - linalg_traits<M>::sub_orientation>::potype()); + return mat_maxnorm + (m, typename principal_orientation_type + <typename linalg_traits<M>::sub_orientation>::potype()); } /* ******************************************************************** */ - /* Clean */ + /* Clean */ /* ******************************************************************** */ /** Clean a vector or matrix (replace near-zero entries with zeroes). */ @@ -909,9 +909,9 @@ namespace gmm { auto it = vect_begin(l), ite = vect_end(l); for (; it != ite; ++it){ if (gmm::abs((*it).real()) < T(threshold)) - *it = std::complex<T>(T(0), (*it).imag()); + *it = std::complex<T>(T(0), (*it).imag()); if (gmm::abs((*it).imag()) < T(threshold)) - *it = std::complex<T>((*it).real(), T(0)); + *it = std::complex<T>((*it).real(), T(0)); } } @@ -935,9 +935,9 @@ namespace gmm { } template <typename L> inline void clean(L &l, double threshold, - abstract_vector) { + abstract_vector) { gmm::clean(l, threshold, typename linalg_traits<L>::storage_type(), - typename linalg_traits<L>::value_type()); + typename linalg_traits<L>::value_type()); } template <typename L> inline void clean(const L &l, double threshold); @@ -953,10 +953,10 @@ namespace gmm { } template <typename L> inline void clean(L &l, double threshold, - abstract_matrix) { + abstract_matrix) { gmm::clean(l, threshold, - typename principal_orientation_type<typename - linalg_traits<L>::sub_orientation>::potype()); + typename principal_orientation_type<typename + linalg_traits<L>::sub_orientation>::potype()); } template <typename L> inline void clean(L &l, double threshold) @@ -966,7 +966,7 @@ namespace gmm { { gmm::clean(linalg_const_cast(l), threshold); } /* ******************************************************************** */ - /* Copy */ + /* Copy */ /* ******************************************************************** */ ///@endcond /** Copy vectors or matrices. @@ -977,10 +977,10 @@ namespace gmm { void copy(const L1& l1, L2& l2) { if ((const void *)(&l1) != (const void *)(&l2)) { if (same_origin(l1,l2)) - GMM_WARNING2("Warning : a conflict is possible in copy\n"); + GMM_WARNING2("Warning : a conflict is possible in copy\n"); copy(l1, l2, typename linalg_traits<L1>::linalg_type(), - typename linalg_traits<L2>::linalg_type()); + typename linalg_traits<L2>::linalg_type()); } } ///@cond DOXY_SHOW_ALL_FUNCTIONS @@ -993,7 +993,7 @@ namespace gmm { GMM_ASSERT2(vect_size(l1) == vect_size(l2), "dimensions mismatch, " << vect_size(l1) << " !=" << vect_size(l2)); copy_vect(l1, l2, typename linalg_traits<L1>::storage_type(), - typename linalg_traits<L2>::storage_type()); + typename linalg_traits<L2>::storage_type()); } template <typename L1, typename L2> inline @@ -1002,7 +1002,7 @@ namespace gmm { if (!m || !n) return; GMM_ASSERT2(n==mat_ncols(l2) && m==mat_nrows(l2), "dimensions mismatch"); copy_mat(l1, l2, typename linalg_traits<L1>::sub_orientation(), - typename linalg_traits<L2>::sub_orientation()); + typename linalg_traits<L2>::sub_orientation()); } template <typename V1, typename V2, typename C1, typename C2> inline @@ -1159,25 +1159,33 @@ namespace gmm { if (ite1 - it1 > 0) { clear(l2); auto it2 = vect_begin(l2), ite2 = vect_end(l2); - while (*(ite1-1) == typename linalg_traits<L1>::value_type(0)) ite1--; + while (*(ite1-1) == typename linalg_traits<L1>::value_type(0)) + ite1--; if (it2 == ite2) { - l2[it1.index()] = *it1; ++it1; - l2[ite1.index()-1] = *(ite1-1); --ite1; - if (it1 < ite1) - { it2 = vect_begin(l2); ++it2; std::copy(it1, ite1, it2); } - } - else { - ptrdiff_t m = it1.index() - it2.index(); - if (m >= 0 && ite1.index() <= ite2.index()) - std::copy(it1, ite1, it2 + m); - else { - if (m < 0) l2[it1.index()] = *it1; - if (ite1.index() > ite2.index()) l2[ite1.index()-1] = *(ite1-1); - it2 = vect_begin(l2); ite2 = vect_end(l2); - m = it1.index() - it2.index(); - std::copy(it1, ite1, it2 + m); - } + l2[it1.index()] = *it1; + ++it1; + l2[ite1.index()-1] = *(ite1-1); + --ite1; + if (it1 < ite1) { + it2 = vect_begin(l2); + ++it2; + std::copy(it1, ite1, it2); + } + } else { + ptrdiff_t m = it1.index() - it2.index(); + if (m >= 0 && ite1.index() <= ite2.index()) + std::copy(it1, ite1, it2 + m); + else { + if (m < 0) + l2[it1.index()] = *it1; + if (ite1.index() > ite2.index()) + l2[ite1.index()-1] = *(ite1-1); + it2 = vect_begin(l2); + ite2 = vect_end(l2); + m = it1.index() - it2.index(); + std::copy(it1, ite1, it2 + m); + } } } } @@ -1221,7 +1229,7 @@ namespace gmm { // cout << "*it = " << *it << endl; // cout << "it.index() = " << it.index() << endl; if (*it != (typename linalg_traits<L1>::value_type)(0)) - l2[it.index()] = *it; + l2[it.index()] = *it; } } @@ -1231,7 +1239,7 @@ namespace gmm { auto it = vect_const_begin(l1), ite = vect_const_end(l1); for (size_type i = 0; it != ite; ++it, ++i) if (*it != (typename linalg_traits<L1>::value_type)(0)) - l2[i] = *it; + l2[i] = *it; } template <typename L1, typename L2> // to be optimised ... @@ -1240,7 +1248,7 @@ namespace gmm { auto it = vect_const_begin(l1), ite = vect_const_end(l1); for (size_type i = 0; it != ite; ++it, ++i) if (*it != (typename linalg_traits<L1>::value_type)(0)) - l2[i] = *it; + l2[i] = *it; } @@ -1250,7 +1258,7 @@ namespace gmm { auto it = vect_const_begin(l1), ite = vect_const_end(l1); for (; it != ite; ++it) if (*it != (typename linalg_traits<L1>::value_type)(0)) - l2[it.index()] = *it; + l2[it.index()] = *it; } /* ******************************************************************** */ @@ -1278,17 +1286,17 @@ namespace gmm { GMM_ASSERT2(vect_size(l1) == vect_size(l2), "dimensions mismatch, " << vect_size(l1) << " !=" << vect_size(l2)); add(l1, l2, typename linalg_traits<L1>::storage_type(), - typename linalg_traits<L2>::storage_type()); + typename linalg_traits<L2>::storage_type()); } template <typename L1, typename L2> inline void add_spec(const L1& l1, L2& l2, abstract_matrix) { GMM_ASSERT2(mat_nrows(l1)==mat_nrows(l2) && mat_ncols(l1)==mat_ncols(l2), "dimensions mismatch l1 is " << mat_nrows(l1) << "x" - << mat_ncols(l1) << " and l2 is " << mat_nrows(l2) - << "x" << mat_ncols(l2)); + << mat_ncols(l1) << " and l2 is " << mat_nrows(l2) + << "x" << mat_ncols(l2)); add(l1, l2, typename linalg_traits<L1>::sub_orientation(), - typename linalg_traits<L2>::sub_orientation()); + typename linalg_traits<L2>::sub_orientation()); } template <typename L1, typename L2> @@ -1447,8 +1455,8 @@ namespace gmm { add(l1, l3); else add(l1, l2, l3, typename linalg_traits<L1>::storage_type(), - typename linalg_traits<L2>::storage_type(), - typename linalg_traits<L3>::storage_type()); + typename linalg_traits<L2>::storage_type(), + typename linalg_traits<L3>::storage_type()); } template <typename IT1, typename IT2, typename IT3> @@ -1466,7 +1474,7 @@ namespace gmm { template <typename IT1, typename IT2, typename IT3> void add_to_full_(IT1 it1, IT1 ite1, IT2 it2, IT2 ite2, - IT3 it3, IT3 ite3) { + IT3 it3, IT3 ite3) { typedef typename std::iterator_traits<IT3>::value_type T; IT3 it = it3; for (; it != ite3; ++it) *it = T(0); @@ -1476,36 +1484,36 @@ namespace gmm { template <typename L1, typename L2, typename L3> inline void add(const L1& l1, const L2& l2, L3& l3, - abstract_dense, abstract_dense, abstract_dense) { + abstract_dense, abstract_dense, abstract_dense) { add_full_(vect_const_begin(l1), vect_const_begin(l2), - vect_begin(l3), vect_end(l3)); + vect_begin(l3), vect_end(l3)); } // generic function for add(v1, v2, v3). // Need to be specialized to optimize particular additions. template <typename L1, typename L2, typename L3, - typename ST1, typename ST2, typename ST3> + typename ST1, typename ST2, typename ST3> inline void add(const L1& l1, const L2& l2, L3& l3, ST1, ST2, ST3) { copy(l2, l3); add(l1, l3, ST1(), ST3()); } template <typename L1, typename L2, typename L3> inline void add(const L1& l1, const L2& l2, L3& l3, - abstract_sparse, abstract_dense, abstract_dense) { + abstract_sparse, abstract_dense, abstract_dense) { add_almost_full_(vect_const_begin(l1), vect_const_end(l1), - vect_const_begin(l2), vect_begin(l3), vect_end(l3)); + vect_const_begin(l2), vect_begin(l3), vect_end(l3)); } template <typename L1, typename L2, typename L3> inline void add(const L1& l1, const L2& l2, L3& l3, - abstract_dense, abstract_sparse, abstract_dense) + abstract_dense, abstract_sparse, abstract_dense) { add(l2, l1, l3, abstract_sparse(), abstract_dense(), abstract_dense()); } template <typename L1, typename L2, typename L3> inline void add(const L1& l1, const L2& l2, L3& l3, - abstract_sparse, abstract_sparse, abstract_dense) { + abstract_sparse, abstract_sparse, abstract_dense) { add_to_full_(vect_const_begin(l1), vect_const_end(l1), - vect_const_begin(l2), vect_const_end(l2), - vect_begin(l3), vect_end(l3)); + vect_const_begin(l2), vect_const_end(l2), + vect_begin(l3), vect_end(l3)); } @@ -1517,11 +1525,11 @@ namespace gmm { while (it1 != ite1 && it2 != ite2) { ptrdiff_t d = it1.index() - it2.index(); if (d < 0) - { l3[it1.index()] += *it1; ++it1; } + { l3[it1.index()] += *it1; ++it1; } else if (d > 0) - { l3[it2.index()] += *it2; ++it2; } + { l3[it2.index()] += *it2; ++it2; } else - { l3[it1.index()] = *it1 + *it2; ++it1; ++it2; } + { l3[it1.index()] = *it1 + *it2; ++it1; ++it2; } } for (; it1 != ite1; ++it1) l3[it1.index()] += *it1; for (; it2 != ite2; ++it2) l3[it2.index()] += *it2; @@ -1533,10 +1541,10 @@ namespace gmm { template <typename L1, typename L2, typename L3> void add(const L1& l1, const L2& l2, L3& l3, - abstract_sparse, abstract_sparse, abstract_sparse) { + abstract_sparse, abstract_sparse, abstract_sparse) { add_spspsp(l1, l2, l3, typename linalg_and<typename - linalg_traits<L1>::index_sorted, - typename linalg_traits<L2>::index_sorted>::bool_type()); + linalg_traits<L1>::index_sorted, + typename linalg_traits<L2>::index_sorted>::bool_type()); } template <typename L1, typename L2> @@ -1558,13 +1566,13 @@ namespace gmm { while (ie1 && *(ite1-1) == T(0)) { ite1--; --ie1; } if (it2 == ite2 || i1 < it2.index()) { - l2[i1] = *it1; ++i1; ++it1; - if (it1 == ite1) return; - it2 = vect_begin(l2); ite2 = vect_end(l2); + l2[i1] = *it1; ++i1; ++it1; + if (it1 == ite1) return; + it2 = vect_begin(l2); ite2 = vect_end(l2); } if (ie1 > ite2.index()) { - --ite1; l2[ie1 - 1] = *ite1; - it2 = vect_begin(l2); + --ite1; l2[ie1 - 1] = *ite1; + it2 = vect_begin(l2); } it2 += i1 - it2.index(); for (; it1 != ite1; ++it1, ++it2) { *it2 += *it1; } @@ -1607,7 +1615,7 @@ namespace gmm { auto it1 = vect_const_begin(l1), ite1 = vect_const_end(l1); for (; it1 != ite1; ++it1) if (*it1 != typename linalg_traits<L1>::value_type(0)) - l2[it1.index()] += *it1; + l2[it1.index()] += *it1; } template <typename L1, typename L2> @@ -1622,12 +1630,12 @@ namespace gmm { auto it2 = vect_begin(l2), ite2 = vect_end(l2); while (*(ite1-1) == T1(0)) ite1--; if (it2 == ite2 || it1.index() < it2.index()) { - l2[it1.index()] = T2(0); - it2 = vect_begin(l2); ite2 = vect_end(l2); + l2[it1.index()] = T2(0); + it2 = vect_begin(l2); ite2 = vect_end(l2); } if (ite1.index() > ite2.index()) { - l2[ite1.index() - 1] = T2(0); - it2 = vect_begin(l2); + l2[ite1.index() - 1] = T2(0); + it2 = vect_begin(l2); } it2 += it1.index() - it2.index(); for (; it1 != ite1; ++it1, ++it2) *it2 += *it1; @@ -1642,7 +1650,7 @@ namespace gmm { } /* ******************************************************************** */ - /* Matrix-vector mult */ + /* Matrix-vector mult */ /* ******************************************************************** */ ///@endcond /** matrix-vector or matrix-matrix product. @@ -1667,12 +1675,12 @@ namespace gmm { GMM_ASSERT2(n==vect_size(l2) && m==vect_size(l3), "dimensions mismatch"); if (!same_origin(l2, l3)) mult_spec(l1, l2, l3, typename principal_orientation_type<typename - linalg_traits<L1>::sub_orientation>::potype()); + linalg_traits<L1>::sub_orientation>::potype()); else { GMM_WARNING2("Warning, A temporary is used for mult\n"); typename temporary_vector<L3>::vector_type temp(vect_size(l3)); mult_spec(l1, l2, temp, typename principal_orientation_type<typename - linalg_traits<L1>::sub_orientation>::potype()); + linalg_traits<L1>::sub_orientation>::potype()); copy(temp, l3); } } @@ -1705,8 +1713,8 @@ namespace gmm { auto itr = mat_row_const_begin(l1); for (; it != ite; ++it, ++itr) *it = vect_sp(linalg_traits<L1>::row(itr), l2, - typename linalg_traits<L1>::storage_type(), - typename linalg_traits<L2>::storage_type()); + typename linalg_traits<L1>::storage_type(), + typename linalg_traits<L2>::storage_type()); } template <typename L1, typename L2, typename L3> @@ -1760,14 +1768,14 @@ namespace gmm { GMM_ASSERT2(n==vect_size(l2) && m==vect_size(l4), "dimensions mismatch"); if (!same_origin(l2, l4)) { mult_add_spec(l1, l2, l4, typename principal_orientation_type<typename - linalg_traits<L1>::sub_orientation>::potype()); + linalg_traits<L1>::sub_orientation>::potype()); } else { GMM_WARNING2("Warning, A temporary is used for mult\n"); typename temporary_vector<L2>::vector_type temp(vect_size(l2)); copy(l2, temp); mult_add_spec(l1,temp, l4, typename principal_orientation_type<typename - linalg_traits<L1>::sub_orientation>::potype()); + linalg_traits<L1>::sub_orientation>::potype()); } } @@ -1784,14 +1792,14 @@ namespace gmm { GMM_ASSERT2(n==vect_size(l2) && m==vect_size(l3), "dimensions mismatch"); if (!same_origin(l2, l3)) { mult_add_spec(l1, l2, l3, typename principal_orientation_type<typename - linalg_traits<L1>::sub_orientation>::potype()); + linalg_traits<L1>::sub_orientation>::potype()); } else { GMM_WARNING2("Warning, A temporary is used for mult\n"); typename temporary_vector<L3>::vector_type temp(vect_size(l2)); copy(l2, temp); mult_add_spec(l1,temp, l3, typename principal_orientation_type<typename - linalg_traits<L1>::sub_orientation>::potype()); + linalg_traits<L1>::sub_orientation>::potype()); } } ///@cond DOXY_SHOW_ALL_FUNCTIONS @@ -1840,7 +1848,7 @@ namespace gmm { auto it = vect_const_begin(l2), ite = vect_const_end(l2); for (; it != ite; ++it) if (*it != typename linalg_traits<L2>::value_type(0)) - add(scaled(mat_const_col(l1, it.index()), *it), l3); + add(scaled(mat_const_col(l1, it.index()), *it), l3); } template <typename L1, typename L2, typename L3> @@ -1848,7 +1856,7 @@ namespace gmm { auto it = vect_const_begin(l2), ite = vect_const_end(l2); for (; it != ite; ++it) if (*it != typename linalg_traits<L2>::value_type(0)) - add(scaled(mat_const_col(l1, it.index()), *it), l3); + add(scaled(mat_const_col(l1, it.index()), *it), l3); } template <typename L1, typename L2, typename L3> inline @@ -1869,7 +1877,7 @@ namespace gmm { /* ******************************************************************** */ - /* Matrix-matrix mult */ + /* Matrix-matrix mult */ /* ******************************************************************** */ @@ -1953,22 +1961,22 @@ namespace gmm { size_type n = mat_ncols(l1); if (n == 0) { gmm::clear(l3); return; } GMM_ASSERT2(n == mat_nrows(l2) && mat_nrows(l1) == mat_nrows(l3) && - mat_ncols(l2) == mat_ncols(l3), "dimensions mismatch"); + mat_ncols(l2) == mat_ncols(l3), "dimensions mismatch"); if (same_origin(l2, l3) || same_origin(l1, l3)) { GMM_WARNING2("A temporary is used for mult"); temp_mat_type temp(mat_nrows(l3), mat_ncols(l3)); mult_spec(l1, l2, temp, typename mult_t< - typename linalg_traits<L1>::sub_orientation, - typename linalg_traits<L2>::sub_orientation, - typename linalg_traits<temp_mat_type>::sub_orientation>::t()); + typename linalg_traits<L1>::sub_orientation, + typename linalg_traits<L2>::sub_orientation, + typename linalg_traits<temp_mat_type>::sub_orientation>::t()); copy(temp, l3); } else mult_spec(l1, l2, l3, typename mult_t< - typename linalg_traits<L1>::sub_orientation, - typename linalg_traits<L2>::sub_orientation, - typename linalg_traits<L3>::sub_orientation>::t()); + typename linalg_traits<L1>::sub_orientation, + typename linalg_traits<L2>::sub_orientation, + typename linalg_traits<L3>::sub_orientation>::t()); } // Completely generic but inefficient @@ -1979,9 +1987,10 @@ namespace gmm { GMM_WARNING2("Inefficient generic matrix-matrix mult is used"); for (size_type i = 0; i < mat_nrows(l3) ; ++i) for (size_type j = 0; j < mat_ncols(l3) ; ++j) { - T a(0); - for (size_type k = 0; k < mat_nrows(l2) ; ++k) a += l1(i, k)*l2(k, j); - l3(i, j) = a; + T a(0); + for (size_type k = 0; k < mat_nrows(l2) ; ++k) + a += l1(i, k)*l2(k, j); + l3(i, j) = a; } } @@ -2007,20 +2016,20 @@ namespace gmm { void mult_spec(const L1& l1, const L2& l2, L3& l3, rcmult) { if (is_sparse(l1) && is_sparse(l2)) { GMM_WARNING3("Inefficient row matrix - col matrix mult for " - "sparse matrices, using temporary"); - mult_row_col_with_temp(l1, l2, l3, - typename principal_orientation_type<typename - linalg_traits<L3>::sub_orientation>::potype()); + "sparse matrices, using temporary"); + mult_row_col_with_temp + (l1, l2, l3, typename principal_orientation_type + <typename linalg_traits<L3>::sub_orientation>::potype()); } else { auto it2b = linalg_traits<L2>::col_begin(l2), it2 = it2b, - ite = linalg_traits<L2>::col_end(l2); + ite = linalg_traits<L2>::col_end(l2); size_type i,j, k = mat_nrows(l1); for (i = 0; i < k; ++i) { - typename linalg_traits<L1>::const_sub_row_type r1=mat_const_row(l1, i); - for (it2 = it2b, j = 0; it2 != ite; ++it2, ++j) - l3(i,j) = vect_sp(r1, linalg_traits<L2>::col(it2)); + typename linalg_traits<L1>::const_sub_row_type r1=mat_const_row(l1, i); + for (it2 = it2b, j = 0; it2 != ite; ++it2, ++j) + l3(i,j) = vect_sp(r1, linalg_traits<L2>::col(it2)); } } } @@ -2039,7 +2048,7 @@ namespace gmm { size_type nn = mat_nrows(l3), mm = mat_nrows(l2); for (size_type i = 0; i < nn; ++i) { for (size_type j = 0; j < mm; ++j) { - add(scaled(mat_const_row(l2, j), l1(i, j)), mat_row(l3, i)); + add(scaled(mat_const_row(l2, j), l1(i, j)), mat_row(l3, i)); } } } @@ -2053,7 +2062,7 @@ namespace gmm { typename linalg_traits<L1>::const_sub_row_type rl1=mat_const_row(l1, i); auto it = vect_const_begin(rl1), ite = vect_const_end(rl1); for (; it != ite; ++it) - add(scaled(mat_const_row(l2, it.index()), *it), mat_row(l3, i)); + add(scaled(mat_const_row(l2, it.index()), *it), mat_row(l3, i)); } } @@ -2065,29 +2074,29 @@ namespace gmm { template <typename L1, typename L2, typename L3> inline void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult) { - mult_spec(l1, l2,l3,c_mult(),typename linalg_traits<L2>::storage_type(), - typename linalg_traits<L2>::sub_orientation()); + mult_spec(l1, l2, l3, c_mult(), typename linalg_traits<L2>::storage_type(), + typename linalg_traits<L2>::sub_orientation()); } template <typename L1, typename L2, typename L3, typename ORIEN> void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult, - abstract_dense, ORIEN) { + abstract_dense, ORIEN) { typedef typename linalg_traits<L2>::value_type T; size_type nn = mat_ncols(l3), mm = mat_ncols(l1); for (size_type i = 0; i < nn; ++i) { clear(mat_col(l3, i)); for (size_type j = 0; j < mm; ++j) { - T b = l2(j, i); - if (b != T(0)) add(scaled(mat_const_col(l1, j), b), mat_col(l3, i)); + T b = l2(j, i); + if (b != T(0)) add(scaled(mat_const_col(l1, j), b), mat_col(l3, i)); } } } template <typename L1, typename L2, typename L3, typename ORIEN> void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult, - abstract_sparse, ORIEN) { + abstract_sparse, ORIEN) { // optimizable clear(l3); size_type nn = mat_ncols(l3); @@ -2095,27 +2104,27 @@ namespace gmm { typename linalg_traits<L2>::const_sub_col_type rc2 = mat_const_col(l2, i); auto it = vect_const_begin(rc2), ite = vect_const_end(rc2); for (; it != ite; ++it) - add(scaled(mat_const_col(l1, it.index()), *it), mat_col(l3, i)); + add(scaled(mat_const_col(l1, it.index()), *it), mat_col(l3, i)); } } template <typename L1, typename L2, typename L3> void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult, - abstract_sparse, row_major) { + abstract_sparse, row_major) { typedef typename linalg_traits<L2>::value_type T; GMM_WARNING3("Inefficient matrix-matrix mult for sparse matrices"); clear(l3); size_type mm = mat_nrows(l2), nn = mat_ncols(l3); for (size_type i = 0; i < nn; ++i) for (size_type j = 0; j < mm; ++j) { - T a = l2(i,j); - if (a != T(0)) add(scaled(mat_const_col(l1, j), a), mat_col(l3, i)); + T a = l2(i,j); + if (a != T(0)) add(scaled(mat_const_col(l1, j), a), mat_col(l3, i)); } } template <typename L1, typename L2, typename L3, typename ORIEN> inline void mult_spec(const L1& l1, const L2& l2, L3& l3, c_mult, - abstract_skyline, ORIEN) + abstract_skyline, ORIEN) { mult_spec(l1, l2, l3, c_mult(), abstract_sparse(), ORIEN()); } @@ -2146,7 +2155,7 @@ namespace gmm { typename linalg_traits<L1>::const_sub_col_type rc1 = mat_const_col(l1, i); auto it = vect_const_begin(rc1), ite = vect_const_end(rc1); for (; it != ite; ++it) - add(scaled(mat_const_row(l2, i), *it), mat_row(l3, it.index())); + add(scaled(mat_const_row(l2, i), *it), mat_row(l3, it.index())); } } @@ -2156,7 +2165,7 @@ namespace gmm { /* ******************************************************************** */ - /* Symmetry test. */ + /* Symmetry test. */ /* ******************************************************************** */ ///@endcond @@ -2165,8 +2174,9 @@ namespace gmm { @param tol a threshold. */ template <typename MAT> inline - bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol - = magnitude_of_linalg(MAT)(-1)) { + bool is_symmetric(const MAT &A, + magnitude_of_linalg(MAT) tol = magnitude_of_linalg(MAT)(-1)) + { typedef magnitude_of_linalg(MAT) R; if (tol < R(0)) tol = default_tol(R()) * mat_maxnorm(A); if (mat_nrows(A) != mat_ncols(A)) return false; @@ -2176,48 +2186,48 @@ namespace gmm { template <typename MAT> bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol, - abstract_dense) { + abstract_dense) { size_type m = mat_nrows(A); for (size_type i = 1; i < m; ++i) for (size_type j = 0; j < i; ++j) - if (gmm::abs(A(i, j)-A(j, i)) > tol) return false; + if (gmm::abs(A(i, j)-A(j, i)) > tol) return false; return true; } template <typename MAT> bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol, - abstract_sparse) { + abstract_sparse) { return is_symmetric(A, tol, typename principal_orientation_type<typename - linalg_traits<MAT>::sub_orientation>::potype()); + linalg_traits<MAT>::sub_orientation>::potype()); } template <typename MAT> bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol, - row_major) { + row_major) { for (size_type i = 0; i < mat_nrows(A); ++i) { typename linalg_traits<MAT>::const_sub_row_type row = mat_const_row(A, i); auto it = vect_const_begin(row), ite = vect_const_end(row); for (; it != ite; ++it) - if (gmm::abs(*it - A(it.index(), i)) > tol) return false; + if (gmm::abs(*it - A(it.index(), i)) > tol) return false; } return true; } template <typename MAT> bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol, - col_major) { + col_major) { for (size_type i = 0; i < mat_ncols(A); ++i) { typename linalg_traits<MAT>::const_sub_col_type col = mat_const_col(A, i); auto it = vect_const_begin(col), ite = vect_const_end(col); for (; it != ite; ++it) - if (gmm::abs(*it - A(i, it.index())) > tol) return false; + if (gmm::abs(*it - A(i, it.index())) > tol) return false; } return true; } template <typename MAT> bool is_symmetric(const MAT &A, magnitude_of_linalg(MAT) tol, - abstract_skyline) + abstract_skyline) { return is_symmetric(A, tol, abstract_sparse()); } ///@endcond @@ -2226,8 +2236,9 @@ namespace gmm { @param tol a threshold. */ template <typename MAT> inline - bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol - = magnitude_of_linalg(MAT)(-1)) { + bool is_hermitian(const MAT &A, + magnitude_of_linalg(MAT) tol = magnitude_of_linalg(MAT)(-1)) + { typedef magnitude_of_linalg(MAT) R; if (tol < R(0)) tol = default_tol(R()) * mat_maxnorm(A); if (mat_nrows(A) != mat_ncols(A)) return false; @@ -2237,48 +2248,48 @@ namespace gmm { template <typename MAT> bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol, - abstract_dense) { + abstract_dense) { size_type m = mat_nrows(A); for (size_type i = 1; i < m; ++i) for (size_type j = 0; j < i; ++j) - if (gmm::abs(A(i, j)-gmm::conj(A(j, i))) > tol) return false; + if (gmm::abs(A(i, j)-gmm::conj(A(j, i))) > tol) return false; return true; } template <typename MAT> bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol, - abstract_sparse) { + abstract_sparse) { return is_hermitian(A, tol, typename principal_orientation_type<typename - linalg_traits<MAT>::sub_orientation>::potype()); + linalg_traits<MAT>::sub_orientation>::potype()); } template <typename MAT> bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol, - row_major) { + row_major) { for (size_type i = 0; i < mat_nrows(A); ++i) { typename linalg_traits<MAT>::const_sub_row_type row = mat_const_row(A, i); auto it = vect_const_begin(row), ite = vect_const_end(row); for (; it != ite; ++it) - if (gmm::abs(gmm::conj(*it) - A(it.index(), i)) > tol) return false; + if (gmm::abs(gmm::conj(*it) - A(it.index(), i)) > tol) return false; } return true; } template <typename MAT> bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol, - col_major) { + col_major) { for (size_type i = 0; i < mat_ncols(A); ++i) { typename linalg_traits<MAT>::const_sub_col_type col = mat_const_col(A, i); auto it = vect_const_begin(col), ite = vect_const_end(col); for (; it != ite; ++it) - if (gmm::abs(gmm::conj(*it) - A(i, it.index())) > tol) return false; + if (gmm::abs(gmm::conj(*it) - A(i, it.index())) > tol) return false; } return true; } template <typename MAT> bool is_hermitian(const MAT &A, magnitude_of_linalg(MAT) tol, - abstract_skyline) + abstract_skyline) { return is_hermitian(A, tol, abstract_sparse()); } ///@endcond } diff --git a/src/gmm/gmm_dense_lu.h b/src/gmm/gmm_dense_lu.h index 33fbbcb1..1a8a26bc 100644 --- a/src/gmm/gmm_dense_lu.h +++ b/src/gmm/gmm_dense_lu.h @@ -132,24 +132,24 @@ namespace gmm { size_type info(0), i, j, jp, M(mat_nrows(A)), N(mat_ncols(A)); size_type NN = std::min(M, N); std::vector<T> c(M), r(N); - + GMM_ASSERT2(ipvt.size()+1 >= NN, "IPVT too small"); for (i = 0; i+1 < NN; ++i) ipvt.set(i, i); - + if (M || N) { for (j = 0; j+1 < NN; ++j) { - R max = gmm::abs(A(j,j)); jp = j; - for (i = j+1; i < M; ++i) /* find pivot. */ - if (gmm::abs(A(i,j)) > max) { jp = i; max = gmm::abs(A(i,j)); } - ipvt.set(j, jp + 1); - - if (max == R(0)) { info = j + 1; break; } + R max = gmm::abs(A(j,j)); jp = j; + for (i = j+1; i < M; ++i) /* find pivot. */ + if (gmm::abs(A(i,j)) > max) { jp = i; max = gmm::abs(A(i,j)); } + ipvt.set(j, jp + 1); + + if (max == R(0)) { info = j + 1; break; } if (jp != j) for (i = 0; i < N; ++i) std::swap(A(jp, i), A(j, i)); - + for (i = j+1; i < M; ++i) { A(i, j) /= A(j,j); c[i-j-1] = -A(i, j); } for (i = j+1; i < N; ++i) r[i-j-1] = A(j, i); // avoid the copy ? - rank_one_update(sub_matrix(A, sub_interval(j+1, M-j-1), - sub_interval(j+1, N-j-1)), c, conjugated(r)); + rank_one_update(sub_matrix(A, sub_interval(j+1, M-j-1), + sub_interval(j+1, N-j-1)), c, conjugated(r)); } ipvt.set(NN-1, NN); } @@ -159,14 +159,14 @@ namespace gmm { /** LU Solve : Solve equation Ax=b, given an LU factored matrix.*/ // Thanks to Valient Gough for this routine! template <typename DenseMatrix, typename VectorB, typename VectorX, - typename Pvector> + typename Pvector> void lu_solve(const DenseMatrix &LU, const Pvector& pvector, - VectorX &x, const VectorB &b) { + VectorX &x, const VectorB &b) { typedef typename linalg_traits<DenseMatrix>::value_type T; copy(b, x); for(size_type i = 0; i < pvector.size(); ++i) { size_type perm = pvector.get(i)-1; // permutations stored in 1's offset - if(i != perm) { T aux = x[i]; x[i] = x[perm]; x[perm] = aux; } + if (i != perm) { T aux = x[i]; x[i] = x[perm]; x[perm] = aux; } } /* solve Ax = b -> LUx = b -> Ux = L^-1 b. */ lower_tri_solve(LU, x, true); @@ -185,16 +185,20 @@ namespace gmm { } template <typename DenseMatrix, typename VectorB, typename VectorX, - typename Pvector> + typename Pvector> void lu_solve_transposed(const DenseMatrix &LU, const Pvector& pvector, - VectorX &x, const VectorB &b) { + VectorX &x, const VectorB &b) { typedef typename linalg_traits<DenseMatrix>::value_type T; copy(b, x); lower_tri_solve(transposed(LU), x, false); upper_tri_solve(transposed(LU), x, true); - for(size_type i = pvector.size(); i > 0; --i) { + for (size_type i = pvector.size(); i > 0; --i) { size_type perm = pvector.get(i-1)-1; // permutations stored in 1's offset - if(i-1 != perm) { T aux = x[i-1]; x[i-1] = x[perm]; x[perm] = aux; } + if (i-1 != perm) { + T aux = x[i-1]; + x[i-1] = x[perm]; + x[perm] = aux; + } } } @@ -202,7 +206,7 @@ namespace gmm { ///@cond DOXY_SHOW_ALL_FUNCTIONS template <typename DenseMatrixLU, typename DenseMatrix, typename Pvector> void lu_inverse(const DenseMatrixLU& LU, const Pvector& pvector, - DenseMatrix& AInv, col_major) { + DenseMatrix& AInv, col_major) { typedef typename linalg_traits<DenseMatrixLU>::value_type T; std::vector<T> tmp(pvector.size(), T(0)); std::vector<T> result(pvector.size()); @@ -216,7 +220,7 @@ namespace gmm { template <typename DenseMatrixLU, typename DenseMatrix, typename Pvector> void lu_inverse(const DenseMatrixLU& LU, const Pvector& pvector, - DenseMatrix& AInv, row_major) { + DenseMatrix& AInv, row_major) { typedef typename linalg_traits<DenseMatrixLU>::value_type T; std::vector<T> tmp(pvector.size(), T(0)); std::vector<T> result(pvector.size()); @@ -236,10 +240,10 @@ namespace gmm { /** Given an LU factored matrix, build the inverse of the matrix. */ template <typename DenseMatrixLU, typename DenseMatrix, typename Pvector> void lu_inverse(const DenseMatrixLU& LU, const Pvector& pvector, - const DenseMatrix& AInv_) { + const DenseMatrix& AInv_) { DenseMatrix& AInv = const_cast<DenseMatrix&>(AInv_); lu_inverse(LU, pvector, AInv, typename principal_orientation_type<typename - linalg_traits<DenseMatrix>::sub_orientation>::potype()); + linalg_traits<DenseMatrix>::sub_orientation>::potype()); } /** Given a dense matrix, build the inverse of the matrix, and diff --git a/src/gmm/gmm_lapack_interface.h b/src/gmm/gmm_lapack_interface.h index fa38addf..724ce0a1 100644 --- a/src/gmm/gmm_lapack_interface.h +++ b/src/gmm/gmm_lapack_interface.h @@ -124,9 +124,9 @@ namespace gmm { /* ********************************************************************* */ # define getrs_interface(f_name, trans1, lapack_name, base_type) inline \ - void f_name(const dense_matrix<base_type > &A, \ - const lapack_ipvt &ipvt, std::vector<base_type > &x, \ - const std::vector<base_type > &b) { \ + void f_name(const dense_matrix<base_type> &A, \ + const lapack_ipvt &ipvt, std::vector<base_type > &x, \ + const std::vector<base_type> &b) { \ GMMLAPACK_TRACE("getrs_interface"); \ BLAS_INT n = BLAS_INT(mat_nrows(A)), info(0), nrhs(1); \ gmm::copy(b, x); trans1; \ @@ -151,8 +151,8 @@ namespace gmm { /* ********************************************************************* */ # define getri_interface(lapack_name, base_type) inline \ - void lu_inverse(const dense_matrix<base_type > &LU, \ - const lapack_ipvt &ipvt, const dense_matrix<base_type > &A_) { \ + void lu_inverse(const dense_matrix<base_type> &LU, \ + const lapack_ipvt &ipvt, const dense_matrix<base_type> &A_) { \ GMMLAPACK_TRACE("getri_interface"); \ dense_matrix<base_type >& \ A = const_cast<dense_matrix<base_type > &>(A_); \ @@ -206,8 +206,8 @@ namespace gmm { BLAS_INT m = BLAS_INT(mat_nrows(A)), n=BLAS_INT(mat_ncols(A)); \ BLAS_INT info(0), lwork(-1); \ base_type work1; \ - if (m && n) { \ - std::copy(A.begin(), A.end(), Q.begin()); \ + if (m && n) { \ + std::copy(A.begin(), A.end(), Q.begin()); \ std::vector<base_type > tau(n); \ lapack_name1(&m, &n, &Q(0,0), &m, &tau[0], &work1 , &lwork, &info); \ lwork = BLAS_INT(gmm::real(work1)); \ @@ -363,11 +363,11 @@ namespace gmm { resize(S, n, n); copy(A, S); \ resize(Q, n, n); \ base_type rconde(0), rcondv(0); \ - BLAS_INT info(0); \ + BLAS_INT info(0); \ lapack_name(&jobvs, &sort, &select, &sense, &n, &S(0,0), &n, \ &sdim, &wr[0], &wi[0], &Q(0,0), &n, &rconde, &rcondv, \ &work[0], &lwork, &iwork[0], &liwork, &bwork[0], &info);\ - GMM_ASSERT1(!info, "SCHUR algorithm failed"); \ + GMM_ASSERT1(!info, "SCHUR algorithm failed"); \ } # define geesx_interface2(lapack_name, base_type) inline \ @@ -386,7 +386,7 @@ namespace gmm { resize(S, n, n); copy(A, S); \ resize(Q, n, n); \ base_type rconde(0), rcondv(0); \ - BLAS_INT info(0); \ + BLAS_INT info(0); \ lapack_name(&jobvs, &sort, &select, &sense, &n, &S(0,0), &n, \ &sdim, &w[0], &Q(0,0), &n, &rconde, &rcondv, \ &work[0], &lwork, &rwork[0], &bwork[0], &info); \ @@ -424,7 +424,7 @@ namespace gmm { resize(U, m, m); \ resize(Vtransposed, n, n); \ char job = 'A'; \ - BLAS_INT info(0); \ + BLAS_INT info(0); \ lapack_name(&job, &job, &m, &n, &X(0,0), &m, &sigma[0], &U(0,0), \ &m, &Vtransposed(0,0), &n, &work[0], &lwork, &info); \ } @@ -444,7 +444,7 @@ namespace gmm { resize(U, m, m); \ resize(Vtransposed, n, n); \ char job = 'A'; \ - BLAS_INT info(0); \ + BLAS_INT info(0); \ lapack_name(&job, &job, &m, &n, &X(0,0), &m, &sigma[0], &U(0,0), \ &m, &Vtransposed(0,0), &n, &work[0], &lwork, \ &rwork[0], &info); \