From 520bd0d3c44ee0f3fbf43dd0fe1d78b7bfe84ffa Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Fri, 20 Oct 2023 09:32:57 +0200
Subject: [PATCH] Added: matrix::augmentRows()

---
 src/LinAlg/Test/TestMatrix.C | 30 +++++++++++++++++++++++++++---
 src/LinAlg/matrix.h          | 25 ++++++++++++++++++++++---
 2 files changed, 49 insertions(+), 6 deletions(-)

diff --git a/src/LinAlg/Test/TestMatrix.C b/src/LinAlg/Test/TestMatrix.C
index cb1c76e7..d49012f4 100644
--- a/src/LinAlg/Test/TestMatrix.C
+++ b/src/LinAlg/Test/TestMatrix.C
@@ -10,7 +10,7 @@
 //!
 //==============================================================================
 
-#include "MatVec.h"
+#include "matrixnd.h"
 
 #include "gtest/gtest.h"
 #include <numeric>
@@ -83,7 +83,29 @@ TEST(TestMatrix, AddRows)
 }
 
 
-TEST(TestMatrix, Augment)
+TEST(TestMatrix, AugmentRows)
+{
+  utl::matrix<int> a(5,3), b(4,3), c(3,2);
+  std::iota(a.begin(),a.end(),1);
+  std::iota(b.begin(),b.end(),16);
+  std::cout <<"A:"<< a;
+  std::cout <<"B:"<< b;
+  size_t nA = a.size();
+  size_t na = a.rows();
+  size_t nb = b.rows();
+  ASSERT_TRUE(a.augmentRows(b));
+  ASSERT_FALSE(a.augmentRows(c));
+  std::cout <<"C:"<< a;
+  for (size_t j = 1; j <= a.cols(); j++)
+    for (size_t i = 1; i <= a.rows(); i++)
+      if (i <= na)
+        EXPECT_EQ(a(i,j), i+na*(j-1));
+      else
+        EXPECT_EQ(a(i,j), nA-na+i+nb*(j-1));
+}
+
+
+TEST(TestMatrix, AugmentCols)
 {
   utl::matrix<int> a(3,5), b(3,4), c(2,3);
   std::iota(a.begin(),a.end(),1);
@@ -128,7 +150,9 @@ TEST(TestMatrix, Multiply)
   EXPECT_FLOAT_EQ(y.sum(),0.0);
 
   u.std::vector<double>::resize(5);
-  v = 0.5*A*u;
+  ASSERT_TRUE(A.multiply(u,v));
+  v *= 0.5;
+
   EXPECT_FLOAT_EQ(v(1),67.5);
   EXPECT_FLOAT_EQ(v(2),75.0);
   EXPECT_FLOAT_EQ(v(3),82.5);
diff --git a/src/LinAlg/matrix.h b/src/LinAlg/matrix.h
index ee1f5a62..fef9127e 100644
--- a/src/LinAlg/matrix.h
+++ b/src/LinAlg/matrix.h
@@ -437,11 +437,31 @@ namespace utl //! General utility classes and functions.
       return *this;
     }
 
+    //! \brief Increase the number of rows by augmenting the given matrix.
+    bool augmentRows(const matrix<T>& B)
+    {
+      if (B.ncol != ncol)
+        return false;
+
+      size_t oldRows = nrow;
+      nrow += B.nrow;
+      this->elem.std::template vector<T>::resize(nrow*ncol,T(0));
+      T* oldMat = this->ptr() + oldRows*(ncol-1);
+      for (size_t c = ncol; c > 0; c--, oldMat -= oldRows)
+      {
+        if (c > 1)
+          memmove(this->ptr(c-1),oldMat,oldRows*sizeof(T));
+        for (size_t r = nrow; r > oldRows; r--)
+          this->elem[r-1+nrow*(c-1)] = B(r-oldRows,c);
+      }
+      return true;
+    }
+
     //! \brief Increase the number of columns by augmenting the given matrix.
     bool augmentCols(const matrix<T>& B)
     {
       if (B.nrow != nrow)
-	return false;
+        return false;
 
       this->elem.insert(this->elem.end(),B.elem.begin(),B.elem.end());
       ncol += B.ncol;
@@ -568,8 +588,7 @@ namespace utl //! General utility classes and functions.
 
     //! \brief Extract a block of the matrix to another matrix.
     void extractBlock(matrix<T>& block, size_t r, size_t c,
-                      bool addTo = false,
-                      bool transposed = false) const
+                      bool addTo = false, bool transposed = false) const
     {
       size_t nr = transposed ? block.cols() : block.rows();
       size_t nc = transposed ? block.rows() : block.cols();