Changeset 4105

Timestamp:

Sep 24, 2021, 6:55:48 AM (2 years ago)

Author:

Peter

Message:

implement +=, -=, and *= operators for DiagonalMatrix?; prefer using arithmatics provided by Vector class rather than own loops.

Location:

trunk

Files:

• test/diagonal_matrix.cc (modified) (2 diffs)
• yat/utility/DiagonalMatrix.cc (modified) (5 diffs)
• yat/utility/DiagonalMatrix.h (modified) (2 diffs)

trunk/test/diagonal_matrix.cc

@@ -2 +2 @@
 
 /*
-  Copyright (C) 2017 Peter Johansson
+  Copyright (C) 2017, 2021 Peter Johansson
 
   This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -118 +118 @@
   suite.out() << "DiagonalMatrix * DiagnoalMatrix\n";
   check(D1 * D1, D2 * D2, suite);
-
+  // test different dimension cases square x square; square x portrait etc
+  for (size_t i=1; i<4; ++i)
+    for (size_t j=1; j<4; ++j) {
+      DiagonalMatrix X1(i,j);
+      Matrix M1(i,j);
+      for (size_t k1=0; k1<X1.rows() && k1<X1.columns(); ++k1) {
+        X1(k1) = 10+k1;
+        M1(k1, k1) = X1(k1, k1);
+      }
+      for (size_t k=1; k<4; ++k) {
+        DiagonalMatrix X2(j,k);
+        Matrix M2(j,k);
+        for (size_t k2=0; k2<X2.rows() && k2<X2.columns(); ++k2) {
+          X2(k2) = 2+k2;
+          M2(k2, k2) = X2(k2, k2);
+        }
+        check(X1*X2, M1*M2, suite);
+      }
+    }
   suite.out() << "DiagonalMatrix * Matrix\n";
   check(D1 * B, D2 * B, suite);

trunk/yat/utility/DiagonalMatrix.cc

@@ -2 +2 @@
 
 /*
-  Copyright (C) 2017 Peter Johansson
+  Copyright (C) 2017, 2021 Peter Johansson
 
   This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -26 +26 @@
 #include "Matrix.h"
 #include "VectorBase.h"
+#include "VectorConstView.h"
+#include "VectorView.h"
 
 #include <algorithm>
@@ -90 +92 @@
 
 
+  DiagonalMatrix& DiagonalMatrix::operator*=(const DiagonalMatrix& rhs)
+  {
+    assert(columns() == rhs.rows());
+    // length of new diagonal
+    size_t n = std::min(row_, rhs.col_);
+    size_t n1 = data_.size();
+    size_t n2 = rhs.data_.size();
+    assert(n==n1 || (n>n1 && n>n2) || (n<n1 && n==n2));
+
+    if (n == n1) {
+      assert(n <= n2);
+      data_.mul(VectorConstView(rhs.data_, 0, n));
+    }
+    else if (n>n1) {
+      assert(n>n2);
+      Vector tmp(n, 0);
+      if (n1 <= n2) {
+        VectorView view(tmp, 0, n1);
+        view = data_;
+        view.mul(VectorConstView(rhs.data_, 0, n1));
+      }
+      else {
+        VectorView view(tmp, 0, n2);
+        view = VectorConstView(data_, 0, n2);
+        view.mul(rhs.data_);
+      }
+      data_ = std::move(tmp);
+    }
+    else {
+      assert(n < n1);
+      assert(n == n2);
+      Vector tmp(rhs.data_);
+      tmp.mul(VectorConstView(data_, 0, n));
+      data_ = std::move(tmp);
+    }
+    col_ = rhs.col_;
+    assert(data_.size() == n);
+    return *this;
+  }
+
+
   DiagonalMatrix operator*(const DiagonalMatrix& lhs, const DiagonalMatrix& rhs)
   {
     assert(lhs.columns() == rhs.rows());
-    DiagonalMatrix res(lhs.rows(), rhs.columns());
-    size_t n = std::min(res.rows(), res.columns());
-    for (size_t i=0; i<n; ++i)
-      res(i) = lhs(i,i) * rhs(i,i);
+    DiagonalMatrix res(lhs);
+    res *= rhs;
     return res;
+  }
+
+
+  DiagonalMatrix& DiagonalMatrix::operator+=(const DiagonalMatrix& rhs)
+  {
+    assert(rows() == rhs.rows());
+    assert(columns() == rhs.columns());
+    assert(data_.size() == rhs.data_.size());
+    data_ += rhs.data_;
+    return *this;
   }
 
@@ -107 +158 @@
     assert(lhs.columns() == rhs.columns());
     DiagonalMatrix res(lhs);
-    size_t n = std::min(res.rows(), res.columns());
-    for (size_t i=0; i<n; ++i)
-      res(i) += rhs(i, i);
+    res += rhs;
     return res;
+  }
+
+
+  DiagonalMatrix& DiagonalMatrix::operator-=(const DiagonalMatrix& rhs)
+  {
+    assert(rows() == rhs.rows());
+    assert(columns() == rhs.columns());
+    assert(data_.size() == rhs.data_.size());
+    data_ -= rhs.data_;
+    return *this;
   }
 
@@ -120 +179 @@
     assert(lhs.columns() == rhs.columns());
     DiagonalMatrix res(lhs);
-    size_t n = std::min(res.rows(), res.columns());
-    for (size_t i=0; i<n; ++i)
-      res(i) -= rhs(i, i);
+    res -= rhs;
     return res;
   }

trunk/yat/utility/DiagonalMatrix.h

@@ -5 +5 @@
 
 /*
-  Copyright (C) 2017 Peter Johansson
+  Copyright (C) 2017, 2021 Peter Johansson
 
   This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -92 +92 @@
      */
     double& operator()(size_t i);
+
+    /**
+       \brief Multiplication and assign operator
+
+       Same as doing *this = *this * rhs.
+
+       \since New in yat 0.20
+     */
+    DiagonalMatrix& operator*=(const DiagonalMatrix& rhs);
+
+    /**
+       \brief Add and assign operator
+
+       Elementwise addition of the diagnonal elements.
+
+       \since New in yat 0.20
+     */
+    DiagonalMatrix& operator+=(const DiagonalMatrix& rhs);
+
+    /**
+       \brief Subtract and assign operator
+
+       Elementwise addition of the diagnonal elements.
+
+       \since New in yat 0.20
+     */
+    DiagonalMatrix& operator-=(const DiagonalMatrix& rhs);
   private:
     Vector data_;