Page MenuHomec4science
Files F90757029

linalg_arnoldi_engine.py
No OneTemporary
Actions

File Metadata

Created
Mon, Nov 4, 11:49

linalg_arnoldi_engine.py
View Options

# 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 <http://www.gnu.org/licenses/>.
#
from copy import copy
import numpy as np
from rrompy.linalg.pod_engine import PODEngine
from rrompy.linalg.linalg_krylov_engine import LinAlgKrylovEngine
from rrompy.utilities.base.types import Np1D, Np2D, List
__all__ = ['LinAlgArnoldiEngine']
class LinAlgArnoldiEngine(LinAlgKrylovEngine):
"""HERE"""
def __init__(self, As:List[Np2D], bs:List[Np1D], energyNormMatrix:Np2D,
mu0 : complex = 0.):
super().__init__(As, bs, mu0)
self.energyNormMatrix = energyNormMatrix
def resetHistory(self):
super().resetHistory()
self.HArnoldi = None
self.RArnoldi = None
self.RHSs = None
self.samplesAug = None
@property
def energyNormMatrix(self):
"""Value of energyNormMatrix. Its assignment resets history."""
return self._energyNormMatrix
@energyNormMatrix.setter
def energyNormMatrix(self, energyNormMatrix):
self._energyNormMatrix = energyNormMatrix
self.PODEngine = PODEngine(self.energyNormMatrix)
def preprocesssamples(self):
ns = self.nsamples
if ns <= 0 or self.nAs() <= 1: return
return self.samplesAug[:, ns - 1].reshape((self.nAs() - 1, - 1)).T
def preprocessb(self, bs : Np1D = None, overwrite : bool = False):
if bs is None: bs = self.bs
ns = self.nsamples
if ns < len(bs):
r = copy(bs[ns])
else:
r = np.zeros_like(bs[0])
if min(len(bs), ns + 1) > 1:
if ns == 1:
r = r / self.RArnoldi[0, 0]
else:
r = ((r - self.RHSs[:, :ns-1].dot(self.RArnoldi[:ns-1, ns-1]))
/ self.RArnoldi[ns-1, ns-1])
if overwrite:
self.RHSs[:, ns - 1] = r
else:
if ns == 1:
self.RHSs = r.reshape((- 1, 1))
else:
self.RHSs = np.hstack((self.RHSs, r[:, None]))
return r
def postprocessu(self, u:Np1D, overwrite : bool = False):
ns = self.nsamples
if ns == 0:
u, h, uAug = self.PODEngine.GS(u, np.empty((0, 0)))
r = h[0]
uAug = np.pad(u, ((self.nAs() - 2) * u.size, 0), "constant")
else:
uAug = np.concatenate((self.samplesAug[u.size :, ns - 1], u),
axis = None)
u, h, uAug = self.PODEngine.GS(u, self.samples[:, : ns], ns, uAug,
self.samplesAug[:, : ns])
if overwrite:
self.HArnoldi[: ns + 1, ns] = h
if ns > 0:
r = self.HArnoldi[: ns + 1, 1 : ns + 1].dot(
self.RArnoldi[: ns, ns - 1])
self.RArnoldi[: ns + 1, ns] = r
self.samplesAug[:, ns] = uAug
else:
if self.nsamples == 0:
self.HArnoldi = h.reshape((1, 1))
self.RArnoldi = r.reshape((1, 1))
self.samplesAug = uAug.reshape((-1, 1))
else:
self.HArnoldi=np.block([[self.HArnoldi, h[:-1, None]],
[np.zeros((1, self.nsamples)), h[-1]]])
if ns > 0:
r = self.HArnoldi[: ns + 1, 1 : ns + 1].dot(
self.RArnoldi[: ns, ns - 1])
self.RArnoldi=np.block([[self.RArnoldi, r[:-1, None]],
[np.zeros((1, self.nsamples)), r[-1]]])
self.samplesAug = np.hstack((self.samplesAug, uAug[:, None]))
return u
def preallocateSamples(self, u:Np1D, n:int):
super().preallocateSamples(u, n)
h = self.HArnoldi
r = self.RArnoldi
saug = self.samplesAug
self.HArnoldi = np.zeros((n, n), dtype = u.dtype)
self.HArnoldi[0, 0] = h[0, 0]
self.RArnoldi = np.zeros((n, n), dtype = u.dtype)
self.RArnoldi[0, 0] = r[0, 0]
self.RHSs = np.empty((u.size, n - 1), dtype = u.dtype)
self.samplesAug = np.empty((saug.shape[0], n), dtype = u.dtype)
self.samplesAug[:, 0] = saug[:, 0]

Event Timeline

c4science · Help