From b8984320717ca2b5e8a51bfe55f3b8d9ede41940 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Thu, 30 Dec 2021 11:00:06 -0700
Subject: [PATCH] Avoid a temp array in the UMFPACK interfaces.

---
 source/lac/sparse_direct.cc | 50 ++++++++++++++++---------------------
 1 file changed, 22 insertions(+), 28 deletions(-)

diff --git a/source/lac/sparse_direct.cc b/source/lac/sparse_direct.cc
index c34239b21a..95a1334087 100644
--- a/source/lac/sparse_direct.cc
+++ b/source/lac/sparse_direct.cc
@@ -258,34 +258,28 @@ SparseDirectUMFPACK::factorize(const Matrix &matrix)
   // iterators are sorted so that they traverse each row from start to end
   // before moving on to the next row. however, this isn't true for block
   // matrices, so we have to do a bit of book keeping
-  {
-    // have an array that for each row points to the first entry not yet
-    // written to
-    std::vector<long int> row_pointers = Ap;
-
-    // loop over the elements of the matrix row by row, as suggested in the
-    // documentation of the sparse matrix iterator class
-    for (size_type row = 0; row < matrix.m(); ++row)
-      {
-        for (typename Matrix::const_iterator p = matrix.begin(row);
-             p != matrix.end(row);
-             ++p)
-          {
-            // write entry into the first free one for this row
-            Ai[row_pointers[row]] = p->column();
-            Ax[row_pointers[row]] = std::real(p->value());
-            if (numbers::NumberTraits<number>::is_complex == true)
-              Az[row_pointers[row]] = std::imag(p->value());
-
-            // then move pointer ahead
-            ++row_pointers[row];
-          }
-      }
-
-    // at the end, we should have written all rows completely
-    for (size_type i = 0; i < Ap.size() - 1; ++i)
-      Assert(row_pointers[i] == Ap[i + 1], ExcInternalError());
-  }
+
+  // loop over the elements of the matrix row by row, as suggested in the
+  // documentation of the sparse matrix iterator class
+  for (size_type row = 0; row < matrix.m(); ++row)
+    {
+      long int index = Ap[row];
+      for (typename Matrix::const_iterator p = matrix.begin(row);
+           p != matrix.end(row);
+           ++p)
+        {
+          // write entry into the first free one for this row
+          Ai[index] = p->column();
+          Ax[index] = std::real(p->value());
+          if (numbers::NumberTraits<number>::is_complex == true)
+            Az[index] = std::imag(p->value());
+
+          // then move pointer ahead
+          ++index;
+        }
+
+      Assert(index == Ap[row + 1], ExcInternalError());
+    }
 
   // make sure that the elements in each row are sorted. we have to be more
   // careful for block sparse matrices, so ship this task out to a
-- 
2.39.5