Creat copy implicitly, possibly avoid copy.

dune/matrix-vector/triangularsolve.hh

@@ -4,11 +4,10 @@
 namespace Dune {
 namespace MatrixVector {
   template <class Matrix, class Vector, class BitVector>
-  static void lowerTriangularSolve(Matrix const& L, Vector const& b, Vector& x,
+  static void lowerTriangularSolve(Matrix const& L, Vector b, Vector& x,
                                    BitVector const* ignore,
                                    bool transpose = false) {
     x = 0;
-    Vector r = b;
     if (transpose) {
       for (auto it = L.begin(); it != L.end(); ++it) {
         const size_t i = it.index();
@@ -16,10 +15,10 @@ namespace MatrixVector {
           continue;
         auto cIt = it->begin();
         assert(cIt.index() == it.index());
-        x[i] = r[i] / *cIt;
+        x[i] = b[i] / *cIt;
         for (; cIt != it->end(); ++cIt) {
           const size_t j = cIt.index();
-          r[j] -= x[i] * *cIt;
+          b[j] -= x[i] * *cIt;
         }
       }
     } else {
@@ -30,22 +29,21 @@ namespace MatrixVector {
         for (auto cIt = it->begin(); cIt != it->end(); ++cIt) {
           size_t j = cIt.index();
           if (i == j) {
-            x[i] = r[i] / *cIt;
+            x[i] = b[i] / *cIt;
             break;
           }
           assert(j < i);
-          r[i] -= *cIt * x[j];
+          b[i] -= *cIt * x[j];
         }
       }
     }
   }
 
   template <class Matrix, class Vector, class BitVector>
-  static void upperTriangularSolve(Matrix const& U, Vector const& b, Vector& x,
+  static void upperTriangularSolve(Matrix const& U, Vector b, Vector& x,
                                    BitVector const* ignore,
                                    bool transpose = false) {
     x = 0;
-    Vector r = b;
     if (transpose) {
       for (auto it = U.beforeEnd(); it != U.beforeBegin(); --it) {
         size_t i = it.index();
@@ -54,13 +52,13 @@ namespace MatrixVector {
         auto diagonal = *cIt;
         if (ignore != nullptr and (*ignore)[i].all())
           continue;
-        x[i] = r[i] / diagonal;
+        x[i] = b[i] / diagonal;
 
         cIt--;
         for (; cIt != it->beforeBegin(); --cIt) {
           size_t j = cIt.index();
           assert(j < i);
-          r[j] -= *cIt * x[i];
+          b[j] -= *cIt * x[i];
         }
       }
     } else {
@@ -75,9 +73,9 @@ namespace MatrixVector {
         for (; cIt != it->end(); ++cIt) {
           size_t j = cIt.index();
           assert(j > i);
-          r[i] -= *cIt * x[j];
+          b[i] -= *cIt * x[j];
         }
-        x[i] = r[i] / diagonal;
+        x[i] = b[i] / diagonal;
       }
     }
   }