Start to reuse LU in a GetDP way...

WARNING: change in the GetDP .pro file interface!

bin/beyn.py

@@ -101,7 +101,7 @@ def simple(operator, origin, radius,
     norm   = np.empty((nFound))
     for i in range(nFound):
         operator.apply(Q[:, i], myLambda[i]) # Apply eigenpair
-        r = operator.solution()              # Get residual
+        r = operator.solution(0)             # Get residual
         norm[i] = np.linalg.norm(r)          # Save the norm
 
     # Done
@@ -174,18 +174,15 @@ def multiSolve(solver, B, w):
     w -- the complex frequency to use
     """
 
-    # Number of solve
-    nSolve = B.shape[1]
-
-    # Initialise X
-    size = solver.size()
-    X    = np.matlib.empty((size, nSolve), dtype=complex)
+    # Solve
+    solver.solve(B, w)
 
-    # Loop and solve
+    # Get solutions
+    nSolve = B.shape[1]
+    size   = solver.size()
+    X      = np.matlib.empty((size, nSolve), dtype=complex)
     for i in range(nSolve):
-        b = B[:, i]
-        solver.solve(b, w)
-        X[:, i] = solver.solution()
+        X[:, i] = solver.solution(i)
 
     # Done
     return X

bin/solver.py

@@ -18,11 +18,13 @@ class GetDPWave:
     The .pro file should use the following variables:
     angularFreqRe -- the real part of the angular frequency to use
     angularFreqIm -- the imaginary part of the angular frequency to use
-    b -- the right hand side (RHS) to use
-    x -- the solution vector computed by GetDP
-    imposeRHS -- should the imposed RHS be taken into account?
-    doPostpro -- should only a view be created for x?
-    doApply -- should only the application of x be computed?
+    nRHS -- the number of right hand side that should be considered (default: 1)
+    b_i / b~{i} -- the ith right hand side to use (first index is 0)
+    x_i / x~{i} -- the ith solution vector computed by GetDP (first index is 0)
+    doInit -- should some initialization be done (default: 0)?
+    doSolve -- should Ax_i = b_i be solved for all i (default: 0)?
+    doPostpro -- should only a view be created for x (default: 0)?
+    doApply -- should only the application of x be computed (default: 0)?
     fileName -- post-processing file name
     """
 
@@ -44,19 +46,20 @@ class GetDPWave:
         self.__resolution = resolution
         self.__optional   = optional
 
-        # Generate DoFs and initialise RHS and solution vectors in GetDP
+        # Generate DoFs and assemble a first system
+        GetDP.GetDPSetNumber("doInit", 1);
         GetDP.GetDP(["getdp",     self.__pro,
                      "-msh",      self.__mesh,
                      "-solve",    self.__resolution,
                      "-v", "2"] + self.__optional)
 
-    def solution(self):
-        """Returns the solution"""
-        return self.__toNumpy(GetDP.GetDPGetNumber("x"))
+    def solution(self, i):
+        """Returns the solution for the ith RHS (first index is 0)"""
+        return self.__toNumpy(GetDP.GetDPGetNumber("x_" + str(i)))
 
     def size(self):
         """Returns the number of degrees of freedom"""
-        return self.solution().shape[0]
+        return self.solution(0).shape[0]
 
     def apply(self, x, w):
         """Applies x to the operator with a pulsation of w
@@ -65,23 +68,34 @@ class GetDPWave:
         """
         GetDP.GetDPSetNumber("angularFreqRe", np.real(w).tolist())
         GetDP.GetDPSetNumber("angularFreqIm", np.imag(w).tolist())
-        GetDP.GetDPSetNumber("x", self.__toGetDP(x))
+        GetDP.GetDPSetNumber("x_0", self.__toGetDP(x))
         GetDP.GetDPSetNumber("doApply", 1)
         GetDP.GetDP(["getdp",  self.__pro,
                      "-msh",   self.__mesh,
                      "-cal"] + self.__optional)
-        GetDP.GetDPSetNumber("doApply", 0)
 
     def solve(self, b, w):
-        """Solves with b as RHS and w as complex angular frequency"""
+        """Solves with b as RHS and w as complex angular frequency
+
+        b is a matrix and each column is a different RHS
+        """
+        # Number of RHS
+        nRHS = b.shape[1]
+
+        # Set variables
+        GetDP.GetDPSetNumber("nRHS", nRHS)
         GetDP.GetDPSetNumber("angularFreqRe", np.real(w).tolist())
         GetDP.GetDPSetNumber("angularFreqIm", np.imag(w).tolist())
-        GetDP.GetDPSetNumber("b", self.__toGetDP(b))
-        GetDP.GetDPSetNumber("imposeRHS", 1)
+
+        # Set RHS
+        for i in range(nRHS):
+            GetDP.GetDPSetNumber("b_" + str(i), self.__toGetDP(b[:, i]))
+
+        # Solve
+        GetDP.GetDPSetNumber("doSolve", 1)
         GetDP.GetDP(["getdp",  self.__pro,
                      "-msh",   self.__mesh,
                      "-cal"] + self.__optional)
-        GetDP.GetDPSetNumber("imposeRHS", 0)
 
     def view(self, x, fileName):
         """Generates a post-processing view
@@ -92,13 +106,12 @@ class GetDPWave:
 
         This method generates a linear system
         """
-        GetDP.GetDPSetNumber("x", self.__toGetDP(x))
+        GetDP.GetDPSetNumber("x_0", self.__toGetDP(x))
         GetDP.GetDPSetNumber("doPostpro", 1)
         GetDP.GetDPSetString("fileName", fileName)
         GetDP.GetDP(["getdp",  self.__pro,
                      "-msh",   self.__mesh,
                      "-cal"] + self.__optional)
-        GetDP.GetDPSetNumber("doPostpro", 0)
 
     @staticmethod
     def __toNumpy(vGetDP):

example/maxwell_linear/maxwell.pro

@@ -16,11 +16,15 @@ Function{
          angularFreqRe, angularFreqIm];
 
   // Algebraic data //
-  DefineConstant[x() = {},  // Solution
-                 b() = {}]; // Right hand side
+  DefineConstant[nRHS = 1];       // Number of RHS for this run
+  For i In {0:nRHS-1}
+    DefineConstant[x~{i}() = {},  // Solution
+                   b~{i}() = {}]; // Right hand side
+  EndFor
 
   // Control data //
-  DefineConstant[imposeRHS = 0,          // Should I use an imposed RHS?
+  DefineConstant[doInit    = 0,          // Should I initialize some stuff?
+                 doSolve   = 0,          // Should I solve Ax = b?
                  doPostpro = 0,          // Should I only create a view for x()?
                  doApply   = 0,          // Should I only apply x(): x <- Ax?
                  fileName  = "eig.pos"]; // Postpro file name
@@ -91,21 +95,32 @@ Resolution{
     Operation{
       Generate[A];
 
-      If(imposeRHS)
-        CopyRightHandSide[b(), A];
+      If(doInit)
+        CopySolution[A, x~{0}()];
       EndIf
 
-      If(!doPostpro && !doApply)
+      If(doSolve)
+        // Full solve for first RHS
+        CopyRightHandSide[b~{0}(), A];
         Solve[A];
-        CopySolution[A, x()];
+        CopySolution[A, x~{0}()];
+
+        // SolveAgain for other RHSs
+        For i In {1:nRHS-1}
+          CopyRightHandSide[b~{i}(), A];
+          SolveAgain[A];
+          CopySolution[A, x~{i}()];
+        EndFor
       EndIf
+
       If(doApply)
-        CopySolution[x(), A];
+        CopySolution[x~{0}(), A];
         Apply[A];
-        CopySolution[A, x()];
+        CopySolution[A, x~{0}()];
       EndIf
+
       If(doPostpro)
-        CopySolution[x(), A];
+        CopySolution[x~{0}(), A];
         PostOperation[Maxwell];
       EndIf
     }