yat/utility/BLAS_level2.h

Code
Comments
Other
Rev Date Author Line
3605 27 Jan 17 peter 1 #ifndef _theplu_yat_utility_blas_level2
3605 27 Jan 17 peter 2 #define _theplu_yat_utility_blas_level2
3605 27 Jan 17 peter 3
3605 27 Jan 17 peter 4 // $Id$
3605 27 Jan 17 peter 5
3605 27 Jan 17 peter 6 /*
4252 18 Nov 22 peter 7   Copyright (C) 2017, 2019, 2020, 2022 Peter Johansson
3605 27 Jan 17 peter 8
3605 27 Jan 17 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
3605 27 Jan 17 peter 10
3605 27 Jan 17 peter 11   The yat library is free software; you can redistribute it and/or
3605 27 Jan 17 peter 12   modify it under the terms of the GNU General Public License as
3605 27 Jan 17 peter 13   published by the Free Software Foundation; either version 3 of the
3605 27 Jan 17 peter 14   License, or (at your option) any later version.
3605 27 Jan 17 peter 15
3605 27 Jan 17 peter 16   The yat library is distributed in the hope that it will be useful,
3605 27 Jan 17 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
3605 27 Jan 17 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3605 27 Jan 17 peter 19   General Public License for more details.
3605 27 Jan 17 peter 20
3605 27 Jan 17 peter 21   You should have received a copy of the GNU General Public License
3605 27 Jan 17 peter 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
3605 27 Jan 17 peter 23 */
3605 27 Jan 17 peter 24
3605 27 Jan 17 peter 25 #include "BasicMatrix.h"
3605 27 Jan 17 peter 26 #include "BasicVector.h"
3908 13 May 20 peter 27 #include "BLAS_level3.h"
3908 13 May 20 peter 28 #include "VectorExpression.h"
3605 27 Jan 17 peter 29 #include "yat_assert.h"
3605 27 Jan 17 peter 30
3908 13 May 20 peter 31 #include "expression/MatrixTraits.h"
3605 27 Jan 17 peter 32
3605 27 Jan 17 peter 33 #include "gsl/gsl_blas.h"
3605 27 Jan 17 peter 34
3605 27 Jan 17 peter 35 #include <cstddef>
3605 27 Jan 17 peter 36
3605 27 Jan 17 peter 37 namespace theplu {
3605 27 Jan 17 peter 38 namespace yat {
3605 27 Jan 17 peter 39 namespace utility {
3605 27 Jan 17 peter 40
3605 27 Jan 17 peter 41   // This file defines operations using both Vector and Matrix
3605 27 Jan 17 peter 42
3605 27 Jan 17 peter 43
3605 27 Jan 17 peter 44   /// \cond IGNORE_DOXYGEN
3605 27 Jan 17 peter 45
3605 27 Jan 17 peter 46   namespace expression {
3605 27 Jan 17 peter 47
3605 27 Jan 17 peter 48     template<typename MATRIX, typename VECTOR>
3605 27 Jan 17 peter 49     class MatrixVector : public VectorExpression<MatrixVector<MATRIX, VECTOR> >
3605 27 Jan 17 peter 50     {
3605 27 Jan 17 peter 51     public:
3908 13 May 20 peter 52
3605 27 Jan 17 peter 53       MatrixVector(const BasicMatrix<MATRIX>& lhs,
3605 27 Jan 17 peter 54                    const BasicVector<VECTOR>& rhs)
3605 27 Jan 17 peter 55       {
3605 27 Jan 17 peter 56         YAT_ASSERT(lhs.columns() == rhs.size());
3605 27 Jan 17 peter 57         YAT_ASSERT(lhs.rows());
3605 27 Jan 17 peter 58         this->allocate_memory(lhs.rows());
3908 13 May 20 peter 59
3908 13 May 20 peter 60         MatrixTraits<BasicMatrix<MATRIX> > traits;
3908 13 May 20 peter 61
3908 13 May 20 peter 62         gsl_blas_dgemv(traits.transpose_type(), traits.factor(lhs),
3908 13 May 20 peter 63                        traits.get_gsl_matrix_p(lhs),
3605 27 Jan 17 peter 64                        rhs.gsl_vector_p(), 0.0, this->v_);
3605 27 Jan 17 peter 65       }
3605 27 Jan 17 peter 66
3605 27 Jan 17 peter 67
3605 27 Jan 17 peter 68       // access an element
3605 27 Jan 17 peter 69       double operator()(size_t i) const
3605 27 Jan 17 peter 70       {
3605 27 Jan 17 peter 71         YAT_ASSERT(this->v_);
3605 27 Jan 17 peter 72         return gsl_vector_get(this->v_, i);
3605 27 Jan 17 peter 73       }
3605 27 Jan 17 peter 74
3605 27 Jan 17 peter 75
3605 27 Jan 17 peter 76       size_t size(void) const
3605 27 Jan 17 peter 77       {
3605 27 Jan 17 peter 78         YAT_ASSERT(this->v_);
3605 27 Jan 17 peter 79         return this->v_->size;
3605 27 Jan 17 peter 80       }
3605 27 Jan 17 peter 81
3605 27 Jan 17 peter 82
3605 27 Jan 17 peter 83       void calculate_gsl_vector_p(void) const
3605 27 Jan 17 peter 84       {
3908 13 May 20 peter 85         // This should never be called since v_ is constructed in
3908 13 May 20 peter 86         // constructor
3605 27 Jan 17 peter 87         YAT_ASSERT(this->v_);
3605 27 Jan 17 peter 88         YAT_ASSERT(0);
3605 27 Jan 17 peter 89       }
3605 27 Jan 17 peter 90     };
3605 27 Jan 17 peter 91
3605 27 Jan 17 peter 92   } // end namespace expression
3605 27 Jan 17 peter 93
3605 27 Jan 17 peter 94   /// \endcond
3605 27 Jan 17 peter 95
3605 27 Jan 17 peter 96
3605 27 Jan 17 peter 97   /**
4252 18 Nov 22 peter 98      \relates MatrixBase
3605 27 Jan 17 peter 99      \relates VectorBase
3605 27 Jan 17 peter 100
3605 27 Jan 17 peter 101      Calculate product of Matrix expression and Vector expression
3605 27 Jan 17 peter 102
3605 27 Jan 17 peter 103      \since new in yat 0.15
3605 27 Jan 17 peter 104    */
3605 27 Jan 17 peter 105   template<class MATRIX, class VECTOR>
3605 27 Jan 17 peter 106   expression::MatrixVector<MATRIX, VECTOR>
3605 27 Jan 17 peter 107   operator*(const BasicMatrix<MATRIX>& lhs, const BasicVector<VECTOR>& rhs)
3605 27 Jan 17 peter 108   {
3605 27 Jan 17 peter 109     YAT_ASSERT(lhs.columns() == rhs.size());
3605 27 Jan 17 peter 110     return expression::MatrixVector<MATRIX, VECTOR>(lhs, rhs);
3605 27 Jan 17 peter 111   }
3605 27 Jan 17 peter 112
3605 27 Jan 17 peter 113
3605 27 Jan 17 peter 114   /**
4252 18 Nov 22 peter 115      \relates MatrixBase
3605 27 Jan 17 peter 116      \relates VectorBase
3605 27 Jan 17 peter 117
3605 27 Jan 17 peter 118      Calculate product of Vector expression and Matrix expression
3605 27 Jan 17 peter 119
3605 27 Jan 17 peter 120      \since new in yat 0.15
3605 27 Jan 17 peter 121    */
3605 27 Jan 17 peter 122   template<class MATRIX, class VECTOR>
3908 13 May 20 peter 123   expression::MatrixVector<MatrixExpression<
3908 13 May 20 peter 124                              expression::TransposedMatrix<MATRIX> >,
3908 13 May 20 peter 125                            VECTOR>
3605 27 Jan 17 peter 126   operator*(const BasicVector<VECTOR>& lhs, const BasicMatrix<MATRIX>& rhs)
3605 27 Jan 17 peter 127   {
3605 27 Jan 17 peter 128     YAT_ASSERT(lhs.size() == rhs.rows());
3908 13 May 20 peter 129     return transpose(rhs) * lhs;
3605 27 Jan 17 peter 130   }
3605 27 Jan 17 peter 131
3605 27 Jan 17 peter 132
3605 27 Jan 17 peter 133 }}} // of namespace utility, yat, and theplu
3605 27 Jan 17 peter 134
3605 27 Jan 17 peter 135 #endif

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