From 8ec7621faa610f647f599739905d23f3254b74b0 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <marsic@temf.tu-darmstadt.de>
Date: Mon, 3 Jul 2017 16:28:10 +0200
Subject: [PATCH] Start to reuse LU in a GetDP way... WARNING: change in the
GetDP .pro file interface!
---
bin/beyn.py | 19 +++++-------
bin/solver.py | 49 +++++++++++++++++++-----------
example/maxwell_linear/maxwell.pro | 35 +++++++++++++++------
3 files changed, 64 insertions(+), 39 deletions(-)
diff --git a/bin/beyn.py b/bin/beyn.py
index c2729a2..ffe0426 100644
--- a/bin/beyn.py
+++ b/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
diff --git a/bin/solver.py b/bin/solver.py
index ed8e82d..e890314 100644
--- a/bin/solver.py
+++ b/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):
diff --git a/example/maxwell_linear/maxwell.pro b/example/maxwell_linear/maxwell.pro
index 0b7a0f2..1e219c8 100644
--- a/example/maxwell_linear/maxwell.pro
+++ b/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
}
--
GitLab