yat/utility/expression/MatrixProduct.h

Code
Comments
Other
Rev Date Author Line
3908 13 May 20 peter 1 #ifndef _theplu_yat_utility_expression_matrix_product
3908 13 May 20 peter 2 #define _theplu_yat_utility_expression_matrix_product
3908 13 May 20 peter 3
3908 13 May 20 peter 4 // $Id$
3908 13 May 20 peter 5
3908 13 May 20 peter 6 /*
3908 13 May 20 peter 7   Copyright (C) 2020 Peter Johansson
3908 13 May 20 peter 8
3908 13 May 20 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
3908 13 May 20 peter 10
3908 13 May 20 peter 11   The yat library is free software; you can redistribute it and/or
3908 13 May 20 peter 12   modify it under the terms of the GNU General Public License as
3908 13 May 20 peter 13   published by the Free Software Foundation; either version 3 of the
3908 13 May 20 peter 14   License, or (at your option) any later version.
3908 13 May 20 peter 15
3908 13 May 20 peter 16   The yat library is distributed in the hope that it will be useful,
3908 13 May 20 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
3908 13 May 20 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3908 13 May 20 peter 19   General Public License for more details.
3908 13 May 20 peter 20
3908 13 May 20 peter 21   You should have received a copy of the GNU General Public License
3908 13 May 20 peter 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
3908 13 May 20 peter 23 */
3908 13 May 20 peter 24
3908 13 May 20 peter 25 #include "MatrixTraits.h"
3908 13 May 20 peter 26 #include "ScaledMatrix.h"
3908 13 May 20 peter 27 #include "TransposedMatrix.h"
3908 13 May 20 peter 28
3908 13 May 20 peter 29 #include "yat/utility/BasicMatrix.h"
3908 13 May 20 peter 30 #include "yat/utility/MatrixExpression.h"
3908 13 May 20 peter 31 #include "yat/utility/yat_assert.h"
3908 13 May 20 peter 32
3908 13 May 20 peter 33 #include <gsl/gsl_blas_types.h>
3908 13 May 20 peter 34 #include <gsl/gsl_blas.h>
3908 13 May 20 peter 35 #include <gsl/gsl_cblas.h>
3908 13 May 20 peter 36 #include <gsl/gsl_matrix.h>
3908 13 May 20 peter 37
3908 13 May 20 peter 38 namespace theplu {
3908 13 May 20 peter 39 namespace yat {
3908 13 May 20 peter 40 namespace utility {
3908 13 May 20 peter 41
3908 13 May 20 peter 42   /// \cond IGNORE_DOXYGEN
3908 13 May 20 peter 43
3908 13 May 20 peter 44 namespace expression {
3908 13 May 20 peter 45
3908 13 May 20 peter 46
3908 13 May 20 peter 47   template<class MATRIX1, class MATRIX2>
3908 13 May 20 peter 48   class MatrixProduct : public MatrixExpression<MatrixProduct<MATRIX1,MATRIX2> >
3908 13 May 20 peter 49   {
3908 13 May 20 peter 50   public:
3908 13 May 20 peter 51     MatrixProduct(const MATRIX1& A, const MATRIX2& B)
3908 13 May 20 peter 52       : A_(A), B_(B) {}
3908 13 May 20 peter 53
3908 13 May 20 peter 54     size_t rows(void) const { return A_.rows(); }
3908 13 May 20 peter 55     size_t columns(void) const { return B_.columns(); }
3908 13 May 20 peter 56
3908 13 May 20 peter 57     double operator()(size_t row, size_t column) const
3908 13 May 20 peter 58     { return gsl_matrix_get(this->gsl_matrix_p(), row, column); }
3908 13 May 20 peter 59
3908 13 May 20 peter 60     void calculate_matrix(gsl_matrix*& result) const
3908 13 May 20 peter 61     {
3908 13 May 20 peter 62       detail::reallocate(result, this->rows(), this->columns());
3908 13 May 20 peter 63       YAT_ASSERT(detail::rows(result) == this->rows());
3908 13 May 20 peter 64       YAT_ASSERT(detail::columns(result) == this->columns());
3908 13 May 20 peter 65       YAT_ASSERT(A_.columns() == B_.rows());
3908 13 May 20 peter 66       MatrixTraits<MATRIX1> traits1;
3908 13 May 20 peter 67       MatrixTraits<MATRIX2> traits2;
3908 13 May 20 peter 68       gsl_blas_dgemm(traits1.transpose_type(), traits2.transpose_type(),
3908 13 May 20 peter 69                      traits1.factor(A_) * traits2.factor(B_),
3908 13 May 20 peter 70                      traits1.get_gsl_matrix_p(A_),
3908 13 May 20 peter 71                      traits2.get_gsl_matrix_p(B_),
3908 13 May 20 peter 72                      0.0, result);
3908 13 May 20 peter 73       YAT_ASSERT(result);
3908 13 May 20 peter 74     }
3908 13 May 20 peter 75   private:
3908 13 May 20 peter 76     const MATRIX1& A_;
3908 13 May 20 peter 77     const MATRIX2& B_;
3908 13 May 20 peter 78   };
3908 13 May 20 peter 79 } // end namespace expression
3908 13 May 20 peter 80
3908 13 May 20 peter 81   /// \endcond
3908 13 May 20 peter 82
3908 13 May 20 peter 83 }}} // of namespace utility, yat, and theplu
3908 13 May 20 peter 84
3908 13 May 20 peter 85 #endif

Generated on Sat Oct 14 02:24:10 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)