Optimise DILU implementation

opm/simulators/linalg/DILU.hpp

@@ -67,33 +67,24 @@ class SeqDilu : public PreconditionerWithUpdate<X, Y>
         // we build the inverseD matrix
         Dinv_.resize(A_.N());
 
-        // is this the correct usage?
         update();
     }
 
     virtual void update() override {
-
+    OPM_TIMEBLOCK(update);
         using block = typename M::block_type;
-
+        for ( auto i = A_.begin(); i != A_.end(); ++i) {
+             Dinv_[i.index()] = A_[i.index()][i.index()];
+        }
         for ( auto i = A_.begin(); i != A_.end(); ++i)
         {
-            block Dinv_temp;
-            for (auto a_ij=i->begin(); a_ij != i->end(); ++a_ij) {
+            block Dinv_temp = Dinv_[i.index()];
+            for (auto a_ij=i->begin(); a_ij.index() < i.index(); ++a_ij) {
                 auto a_ji = A_[a_ij.index()].find(i.index());
-
-                if (a_ij.index() == i.index()) {
-                    Dinv_temp += A_[i.index()][i.index()];
-                }
-
-                // if A[i, j] != 0 and A[j, i] != 0:
-                else if (a_ji != A_[a_ij.index()].end() && (a_ij.index() < i.index())) {
-
-                    auto d_i = Dinv_[a_ij.index()];
-                    d_i.rightmultiply(*a_ji);
-                    d_i.leftmultiply(*a_ij);
-
-                    // d[j] -= A[j, i] * d[i] * A[i, j]
-                    Dinv_temp -= d_i;
+                // if A[i, j] != 0 and A[j, i] != 0
+                if (a_ji != A_[a_ij.index()].end()) {
+                    // Dinv_temp -= A[i, j] * d[j] * A[j, i]
+                    Dinv_temp -= (*a_ij)*Dinv_[a_ij.index()]*(*a_ji);
                 }
             }
             Dinv_temp.invert();
@@ -119,63 +110,49 @@ class SeqDilu : public PreconditionerWithUpdate<X, Y>
  virtual void apply(X& v, const Y& d) override
     {
 
-    // M = (D + L_A) D^-1 (D + U_A)   (a LU decomposition of M)
-    // where L_A and U_A are the strictly lower and upper parts of A and M has the properties:
-    // diag(A) = diag(M)
-    // solving the product M^-1(b-Ax) using upper and lower triangular solve
-    // z = x_k+1 - x_k = M^-1(b - Ax) = (D + L_A)^-1 D (D + U_A)^-1 (b - Ax)
-
-    typedef typename Y::block_type dblock;
-
-    // copy current v
-    X x(v);
-    X z;
-    Y y;
-    z.resize(A_.N());
-    y.resize(A_.N());
-
-    // lower triangular solve: (D + L) y = b - Ax
-    auto endi=A_.end();
-    for (auto i=A_.begin(); i != endi; ++i)
-      {
-        dblock rhsValue(d[i.index()]);
-        auto&& rhs = Impl::asVector(rhsValue);
-        for (auto j=(*i).begin(); j != i->end(); ++j) {
-            // if  A[i][j] != 0
-            // rhs -= A[i][j]* x[j];
-            Impl::asMatrix(*j).mmv(Impl::asVector(x[j.index()]), rhs);
-
-            // rhs -= A[i][j]* y[j];
-            if (j.index() < i.index()) {
-                Impl::asMatrix(*j).mmv(Impl::asVector(y[j.index()]), rhs);
+        // M = (D + L_A) D^-1 (D + U_A)   (a LU decomposition of M)
+        // where L_A and U_A are the strictly lower and upper parts of A and M has the properties:
+        // diag(A) = diag(M)
+        // solving the product M^-1(b-Ax) using upper and lower triangular solve
+        // z = x_k+1 - x_k = M^-1(b - Ax) = (D + L_A)^-1 D (D + U_A)^-1 (b - Ax)
+        OPM_TIMEBLOCK(apply);
+        typedef typename Y::block_type dblock;
+        typedef typename X::block_type vblock;
+
+        // lower triangular solve: (D + L) y = b - Ax
+        auto endi=A_.end();
+        for (auto i=A_.begin(); i != endi; ++i)
+            {
+            dblock rhsValue(d[i.index()]);
+            auto&& rhs = Impl::asVector(rhsValue);
+            for (auto j=(*i).begin(); j.index()<i.index(); ++j) {
+                // if  A[i][j] != 0
+                // rhs -= A[i][j]* y[j], where v_j stores y_i
+                j->mmv(v[j.index()], rhs);
             }
+            // y_i = Dinv_i * rhs
+            // storing y_i in v_i
+            Dinv_[i.index()].mv(rhs, v[i.index()]);  // (D + L)_ii = D_i
         }
-        // y = Dinv * rhs
-        Impl::asMatrix(Dinv_[i.index()]).mv(rhs, y[i.index()]);
-    }
-
 
-    // upper triangular solve: (D + U) z = Dy 
-    auto rendi=A_.beforeBegin();
-    for (auto i=A_.beforeEnd(); i!=rendi; --i)
-    {
-        // rhs = 0
-        dblock rhs;
 
-        for (auto j=(*i).beforeEnd(); j.index()>i.index(); --j) {
-            // if A [i][j] != 0
-            //rhs += A[i][j]*z[j]
-            Impl::asMatrix(*j).umv(Impl::asVector(z[j.index()]), rhs);
+        // upper triangular solve: (D + U) z = Dy 
+        auto rendi=A_.beforeBegin();
+        for (auto i=A_.beforeEnd(); i!=rendi; --i)
+        {
+            // rhs = 0
+            vblock rhs;
 
+            for (auto j=(*i).beforeEnd(); j.index()>i.index(); --j) {
+                // if A[i][j] != 0
+                //rhs += A[i][j]*z[j], where v_j stores z_j
+                j->umv(v[j.index()], rhs);
+            }
+            // calculate update z = M^-1(b - Ax)
+            // z_i = y_i - Dinv_i*rhs
+            // before update v_i is y_i and after update v_i is z_i
+            Dinv_[i.index()].mmv(rhs, v[i.index()]);
         }
-        // calculate update z = M^-1(b - Ax)
-        // z_i = y_i - Dinv*temp
-        dblock temp;
-        Impl::asMatrix(Dinv_[i.index()]).mv(rhs, temp); 
-        z[i.index()] = y[i.index()] - temp;
-    }
-    // update v
-    v += z;
     }
 
     /*!