mirror of
https://github.com/triqs/dft_tools
synced 2024-11-01 11:43:47 +01:00
97 lines
4.4 KiB
C++
97 lines
4.4 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_CONCEPT_TAG_NAME(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( first_dim(P.a), second_dim(P.b)) { 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 (second_dim(a) != first_dim(b)) TRIQS_RUNTIME_ERROR<< "Matrix product : dimension mismatch in A*B "<< a<<" "<< b;
|
|
}
|
|
|
|
domain_type domain() const { return mini_vector<size_t,2>(first_dim(a), second_dim(b));}
|
|
|
|
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(first_dim(rhs),second_dim(rhs)));
|
|
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 (first_dim(lhs) != first_dim(*this))
|
|
TRIQS_RUNTIME_ERROR<< "Matmul : +=/-= operator : first dimension mismatch in A*B "<< first_dim(lhs)<<" vs "<< first_dim(*this);
|
|
if (second_dim(lhs) != second_dim(*this))
|
|
TRIQS_RUNTIME_ERROR<< "Matmul : +=/-= operator : first dimension mismatch in A*B "<< second_dim(lhs)<<" vs "<< second_dim(*this);
|
|
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
|