Add the invertWithUMFPack() utility.

opm/simulators/wells/MSWellHelpers.hpp

@@ -72,7 +72,50 @@ namespace mswellhelpers
 #else
         // this is not thread safe
         OPM_THROW(std::runtime_error, "Cannot use applyUMFPack() without UMFPACK. "
-                  "Reconfigure opm-simulator with SuiteSparse/UMFPACK support and recompile.");
+                  "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile.");
+#endif // HAVE_UMFPACK
+    }
+
+
+
+    /// Applies umfpack and checks for singularity
+    template <typename MatrixType, typename VectorType>
+    MatrixType
+    invertWithUMFPack(const MatrixType& D, std::shared_ptr<Dune::UMFPack<MatrixType> >& linsolver)
+    {
+#if HAVE_UMFPACK
+        const int sz = D.M();
+        const int bsz = D[0][0].M();
+        VectorType e(sz);
+        e = 0.0;
+
+        // Make a block matrix with a full pattern.
+        MatrixType inv(sz, sz, sz*sz, MatrixType::row_wise);
+        for (auto row = inv.createbegin(); row != inv.createend(); ++row) {
+            for (int c = 0; c < sz; ++c) {
+                row.insert(c);
+            }
+        }
+
+        // Create inverse by passing basis vectors to the solver.
+        for (int ii = 0; ii < sz; ++ii) {
+            for (int jj = 0; jj < bsz; ++jj) {
+                e[ii][jj] = 1.0;
+                auto col = applyUMFPack(D, linsolver, e);
+                for (int cc = 0; cc < sz; ++cc) {
+                    for (int dd = 0; dd < bsz; ++dd) {
+                        inv[cc][ii][dd][jj] = col[cc][dd];
+                    }
+                }
+                e[ii][jj] = 0.0;
+            }
+        }
+
+        return inv;
+#else
+        // this is not thread safe
+        OPM_THROW(std::runtime_error, "Cannot use invertWithUMFPack() without UMFPACK. "
+                  "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile.");
 #endif // HAVE_UMFPACK
     }