diff --git a/test/triqs/arrays/inverse2.cpp b/test/triqs/arrays/inverse2.cpp index 2e6a27cf..eaba9d45 100644 --- a/test/triqs/arrays/inverse2.cpp +++ b/test/triqs/arrays/inverse2.cpp @@ -76,7 +76,7 @@ int main(int argc, char **argv) { assert ( (abs(A(i,j) - Wkeep(i,j))) <1.e-10 ); double det = determinant(W); - std::cout<<"determinant "<>(inverse(D.matrix()))<< std::endl; std::cerr << "det_old = " << det_old << "detratio = "<< detratio<< " determin "<< D.determinant() < t2 = t + 2.4; + // tail + auto t = G(inf); TEST(t.order_min()); TEST( t( 2) ) ; diff --git a/triqs/arrays/linalg/det_and_inverse.hpp b/triqs/arrays/linalg/det_and_inverse.hpp index b4f41475..5f45f2ea 100644 --- a/triqs/arrays/linalg/det_and_inverse.hpp +++ b/triqs/arrays/linalg/det_and_inverse.hpp @@ -20,176 +20,182 @@ ******************************************************************************/ #ifndef TRIQS_ARRAYS_LINALG_DET_INV_H #define TRIQS_ARRAYS_LINALG_DET_INV_H - -#include -#include -#include #include "../impl/common.hpp" #include "../matrix.hpp" #include "../blas_lapack/getrf.hpp" #include "../blas_lapack/getri.hpp" -namespace triqs { namespace arrays { +namespace triqs { +namespace arrays { - /// Error which occurs during the matrix inversion + /// Error which occurs during the matrix inversion class matrix_inverse_exception : public triqs::runtime_error {}; /** - * Lazy result of inverse(M) where M can be : + * Lazy result of inverse(M) where M can be : * * a matrix, a matrix_view * * any matrix expression * The object is lazy, it does not do any computation. * It can be copied at no cost * It keeps view of the object A if it a matrix, a copy if it is a formal expression. - */ - template struct inverse_lazy; - template struct inverse_lazy_impl;//debug ony - - /// - template struct determinant_lazy; - - /// Lazy inversion - template inverse_lazy::type> inverse (A && a) { return {std::forward(a)}; } - - /// Lazy computation of det - template determinant_lazy determinant (A const & a) { return determinant_lazy(a); } + */ + template struct inverse_lazy; + + template inverse_lazy::type> inverse(A &&a) { + return {std::forward(a)}; + } // ----------------- implementation ----------------------------------------- - //worker takes a contiguous view and compute the det and inverse in two steps. - //it is separated in case of multiple use (no reallocation of ipvi, etc...) - //NB a view does not resize, only its elements can be changed - template class det_and_inverse_worker { - static_assert ( (is_matrix_view::value),"class must have be a view"); - typedef typename ViewType::value_type VT; - typedef matrix_view V_type; - ViewType V; - const size_t dim; - triqs::arrays::vector ipiv; - short step; + // worker takes a contiguous view and compute the det and inverse in two steps. + // it is separated in case of multiple use (no reallocation of ipvi, etc...) + // A can be a matrix, a matrix_view + template class det_and_inverse_worker { + typedef typename A::value_type value_type; + typedef matrix_view V_type; + A a; + int dim; + triqs::arrays::vector ipiv; + int step, info; + value_type _det; public: - det_and_inverse_worker (ViewType const & a): V(a), dim(first_dim(a)), ipiv(dim), step(0) { - if (first_dim(a)!=second_dim(a)) - TRIQS_RUNTIME_ERROR<<"Inverse/Det error : non-square matrix. Dimensions are : ("< V_type fortran_view (MT const &x) { return (x.indexmap().memory_layout_is_c() ? x.transpose() : x);} + template V_type fortran_view(MT const &x) { + if (x.indexmap().memory_layout_is_c()) + return x.transpose(); + else + return x; + } - void _step1(V_type & W) { - if (step >0) return; - step=1; + void _step1(V_type &W) { + if (step > 0) return; + step = 1; info = lapack::getrf(W, ipiv); - if (info<0) throw matrix_inverse_exception() << "Inverse/Det error : failure of getrf lapack routine "; + if (info < 0) throw matrix_inverse_exception() << "Inverse/Det error : failure of getrf lapack routine "; } - void _compute_det(V_type const & W) { - if (step>1) return; - _det =1; - for (size_t i =0; i 1) return; + _det = 1; + for (size_t i = 0; i < dim; i++) _det *= W(i, i); + bool flip = false; // compute the sign of the permutation + for (size_t i = 0; i < dim; i++) { + if (ipiv(i) != int(i) + 1) flip = !(flip); + } + _det = (flip ? -_det : _det); } - void _step2(V_type & W) { - assert(step==1); //if (step==1) return; - step=2; - _compute_det(W); + void _step2(V_type &W) { + assert(step == 1); // if (step==1) return; + step = 2; + _compute_det(W); info = lapack::getri(W, ipiv); - if (info!=0) throw matrix_inverse_exception() << "Inverse/Det error : matrix is not invertible"; + if (info != 0) throw matrix_inverse_exception() << "Inverse/Det error : matrix is not invertible"; } }; //----------------------------------------------------------- - // an implementation class to gather the common part to matrix and expression.... - template struct inverse_lazy_impl : TRIQS_CONCEPT_TAG_NAME(ImmutableMatrix) { - public: - + template struct inverse_lazy : TRIQS_CONCEPT_TAG_NAME(ImmutableMatrix) { typedef typename std::remove_reference::type A_t; typedef typename std::remove_const::type value_type; typedef typename A_t::domain_type domain_type; - A a; - template inverse_lazy_impl(AA && a_):a (std::forward(a_)) { - if (first_dim(a) != second_dim(a)) TRIQS_RUNTIME_ERROR<< "Inverse : matrix is not square but of size "<< first_dim(a)<<" x "<< second_dim(a); + typedef matrix M_type; + typedef matrix_view M_view_type; + + template inverse_lazy(AA &&a_) : a(std::forward(a_)), M{}, computed{false} { + if (first_dim(a) != second_dim(a)) + TRIQS_RUNTIME_ERROR << "Inverse : matrix is not square but of size " << first_dim(a) << " x " << second_dim(a); + } + + domain_type domain() const { return a.domain(); } + A const & input() const { return a;} + + template value_type const &operator()(K0 const &k0, K1 const &k1) const { + activate(); + return M(k0, k1); + } + + M_type const &operator()() const { + activate(); + return M; + } + + friend std::ostream &operator<<(std::ostream &out, inverse_lazy const &x) { return out << "inverse(" << x.a << ")"; } + + private: + A a; + mutable M_type M; + mutable bool computed; + + void activate() const { + if (computed) return; + M = a; + auto worker = det_and_inverse_worker {M}; + worker.inverse(); + computed = true; } - domain_type domain() const { return a.domain(); } - template value_type operator() (K0 const & k0, K1 const & k1) const { activate(); return _id->M(k0,k1); } - friend std::ostream & operator<<(std::ostream & out,inverse_lazy_impl const&x){return out<<"inverse("< M_type; - typedef matrix_view M_view_type; - M_type M; - internal_data(inverse_lazy_impl const & P):M(P.a){ - det_and_inverse_worker worker(M); - worker.inverse(); - } - }; - friend struct internal_data; - mutable std::shared_ptr _id; - void activate() const { if (!_id) _id= std::make_shared(*this);} }; - // The general case - template - struct inverse_lazy::type>) > : inverse_lazy_impl { - template inverse_lazy(AA && a_):inverse_lazy_impl(std::forward(a_)) {} - }; - - // for matrix and matrix_views, we have more optimisation possible .... - template - struct inverse_lazy::type>) >:inverse_lazy_impl{ - template inverse_lazy(AA && a_):inverse_lazy_impl(std::forward(a_)) {} - - template // Optimized implementation of = - friend void triqs_arrays_assign_delegation (MT & lhs, inverse_lazy const & rhs) { - //std::cerr << " DELEGATING1"<< lhs <::value, "Internal error"); - //std::cerr << " DELEGATING"<< lhs << std::endl; - if ((lhs.indexmap().memory_indices_layout() !=rhs.a.indexmap().memory_indices_layout())|| - (lhs.data_start() != rhs.a.data_start()) || !(has_contiguous_data(lhs))) { rhs.activate(); lhs = rhs._id->M;} - else {// if M = inverse(M) with the SAME object, then we do not need to copy the data - //std::cerr << " DELEGATING"<< lhs << std::endl; - blas_lapack_tools::reflexive_qcache C(lhs);// a reflexive cache will use a temporary "regrouping" copy if and only if needed - det_and_inverse_worker W(C());// the worker to make the inversion of the lhs... - W.inverse(); // worker is working ... - //std::cerr << " DELEGATING"<< lhs << std::endl; - } - //std::cerr << " okok "<< lhs << std::endl; - } - friend std::ostream & operator<<(std::ostream & out,inverse_lazy const&x){return out<<"inverse("< + ENABLE_IF(is_matrix_or_view::type>) + triqs_arrays_assign_delegation(MT &lhs, inverse_lazy const &rhs) { + static_assert(is_matrix_or_view::value, "Can only assign an inverse matrix to a matrix or a matrix_view"); + bool M_eq_inverse_M = ((lhs.indexmap().memory_indices_layout() == rhs.input().indexmap().memory_indices_layout()) && + (lhs.data_start() == rhs.input().data_start()) && (has_contiguous_data(lhs))); + if (!M_eq_inverse_M) { + lhs = rhs(); + } else { + blas_lapack_tools::reflexive_qcache C(lhs); + // a reflexive cache will use a temporary "regrouping" copy if and only if needed + det_and_inverse_worker W(C()); // the worker to make the inversion of the lhs... + W.inverse(); // worker is working ... + } + } //------------------- det ---------------------------------------- - template struct determinant_lazy { - typedef typename A::value_type value_type; - typedef typename const_view_type_if_exists_else_type::type A_type; - A_type a; - determinant_lazy(A const & a_):a(a_){} - operator value_type() { activate(); return _id->det; } - value_type const & operator()() { activate(); return _id->det; } - friend std::ostream & operator<<(std::ostream & out, determinant_lazy const & x){ return out<<"determinant("< worker(M); det = worker.det();} - }; - friend struct internal_data; - mutable std::shared_ptr _id; - void activate() const { if (!_id) _id= std::make_shared(*this);} - }; + template typename std::remove_reference::type::value_type determinant(A &&a) { + // makes a temporary copy of A if A is a const & + // If a is a matrix &&, it is moved into the worker. + auto worker = det_and_inverse_worker::type::value_type>>(std::forward(a)); + return worker.det(); + } +} - }} // namespace triqs::arrays +#ifndef TRIQS_COMPILER_OBSOLETE_GCC +namespace clef { + TRIQS_CLEF_MAKE_FNT_LAZY(determinant); +} +#endif + +} // namespace triqs::arrays #endif