From c8fcb40164ed8767b3a2b9337bd32967ed72f27e Mon Sep 17 00:00:00 2001 From: Olivier Parcollet Date: Thu, 22 Aug 2013 11:41:17 +0200 Subject: [PATCH] arrays : rm linalg/determinant and inverse (links) It is clearer to have only one file : det_and_inverse.hpp changed the included in tests. --- test/speed/array1.cpp | 4 +- test/speed/array_short1.cpp | 4 +- test/triqs/arrays/a_x_ty.cpp | 4 +- test/triqs/arrays/blas_lapack.cpp | 4 +- test/triqs/arrays/mat_vec_mul.cpp | 4 +- test/triqs/arrays/mat_vec_mul2.cpp | 4 +- triqs/arrays.hpp | 3 +- .../expression_template/matrix_algebra.hpp | 2 +- triqs/arrays/linalg/det_and_inverse.hpp | 2 +- triqs/arrays/linalg/determinant.hpp | 187 ------------------ triqs/arrays/linalg/inverse.hpp | 187 ------------------ triqs/arrays/matrix_stack_view.hpp | 2 +- triqs/det_manip/det_manip.hpp | 3 +- .../bravais_lattice_and_brillouin_zone.cpp | 2 +- 14 files changed, 18 insertions(+), 394 deletions(-) delete mode 100644 triqs/arrays/linalg/determinant.hpp delete mode 100644 triqs/arrays/linalg/inverse.hpp diff --git a/test/speed/array1.cpp b/test/speed/array1.cpp index 263ad379..a6bd5629 100644 --- a/test/speed/array1.cpp +++ b/test/speed/array1.cpp @@ -22,8 +22,8 @@ #include #include #include -#include -#include +#include +#include using namespace triqs::arrays; using namespace triqs; diff --git a/test/speed/array_short1.cpp b/test/speed/array_short1.cpp index 640c7a23..d7c088a3 100644 --- a/test/speed/array_short1.cpp +++ b/test/speed/array_short1.cpp @@ -22,8 +22,8 @@ #include #include #include -#include -#include +#include +#include using namespace triqs::arrays; using namespace triqs; diff --git a/test/triqs/arrays/a_x_ty.cpp b/test/triqs/arrays/a_x_ty.cpp index 49e37328..56bbda2a 100644 --- a/test/triqs/arrays/a_x_ty.cpp +++ b/test/triqs/arrays/a_x_ty.cpp @@ -26,8 +26,8 @@ #include "./src/matrix.hpp" #include "./src/linalg/matmul.hpp" #include "./src/linalg/mat_vec_mul.hpp" -#include "./src/linalg/inverse.hpp" -#include "./src/linalg/determinant.hpp" +#include "./src/linalg/det_and_inverse.hpp" +#include "./src/linalg/det_and_inverse.hpp" #include "./src/linalg/a_x_ty.hpp" #include "./src/blas_lapack/dot.hpp" diff --git a/test/triqs/arrays/blas_lapack.cpp b/test/triqs/arrays/blas_lapack.cpp index d1e7dd44..c3ee4eba 100644 --- a/test/triqs/arrays/blas_lapack.cpp +++ b/test/triqs/arrays/blas_lapack.cpp @@ -25,8 +25,8 @@ #include "./src/matrix.hpp" #include "./src/linalg/matmul.hpp" #include "./src/linalg/mat_vec_mul.hpp" -#include "./src/linalg/inverse.hpp" -#include "./src/linalg/determinant.hpp" +#include "./src/linalg/det_and_inverse.hpp" +#include "./src/linalg/det_and_inverse.hpp" #include "./src/blas_lapack/gemm.hpp" #include "./src/blas_lapack/gemv.hpp" diff --git a/test/triqs/arrays/mat_vec_mul.cpp b/test/triqs/arrays/mat_vec_mul.cpp index 922ff88a..db02659b 100644 --- a/test/triqs/arrays/mat_vec_mul.cpp +++ b/test/triqs/arrays/mat_vec_mul.cpp @@ -23,8 +23,8 @@ #include "./src/array.hpp" #include "./src/vector.hpp" #include "./src/matrix.hpp" -#include "./src/linalg/inverse.hpp" -#include "./src/linalg/determinant.hpp" +#include "./src/linalg/det_and_inverse.hpp" +#include "./src/linalg/det_and_inverse.hpp" #include using std::cout; using std::endl; diff --git a/test/triqs/arrays/mat_vec_mul2.cpp b/test/triqs/arrays/mat_vec_mul2.cpp index e0ff8b11..4c71270b 100644 --- a/test/triqs/arrays/mat_vec_mul2.cpp +++ b/test/triqs/arrays/mat_vec_mul2.cpp @@ -24,8 +24,8 @@ #include "./src/vector.hpp" #include "./src/matrix.hpp" #include "./src/asserts.hpp" -#include "./src/linalg/inverse.hpp" -#include "./src/linalg/determinant.hpp" +#include "./src/linalg/det_and_inverse.hpp" +#include "./src/linalg/det_and_inverse.hpp" #include using std::cout; using std::endl; diff --git a/triqs/arrays.hpp b/triqs/arrays.hpp index fcf27b4e..f9bd4b13 100644 --- a/triqs/arrays.hpp +++ b/triqs/arrays.hpp @@ -34,8 +34,7 @@ #include // Linear algebra ?? Keep here ? -//#include -//#include +//#include #endif diff --git a/triqs/arrays/expression_template/matrix_algebra.hpp b/triqs/arrays/expression_template/matrix_algebra.hpp index b66779ea..bd21e673 100644 --- a/triqs/arrays/expression_template/matrix_algebra.hpp +++ b/triqs/arrays/expression_template/matrix_algebra.hpp @@ -24,7 +24,7 @@ #include "../matrix.hpp" #include "../linalg/matmul.hpp" #include "../linalg/mat_vec_mul.hpp" -#include "../linalg/inverse.hpp" +#include "../linalg/det_and_inverse.hpp" namespace triqs { namespace arrays { template diff --git a/triqs/arrays/linalg/det_and_inverse.hpp b/triqs/arrays/linalg/det_and_inverse.hpp index d25dc4c4..e8856ff2 100644 --- a/triqs/arrays/linalg/det_and_inverse.hpp +++ b/triqs/arrays/linalg/det_and_inverse.hpp @@ -111,7 +111,7 @@ namespace triqs { namespace arrays { // an implementation class to gather the common part to matrix and expression.... template struct inverse_lazy_impl : TRIQS_MODEL_CONCEPT(ImmutableMatrix) { - typedef typename boost::remove_const::type value_type; + typedef typename std::remove_const::type value_type; typedef typename A::domain_type domain_type; typedef typename domain_type::index_value_type index_value_type; typedef typename const_view_type_if_exists_else_type::type A_type; diff --git a/triqs/arrays/linalg/determinant.hpp b/triqs/arrays/linalg/determinant.hpp deleted file mode 100644 index 54d587be..00000000 --- a/triqs/arrays/linalg/determinant.hpp +++ /dev/null @@ -1,187 +0,0 @@ -/******************************************************************************* - * - * TRIQS: a Toolbox for Research in Interacting Quantum Systems - * - * Copyright (C) 2011 by O. Parcollet - * - * TRIQS is free software: you can redistribute it and/or modify it under the - * terms of the GNU General Public License as published by the Free Software - * Foundation, either version 3 of the License, or (at your option) any later - * version. - * - * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - * details. - * - * You should have received a copy of the GNU General Public License along with - * TRIQS. If not, see . - * - ******************************************************************************/ -#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 { - - /// Error which occurs during the matrix inversion - class matrix_inverse_exception : public triqs::runtime_error {}; - - /** - * 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 determinant_lazy; - - /// Lazy inversion - template inverse_lazy inverse (A const & a) { return inverse_lazy(a); } - - /// Lazy computation of det - template determinant_lazy determinant (A const & a) { return determinant_lazy(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; - - public: - det_and_inverse_worker (ViewType const & a): V(a), dim(a.dim0()), ipiv(dim), step(0) { - if (a.dim0()!=a.dim1()) - 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);} - - 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 "; - } - - void _compute_det(V_type const & W) { - if (step>1) return; - _det =1; - for (size_t i =0; i struct inverse_lazy_impl : TRIQS_MODEL_CONCEPT(ImmutableMatrix) { - typedef typename boost::remove_const::type value_type; - typedef typename A::domain_type domain_type; - typedef typename domain_type::index_value_type index_value_type; - typedef typename const_view_type_if_exists_else_type::type A_type; - const A_type a; - inverse_lazy_impl(A const & a_):a (a_) { - if (a.dim0() != a.dim1()) TRIQS_RUNTIME_ERROR<< "Inverse : matrix is not square but of size "<< a.dim0()<<" x "<< a.dim1(); - } - domain_type domain() const { return a.domain(); } - size_t dim0() const { return a.dim0();} - size_t dim1() const { return a.dim1();} - 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 boost::shared_ptr _id; - void activate() const { if (!_id) _id= boost::make_shared(*this);} - }; - - // The general case - template - struct inverse_lazy >::type > : inverse_lazy_impl { - inverse_lazy(A const & a):inverse_lazy_impl(a) { } - }; - - // for matrix and matrix_views, we have more optimisation possible .... - template - struct inverse_lazy >::type >:inverse_lazy_impl{ - inverse_lazy(A const & a):inverse_lazy_impl(a) { } - - template // Optimized implementation of = - friend void triqs_arrays_assign_delegation (MT & lhs, inverse_lazy const & rhs) { - static_assert(is_matrix_or_view::value, "Internal error"); - 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 - 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 ... - } - } - friend std::ostream & operator<<(std::ostream & out,inverse_lazy const&x){return out<<"inverse("< struct determinant_lazy { // : { Tag::expression_terminal, Tag::scalar_expression_terminal { - 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 boost::shared_ptr _id; - void activate() const { if (!_id) _id= boost::make_shared(*this);} - }; - - }} // namespace triqs::arrays -#endif diff --git a/triqs/arrays/linalg/inverse.hpp b/triqs/arrays/linalg/inverse.hpp deleted file mode 100644 index e8856ff2..00000000 --- a/triqs/arrays/linalg/inverse.hpp +++ /dev/null @@ -1,187 +0,0 @@ -/******************************************************************************* - * - * TRIQS: a Toolbox for Research in Interacting Quantum Systems - * - * Copyright (C) 2011 by O. Parcollet - * - * TRIQS is free software: you can redistribute it and/or modify it under the - * terms of the GNU General Public License as published by the Free Software - * Foundation, either version 3 of the License, or (at your option) any later - * version. - * - * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - * details. - * - * You should have received a copy of the GNU General Public License along with - * TRIQS. If not, see . - * - ******************************************************************************/ -#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 { - - /// Error which occurs during the matrix inversion - class matrix_inverse_exception : public triqs::runtime_error {}; - - /** - * 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 determinant_lazy; - - /// Lazy inversion - template inverse_lazy inverse (A const & a) { return inverse_lazy(a); } - - /// Lazy computation of det - template determinant_lazy determinant (A const & a) { return determinant_lazy(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; - - public: - det_and_inverse_worker (ViewType const & a): V(a), dim(a.dim0()), ipiv(dim), step(0) { - if (a.dim0()!=a.dim1()) - 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);} - - 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 "; - } - - void _compute_det(V_type const & W) { - if (step>1) return; - _det =1; - for (size_t i =0; i struct inverse_lazy_impl : TRIQS_MODEL_CONCEPT(ImmutableMatrix) { - typedef typename std::remove_const::type value_type; - typedef typename A::domain_type domain_type; - typedef typename domain_type::index_value_type index_value_type; - typedef typename const_view_type_if_exists_else_type::type A_type; - const A_type a; - inverse_lazy_impl(A const & a_):a (a_) { - if (a.dim0() != a.dim1()) TRIQS_RUNTIME_ERROR<< "Inverse : matrix is not square but of size "<< a.dim0()<<" x "<< a.dim1(); - } - domain_type domain() const { return a.domain(); } - size_t dim0() const { return a.dim0();} - size_t dim1() const { return a.dim1();} - 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 { - inverse_lazy(A const & a):inverse_lazy_impl(a) { } - }; - - // for matrix and matrix_views, we have more optimisation possible .... - template - struct inverse_lazy >::type >:inverse_lazy_impl{ - inverse_lazy(A const & a):inverse_lazy_impl(a) { } - - template // Optimized implementation of = - friend void triqs_arrays_assign_delegation (MT & lhs, inverse_lazy const & rhs) { - static_assert(is_matrix_or_view::value, "Internal error"); - 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 - 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 ... - } - } - friend std::ostream & operator<<(std::ostream & out,inverse_lazy const&x){return out<<"inverse("< struct determinant_lazy { // : { Tag::expression_terminal, Tag::scalar_expression_terminal { - 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);} - }; - - }} // namespace triqs::arrays -#endif diff --git a/triqs/arrays/matrix_stack_view.hpp b/triqs/arrays/matrix_stack_view.hpp index d5c0242f..b8474c8e 100644 --- a/triqs/arrays/matrix_stack_view.hpp +++ b/triqs/arrays/matrix_stack_view.hpp @@ -23,7 +23,7 @@ #include "./array.hpp" #include "./matrix.hpp" #include "./matrix_view_proxy.hpp" -#include +#include namespace triqs { namespace arrays { diff --git a/triqs/det_manip/det_manip.hpp b/triqs/det_manip/det_manip.hpp index c0dea051..4929bbcf 100644 --- a/triqs/det_manip/det_manip.hpp +++ b/triqs/det_manip/det_manip.hpp @@ -26,8 +26,7 @@ #include //#include #include -#include -#include +#include #include #include #include diff --git a/triqs/lattice/bravais_lattice_and_brillouin_zone.cpp b/triqs/lattice/bravais_lattice_and_brillouin_zone.cpp index f93df772..dfc00636 100644 --- a/triqs/lattice/bravais_lattice_and_brillouin_zone.cpp +++ b/triqs/lattice/bravais_lattice_and_brillouin_zone.cpp @@ -20,7 +20,7 @@ ******************************************************************************/ #include "bravais_lattice_and_brillouin_zone.hpp" #include -#include +#include #include namespace triqs { namespace lattice_tools {