Lazy evaluation of operations on global matrices

CHANGELOG.md

@@ -15,6 +15,7 @@
 4. update `SuperLU` to version `6.0.1`
 5. better `RCM` algorithm that may result in a smaller bandwidth, thus, poentially more efficient solving [#175](https://github.com/TLCFEM/suanPan/pull/175)
 6. update `Armadillo` to version `12.6.7` [#180](https://github.com/TLCFEM/suanPan/pull/180)
+7. enable lazy evalution and avoid temporary global matrices, faster dynamic analysis [#183](https://github.com/TLCFEM/suanPan/pull/183)
 
 ## version 3.2
 

Domain/MetaMat/DenseMat.hpp

@@ -92,40 +92,26 @@ template<sp_d T> class DenseMat : public MetaMat<T> {
 
     T* memptr() override { return memory.get(); }
 
-    void operator+=(const shared_ptr<MetaMat<T>>& M) override {
+    void scale_accu(const T scalar, const shared_ptr<MetaMat<T>>& M) override {
         if(nullptr == M) return;
-        if(!M->triplet_mat.is_empty()) return this->operator+=(M->triplet_mat);
+        if(!M->triplet_mat.is_empty()) return this->scale_accu(scalar, M->triplet_mat);
         if(this->n_rows != M->n_rows || this->n_cols != M->n_cols || this->n_elem != M->n_elem) throw invalid_argument("size mismatch");
         if(nullptr == M->memptr()) return;
         this->factored = false;
-        arrayops::inplace_plus(memptr(), M->memptr(), this->n_elem);
+        if(1. == scalar) arrayops::inplace_plus(memptr(), M->memptr(), this->n_elem);
+        else if(-1. == scalar) arrayops::inplace_minus(memptr(), M->memptr(), this->n_elem);
+        else suanpan::for_each(this->n_elem, [&](const uword I) { memptr()[I] += scalar * M->memptr()[I]; });
     }
 
-    void operator-=(const shared_ptr<MetaMat<T>>& M) override {
-        if(nullptr == M) return;
-        if(!M->triplet_mat.is_empty()) return this->operator-=(M->triplet_mat);
-        if(this->n_rows != M->n_rows || this->n_cols != M->n_cols || this->n_elem != M->n_elem) throw invalid_argument("size mismatch");
-        if(nullptr == M->memptr()) return;
-        this->factored = false;
-        arrayops::inplace_minus(memptr(), M->memptr(), this->n_elem);
-    }
-
-    void operator+=(const triplet_form<T, uword>& M) override {
-        if(this->n_rows != M.n_rows || this->n_cols != M.n_cols) throw invalid_argument("size mismatch");
-        this->factored = false;
-        const auto row = M.row_mem();
-        const auto col = M.col_mem();
-        const auto val = M.val_mem();
-        for(uword I = 0llu; I < M.n_elem; ++I) this->at(row[I], col[I]) += val[I];
-    }
-
-    void operator-=(const triplet_form<T, uword>& M) override {
+    void scale_accu(const T scalar, const triplet_form<T, uword>& M) override {
         if(this->n_rows != M.n_rows || this->n_cols != M.n_cols) throw invalid_argument("size mismatch");
         this->factored = false;
         const auto row = M.row_mem();
         const auto col = M.col_mem();
         const auto val = M.val_mem();
-        for(uword I = 0llu; I < M.n_elem; ++I) this->at(row[I], col[I]) -= val[I];
+        if(1. == scalar) for(auto I = 0llu; I < M.n_elem; ++I) this->at(row[I], col[I]) += val[I];
+        else if(-1. == scalar) for(auto I = 0llu; I < M.n_elem; ++I) this->at(row[I], col[I]) -= val[I];
+        else for(auto I = 0llu; I < M.n_elem; ++I) this->at(row[I], col[I]) += scalar * val[I];
     }
 
     void operator*=(const T value) override {

Domain/MetaMat/MetaMat.hpp

@@ -36,6 +36,39 @@
 
 template<typename T, typename U> concept ArmaContainer = std::is_floating_point_v<U> && (std::is_convertible_v<T, Mat<U>> || std::is_convertible_v<T, SpMat<U>>) ;
 
+template<sp_d T> class MetaMat;
+
+template<sp_d T> class op_add {
+    friend MetaMat<T>;
+
+    shared_ptr<MetaMat<T>> mat_a, mat_b;
+
+public:
+    explicit op_add(const shared_ptr<MetaMat<T>>& A)
+        : mat_a(A)
+        , mat_b(nullptr) {}
+
+    op_add(const shared_ptr<MetaMat<T>>& A, const shared_ptr<MetaMat<T>>& B)
+        : mat_a(A)
+        , mat_b(B) {}
+};
+
+template<sp_d T> class op_scale {
+    friend MetaMat<T>;
+
+    T scalar;
+    op_add<T> bracket;
+
+public:
+    op_scale(const T A, const shared_ptr<MetaMat<T>>& B)
+        : scalar(A)
+        , bracket(B) {}
+
+    op_scale(const T A, op_add<T>&& B)
+        : scalar(A)
+        , bracket(std::forward<op_add<T>>(B)) {}
+};
+
 template<sp_d T> class MetaMat {
 protected:
     bool factored = false;
@@ -139,11 +172,28 @@ template<sp_d T> class MetaMat {
     [[nodiscard]] virtual const T* memptr() const = 0;
     virtual T* memptr() = 0;
 
-    virtual void operator+=(const shared_ptr<MetaMat>&) = 0;
-    virtual void operator-=(const shared_ptr<MetaMat>&) = 0;
+    virtual void scale_accu(T, const shared_ptr<MetaMat>&) = 0;
+    virtual void scale_accu(T, const triplet_form<T, uword>&) = 0;
+
+    void operator+=(const shared_ptr<MetaMat>& M) { return this->scale_accu(1., M); }
+
+    void operator-=(const shared_ptr<MetaMat>& M) { return this->scale_accu(-1., M); }
+
+    void operator+=(const op_scale<T>& M) {
+        const auto& bracket = M.bracket;
+        if(nullptr != bracket.mat_a) this->scale_accu(M.scalar, bracket.mat_a);
+        if(nullptr != bracket.mat_b) this->scale_accu(M.scalar, bracket.mat_b);
+    }
+
+    void operator-=(const op_scale<T>& M) {
+        const auto& bracket = M.bracket;
+        if(nullptr != bracket.mat_a) this->scale_accu(-M.scalar, bracket.mat_a);
+        if(nullptr != bracket.mat_b) this->scale_accu(-M.scalar, bracket.mat_b);
+    }
+
+    void operator+=(const triplet_form<T, uword>& M) { return this->scale_accu(1., M); }
 
-    virtual void operator+=(const triplet_form<T, uword>&) = 0;
-    virtual void operator-=(const triplet_form<T, uword>&) = 0;
+    void operator-=(const triplet_form<T, uword>& M) { return this->scale_accu(-1., M); }
 
     virtual Mat<T> operator*(const Mat<T>&) const = 0;
 

Domain/MetaMat/SparseMat.hpp

@@ -68,28 +68,18 @@ template<sp_d T> class SparseMat : public MetaMat<T> {
 
     T* memptr() override { throw invalid_argument("not supported"); }
 
-    void operator+=(const shared_ptr<MetaMat<T>>& in_mat) override {
+    void scale_accu(const T scalar, const shared_ptr<MetaMat<T>>& in_mat) override {
         if(nullptr == in_mat) return;
-        if(!in_mat->triplet_mat.is_empty()) return this->operator+=(in_mat->triplet_mat);
+        if(!in_mat->triplet_mat.is_empty()) return this->scale_accu(scalar, in_mat->triplet_mat);
         this->factored = false;
-        for(uword I = 0llu; I < in_mat->n_rows; ++I) for(uword J = 0llu; J < in_mat->n_cols; ++J) if(const auto t_val = in_mat->operator()(I, J); !suanpan::approx_equal(T(0), t_val)) at(I, J) = t_val;
+        for(auto I = 0llu; I < in_mat->n_rows; ++I) for(auto J = 0llu; J < in_mat->n_cols; ++J) if(const auto t_val = in_mat->operator()(I, J); !suanpan::approx_equal(T(0), t_val)) at(I, J) = scalar * t_val;
     }
 
-    void operator-=(const shared_ptr<MetaMat<T>>& in_mat) override {
-        if(nullptr == in_mat) return;
-        if(!in_mat->triplet_mat.is_empty()) return this->operator-=(in_mat->triplet_mat);
-        this->factored = false;
-        for(uword I = 0llu; I < in_mat->n_rows; ++I) for(uword J = 0llu; J < in_mat->n_cols; ++J) if(const auto t_val = in_mat->operator()(I, J); !suanpan::approx_equal(T(0), t_val)) at(I, J) = -t_val;
-    }
-
-    void operator+=(const triplet_form<T, uword>& in_mat) override {
-        this->factored = false;
-        this->triplet_mat += in_mat;
-    }
-
-    void operator-=(const triplet_form<T, uword>& in_mat) override {
+    void scale_accu(const T scalar, const triplet_form<T, uword>& in_mat) override {
         this->factored = false;
-        this->triplet_mat -= in_mat;
+        if(1. == scalar) this->triplet_mat += in_mat;
+        else if(-1. == scalar) this->triplet_mat -= in_mat;
+        else this->triplet_mat.assemble(in_mat, 0, 0, scalar);
     }
 
     Mat<T> operator*(const Mat<T>& in_mat) const override { return this->triplet_mat * in_mat; }

Domain/MetaMat/operator_times.hpp

@@ -18,20 +18,31 @@
 #ifndef OPERATOR_TIMES_HPP
 #define OPERATOR_TIMES_HPP
 
+template<sp_d T> op_add<T> operator+(const shared_ptr<MetaMat<T>>& A, const shared_ptr<MetaMat<T>>& B) { return op_add<T>(A, B); }
+
+template<sp_d T> op_scale<T> operator*(const T value, const shared_ptr<MetaMat<T>>& M) { return op_scale<T>(value, M); }
+
+template<sp_d T> op_scale<T> operator*(const T value, op_add<T>&& M) { return op_scale<T>(value, std::forward<op_add<T>>(M)); }
+
+template<sp_d T> const shared_ptr<MetaMat<T>>& operator+=(const shared_ptr<MetaMat<T>>& M, const op_scale<T>& A) {
+    M->operator+=(A);
+    return M;
+}
+
 template<sp_d T> unique_ptr<MetaMat<T>> operator*(const T value, unique_ptr<MetaMat<T>>&& M) {
     if(nullptr == M) return nullptr;
 
     M->operator*=(value);
     return std::forward<unique_ptr<MetaMat<T>>>(M);
 }
 
-template<sp_d T> unique_ptr<MetaMat<T>> operator*(const T value, const shared_ptr<MetaMat<T>>& M) {
-    if(nullptr == M) return nullptr;
-
-    auto N = M->make_copy();
-    N->operator*=(value);
-    return N;
-}
+//template<sp_d T> unique_ptr<MetaMat<T>> operator*(const T value, const shared_ptr<MetaMat<T>>& M) {
+//    if(nullptr == M) return nullptr;
+//
+//    auto N = M->make_copy();
+//    N->operator*=(value);
+//    return N;
+//}
 
 template<sp_d T> Mat<T> operator*(const shared_ptr<MetaMat<T>>& M, const Mat<T>& A) { return M->operator*(A); }
 
@@ -42,11 +53,11 @@ template<sp_d T> const shared_ptr<MetaMat<T>>& operator*=(const shared_ptr<MetaM
     return M;
 }
 
-template<sp_d T> unique_ptr<MetaMat<T>> operator+(const shared_ptr<MetaMat<T>>& A, const shared_ptr<MetaMat<T>>& B) {
-    auto M = B->make_copy();
-    M->operator+=(A);
-    return M;
-}
+//template<sp_d T> unique_ptr<MetaMat<T>> operator+(const shared_ptr<MetaMat<T>>& A, const shared_ptr<MetaMat<T>>& B) {
+//    auto M = B->make_copy();
+//    M->operator+=(A);
+//    return M;
+//}
 
 template<sp_d T> unique_ptr<MetaMat<T>> operator+(const shared_ptr<MetaMat<T>>& A, unique_ptr<MetaMat<T>>&& B) {
     B->operator+=(A);