*** Wartungsfenster jeden ersten Mittwoch vormittag im Monat ***

Skip to content
Snippets Groups Projects
Commit 1d449bc9 authored by Innerberger, Michael's avatar Innerberger, Michael
Browse files

Merge branch 'refactor/solver-interface' into feature-iterative-solvers

parents f79d4b6d 15ccb501
Branches feature-iterative-solvers
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 461 additions and 999 deletions
examples/algebraicSolver.m
+ 1
− 1
View file @ 1d449bc9
......@@ -33,7 +33,7 @@ for p = 1:pmax
% choose the iterative solver of choice: multigrid (with the variants
% lowOrderVcycle and highOrderVcycle), pcg (with jacobi and additive
% Schwarz preconditioner), and the cg solver
[solver, P] = chooseIterativeSolver(fes, blf, 'multigrid', 'lowOrderVcycle');
[solver, P] = IterativeSolver.choose(fes, blf, 'multigrid', 'lowOrderVcycle');
solver.tol = 1e-6;
solver.maxIter = 100;
......
examples/goafemIterativeSolver.m
+ 1
− 1
View file @ 1d449bc9
......@@ -37,7 +37,7 @@ for p = 1:pmax
lfG.qrf = QuadratureRule.ofOrder(2*p);
%% set up solver and lifting operator for nested iteration
[solver, P] = chooseIterativeSolver(fes, blf, 'multigrid', 'lowOrderVcycle');
[solver, P] = IterativeSolver.choose(fes, blf, 'multigrid', 'lowOrderVcycle');
solver.tol = 1e-6;
solver.maxIter = 1000;
......
lib/solvers/IterativeSolver.m → lib/solvers/@IterativeSolver/IterativeSolver.m
+ 9
− 0
View file @ 1d449bc9
......@@ -25,6 +25,10 @@
%
% tf = solver.isConverged(solver) returns state of convergence of the solver
% for each column of the right-hand side.
%
% [solver, P] = IterativeSolver.choose(fes, blf, class, variant) convenience
% factory method. Returns instances of IterativeSolver and Prolongation
% accoding to the given parameters.
classdef IterativeSolver < handle
......@@ -102,6 +106,11 @@ classdef IterativeSolver < handle
obj.activeComponents = obj.activeComponents & ~tf;
end
end
%% convenience factory function
methods (Static, Access=public)
[solver, P] = choose(fes, blf, class, variant)
end
%% validation functions to use within this class
methods (Static, Access=protected)
......
lib/solvers/chooseIterativeSolver.m → lib/solvers/@IterativeSolver/choose.m
+ 91
− 0
View file @ 1d449bc9
% chooseIterativeSolver Return suitable solver (instance of subclass of
% IterativeSolver) and suitable Prolongation P.
% choose Return suitable solver (instance of subclass of IterativeSolver) and
% suitable Prolongation P.
%
% solver = choose(fes, blf, class) returns a solver of class 'class' for given
% finite element space fes and bilinear form blf. For each class, sensible
% default variants are chosen.
%
% solver = choose(fes, blf, class, variant) returns a solver of class 'class'
% of specific variant 'variant'.
%
% [solver, P] = choose(__) additionally to the solver, returns a suitable
% Prolongation P, that can be used to prolongate the solution of the
% solver to the next finer mesh.
function [solver, P] = chooseIterativeSolver(fes, blf, class, variant)
function [solver, P] = choose(fes, blf, class, variant)
arguments
fes FeSpace
blf BilinearForm
......@@ -18,39 +29,55 @@ P = Prolongation.chooseFor(fes);
switch class
% non-preconditioned CG
case "cg"
solver = CgSolver();
solver = PcgSolver(NoPreconditioner());
% preconditioned CG family
case "pcg"
switch variant
case "iChol"
solver = ICholPcgSolver(fes);
preconditioner = IncompleteCholesky();
case "jacobi"
if order == 1, solver = JacobiPcgSolver(fes);
else, solver = BlockJacobiPcgSolver(fes, blf); end
if order == 1
preconditioner = DiagonalJacobi();
else
preconditioner = BlockJacobi(fes, blf);
end
case {"", "additiveSchwarzLowOrder"}
if order == 1, solver = LowestOrderAdditiveSchwarzPcg(fes, blf, P);
if order == 1
preconditioner = OptimalMLAdditiveSchwarz(P1JacobiCascade(fes, blf, P));
else
solver = AdditiveSchwarzLowOrderPcg(fes, blf);
preconditioner = OptimalMLAdditiveSchwarz(JacobiLowOrderCascade(fes, blf));
end
case "additiveSchwarzHighOrder"
solver = AdditiveSchwarzHighOrderPcg(fes, blf, P);
if order == 1
preconditioner = OptimalMLAdditiveSchwarz(P1JacobiCascade(fes, blf, P));
else
preconditioner = OptimalMLAdditiveSchwarz(JacobiHighOrderCascade(fes, blf, P));
end
otherwise
error('No PCG variant %s!', variant)
end
solver = PcgSolver(preconditioner);
% geometric multigrid family
% geometric , Pmultigrid family
case "multigrid"
switch variant
case {"", "lowOrderVcycle"}
if order == 1, solver = LowestOrderLocalMg(fes, blf, P);
else, solver = LocalMgLowOrderVcycle(fes, blf); end
if order == 1
smoother = P1JacobiCascade(fes, blf, P);
else
smoother = JacobiLowOrderCascade(fes, blf);
end
case "highOrderVcycle"
if order == 1, solver = LowestOrderLocalMg(fes, blf, P);
else, solver = LocalMgHighOrderVcycle(fes, blf, P); end
if order == 1
smoother = P1JacobiCascade(fes, blf, P);
else
smoother = JacobiHighOrderCascade(fes, blf, P);
end
otherwise
error('No multigrid variant %s!', variant)
end
solver = OptimalVcycleMultigridSolver(smoother);
% direct solve (for testing purposes)
case "direct"
......@@ -62,5 +89,3 @@ switch class
end
end
lib/solvers/AdditiveSchwarzHighOrderPcg.m deleted 100644 → 0
+ 0
− 161
View file @ f79d4b6d
% AdditiveSchwarzHighOrderPcg (subclass of PcgSolver) Solves linear equations
% iteratively using the CG method with optimal multilevel additive
% Schwartz preconditioner.
classdef AdditiveSchwarzHighOrderPcg < PcgSolver
%% properties
properties (GetAccess=public,SetAccess=protected)
nLevels
end
properties (Access=protected)
blf
hoFes
loFes
coarseLoFes
coarseHoFes
p1Acoarse
smoother
P
intergridMatrix
freeVertices
freeVerticesOld
freeDofs
freeDofsOld
changedPatches
patchwiseA
end
properties (Access=private)
listenerHandle
end
%% methods
methods (Access=public)
function obj = AdditiveSchwarzHighOrderPcg(fes, blf, P)
arguments
fes FeSpace
blf BilinearForm
P Prolongation
end
obj = obj@PcgSolver();
assert(fes.finiteElement.order > 1, ...
'AdditiveSchwartzPcg only works for higher order finite elements.')
assert(isempty(blf.b), ...
'Additive Schwarz PCG solvers only tested for symmetric problems.')
obj.nLevels = 0;
mesh = fes.mesh;
obj.loFes = FeSpace(mesh, LowestOrderH1Fe, 'dirichlet', ':');
obj.hoFes = fes;
obj.P = P;
obj.blf = blf;
obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
end
function setupSystemMatrix(obj, A)
obj.nLevels = obj.nLevels + 1;
L = obj.nLevels;
obj.freeVerticesOld = obj.freeVertices;
obj.freeVertices = getFreeDofs(obj.loFes);
obj.freeDofsOld = obj.freeDofs;
obj.freeDofs = getFreeDofs(obj.hoFes);
ppMatrix = assemble(obj.blf, obj.hoFes);
ppMatrix = ppMatrix(getFreeDofs(obj.hoFes), getFreeDofs(obj.hoFes));
p1Matrix = assemble(obj.blf, obj.loFes);
p1Matrix = p1Matrix(obj.freeVertices, obj.freeVertices);
if L == 1
obj.p1Acoarse = p1Matrix;
fes = obj.hoFes;
coarseMesh = clone(fes.mesh);
obj.coarseLoFes = FeSpace(coarseMesh, LowestOrderH1Fe, 'dirichlet', fes.bnd.dirichlet);
obj.coarseHoFes = FeSpace(coarseMesh, HigherOrderH1Fe(fes.finiteElement.order), 'dirichlet', fes.bnd.dirichlet);
else
obj.smoother{L} = full(diag(ppMatrix)).^(-1);
obj.intergridMatrix{L} = obj.P.matrix(obj.freeDofs, obj.freeDofsOld);
obj.changedPatches{L} = obj.changedPatches{L}(obj.freeVertices);
obj.patchwiseA{L} = assemblePatchwise(obj.blf, obj.hoFes, obj.changedPatches{L});
end
setupSystemMatrix@PcgSolver(obj, A);
end
function setupRhs(obj, b, varargin)
setupRhs@PcgSolver(obj, b, varargin{:});
end
% preconditioner: inverse of diagonal on each level
function Cx = preconditionAction(obj, residual)
assert(~isempty(obj.nLevels), 'Data for multilevel iteration not given!')
L = obj.nLevels;
rho = cell(L, 1);
if L == 1
Cx = obj.A \ residual;
return
end
% descending cascade
for k = L:-1:2
rho{k} = localSmoothing(obj, k, residual);
residual = obj.intergridMatrix{k}'*residual;
end
% exact solve on coarsest level
residual = interpolateFreeData(residual, obj.coarseHoFes, obj.coarseLoFes);
sigma = obj.p1Acoarse \ residual;
sigma = interpolateFreeData(sigma, obj.coarseLoFes, obj.coarseHoFes);
% ascending cascade
for k = 2:L
sigma = obj.intergridMatrix{k}*sigma;
sigma = sigma + rho{k};
end
Cx = sigma;
end
end
methods (Access=protected)
% callback function to be executed before mesh refinement:
% -> patches change for all vertices adjacent to refined edges and
% all new vertices
function getChangedPatches(obj, mesh, bisecData)
% Selects vertices of bisected edges and all newly created
% vertices
nCOld = mesh.nCoordinates;
nCNew = mesh.nCoordinates + nnz(bisecData.bisectedEdges) + ...
bisecData.nRefinedElements'*bisecData.nInnerNodes;
bisectedEdgeNodes = unique(mesh.edges(:,bisecData.bisectedEdges));
obj.changedPatches{obj.nLevels+1} = false(nCNew, 1);
idx = [bisectedEdgeNodes; ((nCOld+1):nCNew)'];
obj.changedPatches{obj.nLevels+1}(idx) = true;
end
function rho = localSmoothing(obj, k, res)
rho = obj.patchwiseA{k} \ res;
end
end
end
%% auxiliary functions
function interpolatedData = interpolateFreeData(data, fromFes, toFes)
freeDofs = getFreeDofs(toFes);
nComponents = size(data, 2);
interpolatedData = zeros(numel(freeDofs), nComponents);
feFunctionWrapper = FeFunction(fromFes);
for k = 1:nComponents
feFunctionWrapper.setFreeData(data(:,k));
wholeData = nodalInterpolation(feFunctionWrapper, toFes);
interpolatedData(:,k) = wholeData(freeDofs);
end
end
\ No newline at end of file
lib/solvers/AdditiveSchwarzLowOrderPcg.m deleted 100644 → 0
+ 0
− 152
View file @ f79d4b6d
% AdditiveSchwartzPcg (subclass of PcgSolver) Solves linear equations
% iteratively using the CG method with optimal multilevel additive
% Schwartz preconditioner.
classdef AdditiveSchwarzLowOrderPcg < PcgSolver
%% properties
properties (GetAccess=public,SetAccess=protected)
nLevels
end
properties (Access=protected)
blf
hoFes
loFes
p1Acoarse
p1Smoother
P
inclusionMatrix
intergridMatrix
freeVertices
freeVerticesOld
changedPatches
patchwiseA
patchwiseP1Matrix
end
properties (Access=private)
listenerHandle
end
%% methods
methods (Access=public)
function obj = AdditiveSchwarzLowOrderPcg(fes, blf)
arguments
fes FeSpace
blf BilinearForm
end
obj = obj@PcgSolver();
assert(fes.finiteElement.order > 1, ...
'AdditiveSchwartzPcg only works for higher order finite elements.')
assert(isempty(blf.b), ...
'Additive Schwarz PCG solvers only tested for symmetric problems.')
obj.nLevels = 0;
mesh = fes.mesh;
obj.loFes = FeSpace(mesh, LowestOrderH1Fe, 'dirichlet', ':');
obj.hoFes = fes;
obj.P = Prolongation.chooseFor(obj.loFes);
obj.blf = blf;
obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
end
function setupSystemMatrix(obj, A)
obj.nLevels = obj.nLevels + 1;
L = obj.nLevels;
obj.freeVerticesOld = obj.freeVertices;
obj.freeVertices = getFreeDofs(obj.loFes);
p1Matrix = assemble(obj.blf, obj.loFes);
p1Matrix = p1Matrix(obj.freeVertices, obj.freeVertices);
if L == 1
obj.p1Acoarse = p1Matrix;
else
obj.p1Smoother{L} = full(diag(p1Matrix)).^(-1);
obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
obj.changedPatches{L} = obj.changedPatches{L}(obj.freeVertices);
obj.patchwiseA = assemblePatchwise(obj.blf, obj.hoFes);
obj.inclusionMatrix = SpaceProlongation(obj.loFes, obj.hoFes) ...
.matrix(getFreeDofs(obj.loFes), getFreeDofs(obj.hoFes));
end
setupSystemMatrix@PcgSolver(obj, A);
end
function setupRhs(obj, b, varargin)
setupRhs@PcgSolver(obj, b, varargin{:});
end
% preconditioner: inverse of diagonal on each level
function Cx = preconditionAction(obj, residual)
assert(~isempty(obj.nLevels), 'Data for multilevel iteration not given!')
L = obj.nLevels;
rho = cell(L, 1);
if L == 1
Cx = obj.A \ residual;
return
end
% descending cascade
rho{L} = hoGlobalSmoothing(obj, residual);
residual = obj.inclusionMatrix * residual;
for k = L:-1:3
residual = obj.intergridMatrix{k}'*residual;
rho{k-1} = p1LocalSmoothing(obj, k-1, residual);
end
% exact solve on coarsest level
residual = obj.intergridMatrix{2}'*residual;
sigma = obj.p1Acoarse \ residual;
sigma = obj.intergridMatrix{2}*sigma;
% ascending cascade
for k = 3:L
sigma = sigma + rho{k-1};
sigma = obj.intergridMatrix{k}*sigma;
end
sigma = obj.inclusionMatrix' * sigma;
sigma = sigma + rho{L};
Cx = sigma;
end
end
methods (Access=protected)
% callback function to be executed before mesh refinement:
% -> patches change for all vertices adjacent to refined edges and
% all new vertices
function getChangedPatches(obj, mesh, bisecData)
% Selects vertices of bisected edges and all newly created
% vertices
nCOld = mesh.nCoordinates;
nCNew = mesh.nCoordinates + nnz(bisecData.bisectedEdges) + ...
bisecData.nRefinedElements'*bisecData.nInnerNodes;
bisectedEdgeNodes = unique(mesh.edges(:,bisecData.bisectedEdges));
obj.changedPatches{obj.nLevels+1} = false(nCNew, 1);
idx = [bisectedEdgeNodes; ((nCOld+1):nCNew)'];
obj.changedPatches{obj.nLevels+1}(idx) = true;
end
function rho = p1LocalSmoothing(obj, k, res)
idx = obj.changedPatches{k};
rho = zeros(size(res));
rho(idx,:) = obj.p1Smoother{k}(idx).*res(idx,:);
% DEBUG: Local patchwise smoothing in P1
% rho = obj.patchwiseP1Matrix{k} \ res;
end
function rho = hoGlobalSmoothing(obj, res)
rho = obj.patchwiseA \ res;
end
end
end
\ No newline at end of file
lib/solvers/ICholPcgSolver.m deleted 100644 → 0
+ 0
− 39
View file @ f79d4b6d
% ICholPcgSolver (subclass of PcgSolver) Solves linear equations iteratively
% using the CG method with incomplete Cholesky decomposition preconditioner.
classdef ICholPcgSolver < PcgSolver
%% properties
properties (GetAccess=public,SetAccess=protected)
fes
end
properties (Access=private)
C
end
%% methods
methods (Access=public)
% preconditioner action: inverse of diagonal
function obj = ICholPcgSolver(fes)
arguments
fes FeSpace
end
obj = obj@PcgSolver();
obj.fes = fes;
end
function setupSystemMatrix(obj, A)
obj.C = ichol(A);
setupSystemMatrix@PcgSolver(obj, A);
end
function setupRhs(obj, b, varargin)
setupRhs@PcgSolver(obj, b, varargin{:});
end
function Cx = preconditionAction(obj, x)
Cx = obj.C' \ (obj.C \ x);
end
end
end
\ No newline at end of file
lib/solvers/JacobiPcgSolver.m deleted 100644 → 0
+ 0
− 39
View file @ f79d4b6d
% JacobiPcgSolver (subclass of PcgSolver) Solves linear equations
% iteratively using the CG method with diagonal preconditioner.
classdef JacobiPcgSolver < PcgSolver
%% properties
properties (GetAccess=public,SetAccess=protected)
fes
end
properties (Access=private)
C
end
%% methods
methods (Access=public)
% preconditioner action: inverse of diagonal
function obj = JacobiPcgSolver(fes)
arguments
fes FeSpace
end
obj = obj@PcgSolver();
obj.fes = fes;
end
function setupSystemMatrix(obj, A)
obj.C = full(diag(A)).^(-1);
setupSystemMatrix@PcgSolver(obj, A);
end
function setupRhs(obj, b, varargin)
setupRhs@PcgSolver(obj, b, varargin{:});
end
function Cx = preconditionAction(obj, x)
Cx = obj.C .* x;
end
end
end
\ No newline at end of file
lib/solvers/LocalMgHighOrderVcycle.m deleted 100644 → 0
+ 0
− 172
View file @ f79d4b6d
% OptimalLocalMGSolver (subclass of MGSolver) Solves linear equations
% iteratively using high-order p-optimal multilevel local smoothing (additive
% Schwarz/block Jacobi).
% One iteration = Vcycle(0,1) (no pre-/one post-smoothing step).
%
% With the solver comes an
classdef LocalMgHighOrderVcycle < MGSolver
%% properties
properties (GetAccess=public,SetAccess=protected)
nLevels
end
properties (Access=protected)
blf
hoFes
loFes
P
intergridMatrix
inclusionMatrix
freeDofs
freeDofsOld
changedPatches
patchwiseA
Acoarse
Afine
coarseLoFes
coarseHoFes
end
%% event data
properties (Access=private)
listenerHandle
end
%% methods
methods (Access=public)
function obj = LocalMgHighOrderVcycle(fes, blf, P)
arguments
fes FeSpace
blf BilinearForm
P Prolongation
end
obj = obj@MGSolver();
assert(fes.finiteElement.order > 1, ...
'LocalMgLowOrderVcycle only works for higher order finite elements.')
assert(isempty(blf.b), ...
'Multigrid solvers only tested for symmetric problems.')
obj.nLevels = 0;
mesh = fes.mesh;
obj.P = P;
obj.loFes = FeSpace(mesh, LowestOrderH1Fe(), 'dirichlet', fes.bnd.dirichlet);
obj.hoFes = fes;
obj.blf = blf;
obj.Afine = cell(1,1);
obj.Acoarse = assemble(obj.blf, obj.loFes);
obj.Acoarse = obj.Acoarse(getFreeDofs(obj.loFes), getFreeDofs(obj.loFes));
% set up FE spaces on coarsest level for projection
coarseMesh = clone(fes.mesh);
obj.coarseLoFes = FeSpace(coarseMesh, LowestOrderH1Fe, 'dirichlet', fes.bnd.dirichlet);
obj.coarseHoFes = FeSpace(coarseMesh, HigherOrderH1Fe(fes.finiteElement.order), 'dirichlet', fes.bnd.dirichlet);
obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
end
function setupSystemMatrix(obj, A)
obj.nLevels = obj.nLevels + 1;
L = obj.nLevels;
freeVertices = getFreeDofs(obj.loFes);
obj.freeDofsOld = obj.freeDofs;
obj.freeDofs = getFreeDofs(obj.hoFes);
if L >= 2
obj.intergridMatrix{L} = obj.P.matrix(obj.freeDofs, obj.freeDofsOld);
obj.changedPatches{L} = freeVertices(obj.changedPatches{L}(freeVertices));
obj.inclusionMatrix = SpaceProlongation(obj.coarseLoFes, obj.coarseHoFes) ...
.matrix(getFreeDofs(obj.coarseLoFes), getFreeDofs(obj.coarseHoFes));
if L == 2
obj.patchwiseA{L} = assemblePatchwise(obj.blf, obj.hoFes, ':');
else
obj.patchwiseA{L} = assemblePatchwise(obj.blf, obj.hoFes, obj.changedPatches{L});
end
end
setupSystemMatrix@MGSolver(obj, A);
obj.Afine{L} = obj.A;
end
function setupRhs(obj, b, varargin)
setupRhs@MGSolver(obj, b, varargin{:});
end
% Geometric MultiGrid
function [Cx, algError2] = Vcycle(obj, res)
assert(isequal(size(res, 1), size(obj.A, 1)), ...
'Setup for multilevel iteration not complete!')
% if there is only one coarse level: exact solve
if obj.nLevels == 1
Cx = obj.A \ res;
algError2 = scalarProduct(Cx, obj.A*Cx);
return
end
L = obj.nLevels;
algError2 = zeros(1, size(res, 2)); % built-in estimator of the algebraic error
% descending cascade in p: no smoothing
residual{L} = res;
for k = L:-1:2
residual{k-1} = obj.intergridMatrix{k}'*residual{k};
end
% coarse level in P1
residual{1} = obj.inclusionMatrix * residual{1};
sigma = obj.Acoarse \ residual{1};
algError2 = algError2 + scalarProduct(sigma, obj.Acoarse*sigma);
% ascending cascade in Pp
sigma = obj.inclusionMatrix' * sigma;
for k = 2:L
sigma = obj.intergridMatrix{k}*sigma;
uptres = residual{k} - obj.Afine{k}*sigma; %updated residual
rho = obj.patchwiseA{k} \ uptres;
[aeUpd, sUpd] = computeOptimalUpdate(obj.Afine{k}, uptres, rho);
sigma = sigma + sUpd;
algError2 = algError2 + aeUpd;
end
Cx = sigma;
end
end
methods (Access=protected)
% callback function to be executed before mesh refinement:
% -> patches change for all vertices adjacent to refined edges and
% all new vertices
function getChangedPatches(obj, mesh, bisecData)
nCOld = mesh.nCoordinates;
nCNew = mesh.nCoordinates + nnz(bisecData.bisectedEdges) + ...
bisecData.nRefinedElements'*bisecData.nInnerNodes;
bisectedEdgeNodes = unique(mesh.edges(:,bisecData.bisectedEdges));
obj.changedPatches{obj.nLevels+1} = false(nCNew, 1);
idx = [unique(bisectedEdgeNodes); ((nCOld+1):nCNew)'];
obj.changedPatches{obj.nLevels+1}(idx) = true;
end
end
end
%% auxiliary functions
% error correction with optimal stepsize
function [etaUpdate, sigmaUpdate] = computeOptimalUpdate(A, res, rho)
rhoArho = scalarProduct(rho, A*rho);
lambda = scalarProduct(res, rho) ./ rhoArho;
sigmaUpdate = lambda.*rho;
etaUpdate = rhoArho.*(lambda.^2);
if any(lambda > 3)
warning('MG step-sizes no longer bound by d+1. Optimality of step size cannot be guaranteed!')
end
end
function a = scalarProduct(x, y)
a = sum(x.*y, 1);
end
\ No newline at end of file
lib/solvers/LocalMgLowOrderVcycle.m deleted 100644 → 0
+ 0
− 179
View file @ f79d4b6d
% OptimalLocalMGSolver (subclass of MGSolver) Solves linear equations
% iteratively using high-order p-optimal multilevel local smoothing (additive
% Schwarz/block Jacobi).
% One iteration = Vcycle(0,1) (no pre-/one post-smoothing step).
%
% With the solver comes an
classdef LocalMgLowOrderVcycle < MGSolver
%% properties
properties (GetAccess=public,SetAccess=protected)
nLevels
end
properties (Access=protected)
blf
hoFes
loFes
p1Matrix
p1Smoother
P
intergridMatrix
inclusionMatrix
freeVertices
freeVerticesOld
changedPatches
patchwiseA
end
%% event data
properties (Access=private)
listenerHandle
end
%% methods
methods (Access=public)
function obj = LocalMgLowOrderVcycle(fes, blf)
arguments
fes FeSpace
blf BilinearForm
end
obj = obj@MGSolver();
assert(fes.finiteElement.order > 1, ...
'LocalMgLowOrderVcycle only works for higher order finite elements.')
assert(isempty(blf.b), ...
'Multigrid solvers only tested for symmetric problems.')
obj.nLevels = 0;
mesh = fes.mesh;
obj.loFes = FeSpace(mesh, LowestOrderH1Fe(), 'dirichlet', ':');
obj.hoFes = fes;
obj.P = Prolongation.chooseFor(obj.loFes);
obj.blf = blf;
obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
end
function setupSystemMatrix(obj, A)
obj.nLevels = obj.nLevels + 1;
L = obj.nLevels;
obj.freeVerticesOld = obj.freeVertices;
obj.freeVertices = getFreeDofs(obj.loFes);
freeVerticesHigher = getFreeDofs(obj.hoFes);
obj.p1Matrix{L} = assemble(obj.blf, obj.loFes);
obj.p1Matrix{L} = obj.p1Matrix{L}(obj.freeVertices, obj.freeVertices);
obj.p1Smoother{L} = full(diag(obj.p1Matrix{L})).^(-1);
if L >= 2
obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeVertices));
obj.patchwiseA = assemblePatchwise(obj.blf, obj.hoFes);
obj.inclusionMatrix = SpaceProlongation(obj.loFes, obj.hoFes) ...
.matrix(getFreeDofs(obj.loFes), getFreeDofs(obj.hoFes));
end
setupSystemMatrix@MGSolver(obj, A);
end
function setupRhs(obj, b, varargin)
setupRhs@MGSolver(obj, b, varargin{:});
end
% Geometric MultiGrid
function [Cx, algError2] = Vcycle(obj, res)
assert(isequal(size(res, 1), size(obj.A, 1)), ...
'Setup for multilevel iteration not complete!')
% if there is only one coarse level: exact solve
if obj.nLevels == 1
Cx = obj.A \ res;
algError2 = scalarProduct(Cx, obj.A*Cx);
return
end
L = obj.nLevels;
algError2 = zeros(1, size(res, 2)); % built-in estimator of the algebraic error
% descending cascade in P1: no smoothing
residual{L} = obj.inclusionMatrix * res; %interpolateFreeData(res, obj.hoFes, obj.loFes);
for k = L:-1:2
residual{k-1} = obj.intergridMatrix{k}'*residual{k};
end
% exact solve on coarsest level to compute accumulative lifting
% of the residual (sigma)
sigma = obj.p1Matrix{1} \ residual{1};
algError2 = algError2 + scalarProduct(sigma, obj.p1Matrix{1}*sigma);
% ascending cascade in P1 AND local smoothing
for k = 2:(L-1)
sigma = obj.intergridMatrix{k}*sigma;
uptres = residual{k} - obj.p1Matrix{k}*sigma; %updated residual
rho = p1LocalSmoothing(obj, k, uptres);
[aeUpd, sUpd] = computeOptimalUpdate(obj.p1Matrix{k}, uptres, rho);
sigma = sigma + sUpd;
algError2 = algError2 + aeUpd;
end
% smoothing on finest level dependent on p
sigma = obj.intergridMatrix{L}*sigma;
sigma = obj.inclusionMatrix' * sigma; %interpolateFreeData(sigma, obj.loFes, obj.hoFes);
uptres = res - obj.A*sigma;
rho = hoGlobalSmoothing(obj, uptres);
[aeUpd, sUpd] = computeOptimalUpdate(obj.A, uptres, rho);
sigma = sigma + sUpd;
algError2 = algError2 + aeUpd;
Cx = sigma;
end
end
methods (Access=protected)
% callback function to be executed before mesh refinement:
% -> patches change for all vertices adjacent to refined edges and
% all new vertices
function getChangedPatches(obj, mesh, bisecData)
nCOld = mesh.nCoordinates;
nCNew = mesh.nCoordinates + nnz(bisecData.bisectedEdges) + ...
bisecData.nRefinedElements'*bisecData.nInnerNodes;
bisectedEdgeNodes = unique(mesh.edges(:,bisecData.bisectedEdges));
obj.changedPatches{obj.nLevels+1} = false(nCNew, 1);
idx = [unique(bisectedEdgeNodes); ((nCOld+1):nCNew)'];
obj.changedPatches{obj.nLevels+1}(idx) = true;
end
function rho = p1LocalSmoothing(obj, k, res)
idx = obj.changedPatches{k};
rho = zeros(size(res));
rho(idx,:) = obj.p1Smoother{k}(idx).*res(idx,:);
end
function rho = hoGlobalSmoothing(obj, res)
rho = obj.patchwiseA \ res;
end
end
end
%% auxiliary functions
% error correction with optimal stepsize
function [etaUpdate, sigmaUpdate] = computeOptimalUpdate(A, res, rho)
rhoArho = scalarProduct(rho, A*rho);
if max(abs(rho)) < eps
lambda = 1;
else
lambda = scalarProduct(res, rho) ./ rhoArho;
end
sigmaUpdate = lambda.*rho;
etaUpdate = rhoArho.*(lambda.^2);
if any(lambda > 3)
warning('MG step-sizes no longer bound by d+1. Optimality of step size cannot be guaranteed!')
end
end
function a = scalarProduct(x, y)
a = sum(x.*y, 1);
end
lib/solvers/LowestOrderLocalMg.m deleted 100644 → 0
+ 0
− 162
View file @ f79d4b6d
% LowestOrderLocalMg (subclass of MGSolver) Solves linear equations
% stemming from lowest order AFEM computations iteratively using a
% multigrid solver where one iteration is characterized by Vcycle(0,1)
% (no pre-/one post-smoothing step).
% On every level, local Jacobi smoothing is employed on changed patches
% and the smoothed update is weighted by an optimal step size.
%
% See also [Innerberger, Miraçi, Praetorius, Streitberger; Optimal
% computational costs of AFEM with optimal local hp-robust multigrid
% solver].
classdef LowestOrderLocalMg < MGSolver
%% properties
properties (GetAccess=public,SetAccess=protected)
nLevels
end
properties (Access=protected)
fes
blf
matrix
smoother
P
intergridMatrix
freeVertices
freeVerticesOld
changedPatches
end
%% event data
properties (Access=private)
listenerHandle
end
%% methods
methods (Access=public)
function obj = LowestOrderLocalMg(fes, blf, P)
arguments
fes FeSpace
blf BilinearForm
P Prolongation
end
obj = obj@MGSolver();
assert(fes.finiteElement.order == 1, ...
'LowestOrderLocalMg only works for lowest order finite elements.')
assert(isempty(blf.b), ...
'Multigrid solvers only tested for symmetric problems.')
obj.nLevels = 0;
mesh = fes.mesh;
obj.fes = fes;
obj.P = P;
obj.blf = blf;
obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
end
function setupSystemMatrix(obj, A)
obj.nLevels = obj.nLevels + 1;
L = obj.nLevels;
obj.freeVerticesOld = obj.freeVertices;
obj.freeVertices = getFreeDofs(obj.fes);
obj.matrix{L} = A;
obj.smoother{L} = full(diag(obj.matrix{L})).^(-1);
if L >= 2
obj.intergridMatrix{L} = obj.P.matrix(obj.freeVertices, obj.freeVerticesOld);
obj.changedPatches{L} = find(obj.changedPatches{L}(obj.freeVertices));
end
setupSystemMatrix@MGSolver(obj, A);
end
function setupRhs(obj, b, varargin)
setupRhs@MGSolver(obj, b, varargin{:});
end
% Geometric MultiGrid
function [Cx, algError2] = Vcycle(obj, res)
assert(isequal(size(res, 1), size(obj.A, 1)), ...
'Setup for multilevel iteration not complete!')
% if there is only one coarse level: exact solve
if obj.nLevels == 1
Cx = obj.A \ res;
algError2 = scalarProduct(Cx, obj.A*Cx);
return
end
L = obj.nLevels;
algError2 = zeros(1, size(res, 2)); % built-in estimator of the algebraic error
% descending cascade in P1: no smoothing
residual{L} = res;
for k = L:-1:2
residual{k-1} = obj.intergridMatrix{k}'*residual{k};
end
% exact solve on coarsest level to compute accumulative lifting
% of the residual (sigma)
sigma = obj.matrix{1} \ residual{1};
algError2 = algError2 + scalarProduct(sigma, obj.matrix{1}*sigma);
% ascending cascade in P1 AND local smoothing
for k = 2:L
sigma = obj.intergridMatrix{k}*sigma;
uptres = residual{k} - obj.matrix{k}*sigma; %updated residual
rho = localSmoothing(obj, k, uptres);
[aeUpd, sUpd] = computeOptimalUpdate(obj.matrix{k}, uptres, rho);
sigma = sigma + sUpd;
algError2 = algError2 + aeUpd;
end
Cx = sigma;
end
end
methods (Access=protected)
% callback function to be executed before mesh refinement:
% -> patches change for all vertices adjacent to refined edges and
% all new vertices
function getChangedPatches(obj, mesh, bisecData)
nCOld = mesh.nCoordinates;
nCNew = mesh.nCoordinates + nnz(bisecData.bisectedEdges) + ...
bisecData.nRefinedElements'*bisecData.nInnerNodes;
bisectedEdgeNodes = unique(mesh.edges(:,bisecData.bisectedEdges));
obj.changedPatches{obj.nLevels+1} = false(nCNew, 1);
idx = [unique(bisectedEdgeNodes); ((nCOld+1):nCNew)'];
obj.changedPatches{obj.nLevels+1}(idx) = true;
end
function rho = localSmoothing(obj, k, res)
idx = obj.changedPatches{k};
rho = zeros(size(res));
rho(idx,:) = obj.smoother{k}(idx).*res(idx,:);
end
end
end
%% auxiliary functions
% error correction with optimal stepsize
function [etaUpdate, sigmaUpdate] = computeOptimalUpdate(A, res, rho)
rhoArho = scalarProduct(rho, A*rho);
if max(abs(rho)) < eps
lambda = 1;
else
lambda = scalarProduct(res, rho) ./ rhoArho;
end
sigmaUpdate = lambda.*rho;
etaUpdate = rhoArho.*(lambda.^2);
if any(lambda > 3)
warning('MG step-sizes no longer bound by d+1. Optimality of step size cannot be guaranteed!')
end
end
function a = scalarProduct(x, y)
a = sum(x.*y, 1);
end
lib/solvers/MGSolver.m deleted 100644 → 0
+ 0
− 74
View file @ f79d4b6d
% MGSolver (abstract subclass of IterativeSolver) Solves linear equations
% iteratively, using a Vcycle geometric multigrid solver.
%
% Abstract methods to be implemented by subclasses:
%
% solver.Vcycle(x) performs one Vcycle iteration.
%
% See also IterativeSolver
classdef MGSolver < IterativeSolver
%% properties
properties (SetAccess=protected, GetAccess=public)
residualCNorm
algEstimator
end
properties (Access=protected)
residual
Cresidual
normb
end
%% methods
methods (Access=public)
% initialize iteration
function setupSystemMatrix(obj, A, varargin)
setupSystemMatrix@IterativeSolver(obj, A);
end
function setupRhs(obj, b, varargin)
setupRhs@IterativeSolver(obj, b, varargin{:});
% initialize residual & hierarchy
obj.residual = b - obj.A*obj.x;
[obj.Cresidual, obj.algEstimator] = obj.Vcycle(obj.residual);
obj.residualCNorm = sum(obj.residual.*obj.Cresidual, 1);
obj.normb = sqrt(sum(b.^2, 1));
end
% stopping criterion
function tf = isConverged(obj)
tf = ((obj.iterationCount >= obj.maxIter) ...
| (sqrt(obj.residualCNorm) < obj.tol) ...
| (sqrt(obj.residualCNorm)./obj.normb < obj.tol));
end
end
methods (Access=protected)
% one step
function computeUpdate(obj)
% only iterate on active components
idx = find(obj.activeComponents);
% update solution & residual
obj.x(:,idx) = obj.x(:,idx) + obj.Cresidual(:,idx);
obj.residual(:,idx) = obj.residual(:,idx) - obj.A * obj.Cresidual(:,idx);
% residual & update error correction
[obj.Cresidual(:,idx), obj.algEstimator(idx)] = obj.Vcycle(obj.residual(:,idx));
obj.residualCNorm(:,idx) = sum(obj.residual(:,idx).*obj.Cresidual(:,idx), 1);
end
end
methods (Abstract,Access=public)
Vcycle(obj, x)
end
end
lib/solvers/OptimalVcycleMultigridSolver.m 0 → 100644
+ 144
− 0
View file @ 1d449bc9
% OptimalVcycleMultigridSolver (subclass of IterativeSolver) Solves linear
% equations iteratively, using a geometric multigrid solver with specified
% MultiLevelSmoother.
% One iteration = Vcycle(0,1) (no pre-/one post-smoothing step).
%
% To obtain the solver described in [Innerberger, Miraçi, Praetorius,
% Streitberger; Optimal computational costs of AFEM with optimal local
% hp-robust multigrid solver], choose the smoother as DiagonalJacobi (p=1),
% LowOrderBlockJacobi (p > 1, 1->1->p iteration), or HighOrderBlockJacobi
% (p > 1, 1->p->p iteration).
% This solver also comes with built-in estimator for the algebraic error,
% which is also returned by this multigrid solver.
%
% See also IterativeSolver, MultiLevelSmoother, DiagonalJacobi,
% LowOrderBlockJacobi, HighOrderBlockJacobi
classdef OptimalVcycleMultigridSolver < IterativeSolver
%% properties
properties (SetAccess=protected, GetAccess=public)
residualCNorm
algEstimator
residual
Cresidual
end
properties (Dependent, Access=public)
nLevels
end
properties (Access=protected)
smoother
projectedMatrix
lhsNorm
end
methods
function obj = OptimalVcycleMultigridSolver(smoother)
obj = obj@IterativeSolver();
obj.smoother = smoother;
end
function n = get.nLevels(obj)
n = obj.smoother.nLevels;
end
function setupSystemMatrix(obj, A, varargin)
% override superclass method, because ALL matrices are needed
setupSystemMatrix@IterativeSolver(obj, A, varargin{:});
projectedA = obj.smoother.setup(A, varargin{:});
obj.projectedMatrix{obj.nLevels} = projectedA;
end
function setupRhs(obj, b, varargin)
setupRhs@IterativeSolver(obj, b, varargin{:});
% initialize residual & hierarchy
obj.residual = b - obj.A * obj.x;
[obj.Cresidual, obj.algEstimator] = obj.Vcycle(obj.residual);
obj.residualCNorm = sqrt(dot(obj.residual, obj.Cresidual, Dim.Vector));
obj.lhsNorm = vecnorm(b, 2, Dim.Vector);
end
% stopping criterion
function tf = isConverged(obj)
tf = ((obj.iterationCount >= obj.maxIter) ...
| (obj.residualCNorm < obj.tol) ...
| (obj.residualCNorm ./ obj.lhsNorm < obj.tol));
end
end
methods (Access=protected)
% one step
function computeUpdate(obj)
% only iterate on active components
idx = find(obj.activeComponents);
% update solution & residual
obj.x(:,idx) = obj.x(:,idx) + obj.Cresidual(:,idx);
obj.residual(:,idx) = obj.residual(:,idx) - obj.A * obj.Cresidual(:,idx);
% residual & update error correction
[obj.Cresidual(:,idx), obj.algEstimator(idx)] = obj.Vcycle(obj.residual(:,idx));
obj.residualCNorm(:,idx) = sqrt(dot(obj.residual(:,idx), obj.Cresidual(:,idx), Dim.Vector));
end
function [Cx, algEstimator] = Vcycle(obj, x)
assert(isequal(size(x, 1), size(obj.A, 1)), ...
'Setup for multilevel iteration not complete!')
% if there is only one coarse level: exact solve
L = obj.nLevels;
if L == 1
Cx = obj.A \ x;
algEstimator = sqrt(dot(Cx, obj.A * Cx, Dim.Vector));
return
end
% descending cascade: no smoothing
res{L} = x;
for k = L:-1:2
res{k-1} = obj.smoother.restrict(res{k}, k);
end
% exact solve on coarsest level to compute accumulative lifting of x
Cx = obj.projectedMatrix{1} \ res{1};
eta2 = dot(Cx, obj.projectedMatrix{1} * Cx, Dim.Vector);
% ascending cascade: (local) smoothing and compute error estimator
for k = 2:(L-1)
Cx = obj.smoother.prolongate(Cx, k);
updatedResidual = res{k} - obj.projectedMatrix{k} * Cx;
rho = obj.smoother.smooth(updatedResidual, k);
[etaUpdate2, xUpdate] = computeOptimalUpdate(obj.projectedMatrix{k}, updatedResidual, rho);
Cx = Cx + xUpdate;
eta2 = eta2 + etaUpdate2;
end
% final smoothing and residual update (with non-projected matrix)
Cx = obj.smoother.prolongate(Cx, L);
updatedResidual = res{L} - obj.A * Cx;
rho = obj.smoother.smooth(updatedResidual, L);
[etaUpdate2, xUpdate] = computeOptimalUpdate(obj.A, updatedResidual, rho);
Cx = Cx + xUpdate;
algEstimator = sqrt(eta2 + etaUpdate2);
end
end
end
%% auxiliary functions
% error correction with optimal stepsize
function [etaUpdate2, resUpdate] = computeOptimalUpdate(A, res, rho)
rhoArho = dot(rho, A*rho, Dim.Vector);
if max(abs(rho)) < eps
lambda = 1;
else
lambda = dot(res, rho, Dim.Vector) ./ rhoArho;
end
resUpdate = lambda .* rho;
etaUpdate2 = rhoArho .* lambda.^2;
if any(lambda > 3)
warning('MG step-sizes no longer bound by d+1. Optimality of step size cannot be guaranteed!')
end
end
lib/solvers/PcgSolver.m
+ 26
− 19
View file @ 1d449bc9
......@@ -10,27 +10,41 @@ classdef PcgSolver < IterativeSolver
end
properties (Access=protected)
C
residual
Cresidual
searchDirection
normb
end
%% constructor
methods (Access=public)
function obj = PcgSolver(preconditioner)
arguments
preconditioner (1,1) Preconditioner
end
obj = obj@IterativeSolver();
obj.C = preconditioner;
end
end
%% extend superclass methods
methods (Access=public)
function setupSystemMatrix(obj, A)
setupSystemMatrix@IterativeSolver(obj, A);
obj.C.setup(A);
end
function setupRhs(obj, b, varargin)
setupRhs@IterativeSolver(obj, b, varargin{:});
% initialize residual & search direction
obj.residual = b - obj.A*obj.x;
obj.Cresidual = obj.preconditionAction(obj.residual);
obj.residual = b - obj.A * obj.x;
obj.Cresidual = obj.C.apply(obj.residual);
obj.searchDirection = obj.Cresidual;
obj.residualCNorm = sum(obj.residual.*obj.Cresidual, 1);
obj.normb = sqrt(sum(b.^2, 1));
obj.residualCNorm = sqrt(dot(obj.residual, obj.Cresidual, 1));
obj.normb = vecnorm(b, 2, 1);
end
end
......@@ -38,8 +52,8 @@ classdef PcgSolver < IterativeSolver
methods (Access=public)
function tf = isConverged(obj)
tf = ((obj.iterationCount >= obj.maxIter) ...
| (sqrt(obj.residualCNorm) < obj.tol) ...
| (sqrt(obj.residualCNorm)./obj.normb < obj.tol));
| (obj.residualCNorm < obj.tol) ...
| (obj.residualCNorm ./ obj.normb < obj.tol));
end
end
......@@ -50,30 +64,23 @@ classdef PcgSolver < IterativeSolver
% update solution
AsearchDirection = obj.A * obj.searchDirection(:,idx);
if sum(obj.searchDirection(:,idx).*AsearchDirection, 1) < eps
dAd = dot(obj.searchDirection(:,idx), AsearchDirection, 1);
if dAd < eps
alpha = 1;
else
alpha = obj.residualCNorm(:,idx) ./ sum(obj.searchDirection(:,idx).*AsearchDirection, 1);
alpha = obj.residualCNorm(:,idx).^2 ./ dAd;
end
obj.x(:,idx) = obj.x(:,idx) + alpha .* obj.searchDirection(:,idx);
% DEBUG:
% disp(['alpha = ', num2str(alpha)])
% update residual
obj.residual(:,idx) = obj.residual(:,idx) - alpha .* AsearchDirection;
obj.Cresidual(:,idx) = obj.preconditionAction(obj.residual(:,idx));
obj.Cresidual(:,idx) = obj.C.apply(obj.residual(:,idx));
residualCNormOld = obj.residualCNorm(:,idx);
obj.residualCNorm(:,idx) = sum(obj.residual(:,idx).*obj.Cresidual(:,idx), 1);
obj.residualCNorm(:,idx) = sqrt(dot(obj.residual(:,idx), obj.Cresidual(:,idx), 1));
% update search direction
beta = obj.residualCNorm(:,idx) ./ residualCNormOld;
beta = (obj.residualCNorm(:,idx) ./ residualCNormOld).^2;
obj.searchDirection(:,idx) = obj.Cresidual(:,idx) + beta .* obj.searchDirection(:,idx);
end
end
%% abstract preconditioning
methods (Abstract, Access=public)
preconditionAction(obj, x)
end
end
\ No newline at end of file
lib/solvers/BlockJacobiPcgSolver.m → lib/solvers/preconditioners/BlockJacobi.m
+ 34
− 0
View file @ 1d449bc9
% BlockJacobiPcgSolver (subclass of PcgSolver) Solves linear equations
% iteratively using the CG method with block-diagonal preconditioner.
% BlockJacobi (subclass of Preconditioner) block-diagonal preconditioner.
%
% See also: Preconditioner, PcgSolver
classdef BlockJacobiPcgSolver < PcgSolver
classdef BlockJacobi < Preconditioner
%% properties
properties (GetAccess=public,SetAccess=protected)
fes
end
properties (Access=private)
C
fes
blf
end
%% methods
methods (Access=public)
% preconditioner action: patchwise inverse of diagonal
function obj = BlockJacobiPcgSolver(fes, blf)
function obj = BlockJacobi(fes, blf)
arguments
fes FeSpace
blf BilinearForm
end
obj = obj@PcgSolver();
obj.fes = fes;
obj.blf = blf;
end
function setupSystemMatrix(obj, A)
function setup(obj, A)
obj.C = assemblePatchwise(obj.blf, obj.fes, ':');
setupSystemMatrix@PcgSolver(obj, A);
end
function setupRhs(obj, b, varargin)
setupRhs@PcgSolver(obj, b, varargin{:});
end
function Cx = preconditionAction(obj, x)
function Cx = apply(obj, x)
Cx = obj.C \ x;
end
end
......
lib/solvers/preconditioners/DiagonalJacobi.m 0 → 100644
+ 22
− 0
View file @ 1d449bc9
% DiagonalJacobi (subclass of Preconditioner) diagonal preconditioner.
%
% See also: Preconditioner, PcgSolver
classdef DiagonalJacobi < Preconditioner
%% properties
properties (Access=private)
C
end
%% methods
methods (Access=public)
% preconditioner action: inverse of diagonal
function setup(obj, A)
obj.C = full(diag(A)).^(-1);
end
function Cx = apply(obj, x)
Cx = obj.C .* x;
end
end
end
\ No newline at end of file
lib/solvers/preconditioners/IncompleteCholesky.m 0 → 100644
+ 22
− 0
View file @ 1d449bc9
% IncompleteCholesky (subclass of Preconditioner) incomplete Cholesky
% decomposition preconditioner.
%
% See also: Preconditioner, PcgSolver
classdef IncompleteCholesky < Preconditioner
%% properties
properties (Access=private)
C
end
%% methods
methods (Access=public)
function setup(obj, A)
obj.C = ichol(A);
end
function Cx = apply(obj, x)
Cx = obj.C' \ (obj.C \ x);
end
end
end
\ No newline at end of file
lib/solvers/preconditioners/NoPreconditioner.m 0 → 100644
+ 12
− 0
View file @ 1d449bc9
% NoPreconditioner (subclass of Preconditioner) identity preconditioner (=no-op).
%
% See also: Preconditioner, PcgSolver
classdef NoPreconditioner < Preconditioner
methods (Access=public)
function Cx = apply(~, x)
Cx = x;
end
end
end
\ No newline at end of file
lib/solvers/preconditioners/OptimalMLAdditiveSchwarz.m 0 → 100644
+ 76
− 0
View file @ 1d449bc9
% OptimalMLAdditiveSchwarz (subclass of Preconditioner) optimal multilevel
% additive Schwarz preconditioner with given smoother.
%
% See also: Preconditioner, PcgSolver, MultilevelSmoother
classdef OptimalMLAdditiveSchwarz < Preconditioner
%% properties
properties (Access=protected)
smoother
Acoarse
projectedAcoarse
end
properties (Dependent)
nLevels
end
properties (Access=private)
listenerHandle
end
%% methods
methods (Access=public)
function obj = OptimalMLAdditiveSchwarz(smoother)
arguments
smoother MultilevelSmoother
end
obj = obj@Preconditioner();
obj.smoother = smoother;
end
function setup(obj, A)
if obj.nLevels == 0
obj.Acoarse = A;
obj.projectedAcoarse = obj.smoother.setup(A);
else
obj.smoother.setup(A);
end
end
function Cx = apply(obj, x)
assert(~isempty(obj.nLevels), 'Data for multilevel iteration not given!')
L = obj.nLevels;
rho = cell(L, 1);
if L == 1
Cx = obj.Acoarse \ x;
return
end
% descending cascade
residual = x;
for k = L:-1:2
rho{k} = obj.smoother.smooth(residual, k);
residual = obj.smoother.restrict(residual, k);
end
% exact solve on coarsest level
Cx = obj.projectedAcoarse \ residual;
% ascending cascade
for k = 2:L
Cx = obj.smoother.prolongate(Cx, k);
Cx = Cx + rho{k};
end
end
end
methods
function n = get.nLevels(obj)
n = obj.smoother.nLevels;
end
end
end
\ No newline at end of file
lib/solvers/preconditioners/Preconditioner.m 0 → 100644
+ 23
− 0
View file @ 1d449bc9
% Preconditioner (abstract handle class) Interface for preconditioners.
%
% pc.setup(A, ...) hook to set up the preconditioner for the matrix A.
% Default implementation does nothing.
%
% Abstract methods to implement:
%
% Cx = pc.apply(x) applies the preconditioner to the vector x.
%
% See also: PcgSolver
classdef Preconditioner < handle
methods (Access=public)
function setup(~, ~)
% Default implementation does nothing.
end
end
methods (Abstract, Access=public)
apply(obj, x)
end
end
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

[ 22.Nov 12h : 4,509 users | 1,084 groups | 6,413 projects ]