Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/rrompy/linalg/linalg_krylov_engine.py
c11471199317develop

Last Change
Raw File

linalg_krylov_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.linalg_base_engine import LinAlgBaseEngine
from rrompy.utilities.base.types import Np1D, Np2D
__all__ = ['LinAlgKrylovEngine']
class LinAlgKrylovEngine(LinAlgBaseEngine):
"""HERE"""
def resetHistory(self):
self.samples = None
self.nsamples = 0
def preprocesssamples(self):
if self.samples is None or self.nAs() <= 1: return
if self.nAs() > self.nsamples:
return self.samples[:, : self.nsamples]
return self.samples[:, self.nsamples + 1 - self.nAs() : self.nsamples]
def preprocessb(self, bs : Np1D = None, overwrite : bool = False):
if bs is None: bs = self.bs
if self.nsamples < len(bs):
return copy(bs[self.nsamples])
return np.zeros_like(bs[0])
def postprocessu(self, u:Np1D, overwrite : bool = False):
return u
def preallocateSamples(self, u:Np1D, n:int):
self.samples = np.empty((u.size, n), dtype = u.dtype)
self.samples[:, 0] = u
def nextSample(self, mu:complex, overwrite : bool = False) -> Np1D:
if not np.isclose(mu, self.mu0):
self.shiftCenter(mu, update = True)
samplesOld = self.preprocesssamples()
RHS = self.preprocessb(overwrite = overwrite)
for i in range(1, min(self.nAs(), self.nsamples + 1)):
RHS -= self.As[i].dot(samplesOld[:, - i])
u = self.postprocessu(self.solveLS(mu, bs = [RHS]),
overwrite = overwrite)
if overwrite:
self.samples[:, self.nsamples] = u
else:
if self.nsamples == 0:
self.samples = u[:, None]
else:
self.samples = np.hstack((self.samples, u[:, None]))
self.nsamples += 1
return u
def iterSample(self, mu:complex, n:int) -> Np2D:
if n <= 0:
raise Exception(("Number of Krylov iterations must be positive."))
self.resetHistory()
u = self.nextSample(mu)
if n > 1:
self.preallocateSamples(u, n)
for _ in range(1, n):
self.nextSample(mu, overwrite = True)
return self.samples
c4science · Help