mirror of
https://github.com/triqs/dft_tools
synced 2024-12-26 06:14:14 +01:00
200 lines
6.6 KiB
C++
200 lines
6.6 KiB
C++
/*******************************************************************************
|
|
*
|
|
* 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 <http://www.gnu.org/licenses/>.
|
|
*
|
|
******************************************************************************/
|
|
#ifndef TRIQS_ARRAYS_LINALG_DET_INV_H
|
|
#define TRIQS_ARRAYS_LINALG_DET_INV_H
|
|
#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 <typename A> struct inverse_lazy;
|
|
|
|
template <class A> inverse_lazy<typename utility::remove_rvalue_ref<A>::type> inverse(A &&a) {
|
|
return {std::forward<A>(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...)
|
|
// A can be a matrix, a matrix_view
|
|
template <typename A> class det_and_inverse_worker {
|
|
typedef typename A::value_type value_type;
|
|
typedef matrix_view<value_type> V_type;
|
|
A a;
|
|
int dim;
|
|
triqs::arrays::vector<int> ipiv;
|
|
int step, info;
|
|
value_type _det;
|
|
|
|
public:
|
|
det_and_inverse_worker(A a_) : a(std::move(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 :(" << first_dim(a) << "," << second_dim(a)
|
|
<< ")\n ";
|
|
if (!(has_contiguous_data(a))) TRIQS_RUNTIME_ERROR << "det_and_inverse_worker only takes a contiguous view";
|
|
}
|
|
|
|
value_type det() {
|
|
V_type W = fortran_view(a);
|
|
_step1(W);
|
|
_compute_det(W);
|
|
return _det;
|
|
}
|
|
|
|
A const &inverse() {
|
|
if (step < 2) {
|
|
V_type W = fortran_view(a);
|
|
_step1(W);
|
|
_step2(W);
|
|
}
|
|
return a;
|
|
}
|
|
|
|
private:
|
|
|
|
// no need of special traversal
|
|
template <typename MT> 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;
|
|
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 < 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);
|
|
info = lapack::getri(W, ipiv);
|
|
if (info != 0) throw matrix_inverse_exception() << "Inverse/Det error : matrix is not invertible";
|
|
}
|
|
};
|
|
|
|
//-----------------------------------------------------------
|
|
|
|
template <typename A> struct inverse_lazy : TRIQS_CONCEPT_TAG_NAME(ImmutableMatrix) {
|
|
typedef typename std::remove_reference<A>::type A_t;
|
|
typedef typename std::remove_const<typename A_t::value_type>::type value_type;
|
|
typedef typename A_t::domain_type domain_type;
|
|
typedef matrix<value_type> M_type;
|
|
typedef matrix_view<value_type> M_view_type;
|
|
|
|
template <typename AA> inverse_lazy(AA &&a_) : a(std::forward<AA>(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 <typename K0, typename K1> 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_view_type> {M};
|
|
worker.inverse();
|
|
computed = true;
|
|
}
|
|
};
|
|
|
|
// Optimized implementation of =
|
|
// if M = inverse(M) with the SAME object, then we do not need to copy the data
|
|
template <typename A, typename MT>
|
|
ENABLE_IF(is_matrix_or_view<typename std::remove_reference<A>::type>)
|
|
triqs_arrays_assign_delegation(MT &lhs, inverse_lazy<A> const &rhs) {
|
|
static_assert(is_matrix_or_view<MT>::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<MT> C(lhs);
|
|
// a reflexive cache will use a temporary "regrouping" copy if and only if needed
|
|
det_and_inverse_worker<typename MT::view_type> W(C()); // the worker to make the inversion of the lhs...
|
|
W.inverse(); // worker is working ...
|
|
}
|
|
}
|
|
|
|
//------------------- det ----------------------------------------
|
|
|
|
template <typename A> typename std::remove_reference<A>::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<matrix<typename std::remove_reference<A>::type::value_type>>(std::forward<A>(a));
|
|
return worker.det();
|
|
}
|
|
}
|
|
|
|
namespace clef {
|
|
TRIQS_CLEF_MAKE_FNT_LAZY(determinant);
|
|
}
|
|
|
|
} // namespace triqs::arrays
|
|
#endif
|