yat/utility/DiagonalMatrix.cc

Code
Comments
Other
Rev Date Author Line
3655 13 Jul 17 peter 1 // $Id$
3655 13 Jul 17 peter 2
3655 13 Jul 17 peter 3 /*
4207 26 Aug 22 peter 4   Copyright (C) 2017, 2021, 2022 Peter Johansson
3655 13 Jul 17 peter 5
3655 13 Jul 17 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
3655 13 Jul 17 peter 7
3655 13 Jul 17 peter 8   The yat library is free software; you can redistribute it and/or
3655 13 Jul 17 peter 9   modify it under the terms of the GNU General Public License as
3655 13 Jul 17 peter 10   published by the Free Software Foundation; either version 3 of the
3655 13 Jul 17 peter 11   License, or (at your option) any later version.
3655 13 Jul 17 peter 12
3655 13 Jul 17 peter 13   The yat library is distributed in the hope that it will be useful,
3655 13 Jul 17 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
3655 13 Jul 17 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3655 13 Jul 17 peter 16   General Public License for more details.
3655 13 Jul 17 peter 17
3655 13 Jul 17 peter 18   You should have received a copy of the GNU General Public License
3655 13 Jul 17 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
3655 13 Jul 17 peter 20 */
3655 13 Jul 17 peter 21
3655 13 Jul 17 peter 22 #include <config.h>
3655 13 Jul 17 peter 23
3655 13 Jul 17 peter 24 #include "DiagonalMatrix.h"
3655 13 Jul 17 peter 25
3655 13 Jul 17 peter 26 #include "Matrix.h"
3655 13 Jul 17 peter 27 #include "VectorBase.h"
4105 24 Sep 21 peter 28 #include "VectorConstView.h"
4105 24 Sep 21 peter 29 #include "VectorView.h"
3655 13 Jul 17 peter 30
3655 13 Jul 17 peter 31 #include <algorithm>
3655 13 Jul 17 peter 32 #include <cassert>
3655 13 Jul 17 peter 33
3655 13 Jul 17 peter 34 namespace theplu {
3655 13 Jul 17 peter 35 namespace yat {
3655 13 Jul 17 peter 36 namespace utility {
3655 13 Jul 17 peter 37
3655 13 Jul 17 peter 38   DiagonalMatrix::DiagonalMatrix(void)
3655 13 Jul 17 peter 39     : row_(0), col_(0)
3655 13 Jul 17 peter 40   {
3655 13 Jul 17 peter 41   }
3655 13 Jul 17 peter 42
3655 13 Jul 17 peter 43
3655 13 Jul 17 peter 44   DiagonalMatrix::DiagonalMatrix(size_t r, size_t c, double x)
3655 13 Jul 17 peter 45     : data_(std::min(r, c), x), row_(r), col_(c)
3655 13 Jul 17 peter 46   {
3655 13 Jul 17 peter 47   }
3655 13 Jul 17 peter 48
3655 13 Jul 17 peter 49
3655 13 Jul 17 peter 50   DiagonalMatrix::DiagonalMatrix(const VectorBase& vec)
3655 13 Jul 17 peter 51     : data_(vec), row_(vec.size()), col_(vec.size())
3655 13 Jul 17 peter 52   {
3655 13 Jul 17 peter 53   }
3655 13 Jul 17 peter 54
3655 13 Jul 17 peter 55
4125 14 Jan 22 peter 56   DiagonalMatrix::DiagonalMatrix(const MatrixBase& m)
3655 13 Jul 17 peter 57     : data_(std::min(m.rows(), m.columns())), row_(m.rows()), col_(m.columns())
3655 13 Jul 17 peter 58   {
3655 13 Jul 17 peter 59     for (size_t i=0; i<data_.size(); ++i)
3655 13 Jul 17 peter 60       data_(i) = m(i,i);
3655 13 Jul 17 peter 61   }
3655 13 Jul 17 peter 62
3655 13 Jul 17 peter 63
3655 13 Jul 17 peter 64   size_t DiagonalMatrix::rows(void) const
3655 13 Jul 17 peter 65   {
3655 13 Jul 17 peter 66     return row_;
3655 13 Jul 17 peter 67   }
3655 13 Jul 17 peter 68
3655 13 Jul 17 peter 69
3655 13 Jul 17 peter 70   size_t DiagonalMatrix::columns(void) const
3655 13 Jul 17 peter 71   {
3655 13 Jul 17 peter 72     return col_;
3655 13 Jul 17 peter 73   }
3655 13 Jul 17 peter 74
3655 13 Jul 17 peter 75
3655 13 Jul 17 peter 76   const double DiagonalMatrix::operator()(size_t row, size_t col) const
3655 13 Jul 17 peter 77   {
3655 13 Jul 17 peter 78     assert(row < rows());
3655 13 Jul 17 peter 79     assert(col < columns());
3655 13 Jul 17 peter 80     if (row!=col)
3655 13 Jul 17 peter 81       return 0.0;
3655 13 Jul 17 peter 82     return data_(row);
3655 13 Jul 17 peter 83   }
3655 13 Jul 17 peter 84
3655 13 Jul 17 peter 85
3655 13 Jul 17 peter 86   double& DiagonalMatrix::operator()(size_t i)
3655 13 Jul 17 peter 87   {
3655 13 Jul 17 peter 88     assert(i < rows());
3655 13 Jul 17 peter 89     assert(i < columns());
3655 13 Jul 17 peter 90     return data_(i);
3655 13 Jul 17 peter 91   }
3655 13 Jul 17 peter 92
3655 13 Jul 17 peter 93
4105 24 Sep 21 peter 94   DiagonalMatrix& DiagonalMatrix::operator*=(const DiagonalMatrix& rhs)
4105 24 Sep 21 peter 95   {
4105 24 Sep 21 peter 96     assert(columns() == rhs.rows());
4105 24 Sep 21 peter 97     // length of new diagonal
4105 24 Sep 21 peter 98     size_t n = std::min(row_, rhs.col_);
4105 24 Sep 21 peter 99     size_t n1 = data_.size();
4105 24 Sep 21 peter 100     size_t n2 = rhs.data_.size();
4105 24 Sep 21 peter 101     assert(n==n1 || (n>n1 && n>n2) || (n<n1 && n==n2));
4105 24 Sep 21 peter 102
4105 24 Sep 21 peter 103     if (n == n1) {
4105 24 Sep 21 peter 104       assert(n <= n2);
4105 24 Sep 21 peter 105       data_.mul(VectorConstView(rhs.data_, 0, n));
4105 24 Sep 21 peter 106     }
4105 24 Sep 21 peter 107     else if (n>n1) {
4105 24 Sep 21 peter 108       assert(n>n2);
4105 24 Sep 21 peter 109       Vector tmp(n, 0);
4105 24 Sep 21 peter 110       if (n1 <= n2) {
4105 24 Sep 21 peter 111         VectorView view(tmp, 0, n1);
4105 24 Sep 21 peter 112         view = data_;
4105 24 Sep 21 peter 113         view.mul(VectorConstView(rhs.data_, 0, n1));
4105 24 Sep 21 peter 114       }
4105 24 Sep 21 peter 115       else {
4105 24 Sep 21 peter 116         VectorView view(tmp, 0, n2);
4105 24 Sep 21 peter 117         view = VectorConstView(data_, 0, n2);
4105 24 Sep 21 peter 118         view.mul(rhs.data_);
4105 24 Sep 21 peter 119       }
4105 24 Sep 21 peter 120       data_ = std::move(tmp);
4105 24 Sep 21 peter 121     }
4105 24 Sep 21 peter 122     else {
4105 24 Sep 21 peter 123       assert(n < n1);
4105 24 Sep 21 peter 124       assert(n == n2);
4105 24 Sep 21 peter 125       Vector tmp(rhs.data_);
4105 24 Sep 21 peter 126       tmp.mul(VectorConstView(data_, 0, n));
4105 24 Sep 21 peter 127       data_ = std::move(tmp);
4105 24 Sep 21 peter 128     }
4105 24 Sep 21 peter 129     col_ = rhs.col_;
4105 24 Sep 21 peter 130     assert(data_.size() == n);
4105 24 Sep 21 peter 131     return *this;
4105 24 Sep 21 peter 132   }
4105 24 Sep 21 peter 133
4105 24 Sep 21 peter 134
3655 13 Jul 17 peter 135   DiagonalMatrix operator*(const DiagonalMatrix& lhs, const DiagonalMatrix& rhs)
3655 13 Jul 17 peter 136   {
3655 13 Jul 17 peter 137     assert(lhs.columns() == rhs.rows());
4105 24 Sep 21 peter 138     DiagonalMatrix res(lhs);
4105 24 Sep 21 peter 139     res *= rhs;
3655 13 Jul 17 peter 140     return res;
3655 13 Jul 17 peter 141   }
3655 13 Jul 17 peter 142
3655 13 Jul 17 peter 143
4105 24 Sep 21 peter 144   DiagonalMatrix& DiagonalMatrix::operator+=(const DiagonalMatrix& rhs)
4105 24 Sep 21 peter 145   {
4105 24 Sep 21 peter 146     assert(rows() == rhs.rows());
4105 24 Sep 21 peter 147     assert(columns() == rhs.columns());
4105 24 Sep 21 peter 148     assert(data_.size() == rhs.data_.size());
4105 24 Sep 21 peter 149     data_ += rhs.data_;
4105 24 Sep 21 peter 150     return *this;
4105 24 Sep 21 peter 151   }
4105 24 Sep 21 peter 152
4105 24 Sep 21 peter 153
3655 13 Jul 17 peter 154   DiagonalMatrix
3655 13 Jul 17 peter 155   operator+(const DiagonalMatrix& lhs, const DiagonalMatrix& rhs)
3655 13 Jul 17 peter 156   {
3655 13 Jul 17 peter 157     assert(lhs.rows() == rhs.rows());
3655 13 Jul 17 peter 158     assert(lhs.columns() == rhs.columns());
3655 13 Jul 17 peter 159     DiagonalMatrix res(lhs);
4105 24 Sep 21 peter 160     res += rhs;
3655 13 Jul 17 peter 161     return res;
3655 13 Jul 17 peter 162   }
3655 13 Jul 17 peter 163
3655 13 Jul 17 peter 164
4105 24 Sep 21 peter 165   DiagonalMatrix& DiagonalMatrix::operator-=(const DiagonalMatrix& rhs)
4105 24 Sep 21 peter 166   {
4105 24 Sep 21 peter 167     assert(rows() == rhs.rows());
4105 24 Sep 21 peter 168     assert(columns() == rhs.columns());
4105 24 Sep 21 peter 169     assert(data_.size() == rhs.data_.size());
4105 24 Sep 21 peter 170     data_ -= rhs.data_;
4105 24 Sep 21 peter 171     return *this;
4105 24 Sep 21 peter 172   }
4105 24 Sep 21 peter 173
4105 24 Sep 21 peter 174
3655 13 Jul 17 peter 175   DiagonalMatrix
3655 13 Jul 17 peter 176   operator-(const DiagonalMatrix& lhs, const DiagonalMatrix& rhs)
3655 13 Jul 17 peter 177   {
3655 13 Jul 17 peter 178     assert(lhs.rows() == rhs.rows());
3655 13 Jul 17 peter 179     assert(lhs.columns() == rhs.columns());
3655 13 Jul 17 peter 180     DiagonalMatrix res(lhs);
4105 24 Sep 21 peter 181     res -= rhs;
3655 13 Jul 17 peter 182     return res;
3655 13 Jul 17 peter 183   }
3655 13 Jul 17 peter 184
3655 13 Jul 17 peter 185
3655 13 Jul 17 peter 186 }}}

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