mirror of
https://github.com/triqs/dft_tools
synced 2024-12-27 06:43:40 +01:00
fc2a620eae
- improve the mem_block and shared_block. - the reference counting is now done in the mem_block and shared_block, removing the need of shared_ptr. - speed tests shows that shared_ptr is very slow (due to thread safety?) the new version is much better, though not perfect. - Hence introducing weak views. - also : -- clean the guard mechanism for python (to allow returning from python without any python ref left). -- clean code, add documentation for mem_block -- remove nan init, which was not working, and corresponding test -- serialisation of view still unchanged (need to forbid serialization of view ??). - tests ok, incl. valgrind tests.
100 lines
4.5 KiB
C++
100 lines
4.5 KiB
C++
/*******************************************************************************
|
|
*
|
|
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
|
*
|
|
* Copyright (C) 2011-2012 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_EXPRESSION_MATMUL_H
|
|
#define TRIQS_ARRAYS_EXPRESSION_MATMUL_H
|
|
#include <boost/type_traits/is_same.hpp>
|
|
#include <boost/typeof/typeof.hpp>
|
|
#include "../blas_lapack/gemm.hpp"
|
|
namespace triqs { namespace arrays {
|
|
|
|
///
|
|
template<typename A, typename B> class matmul_lazy;
|
|
///
|
|
template<typename A, typename B> matmul_lazy<A,B> matmul (A const & a, B const & b) { return matmul_lazy<A,B>(a,b); }
|
|
|
|
// ----------------- implementation -----------------------------------------
|
|
|
|
template<typename A, typename B> class matmul_lazy : TRIQS_MODEL_CONCEPT(ImmutableMatrix) {
|
|
typedef typename boost::remove_const<typename A::value_type>::type V1;
|
|
typedef typename boost::remove_const<typename B::value_type>::type V2;
|
|
//static_assert((boost::is_same<V1,V2>::value),"Different values : not implemented");
|
|
|
|
public:
|
|
typedef BOOST_TYPEOF_TPL( V1() * V2()) value_type; // what is the result of multiplying a V1 by a V2 ?
|
|
typedef typename A::domain_type domain_type;
|
|
typedef typename const_view_type_if_exists_else_type<A>::type A_type;
|
|
typedef typename const_view_type_if_exists_else_type<B>::type B_type;
|
|
const A_type a; const B_type b;
|
|
|
|
private:
|
|
typedef matrix<value_type> matrix_type;
|
|
|
|
struct internal_data { // implementing the pattern LazyPreCompute
|
|
matrix_type R;
|
|
internal_data(matmul_lazy const & P): R( P.a.dim0(), P.b.dim1()) { blas::gemm(1.0,P.a, P.b, 0.0, R); }
|
|
};
|
|
friend struct internal_data;
|
|
mutable std::shared_ptr<internal_data> _id;
|
|
void activate() const { if (!_id) _id= std::make_shared<internal_data>(*this);}
|
|
|
|
public:
|
|
matmul_lazy( A const & a_, B const & b_):a(a_),b(b_){
|
|
if (a.dim1() != b.dim0()) TRIQS_RUNTIME_ERROR<< "Matrix product : dimension mismatch in A*B "<< a<<" "<< b;
|
|
}
|
|
|
|
domain_type domain() const { return mini_vector<size_t,2>(a.dim0(), b.dim1());}
|
|
//domain_type domain() const { return indexmaps::cuboid::domain_t<2>(mini_vector<size_t,2>(a.dim0(), b.dim1()));}
|
|
size_t dim0() const { return a.dim0();}
|
|
size_t dim1() const { return b.dim1();}
|
|
|
|
template<typename K0, typename K1> value_type operator() (K0 const & k0, K1 const & k1) const { activate(); return _id->R(k0,k1); }
|
|
|
|
// TO BE REMOVED because of the aliasing question
|
|
// Optimized implementation of =, +=, -=
|
|
template<typename LHS>
|
|
friend void triqs_arrays_assign_delegation (LHS & lhs, matmul_lazy const & rhs) {
|
|
static_assert((is_matrix_or_view<LHS>::value), "LHS is not a matrix");
|
|
resize_or_check_if_view(lhs,make_shape(rhs.dim0(),rhs.dim1()));
|
|
blas::gemm(1.0,rhs.a, rhs.b, 0.0, lhs);
|
|
}
|
|
|
|
template<typename LHS>
|
|
friend void triqs_arrays_compound_assign_delegation (LHS & lhs, matmul_lazy const & rhs, mpl::char_<'A'>) { rhs.assign_comp_impl(lhs,1.0);}
|
|
template<typename LHS>
|
|
friend void triqs_arrays_compound_assign_delegation (LHS & lhs, matmul_lazy const & rhs, mpl::char_<'S'>) { rhs.assign_comp_impl(lhs,-1.0);}
|
|
|
|
private:
|
|
template<typename LHS> void assign_comp_impl (LHS & lhs, double S) const {
|
|
static_assert((is_matrix_or_view<LHS>::value), "LHS is not a matrix");
|
|
if (lhs.dim0() != dim0())
|
|
TRIQS_RUNTIME_ERROR<< "Matmul : +=/-= operator : first dimension mismatch in A*B "<< lhs.dim0()<<" vs "<< dim0();
|
|
if (lhs.dim1() != dim1())
|
|
TRIQS_RUNTIME_ERROR<< "Matmul : +=/-= operator : first dimension mismatch in A*B "<< lhs.dim1()<<" vs "<< dim1();
|
|
blas::gemm(S,a, b, 1.0, lhs);
|
|
}
|
|
|
|
friend std::ostream & operator<<(std::ostream & out, matmul_lazy<A,B> const & x){return out<<x.a<<" * "<<x.b;}
|
|
|
|
};// class matmul_lazy
|
|
|
|
}}//namespace triqs::arrays
|
|
#endif
|