diff --git a/rrompy/hfengines/base/hfengine_base.py b/rrompy/hfengines/base/hfengine_base.py
index 93d1958..20eeee1 100644
--- a/rrompy/hfengines/base/hfengine_base.py
+++ b/rrompy/hfengines/base/hfengine_base.py
@@ -1,296 +1,295 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from abc import abstractmethod
import numpy as np
import scipy.sparse as scsp
from numbers import Number
from copy import copy as softcopy, deepcopy as copy
from rrompy.utilities.base.decorators import nonaffine_construct
-from rrompy.utilities.base.types import (Np1D, Np2D, List, DictAny, paramVal,
- paramList, sampList)
+from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, ListAny,
+ DictAny, paramVal, paramList,
+ sampList)
from rrompy.utilities.numerical import solve as tsolve, dot, customPInv
from rrompy.utilities.numerical.rayleigh_quotient_iteration import (
rayleighQuotientIteration)
from rrompy.utilities.exception_manager import RROMPyAssert
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.solver.linear_solver import setupSolver
try:
from mpi4py import MPI
_MPI_IS_LOADED = True
- from rrompy.utilities.base.partitioning import partition
except ImportError:
_MPI_IS_LOADED = False
__all__ = ['HFEngineBase']
class HFEngineBase:
"""Generic solver for parametric problems."""
_forceSerial = False
def __init__(self, verbosity : int = 10, timestamp : bool = True):
self.verbosity = verbosity
self.timestamp = timestamp
self.setSolver("SPSOLVE", {"use_umfpack" : False})
self.npar = 0
self._C = None
self.outputNormMatrix = 1.
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
def __dir_base__(self):
return [x for x in self.__dir__() if x[:2] != "__"]
def __deepcopy__(self, memo):
return softcopy(self)
@property
def npar(self):
"""Value of npar."""
return self._npar
@npar.setter
def npar(self, npar):
nparOld = self._npar if hasattr(self, "_npar") else -1
if npar != nparOld:
self.rescalingExp = [1.] * npar
self._npar = npar
@property
def spacedim(self):
return 1
def checkParameter(self, mu:paramVal):
muP = checkParameter(mu, self.npar)
if self.npar == 0: muP.reset((1, 0), muP.dtype)
return muP
def checkParameterList(self, mu:paramList):
muL = checkParameterList(mu, self.npar)
if self.npar == 0: muL[0].reset((1, 0), muL[0].dtype)
return muL
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
self.energyNormMatrix = 1.
def buildEnergyNormDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
self.energyNormDualMatrix = 1.
def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False,
dual : bool = False, is_state : bool = True) -> Np2D:
"""Scalar product."""
if is_state or self.isCEye:
if dual:
if not hasattr(self, "energyNormDualMatrix"):
self.buildEnergyNormDualForm()
energyMat = self.energyNormDualMatrix
else:
if not hasattr(self, "energyNormMatrix"):
self.buildEnergyNormForm()
energyMat = self.energyNormMatrix
else:
energyMat = self.outputNormMatrix
if not isinstance(u, (np.ndarray,)): u = u.data
if not isinstance(v, (np.ndarray,)): v = v.data
if onlyDiag:
return np.sum(dot(energyMat, u) * v.conj(), axis = 0)
return dot(dot(energyMat, u).T, v.conj()).T
def norm(self, u:Np2D, dual : bool = False,
is_state : bool = True) -> Np1D:
return np.abs(self.innerProduct(u, u, onlyDiag = True, dual = dual,
is_state = is_state)) ** .5
def baselineA(self):
"""Return 0 of shape consistent with operator of linear system."""
if (hasattr(self, "As") and hasattr(self.As, "__len__")
and self.As[0] is not None):
d = self.As[0].shape[0]
else:
d = self.spacedim
return scsp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(d + 1)),
shape = (d, d), dtype = np.complex)
def baselineb(self):
"""Return 0 of shape consistent with RHS of linear system."""
return np.zeros(self.spacedim, dtype = np.complex)
@nonaffine_construct
@abstractmethod
def A(self, mu : paramVal = [], der : List[int] = 0) -> Np2D:
"""
Assemble terms of operator of linear system and return it (or its
derivative) at a given parameter.
"""
return
@nonaffine_construct
@abstractmethod
def b(self, mu : paramVal = [], der : List[int] = 0) -> Np1D:
"""
Assemble terms of RHS of linear system and return it (or its
derivative) at a given parameter.
"""
return
@property
def C(self):
"""Value of C."""
if self._C is None: self._C = 1.
return self._C
@property
def isCEye(self):
return isinstance(self.C, Number)
def applyC(self, u:sampList):
"""Apply LHS of linear system."""
return dot(self.C, u)
def applyCpInv(self, u:sampList):
"""Apply pseudoinverse of LHS of linear system."""
return dot(customPInv(self.C), u)
def setSolver(self, solverType:str, solverArgs : DictAny = {}):
"""Choose solver type and parameters."""
self._solver, self._solverArgs = setupSolver(solverType, solverArgs)
def solve(self, mu : paramList = [], RHS : sampList = None,
return_state : bool = False) -> sampList:
"""
Find solution of linear system.
Args:
mu: parameter value.
RHS: RHS of linear system. If None, defaults to that of parametric
system. Defaults to None.
return_state: whether to return state before multiplication by c.
Defaults to False.
"""
from rrompy.sampling import sampleList, emptySampleList
+ from rrompy.utilities.base.parallel import (partitionScatter,
+ sampleGather)
if mu == []: mu = self.mu0
mu = self.checkParameterList(mu)[0]
- if len(mu) == 0: return emptySampleList()
- if _MPI_IS_LOADED and not self._forceSerial:
- comm = MPI.COMM_WORLD
- crank = comm.Get_rank()
- if crank == 0:
- mupart, idxpart = partition(mu, comm.Get_size(), True)
- else:
- mupart, idxpart = None, None
- mu0 = self.checkParameterList(mu)[0]
- mu = self.checkParameterList(comm.scatter(mupart, root = 0))[0]
- idxpart = comm.bcast(idxpart, root = 0)
+ nmu = len(mu)
+ if nmu == 0: return emptySampleList()
+ mu, idx = partitionScatter(mu, self._forceSerial)
+ mu = self.checkParameterList(mu)[0]
if len(mu) == 0:
sol = emptySampleList()
else:
if RHS is None:
RHS = sampleList([self.b(m) for m in mu])
else:
RHS = sampleList(RHS)
- if len(RHS) > 1: RHS = sampleList([RHS[x] for x in idxpart])
+ if len(RHS) > 1: RHS = sampleList([RHS[i] for i in idx])
mult = 0 if len(RHS) == 1 else 1
RROMPyAssert(mult * (len(mu) - 1) + 1, len(RHS), "Sample size")
for j in range(len(mu)):
u = tsolve(self.A(mu[j]), RHS[mult * j], self._solver,
self._solverArgs)
if j == 0:
sol = emptySampleList()
sol.reset((len(u), len(mu)), dtype = u.dtype)
sol[j] = u
if not return_state:
sol = sampleList(self.applyC(sol.data))
- if _MPI_IS_LOADED and not self._forceSerial:
- solsdata = comm.allgather(sol.data.T)
- sol = emptySampleList()
- sol.reset((solsdata[0].shape[1], len(mu0)),
- dtype = solsdata[0].dtype)
- for idx, dat in zip(idxpart, solsdata):
- if len(idx) > 0: sol[idx] = dat
- return sol
+ return sampleGather(sol, nmu, self._forceSerial)
def residual(self, mu : paramList = [], u : sampList = None,
post_c : bool = True) -> sampList:
"""
Find residual of linear system for given approximate solution.
Args:
mu: parameter value.
u: numpy complex array with function dofs. If None, set to 0.
post_c: whether to post-process using c. Defaults to True.
"""
from rrompy.sampling import sampleList, emptySampleList
+ from rrompy.utilities.base.parallel import (partitionScatter,
+ sampleGather)
if mu == []: mu = self.mu0
mu = self.checkParameterList(mu)[0]
- if len(mu) == 0: return emptySampleList()
- v = sampleList(np.zeros((self.spacedim, len(mu))))
- if u is not None: v = v + sampleList(u)
- for j in range(len(mu)):
- r = self.b(mu[j]) - dot(self.A(mu[j]), v[j])
- if post_c:
- if j == 0: res = np.empty((len(r), len(mu)), dtype = r.dtype)
- res[:, j] = r
- else:
- if j == 0:
- res = emptySampleList()
- res.reset((len(r), len(mu)), dtype = r.dtype)
- res[j] = r
- if post_c: res = sampleList(self.applyC(res))
- return res
+ nmu = len(mu)
+ if nmu == 0: return emptySampleList()
+ mu, idx = partitionScatter(mu, self._forceSerial)
+ mu = self.checkParameterList(mu)[0]
+ if len(mu) == 0:
+ res = emptySampleList()
+ else:
+ v = sampleList(np.zeros((self.spacedim, len(mu))))
+ if u is not None:
+ u = sampleList(u)
+ v = v + sampleList([u[i] for i in idx])
+ for j in range(len(mu)):
+ r = self.b(mu[j]) - dot(self.A(mu[j]), v[j])
+ if post_c:
+ if j == 0: res = np.empty((len(r), len(mu)),
+ dtype = r.dtype)
+ res[:, j] = r
+ else:
+ if j == 0:
+ res = emptySampleList()
+ res.reset((len(r), len(mu)), dtype = r.dtype)
+ res[j] = r
+ if post_c: res = sampleList(self.applyC(res))
+ return sampleGather(res, nmu, self._forceSerial)
def stabilityFactor(self, mu : paramList = [], u : sampList = None,
nIterP : int = 10, nIterR : int = 10) -> sampList:
"""
Find stability factor of matrix of linear system using iterative
inverse power iteration- and Rayleigh quotient-based procedure.
Args:
mu: parameter values.
u: numpy complex arrays with function dofs.
nIterP: number of iterations of power method.
nIterR: number of iterations of Rayleigh quotient method.
"""
from rrompy.sampling import sampleList
if mu == []: mu = self.mu0
mu = self.checkParameterList(mu)[0]
u = sampleList(u)
if len(u) != len(mu):
u = sampleList(np.tile(u.data.reshape(-1, 1), (1, len(mu))))
stabFact = np.empty(len(mu), dtype = float)
if not hasattr(self, "energyNormMatrix"):
self.buildEnergyNormForm()
for j in range(len(mu)):
stabFact[j] = rayleighQuotientIteration(self.A(mu[j]), u[j],
self.energyNormMatrix,
0., nIterP, nIterR,
self._solver,
self._solverArgs)
return stabFact
diff --git a/rrompy/utilities/base/parallel.py b/rrompy/utilities/base/parallel.py
new file mode 100644
index 0000000..5f6585b
--- /dev/null
+++ b/rrompy/utilities/base/parallel.py
@@ -0,0 +1,77 @@
+# Copyright (C) 2018 by the RROMPy authors
+#
+# This file is part of RROMPy.
+#
+# RROMPy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# RROMPy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with RROMPy. If not, see .
+#
+
+import numpy as np
+from rrompy.utilities.base.types import Tuple, ListAny, paramList, sampList
+from rrompy.sampling import emptySampleList
+try:
+ from mpi4py import MPI
+ _MPI_IS_LOADED = True
+except ImportError:
+ _MPI_IS_LOADED = False
+
+__all__ = ['partitionIdx', 'partition', 'partitionScatter', 'sampleGather']
+
+def partitionIdx(a:int, n:int):
+ x, idxs = 0, []
+ size = 0
+ for j in range(n):
+ if size * (n - j) != a - x: size = int(np.ceil((a - x) / (n - j)))
+ idxs += [list(range(x, x + size))]
+ x += size
+ return idxs
+
+def partition(obj:ListAny, n:int, return_idxs : bool = False):
+ idxs = partitionIdx(len(obj), n)
+ try:
+ res = [obj[idx] for idx in idxs]
+ if return_idxs: return res, idxs
+ return res
+ except:
+ res = [[obj[i] for i in idx] for idx in idxs]
+ if return_idxs: return res, idxs
+ return res
+
+def partitionScatter(mu:paramList, force_serial : bool = False) -> Tuple[
+ paramList, ListAny]:
+ """Split parameters among parallel cores."""
+ if not force_serial and _MPI_IS_LOADED and MPI.COMM_WORLD.Get_size() > 1:
+ comm = MPI.COMM_WORLD
+ muidx = None
+ if comm.Get_rank() == 0:
+ mupart, idxpart = partition(mu, comm.Get_size(), True)
+ muidx = zip(mupart, idxpart)
+ return comm.scatter(muidx, root = 0)
+ return mu, list(range(len(mu)))
+
+def sampleGather(obj:sampList, nmu:int,
+ force_serial : bool = False) -> sampList:
+ """Gather samples from parallel cores."""
+ if not force_serial and _MPI_IS_LOADED and MPI.COMM_WORLD.Get_size() > 1:
+ objList = MPI.COMM_WORLD.allgather(obj.data.T)
+ res = emptySampleList()
+ res.reset((objList[0].shape[1], nmu), dtype = objList[0].dtype)
+ idx = 0
+ for dat in objList:
+ ldat = len(dat)
+ if ldat > 0:
+ res[list(range(idx, idx + ldat))] = dat
+ idx += ldat
+ return res
+ return obj
+
diff --git a/rrompy/utilities/base/partitioning.py b/rrompy/utilities/base/partitioning.py
deleted file mode 100644
index 05bd0cf..0000000
--- a/rrompy/utilities/base/partitioning.py
+++ /dev/null
@@ -1,42 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import numpy as np
-from rrompy.utilities.base.types import ListAny
-
-__all__ = ['partitionIdx', 'partition']
-
-def partitionIdx(a:int, n:int):
- x, idxs = 0, []
- size = 0
- for j in range(n):
- if size * (n - j) != a - x: size = int(np.ceil((a - x) / (n - j)))
- idxs += [list(range(x, x + size))]
- x += size
- return idxs
-
-def partition(obj:ListAny, n:int, return_idxs : bool = False):
- idxs = partitionIdx(len(obj), n)
- try:
- res = [obj[idx] for idx in idxs]
- if return_idxs: return res, idxs
- return res
- except:
- res = [[obj[i] for i in idx] for idx in idxs]
- if return_idxs: return res, idxs
- return res