yat/utility/DiagonalMatrix.h

Code
Comments
Other
Rev Date Author Line
3655 13 Jul 17 peter 1 #ifndef theplu_yat_utility_diagonal_matrix
3655 13 Jul 17 peter 2 #define theplu_yat_utility_diagonal_matrix
3655 13 Jul 17 peter 3
3655 13 Jul 17 peter 4 // $Id$
3655 13 Jul 17 peter 5
3655 13 Jul 17 peter 6 /*
4207 26 Aug 22 peter 7   Copyright (C) 2017, 2021, 2022 Peter Johansson
3655 13 Jul 17 peter 8
3655 13 Jul 17 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
3655 13 Jul 17 peter 10
3655 13 Jul 17 peter 11   The yat library is free software; you can redistribute it and/or
3655 13 Jul 17 peter 12   modify it under the terms of the GNU General Public License as
3655 13 Jul 17 peter 13   published by the Free Software Foundation; either version 3 of the
3655 13 Jul 17 peter 14   License, or (at your option) any later version.
3655 13 Jul 17 peter 15
3655 13 Jul 17 peter 16   The yat library is distributed in the hope that it will be useful,
3655 13 Jul 17 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
3655 13 Jul 17 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3655 13 Jul 17 peter 19   General Public License for more details.
3655 13 Jul 17 peter 20
3655 13 Jul 17 peter 21   You should have received a copy of the GNU General Public License
3655 13 Jul 17 peter 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
3655 13 Jul 17 peter 23 */
3655 13 Jul 17 peter 24
3655 13 Jul 17 peter 25 #include "BasicMatrix.h"
3655 13 Jul 17 peter 26 #include "BasicVector.h"
3655 13 Jul 17 peter 27 #include "BLAS_utility.h"
3655 13 Jul 17 peter 28
3655 13 Jul 17 peter 29 #include "MatrixExpression.h"
3655 13 Jul 17 peter 30 #include "VectorExpression.h"
3655 13 Jul 17 peter 31
3655 13 Jul 17 peter 32 #include "Vector.h"
3655 13 Jul 17 peter 33 #include "yat_assert.h"
3655 13 Jul 17 peter 34
3655 13 Jul 17 peter 35 namespace theplu {
3655 13 Jul 17 peter 36 namespace yat {
3655 13 Jul 17 peter 37 namespace utility {
3655 13 Jul 17 peter 38
4125 14 Jan 22 peter 39   class MatrixBase;
3655 13 Jul 17 peter 40
3655 13 Jul 17 peter 41   /**
3655 13 Jul 17 peter 42      \brief Diagonal Matrix
3655 13 Jul 17 peter 43
3655 13 Jul 17 peter 44      All off-diagonal elements are zero by definition, which can be
3655 13 Jul 17 peter 45      used to speed up some calculation.
3655 13 Jul 17 peter 46
3655 13 Jul 17 peter 47      \since New in yat 0.15
3655 13 Jul 17 peter 48    */
3655 13 Jul 17 peter 49   class DiagonalMatrix
3655 13 Jul 17 peter 50   {
3655 13 Jul 17 peter 51   public:
3655 13 Jul 17 peter 52     /**
3655 13 Jul 17 peter 53        \brief Default constructor
3655 13 Jul 17 peter 54     */
3655 13 Jul 17 peter 55     DiagonalMatrix(void);
3655 13 Jul 17 peter 56
3655 13 Jul 17 peter 57     /**
3655 13 Jul 17 peter 58        Construct a \a row x \a column Matrix with all elements set to x
3655 13 Jul 17 peter 59      */
3655 13 Jul 17 peter 60     DiagonalMatrix(size_t row, size_t column, double x=0.0);
3655 13 Jul 17 peter 61
3655 13 Jul 17 peter 62     /**
3655 13 Jul 17 peter 63        Construct a square matrix in which diagonal is set to \a diagonal
3655 13 Jul 17 peter 64      */
3655 13 Jul 17 peter 65     explicit DiagonalMatrix(const VectorBase& diagonal);
3655 13 Jul 17 peter 66
3655 13 Jul 17 peter 67     /**
3655 13 Jul 17 peter 68        Construct a diagonal matrix from \a m. All non-diagonal
3655 13 Jul 17 peter 69        elements are set to zero.
3655 13 Jul 17 peter 70      */
4125 14 Jan 22 peter 71     explicit DiagonalMatrix(const MatrixBase& m);
3655 13 Jul 17 peter 72
3655 13 Jul 17 peter 73     /**
3655 13 Jul 17 peter 74        \return number of rows
3655 13 Jul 17 peter 75      */
3655 13 Jul 17 peter 76     size_t rows(void) const;
3655 13 Jul 17 peter 77
3655 13 Jul 17 peter 78     /**
3655 13 Jul 17 peter 79        \return number of columns
3655 13 Jul 17 peter 80      */
3655 13 Jul 17 peter 81     size_t columns(void) const;
3655 13 Jul 17 peter 82
3655 13 Jul 17 peter 83     /**
3655 13 Jul 17 peter 84        access element at row \a row and column \a column
3655 13 Jul 17 peter 85      */
3655 13 Jul 17 peter 86     const double operator()(size_t row, size_t col) const;
3655 13 Jul 17 peter 87
3655 13 Jul 17 peter 88     /**
3655 13 Jul 17 peter 89        access element at row \a i and column \a i
3655 13 Jul 17 peter 90
3655 13 Jul 17 peter 91        Elements outside diagnonal are only accessible as const
3655 13 Jul 17 peter 92      */
3655 13 Jul 17 peter 93     double& operator()(size_t i);
4105 24 Sep 21 peter 94
4105 24 Sep 21 peter 95     /**
4105 24 Sep 21 peter 96        \brief Multiplication and assign operator
4105 24 Sep 21 peter 97
4105 24 Sep 21 peter 98        Same as doing *this = *this * rhs.
4105 24 Sep 21 peter 99
4105 24 Sep 21 peter 100        \since New in yat 0.20
4105 24 Sep 21 peter 101      */
4105 24 Sep 21 peter 102     DiagonalMatrix& operator*=(const DiagonalMatrix& rhs);
4105 24 Sep 21 peter 103
4105 24 Sep 21 peter 104     /**
4105 24 Sep 21 peter 105        \brief Add and assign operator
4105 24 Sep 21 peter 106
4105 24 Sep 21 peter 107        Elementwise addition of the diagnonal elements.
4105 24 Sep 21 peter 108
4105 24 Sep 21 peter 109        \since New in yat 0.20
4105 24 Sep 21 peter 110      */
4105 24 Sep 21 peter 111     DiagonalMatrix& operator+=(const DiagonalMatrix& rhs);
4105 24 Sep 21 peter 112
4105 24 Sep 21 peter 113     /**
4105 24 Sep 21 peter 114        \brief Subtract and assign operator
4105 24 Sep 21 peter 115
4105 24 Sep 21 peter 116        Elementwise addition of the diagnonal elements.
4105 24 Sep 21 peter 117
4105 24 Sep 21 peter 118        \since New in yat 0.20
4105 24 Sep 21 peter 119      */
4105 24 Sep 21 peter 120     DiagonalMatrix& operator-=(const DiagonalMatrix& rhs);
3655 13 Jul 17 peter 121   private:
3655 13 Jul 17 peter 122     Vector data_;
3655 13 Jul 17 peter 123     size_t row_;
3655 13 Jul 17 peter 124     size_t col_;
3655 13 Jul 17 peter 125   };
3655 13 Jul 17 peter 126
4128 18 Jan 22 peter 127   /// \cond IGNORE_DOXYGEN
3655 13 Jul 17 peter 128
3655 13 Jul 17 peter 129   // expression classes
3655 13 Jul 17 peter 130   namespace expression
3655 13 Jul 17 peter 131   {
3655 13 Jul 17 peter 132     template<typename T>
3655 13 Jul 17 peter 133     class DiagonalMatrixMatrix :
3655 13 Jul 17 peter 134       public MatrixExpression<DiagonalMatrixMatrix<T> >
3655 13 Jul 17 peter 135     {
3655 13 Jul 17 peter 136     public:
3655 13 Jul 17 peter 137       DiagonalMatrixMatrix(const DiagonalMatrix& lhs,const BasicMatrix<T>& rhs)
3655 13 Jul 17 peter 138         : lhs_(lhs), rhs_(rhs) {}
3655 13 Jul 17 peter 139       size_t rows(void) const { return lhs_.rows(); }
3655 13 Jul 17 peter 140       size_t columns(void) const { return rhs_.columns(); }
3655 13 Jul 17 peter 141       double operator()(size_t i, size_t j) const
3655 13 Jul 17 peter 142       { return lhs_(i, i) * rhs_(i, j); }
3655 13 Jul 17 peter 143
3655 13 Jul 17 peter 144       void calculate_matrix(gsl_matrix*& result) const
3655 13 Jul 17 peter 145       {
3655 13 Jul 17 peter 146         detail::reallocate(result, rows(), columns());
3655 13 Jul 17 peter 147         for (size_t i=0; i<rows(); ++i)
3655 13 Jul 17 peter 148           for (size_t j=0; j<columns(); ++j)
3655 13 Jul 17 peter 149             gsl_matrix_set(result, i, j, (*this)(i, j));
3655 13 Jul 17 peter 150       }
3655 13 Jul 17 peter 151
3655 13 Jul 17 peter 152     private:
3655 13 Jul 17 peter 153       const DiagonalMatrix& lhs_;
3655 13 Jul 17 peter 154       const BasicMatrix<T>& rhs_;
3655 13 Jul 17 peter 155     };
3655 13 Jul 17 peter 156
3655 13 Jul 17 peter 157
3655 13 Jul 17 peter 158     template<class T>
3655 13 Jul 17 peter 159     class DiagonalMatrixVector
3655 13 Jul 17 peter 160       : public VectorExpression<DiagonalMatrixVector<T> >
3655 13 Jul 17 peter 161     {
3655 13 Jul 17 peter 162     public:
3655 13 Jul 17 peter 163       DiagonalMatrixVector(const DiagonalMatrix& lhs, const BasicVector<T>& rhs)
3655 13 Jul 17 peter 164         : m_(lhs), vec_(rhs), size_(lhs.rows())
3655 13 Jul 17 peter 165       {
3655 13 Jul 17 peter 166       }
3655 13 Jul 17 peter 167
3655 13 Jul 17 peter 168       DiagonalMatrixVector(const BasicVector<T>& lhs, const DiagonalMatrix& rhs)
3655 13 Jul 17 peter 169         : m_(rhs), vec_(lhs), size_(rhs.columns())
3655 13 Jul 17 peter 170       {
3655 13 Jul 17 peter 171       }
3655 13 Jul 17 peter 172
3655 13 Jul 17 peter 173       double operator()(size_t i) const
3655 13 Jul 17 peter 174       {
3655 13 Jul 17 peter 175         YAT_ASSERT(i < size());
3655 13 Jul 17 peter 176         return m_(i, i) * vec_(i);
3655 13 Jul 17 peter 177       }
3655 13 Jul 17 peter 178
3655 13 Jul 17 peter 179
3655 13 Jul 17 peter 180       size_t size(void) const
3655 13 Jul 17 peter 181       {
3655 13 Jul 17 peter 182         return size_;
3655 13 Jul 17 peter 183       }
3655 13 Jul 17 peter 184
3655 13 Jul 17 peter 185
3655 13 Jul 17 peter 186       void calculate_gsl_vector_p(void) const
3655 13 Jul 17 peter 187       {
3655 13 Jul 17 peter 188         this->allocate_memory(size_);
3655 13 Jul 17 peter 189         for (size_t i=0; i<size_; ++i)
3655 13 Jul 17 peter 190           gsl_vector_set(this->v_, i, (*this)(i));
3655 13 Jul 17 peter 191       }
3655 13 Jul 17 peter 192     private:
3655 13 Jul 17 peter 193       const DiagonalMatrix& m_;
3655 13 Jul 17 peter 194       const BasicVector<T>& vec_;
3655 13 Jul 17 peter 195       size_t size_;
3655 13 Jul 17 peter 196     };
3655 13 Jul 17 peter 197
3655 13 Jul 17 peter 198
3655 13 Jul 17 peter 199
3655 13 Jul 17 peter 200     template<typename T>
3655 13 Jul 17 peter 201     class MatrixDiagonalMatrix :
3655 13 Jul 17 peter 202       public MatrixExpression<MatrixDiagonalMatrix<T> >
3655 13 Jul 17 peter 203     {
3655 13 Jul 17 peter 204     public:
3655 13 Jul 17 peter 205       MatrixDiagonalMatrix(const BasicMatrix<T>& lhs, const DiagonalMatrix& rhs)
3655 13 Jul 17 peter 206         : lhs_(lhs), rhs_(rhs) {}
3655 13 Jul 17 peter 207       size_t rows(void) const { return lhs_.rows(); }
3655 13 Jul 17 peter 208       size_t columns(void) const { return rhs_.columns(); }
3655 13 Jul 17 peter 209       double operator()(size_t i, size_t j) const
3655 13 Jul 17 peter 210       { return lhs_(i, j) * rhs_(j, j); }
3655 13 Jul 17 peter 211
3655 13 Jul 17 peter 212       void calculate_matrix(gsl_matrix*& result) const
3655 13 Jul 17 peter 213       {
3655 13 Jul 17 peter 214         detail::reallocate(result, rows(), columns());
3655 13 Jul 17 peter 215         for (size_t i=0; i<rows(); ++i)
3655 13 Jul 17 peter 216           for (size_t j=0; j<columns(); ++j)
3655 13 Jul 17 peter 217             gsl_matrix_set(result, i, j, (*this)(i, j));
3655 13 Jul 17 peter 218       }
3655 13 Jul 17 peter 219
3655 13 Jul 17 peter 220     private:
3655 13 Jul 17 peter 221       const BasicMatrix<T>& lhs_;
3655 13 Jul 17 peter 222       const DiagonalMatrix& rhs_;
3655 13 Jul 17 peter 223     };
3655 13 Jul 17 peter 224
3655 13 Jul 17 peter 225
3655 13 Jul 17 peter 226     template<typename T>
3655 13 Jul 17 peter 227     class VectorDiagonalMatrix :
3655 13 Jul 17 peter 228       public BasicVector<VectorDiagonalMatrix<T> >
3655 13 Jul 17 peter 229     {
3655 13 Jul 17 peter 230     public:
3655 13 Jul 17 peter 231       VectorDiagonalMatrix(const BasicVector<T>& rhs, const DiagonalMatrix& lhs)
3655 13 Jul 17 peter 232         : lhs_(lhs), rhs_(rhs) {}
3655 13 Jul 17 peter 233
3655 13 Jul 17 peter 234       double operator()(size_t i) const
3655 13 Jul 17 peter 235       {
3655 13 Jul 17 peter 236         YAT_ASSERT(i < size());
3655 13 Jul 17 peter 237         return lhs_(i) * rhs_(i, i);
3655 13 Jul 17 peter 238       }
3655 13 Jul 17 peter 239
3655 13 Jul 17 peter 240       size_t size(void) const
3655 13 Jul 17 peter 241       {
3655 13 Jul 17 peter 242         return rhs_.columns();
3655 13 Jul 17 peter 243       }
3655 13 Jul 17 peter 244
3655 13 Jul 17 peter 245
3655 13 Jul 17 peter 246       const gsl_vector* gsl_vector_p(void) const
3655 13 Jul 17 peter 247       {
3655 13 Jul 17 peter 248         size_t n = size();
3655 13 Jul 17 peter 249         gsl_vector* vec = detail::create_gsl_vector(n);
3655 13 Jul 17 peter 250         for (size_t i=0; i<n; ++i)
3655 13 Jul 17 peter 251           gsl_vector_set(vec, i, (*this)(i));
3655 13 Jul 17 peter 252         return vec;
3655 13 Jul 17 peter 253       }
3655 13 Jul 17 peter 254
3655 13 Jul 17 peter 255     private:
3655 13 Jul 17 peter 256       const BasicVector<T>& lhs_;
3655 13 Jul 17 peter 257       const DiagonalMatrix& rhs_;
3655 13 Jul 17 peter 258     };
3655 13 Jul 17 peter 259
3655 13 Jul 17 peter 260
3655 13 Jul 17 peter 261     template<class M1, class M2, class OP>
3655 13 Jul 17 peter 262     class DiagonalMatrix
3655 13 Jul 17 peter 263       : public MatrixExpression<DiagonalMatrix<M1, M2, OP> >
3655 13 Jul 17 peter 264     {
3655 13 Jul 17 peter 265     public:
3655 13 Jul 17 peter 266       DiagonalMatrix(const M1& lhs, const M2& rhs)
3655 13 Jul 17 peter 267         : lhs_(lhs), rhs_(rhs)
3655 13 Jul 17 peter 268       {
3655 13 Jul 17 peter 269         YAT_ASSERT(lhs.rows() == rhs.rows());
3655 13 Jul 17 peter 270         YAT_ASSERT(lhs.columns() == rhs.columns());
3655 13 Jul 17 peter 271       }
3655 13 Jul 17 peter 272       double operator()(size_t i, size_t j) const
3655 13 Jul 17 peter 273       {
3655 13 Jul 17 peter 274         return get(i, j, OP());
3655 13 Jul 17 peter 275       }
3655 13 Jul 17 peter 276       size_t rows(void) const { return lhs_.rows(); }
3655 13 Jul 17 peter 277       size_t columns(void) const { return rhs_.columns(); }
3655 13 Jul 17 peter 278       void calculate_matrix(gsl_matrix*& m) const
3655 13 Jul 17 peter 279       {
3655 13 Jul 17 peter 280         detail::reallocate(m, rows(), columns());
3655 13 Jul 17 peter 281         for (size_t i=0; i<rows(); ++i)
3655 13 Jul 17 peter 282           for (size_t j=0; j<columns(); ++j)
3655 13 Jul 17 peter 283             gsl_matrix_set(m, i, j, (*this)(i,j));
3655 13 Jul 17 peter 284       }
3655 13 Jul 17 peter 285
3655 13 Jul 17 peter 286     private:
3655 13 Jul 17 peter 287       const M1& lhs_;
3655 13 Jul 17 peter 288       const M2& rhs_;
3655 13 Jul 17 peter 289
3655 13 Jul 17 peter 290       double get(size_t i, size_t j, Plus) const
3655 13 Jul 17 peter 291       { return lhs_(i, j) + rhs_(i, j); }
3655 13 Jul 17 peter 292
3655 13 Jul 17 peter 293       double get(size_t i, size_t j, Minus) const
3655 13 Jul 17 peter 294       { return lhs_(i, j) - rhs_(i, j); }
3655 13 Jul 17 peter 295     };
3655 13 Jul 17 peter 296
3655 13 Jul 17 peter 297   } // end of namespace expression
3655 13 Jul 17 peter 298
3655 13 Jul 17 peter 299   /// \endcond
3655 13 Jul 17 peter 300
3655 13 Jul 17 peter 301
3655 13 Jul 17 peter 302   /**
3655 13 Jul 17 peter 303      Multiplication of two diagonal matrices
3655 13 Jul 17 peter 304
3655 13 Jul 17 peter 305      Diagonal elements are calculated as lhs(i,i) * rhs(i,i)
3655 13 Jul 17 peter 306
3655 13 Jul 17 peter 307      \since new in yat 0.15
3655 13 Jul 17 peter 308
3655 13 Jul 17 peter 309      \relates DiagonalMatrix
3655 13 Jul 17 peter 310    */
3655 13 Jul 17 peter 311   DiagonalMatrix operator*(const DiagonalMatrix& lhs,
3655 13 Jul 17 peter 312                            const DiagonalMatrix& rhs);
3655 13 Jul 17 peter 313
3655 13 Jul 17 peter 314   /**
3655 13 Jul 17 peter 315      \brief matrix matrix multiplication
3655 13 Jul 17 peter 316
3655 13 Jul 17 peter 317      \since new in yat 0.15
3655 13 Jul 17 peter 318
3655 13 Jul 17 peter 319      \relates DiagonalMatrix
3655 13 Jul 17 peter 320    */
3655 13 Jul 17 peter 321   template<class T>
3655 13 Jul 17 peter 322   expression::DiagonalMatrixMatrix<T>
3655 13 Jul 17 peter 323   operator*(const DiagonalMatrix& lhs, const BasicMatrix<T>& rhs)
3655 13 Jul 17 peter 324   {
3655 13 Jul 17 peter 325     YAT_ASSERT(lhs.columns() == rhs.rows());
3655 13 Jul 17 peter 326     return expression::DiagonalMatrixMatrix<T>(lhs, rhs);
3655 13 Jul 17 peter 327   }
3655 13 Jul 17 peter 328
3655 13 Jul 17 peter 329
3655 13 Jul 17 peter 330   /**
3655 13 Jul 17 peter 331      \brief matrix matrix multiplication
3655 13 Jul 17 peter 332
3655 13 Jul 17 peter 333      \since new in yat 0.15
3655 13 Jul 17 peter 334
3655 13 Jul 17 peter 335      \relates DiagonalMatrix
3655 13 Jul 17 peter 336    */
3655 13 Jul 17 peter 337   template<class T>
3655 13 Jul 17 peter 338   expression::MatrixDiagonalMatrix<T>
3655 13 Jul 17 peter 339   operator*(const BasicMatrix<T>& lhs, const DiagonalMatrix& rhs)
3655 13 Jul 17 peter 340   {
3655 13 Jul 17 peter 341     YAT_ASSERT(lhs.columns() == rhs.rows());
3655 13 Jul 17 peter 342     return expression::MatrixDiagonalMatrix<T>(lhs, rhs);
3655 13 Jul 17 peter 343   }
3655 13 Jul 17 peter 344
3655 13 Jul 17 peter 345
3655 13 Jul 17 peter 346   /**
3655 13 Jul 17 peter 347      \brief matrix matrix addition
3655 13 Jul 17 peter 348
3655 13 Jul 17 peter 349      \since new in yat 0.15
3655 13 Jul 17 peter 350
3655 13 Jul 17 peter 351      \relates DiagonalMatrix
3655 13 Jul 17 peter 352    */
3655 13 Jul 17 peter 353   DiagonalMatrix
3655 13 Jul 17 peter 354   operator+(const DiagonalMatrix& lhs, const DiagonalMatrix& rhs);
3655 13 Jul 17 peter 355
3655 13 Jul 17 peter 356
3655 13 Jul 17 peter 357   /**
3655 13 Jul 17 peter 358      \brief matrix matrix addition
3655 13 Jul 17 peter 359
3655 13 Jul 17 peter 360      \since new in yat 0.15
3655 13 Jul 17 peter 361
3655 13 Jul 17 peter 362      \relates DiagonalMatrix
3655 13 Jul 17 peter 363    */
3655 13 Jul 17 peter 364   template<class T>
3655 13 Jul 17 peter 365   expression::DiagonalMatrix<DiagonalMatrix, BasicMatrix<T>, expression::Plus>
3655 13 Jul 17 peter 366   operator+(const DiagonalMatrix& lhs, const BasicMatrix<T>& rhs)
3655 13 Jul 17 peter 367   {
3655 13 Jul 17 peter 368     YAT_ASSERT(lhs.rows() == rhs.rows());
3655 13 Jul 17 peter 369     YAT_ASSERT(lhs.columns() == rhs.columns());
3655 13 Jul 17 peter 370     return expression::DiagonalMatrix<DiagonalMatrix, BasicMatrix<T>,
3655 13 Jul 17 peter 371                                       expression::Plus>(lhs, rhs);
3655 13 Jul 17 peter 372   }
3655 13 Jul 17 peter 373
3655 13 Jul 17 peter 374
3655 13 Jul 17 peter 375   /**
3655 13 Jul 17 peter 376      \brief matrix matrix addition
3655 13 Jul 17 peter 377
3655 13 Jul 17 peter 378      \since new in yat 0.15
3655 13 Jul 17 peter 379
3655 13 Jul 17 peter 380      \relates DiagonalMatrix
3655 13 Jul 17 peter 381    */
3655 13 Jul 17 peter 382   template<class T>
3655 13 Jul 17 peter 383   expression::DiagonalMatrix<BasicMatrix<T>, DiagonalMatrix, expression::Plus>
3655 13 Jul 17 peter 384   operator+(const BasicMatrix<T>& lhs, const DiagonalMatrix& rhs)
3655 13 Jul 17 peter 385   {
3655 13 Jul 17 peter 386     YAT_ASSERT(lhs.rows() == rhs.rows());
3655 13 Jul 17 peter 387     YAT_ASSERT(lhs.columns() == rhs.columns());
3655 13 Jul 17 peter 388     return expression::DiagonalMatrix<BasicMatrix<T>, DiagonalMatrix,
3655 13 Jul 17 peter 389                                       expression::Plus>(lhs, rhs);
3655 13 Jul 17 peter 390   }
3655 13 Jul 17 peter 391
3655 13 Jul 17 peter 392
3655 13 Jul 17 peter 393   /**
3655 13 Jul 17 peter 394      \brief matrix matrix subtraction
3655 13 Jul 17 peter 395
3655 13 Jul 17 peter 396      \since new in yat 0.15
3655 13 Jul 17 peter 397
3655 13 Jul 17 peter 398      \relates DiagonalMatrix
3655 13 Jul 17 peter 399    */
3655 13 Jul 17 peter 400   DiagonalMatrix
3655 13 Jul 17 peter 401   operator-(const DiagonalMatrix& lhs, const DiagonalMatrix& rhs);
3655 13 Jul 17 peter 402
3655 13 Jul 17 peter 403
3655 13 Jul 17 peter 404   /**
3655 13 Jul 17 peter 405      \brief matrix matrix subtraction
3655 13 Jul 17 peter 406
3655 13 Jul 17 peter 407      \since new in yat 0.15
3655 13 Jul 17 peter 408
3655 13 Jul 17 peter 409      \relates DiagonalMatrix
3655 13 Jul 17 peter 410    */
3655 13 Jul 17 peter 411   template<class T>
3655 13 Jul 17 peter 412   expression::DiagonalMatrix<DiagonalMatrix, BasicMatrix<T>, expression::Minus>
3655 13 Jul 17 peter 413   operator-(const DiagonalMatrix& lhs, const BasicMatrix<T>& rhs)
3655 13 Jul 17 peter 414   {
3655 13 Jul 17 peter 415     YAT_ASSERT(lhs.rows() == rhs.rows());
3655 13 Jul 17 peter 416     YAT_ASSERT(lhs.columns() == rhs.columns());
3655 13 Jul 17 peter 417     return expression::DiagonalMatrix<DiagonalMatrix, BasicMatrix<T>,
3655 13 Jul 17 peter 418                                       expression::Minus>(lhs, rhs);
3655 13 Jul 17 peter 419   }
3655 13 Jul 17 peter 420
3655 13 Jul 17 peter 421
3655 13 Jul 17 peter 422   /**
3655 13 Jul 17 peter 423      \brief matrix matrix subtraction
3655 13 Jul 17 peter 424
3655 13 Jul 17 peter 425      \since new in yat 0.15
3655 13 Jul 17 peter 426
3655 13 Jul 17 peter 427      \relates DiagonalMatrix
3655 13 Jul 17 peter 428    */
3655 13 Jul 17 peter 429   template<class T>
3655 13 Jul 17 peter 430   expression::DiagonalMatrix<BasicMatrix<T>, DiagonalMatrix, expression::Minus>
3655 13 Jul 17 peter 431   operator-(const BasicMatrix<T>& lhs, const DiagonalMatrix& rhs)
3655 13 Jul 17 peter 432   {
3655 13 Jul 17 peter 433     YAT_ASSERT(lhs.rows() == rhs.rows());
3655 13 Jul 17 peter 434     YAT_ASSERT(lhs.columns() == rhs.columns());
3655 13 Jul 17 peter 435     return expression::DiagonalMatrix<BasicMatrix<T>, DiagonalMatrix,
3655 13 Jul 17 peter 436                                       expression::Minus>(lhs, rhs);
3655 13 Jul 17 peter 437   }
3655 13 Jul 17 peter 438
3655 13 Jul 17 peter 439
3655 13 Jul 17 peter 440   /**
3655 13 Jul 17 peter 441      \brief matrix vector multiplication
3655 13 Jul 17 peter 442
3655 13 Jul 17 peter 443      \since new in yat 0.15
3655 13 Jul 17 peter 444
3655 13 Jul 17 peter 445      \relates DiagonalMatrix
3655 13 Jul 17 peter 446    */
3655 13 Jul 17 peter 447   template<class T>
3655 13 Jul 17 peter 448   expression::DiagonalMatrixVector<T>
3655 13 Jul 17 peter 449   operator*(const DiagonalMatrix& lhs, const BasicVector<T>& rhs)
3655 13 Jul 17 peter 450   {
3655 13 Jul 17 peter 451     YAT_ASSERT(lhs.columns() == rhs.size());
3655 13 Jul 17 peter 452     return expression::DiagonalMatrixVector<T>(lhs, rhs);
3655 13 Jul 17 peter 453   }
3655 13 Jul 17 peter 454
3655 13 Jul 17 peter 455
3655 13 Jul 17 peter 456   /**
3655 13 Jul 17 peter 457      \brief vector matrix multiplication
3655 13 Jul 17 peter 458
3655 13 Jul 17 peter 459      \since new in yat 0.15
3655 13 Jul 17 peter 460
3655 13 Jul 17 peter 461      \relates DiagonalMatrix
3655 13 Jul 17 peter 462    */
3655 13 Jul 17 peter 463   template<class T>
3655 13 Jul 17 peter 464   expression::DiagonalMatrixVector<T>
3655 13 Jul 17 peter 465   operator*(const BasicVector<T>& lhs, const DiagonalMatrix& rhs)
3655 13 Jul 17 peter 466   {
3655 13 Jul 17 peter 467     YAT_ASSERT(lhs.size() == rhs.rows());
3655 13 Jul 17 peter 468     return   expression::DiagonalMatrixVector<T>(lhs, rhs);
3655 13 Jul 17 peter 469   }
3655 13 Jul 17 peter 470
3655 13 Jul 17 peter 471
3655 13 Jul 17 peter 472 }}}
3655 13 Jul 17 peter 473 #endif

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