getdp: it is now possible to compute only the application of a vector to the operator.

getdp: beyn.simple() computes now the *absolute* residual norm of an eigenpair.

cavity/maxwell.pro

@@ -22,7 +22,8 @@ Function{
 
   // Control data //
   DefineConstant[imposeRHS = 0,          // Should I use an imposed RHS?
-                 viewOnly  = 0,          // Should I use X for postpro only?
+                 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
 }
 
@@ -95,11 +96,16 @@ Resolution{
         CopyRightHandSide[b(), A];
       EndIf
 
-      If(!viewOnly)
+      If(!doPostpro && !doApply)
         Solve[A];
         CopySolution[A, x()];
       EndIf
-      If(viewOnly)
+      If(doApply)
+        CopySolution[x(), A];
+        Apply[A];
+        CopySolution[A, x()];
+      EndIf
+      If(doPostpro)
         CopySolution[x(), A];
       EndIf
     }

src/beyn.py

@@ -17,7 +17,8 @@ def simple(operator, origin, radius,
     rankTol -- the tolerance for the rank test (optional)
     verbose -- should I be verbose? (optional)
 
-    Returns the computed eigenvalues and eigenvectors
+    Returns the computed eigenvalues, eigenvectors
+    and the associate *absolute* residual norm
     """
 
     # Display the parameter used
@@ -71,14 +72,22 @@ def simple(operator, origin, radius,
     A1 = integrate(operator, myPath, 1, vHat)
     B  = V0.H * A1 * W0 * S0Inv
 
-    # Eigenvalues & eigenvectors
+    # Eigenvalues & eigenvectors (with projection onto V0)
     if(verbose): print "Solving linear EVP..."
     myLambda, QHat = numpy.linalg.eig(B)
     Q = V0 * QHat;
 
+    # Absolute residual norm
+    nFound = myLambda.shape[0]
+    norm   = np.empty((nFound))
+    for i in range(nFound):
+        operator.apply(Q[:, i], myLambda[i]) # Apply eigenpair
+        r = operator.solution()              # Get residual
+        norm[i] = np.linalg.norm(r)          # Save the norm
+
     # Done
     if(verbose): print "Done!"
-    return myLambda, Q
+    return myLambda, Q, norm
 
 
 ## Import only simple (other functions are just helpers)

src/cim.py

@@ -12,9 +12,9 @@ arg = args.parse()
 operator = Solver.GetDPWave(arg.pro, arg.mesh, arg.resolution)
 
 # Compute
-l, V = Beyn.simple(operator, arg.origin, arg.radius,
-                   arg.nodes, arg.maxIt, arg.lStart, arg.lStep, arg.rankTol,
-                   not arg.quiet)
+l, V, r = Beyn.simple(operator, arg.origin, arg.radius,
+                      arg.nodes, arg.maxIt, arg.lStart, arg.lStep, arg.rankTol,
+                      not arg.quiet)
 
 # Display the computed eigenvalues
 if(not arg.quiet):
@@ -26,7 +26,8 @@ for i in range(l.shape[0]):
     if(not arg.quiet): print " # " + str(i) + ": ",
     print "(" + format(np.real(l[i]).tolist(), '+e') + ")",
     print "+",
-    print "(" + format(np.imag(l[i]).tolist(), '+e') + ")*j"
+    print "(" + format(np.imag(l[i]).tolist(), '+e') + ")*j",
+    print "| " + format(r[i], 'e')
 
 if(not arg.quiet):
     print "----------- "

src/solver.py

@@ -12,7 +12,8 @@ class GetDPWave:
     b -- the right hand side (RHS) to use
     x -- the solution vector computed by GetDP
     imposeRHS -- should the imposed RHS be taken into account?
-    viewOnly -- should x be used for post-processing only
+    doPostpro -- should only a view be created for x?
+    doApply -- should only the application of x be computed?
     fileName -- post-processing file name
     """
 
@@ -43,6 +44,18 @@ class GetDPWave:
         """Returns the number of degrees of freedom"""
         return self.solution().shape[0]
 
+    def apply(self, x, w):
+        """Applies x to the operator with a pulsation of w
+
+        This method updates self.solution()
+        """
+        GetDP.GetDPSetNumber("x", self.__toGetDP(x))
+        GetDP.GetDPSetNumber("doApply", 1)
+        GetDP.GetDP(["getdp", self.__pro,
+                     "-msh",  self.__mesh,
+                     "-cal"])
+        GetDP.GetDPSetNumber("doApply", 0)
+
     def solve(self, b, w):
         """Solves with b as RHS and w as complex angular frequency"""
         GetDP.GetDPSetNumber("angularFreqRe", np.real(w).tolist())
@@ -65,13 +78,13 @@ class GetDPWave:
         This method generates a linear system
         """
         GetDP.GetDPSetNumber("x", self.__toGetDP(x))
-        GetDP.GetDPSetNumber("viewOnly", 1)
+        GetDP.GetDPSetNumber("doPostpro", 1)
         GetDP.GetDPSetString("fileName", fileName)
         GetDP.GetDP(["getdp", self.__pro,
                      "-msh",  self.__mesh,
                      "-cal",
                      "-pos",  posName])
-        GetDP.GetDPSetNumber("viewOnly", 0)
+        GetDP.GetDPSetNumber("doPostpro", 0)
 
     @staticmethod
     def __toNumpy(vGetDP):