Cleanup sparse matrix iteration patterns

Improves readability and consistency of sparse matrix element access. This is as a preparation work for a new matrix from dune-istl.
* Use begin and end creator iterators accounting that they may be different types.
* Use iterator arithmetic instead of row implementation details.
* Use sparseRange idiom to loop over the row sparsity pattern.
* Avoid taking mutable references of rows, given that a matrix may return proxy objects.
* Do not assign an non-zero numeric value to the jacobian during its pattern creation.
* Include missing header for block vectors

flowexperimental/comp/wells/CompWellEquations_impl.hpp

@@ -41,14 +41,14 @@ init(const int num_conn,  const std::vector<std::size_t>& cells)
     duneB_.setSize(1, num_conn, num_conn);
     duneC_.setSize(1, num_conn, num_conn);
 
-    for (auto row = duneD_.createbegin(),
-              end = duneD_.createend(); row != end; ++row) {
+    auto endD = duneD_.createend();
+    for (auto row = duneD_.createbegin(); row != endD; ++row) {
         // Add nonzeros for diagonal
         row.insert(row.index());
     }
 
-    for (auto row = duneB_.createbegin(),
-                 end = duneB_.createend(); row != end; ++row) {
+    auto endB = duneB_.createend();
+    for (auto row = duneB_.createbegin(); row != endB; ++row) {
         for (int con = 0 ; con < num_conn; ++con) {
             row.insert(con);
         }
@@ -56,8 +56,8 @@ init(const int num_conn,  const std::vector<std::size_t>& cells)
 
     // make the C^T matrix
     // TODO: let us see whether we should change the naming of DuneC_
-    for (auto row = duneC_.createbegin(),
-                 end = duneC_.createend(); row != end; ++row) {
+    auto endC = duneC_.createend();
+    for (auto row = duneC_.createbegin(); row != endC; ++row) {
         for (int con = 0; con < num_conn; ++con) {
             row.insert(con);
         }

opm/simulators/linalg/MILU.cpp

@@ -167,23 +167,21 @@ void milun_decomposition(const M& A, int n, MILU_VARIANT milu, M& ILU,
 
     auto iluRow = ILU.createbegin();
 
-    for(std::size_t i = 0, iend = A.N(); i < iend; ++i)
+    for (auto&& [_, i] : sparseRange(A))
     {
-        auto& orow = A[inverseOrdering[i]];
-
         Map rowPattern;
-        for ( auto col = orow.begin(), cend = orow.end(); col != cend; ++col)
+        for (auto&& [_, j] : sparseRange(A[inverseOrdering[i]]))
         {
-            rowPattern[ordering[col.index()]] = 0;
+            rowPattern[ordering[j]] = 0;
         }
 
         for(auto ik = rowPattern.begin(); ik->first < i; ++ik)
         {
             if ( ik->second < n )
             {
-                auto& rowk = ILU[ik->first];
+                auto rowk = ILU.begin() + ik->first;
 
-                for ( auto kj = rowk.find(ik->first), endk = rowk.end();
+                for ( auto kj = rowk->find(ik->first), endk = rowk->end();
                       kj != endk; ++kj)
                 {
                     // Assume double and block_type FieldMatrix
@@ -220,14 +218,14 @@ void milun_decomposition(const M& A, int n, MILU_VARIANT milu, M& ILU,
     // copy Entries from A
     for(auto iter=A.begin(), iend = A.end(); iter != iend; ++iter)
     {
-        auto& newRow = ILU[ordering[iter.index()]];
         // reset stored generation
-        std::ranges::fill(newRow, 0);
+        auto newRowIt = ILU.begin() + ordering[iter.index()];
+        std::ranges::fill(*newRowIt, 0);
 
         // copy row.
         for(auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
         {
-            newRow[ordering[col.index()]] = *col;
+            (*newRowIt)[ordering[col.index()]] = *col;
         }
     }
     // call decomposition on pattern

opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp

@@ -465,9 +465,8 @@ update()
                     // have changed, so we must copy all elements.
                     for (std::size_t row = 0; row < A_->N(); ++row) {
                         const auto& Arow = (*A_)[row];
-                        auto& ILUrow = (*ILU_)[row];
                         auto Ait = Arow.begin();
-                        auto Iit = ILUrow.begin();
+                        auto Iit = (*ILU_)[row].begin();
                         for (; Ait != Arow.end(); ++Ait, ++Iit) {
                             *Iit = *Ait;
                         }
@@ -484,7 +483,8 @@ update()
                                                 A_->nonzeroes(), Matrix::row_wise);
                 auto& newA = *ILU_;
                 // Create sparsity pattern
-                for (auto iter = newA.createbegin(), iend = newA.createend(); iter != iend; ++iter)
+                auto endcreateA = newA.createend();
+                for (auto iter = newA.createbegin(); iter != endcreateA; ++iter)
                 {
                     const auto& row = (*A_)[inverseOrdering[iter.index()]];
                     for (auto col = row.begin(), cend = row.end(); col != cend; ++col)
@@ -493,12 +493,12 @@ update()
                     }
                 }
                 // Copy values
-                for (auto iter = A_->begin(), iend = A_->end(); iter != iend; ++iter)
+                for (auto iter = A_->begin(); iter != A_->end(); ++iter)
                 {
-                    auto& newRow = newA[ordering_[iter.index()]];
-                    for (auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
+                    auto newRow = newA.begin() + ordering_[iter.index()];
+                    for (auto&& [A_ij, j] : sparseRange(*newRow))
                     {
-                        newRow[ordering_[col.index()]] = *col;
+                        (*newRow)[ordering_[j]] = A_ij;
                     }
                 }
             }

opm/simulators/linalg/gpubridge/GpuBridge.cpp

@@ -160,7 +160,7 @@ int replaceZeroDiagonal(BridgeMatrix& mat,
     } else {
         for (typename BridgeMatrix::iterator r = mat.begin(); r != mat.end(); ++r) {
             typename BridgeMatrix::size_type offset = diag_indices[r.index()];
-            auto& diag_block = r->getptr()[offset]; // diag_block is a reference to MatrixBlock, located on column r of row r
+            auto& diag_block = *std::next(r->begin(), offset); // diag_block is a reference to MatrixBlock, located on column r of row r
             for (int rr = 0; rr < dim; ++rr) {
                 auto& val = diag_block[rr][rr];
                 if (val == 0.0) {                     // could be replaced by '< 1e-30' or similar

opm/simulators/linalg/istlsparsematrixadapter.hh

@@ -138,15 +138,11 @@ public:
         for (int i = 0; i < diagBlock.rows; ++i)
             diagBlock[i][i] = diag;
 
-        auto& matRow = (*istlMatrix_)[row];
-        auto colIt = matRow.begin();
-        const auto& colEndIt = matRow.end();
-        for (; colIt != colEndIt; ++colIt) {
-            if (colIt.index() == row)
-                *colIt = diagBlock;
+        for (auto&& [a_ij, j] : sparseRange((*istlMatrix_)[row]))
+            if (j == row)
+                a_ij = diagBlock;
             else
-                *colIt = Scalar(0.0);
-        }
+                a_ij = Scalar(0.0);
     }
 
     /*!

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1493,15 +1493,12 @@ namespace Opm {
     {
         int nw =  this->numLocalWellsEnd();
         int rdofs = local_num_cells_;
-        for (int i = 0; i < nw; ++i) {
-            int wdof = rdofs + i;
-            jacobian.entry(wdof,wdof) = 1.0;// better scaling ?
-        }
         const auto wellconnections = this->getMaxWellConnections();
         for (int i = 0; i < nw; ++i) {
+            int wdof = rdofs + i;
+            jacobian.entry(wdof,wdof) = 0.0;
             const auto& perfcells = wellconnections[i];
             for (int perfcell : perfcells) {
-                int wdof = rdofs + i;
                 jacobian.entry(wdof, perfcell) = 0.0;
                 jacobian.entry(perfcell, wdof) = 0.0;
             }

opm/simulators/wells/MultisegmentWellEquations.cpp

@@ -86,8 +86,8 @@ init(const int numPerfs,
     duneC_.setSize(well_.numberOfSegments(), numPerfs, numPerfs);
 
     // we need to add the off diagonal ones
-    for (auto row = duneD_.createbegin(),
-              end = duneD_.createend(); row != end; ++row) {
+    auto endD = duneD_.createend();
+    for (auto row = duneD_.createbegin(); row != endD; ++row) {
         // the number of the row corrspnds to the segment now
         const int seg = row.index();
         // adding the item related to outlet relation
@@ -108,8 +108,8 @@ init(const int numPerfs,
     }
 
     // make the C matrix
-    for (auto row = duneC_.createbegin(),
-              end = duneC_.createend(); row != end; ++row) {
+    auto endC = duneC_.createend();
+    for (auto row = duneC_.createbegin(); row != endC; ++row) {
         // the number of the row corresponds to the segment number now.
         for (const int& perf : perforations[row.index()]) {
             const int local_perf_index = parallel_well_info.activePerfToLocalPerf(perf);
@@ -120,8 +120,8 @@ init(const int numPerfs,
     }
 
     // make the B^T matrix
-    for (auto row = duneB_.createbegin(),
-              end = duneB_.createend(); row != end; ++row) {
+    auto endB = duneB_.createend();
+    for (auto row = duneB_.createbegin(); row != endB; ++row) {
         // the number of the row corresponds to the segment number now.
         for (const int& perf : perforations[row.index()]) {
             const int local_perf_index = parallel_well_info.activePerfToLocalPerf(perf);

opm/simulators/wells/StandardWellEquations.cpp

@@ -66,16 +66,16 @@ init(const int numWellEq,
     duneB_.setSize(1, numPerfs, numPerfs);
     duneC_.setSize(1, numPerfs, numPerfs);
 
-    for (auto row = duneD_.createbegin(),
-              end = duneD_.createend(); row != end; ++row) {
+    auto endD = duneD_.createend();
+    for (auto row = duneD_.createbegin(); row != endD; ++row) {
         // Add nonzeros for diagonal
         row.insert(row.index());
     }
       // the block size is run-time determined now
     duneD_[0][0].resize(numWellEq, numWellEq);
 
-    for (auto row = duneB_.createbegin(),
-              end = duneB_.createend(); row != end; ++row) {
+    auto endB = duneB_.createend();
+    for (auto row = duneB_.createbegin(); row != endB; ++row) {
         for (int perf = 0 ; perf < numPerfs; ++perf) {
             row.insert(perf);
         }
@@ -87,8 +87,8 @@ init(const int numWellEq,
     }
 
          // make the C^T matrix
-    for (auto row = duneC_.createbegin(),
-              end = duneC_.createend(); row != end; ++row) {
+    auto endC = duneC_.createend();
+    for (auto row = duneC_.createbegin(); row != endC; ++row) {
         for (int perf = 0; perf < numPerfs; ++perf) {
             row.insert(perf);
         }

opm/simulators/wells/WellHelpers.cpp

@@ -31,6 +31,8 @@
 
 #include <opm/simulators/wells/ParallelWellInfo.hpp>
 
+#include <dune/istl/bvector.hh>
+
 #include <fmt/format.h>
 
 #include <cstddef>