Page MenuHomec4science
Files F102565282

linalg_base_engine.py
No OneTemporary
Actions

File Metadata

Created
Sat, Feb 22, 01:26

linalg_base_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/>.
#
import numpy as np
from copy import copy
import scipy.sparse.linalg as spla
from rrompy.utilities.base.types import Np1D, Np2D, List, Tuple
__all__ = ['LinAlgBaseEngine']
class LinAlgBaseEngine:
"""HERE"""
def __init__(self, As:List[Np2D], bs:List[Np1D], mu0 : complex = 0.):
self.As = As
self.bs = bs
self.mu0 = mu0
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
def resetHistory(self):
pass
@property
def As(self):
"""Value of As. Its assignment resets history."""
return self._As
@As.setter
def As(self, As):
self._As = As
self.resetHistory()
@property
def bs(self):
"""Value of bs. Its assignment resets history."""
return self._bs
@bs.setter
def bs(self, bs):
self._bs = bs
self.resetHistory()
def nAs(self):
"""Number of terms in As."""
return len(self._As)
def nbs(self):
"""Number of terms in bs."""
return len(self._bs)
def shiftCenter(self, mu:complex, As : List[Np2D] = None,
bs : List[Np1D] = None, update : bool = False,
which : Np1D = None)-> Tuple[List[Np2D], List[Np1D]]:
"""
Build linear system from polynomial blocks.
Args:
mu: Parameter value.
As: Matrix blocks. If None, set to self.As. Defaults to None.
bs: RHS blocks. If None, set to self.bs. Defaults to None.
which: What terms to compute. If None, all are computed. Defaults
to None.
Returns:
Shifted matrix and RHS coefficients of system.
"""
if update and not (As is None and bs is None and which is None):
raise Exception(("Default values for As, bs, and which are "
"required for update."))
if As is None: As = self.As
if bs is None: bs = self.bs
lAs, lbs = len(As), len(bs)
maxL = max(lAs, lbs)
if which is None: which = np.arange(maxL)
mu0 = self.mu0
if np.isclose(mu, mu0):
return ([copy(As[term]) for term in which],
[copy(bs[term]) for term in which])
AsNew = []
bsNew = []
for term in which:
rho = 1.
if term < lAs:
AsNew += [copy(As[term])]
if term < lbs:
bsNew += [copy(bs[term])]
for j in range(1, maxL - term):
rho *= (term + j) / j * (mu - mu0)
if j < lAs:
AsNew[-1] += rho * As[term + j]
if j < lbs:
bsNew[-1] += rho * bs[term + j]
if update:
self.As = AsNew
self.bs = bsNew
self.mu0 = mu
return AsNew, bsNew
def buildLS(self, mu:complex, As : List[Np2D] = None,
bs : List[Np1D] = None) -> Tuple[Np2D, Np1D]:
"""
Build linear system from polynomial blocks.
Args:
mu: Parameter value.
As: Matrix blocks. If None, set to self.As. Defaults to None.
bs: RHS blocks. If None, set to self.bs. Defaults to None.
Returns:
Matrix and RHS of system.
"""
As, bs = self.shiftCenter(mu, As, bs, which = [0])
return As[0], bs[0]
def solveLS(self, mu:complex, As : List[Np2D] = None,
bs : List[Np1D] = None) -> Np1D:
"""
Solve linear system.
Args:
mu: Parameter value.
Returns:
Solution of system.
"""
A, b = self.buildLS(mu, As, bs)
u = spla.spsolve(A, b)
return u

Event Timeline

c4science · Help