diff --git a/examples/afemWithGmres.m b/examples/afemWithGmres.m
index 0a1933189f42cfbe925e28dd59d5a14406ca0170..52eeb676c4aa2031f12d81a9feb7bbb2a189fbc4 100644
--- a/examples/afemWithGmres.m
+++ b/examples/afemWithGmres.m
@@ -4,14 +4,14 @@
% ******************************************************************************
%% parameters
-nDofsMax = 5e5;
+nDofsMax = 1e6;
theta = 0.5;
-pmax = 4;
+pmax = 2;
lamalg = 0.1;
%% initialization for every polynomial degree
-[nElem, nDofs, errEst, iterations, time, directSolveTime] = deal(zeros(pmax, 1000));
-for p = 4:pmax
+%[nElem, nDofs, errEst, iterations, time, directSolveTime] = deal(zeros(pmax, 1000));
+for p = 2:pmax
%% setup geometry & spaces
printLogMessage('*** p = %d (of %d) ***', p, pmax)
mesh = Mesh.loadFromGeometry('Lshape');
@@ -36,16 +36,17 @@ for p = 4:pmax
lf = LinearForm();
blf.a = Constant(mesh, 1);
- blf.b = Constant(mesh, [10;10]);
+ blf.b = Constant(mesh, [10; 10]);
blf.c = Constant(mesh, 0);
symmetricBlf.a = blf.a;
lf.f = Constant(mesh, 1);
%% set up solver and operator for nested iteration
- [solver, P] = IterativeSolver.choose(fes, symmetricBlf, "gmres", "additiveSchwarzLowOrder");
+ [solver, P] = IterativeSolver.choose(fes, symmetricBlf, "gmres", ...
+ "multigrid");
solver.tol = 1e-6;
- solver.maxIter = 10;
+ solver.maxIter = 50;
%% adaptive loop
ell = 0;
diff --git a/lib/solvers/@IterativeSolver/IterativeSolver.m b/lib/solvers/@IterativeSolver/IterativeSolver.m
index 496a7d658dc2db4563cd5f03857037957d5d2329..c7f4889e7aaa7f9b1fbdd36a6120fcb83f88ac4c 100644
--- a/lib/solvers/@IterativeSolver/IterativeSolver.m
+++ b/lib/solvers/@IterativeSolver/IterativeSolver.m
@@ -104,7 +104,6 @@ classdef IterativeSolver < handle
function step(obj)
obj.computeUpdate();
- %obj.PgmresIteration;
obj.iterationCount(obj.activeComponents) = obj.iterationCount(obj.activeComponents) + 1;
end
diff --git a/lib/solvers/@IterativeSolver/choose.m b/lib/solvers/@IterativeSolver/choose.m
index c5aa8e5c20bb40852ada55105b6f905f09215fa9..5f29b28b32be89359196f239fe499976564c1f99 100644
--- a/lib/solvers/@IterativeSolver/choose.m
+++ b/lib/solvers/@IterativeSolver/choose.m
@@ -19,7 +19,8 @@ arguments
class (1,1) string {mustBeMember(class, ["cg", "pcg", "gmres", "multigrid", "direct"])}
variant (1,1) string {mustBeMember(variant, [ "", ...
"iChol", "jacobi", ...
- "additiveSchwarzLowOrder", "additiveSchwarzHighOrder" ...
+ "additiveSchwarzLowOrder", "additiveSchwarzHighOrder", ...
+ "multigrid",...
"lowOrderVcycle", "highOrderVcycle"])} = ""
end
@@ -67,16 +68,25 @@ switch class
else
preconditioner = OptimalMLAdditiveSchwarz(JacobiLowOrderCascade(fes, blf));
end
+ solver = PgmresSolver(preconditioner);
case "additiveSchwarzHighOrder"
if order == 1
preconditioner = OptimalMLAdditiveSchwarz(P1JacobiCascade(fes, blf, P));
else
preconditioner = OptimalMLAdditiveSchwarz(JacobiHighOrderCascade(fes, blf, P));
end
+ solver = PgmresSolver(preconditioner);
+ case "multigrid"
+ if order == 1
+ smoother = P1JacobiCascade(fes, blf, P);
+ else
+ smoother = JacobiLowOrderCascade(fes, blf);
+ end
+ solver = PgmresMultigrid(smoother);
otherwise
error('No PCG variant %s!', variant)
end
- solver = PgmresSolver(preconditioner);
+
% geometric , Pmultigrid family
case "multigrid"
diff --git a/lib/solvers/@PgmresSolverSetup/PgmresSolverSetup.m b/lib/solvers/@PgmresSolverSetup/PgmresSolverSetup.m
deleted file mode 100644
index b61e09f5a2d41b01d89a46154528d8775b495576..0000000000000000000000000000000000000000
--- a/lib/solvers/@PgmresSolverSetup/PgmresSolverSetup.m
+++ /dev/null
@@ -1,111 +0,0 @@
-% PgmresSolverSetup (abstract handle class) Interface for a
-% preconditioned GMRES solver setup.
-%% NOTE: This is still a work in progress and not yet fully functional.
-
-
-classdef PgmresSolverSetup < handle
- %% properties
- properties (Access=public)
- maxIter (1,:) double {mustBePositive}
- tol (1,:) double {mustBePositive}
- end
-
- properties (SetAccess=protected, GetAccess=public)
- A
- symA
- b
- x
- iterationCount
- activeComponents
- end
-
- methods (Abstract, Access=public)
- isConverged(obj)
- end
-
- methods (Abstract, Access=protected)
- PgmresIteration(obj)
- end
-
- %% generic routines for iterative solvers
- % overwrite to include additional steps
- methods (Access=public)
- function obj = PgmresSolverSetup()
- obj.maxIter = 100;
- obj.tol = 1e-8;
- end
-
- function setupSystemMatrix(obj, A, symA)
- arguments
- obj
- A (:,:) double {PgmresSolverSetup.mustBeQuadratic}
- symA (:,:) double {PgmresSolverSetup.mustBeAdmissible}
- end
- obj.A = A;
- obj.symA = symA;
- end
-
- function setupRhs(obj, b, x0)
- arguments
- obj
- b (:,:) double
- x0 (:,:) double {PgmresSolverSetup.mustBeCompatible(b, x0)} = zeros(size(b));
- end
- assert(isequal(size(obj.A,2), size(b,1)), ...
- 'Right-hand size and matrix must be of compatible size.')
- obj.b = b;
- obj.x = x0;
- obj.iterationCount = zeros(1, size(b,2));
- obj.activeComponents = true(1, size(b,2));
- end
-
- function x = solve(obj)
- obj.PgmresIteration() %do the full gmres iteration instead of a stepwise approach
-
- x = obj.x;
- end
-
- % function step(obj)
- % obj.computeUpdate();
- % obj.iterationCount(obj.activeComponents) = obj.iterationCount(obj.activeComponents) + 1;
- % end
-
- function tf = applyStoppingCriterion(obj)
- tf = isConverged(obj);
- 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)
- function mustBeQuadratic(A)
- if ~(size(A,1) == size(A,2))
- eidType = 'Solver:matrixNotQuadratic';
- msgType = 'Matrix must be quadratic.';
- throwAsCaller(MException(eidType,msgType))
- end
- end
-
- function mustBeAdmissible(A)
- % check if matrix for preconditioner is symmetric
- if ~issymmetric(A)
- eidType = 'Solver:matrixNotSymmetric';
- msgType = 'Matrix must be symmetric.';
- throwAsCaller(MException(eidType,msgType))
- end
- end
-
- function mustBeCompatible(b, x0)
- if ~isempty(x0) && ~isequal(size(b), size(x0))
- eidType = 'Solver:dataNotCompatible';
- msgType = 'Size of RHS and initial guess must be compatible.';
- throwAsCaller(MException(eidType,msgType))
- end
- end
- end
-end
\ No newline at end of file
diff --git a/lib/solvers/@PgmresSolverSetup/choose.m b/lib/solvers/@PgmresSolverSetup/choose.m
deleted file mode 100644
index 77900eda3bea667edf868af8842e621072ded5e6..0000000000000000000000000000000000000000
--- a/lib/solvers/@PgmresSolverSetup/choose.m
+++ /dev/null
@@ -1,54 +0,0 @@
-% 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] = choose(fes, blf, variant)
-arguments
- fes FeSpace
- blf BilinearForm
- variant (1,1) string {mustBeMember(variant, [ "", ...
- "iChol", "jacobi", ...
- "additiveSchwarzLowOrder", "additiveSchwarzHighOrder" ...
- "lowOrderVcycle", "highOrderVcycle"])} = ""
-end
-
-order = fes.finiteElement.order;
-P = Prolongation.chooseFor(fes);
-
-switch variant
- case "iChol"
- preconditioner = IncompleteCholesky();
- case "jacobi"
- if order == 1
- preconditioner = DiagonalJacobi();
- else
- preconditioner = BlockJacobi(fes, blf);
- end
- case {"", "additiveSchwarzLowOrder"}
- if order == 1
- preconditioner = OptimalMLAdditiveSchwarz(P1JacobiCascade(fes, blf, P));
- else
- preconditioner = OptimalMLAdditiveSchwarz(JacobiLowOrderCascade(fes, blf));
- end
- case "additiveSchwarzHighOrder"
- 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 = PgmresSolver(preconditioner);
-
-end
\ No newline at end of file
diff --git a/lib/solvers/PgmresMultigrid.m b/lib/solvers/PgmresMultigrid.m
new file mode 100644
index 0000000000000000000000000000000000000000..c8343d6608704eea59aa24e7d82ab87f678f84a9
--- /dev/null
+++ b/lib/solvers/PgmresMultigrid.m
@@ -0,0 +1,218 @@
+% PgmresMG01 (abstract subclass of Iterative) Solves nonsymmetric linear equations
+% iteratively by an optimally preconditioned generalized minimal residual (GMRES) method.
+%
+% See also: IterativeSolver
+%% NOTE: This is still a work in progress and not yet fully functional.
+
+classdef PgmresMultigrid < IterativeSolver
+
+ properties (SetAccess=protected, GetAccess=public)
+ residualCNorm % norm of the preconditioned residual
+ end
+
+ properties (Dependent, Access=public)
+ nLevels
+ end
+
+ properties (Access=protected)
+ smoother
+ projectedMatrix
+ lhsNorm
+ end
+
+ properties (Access=protected)
+ C % preconditioner
+ residual % residual vector
+ Cresidual % preconditioned residual vector
+ searchDirection % search direction
+ normb % norm of right-hand side
+ H % Hessenberg matrix for Arnoldi process
+ V % Krylov subspace basis
+ PV % preconditioned Krylov subspace basis
+ g % modified residual vector
+ sine
+ cosine
+ end
+
+ %% constructor
+ methods
+ function obj = PgmresMultigrid (smoother)
+ arguments
+ smoother (1,1)
+ end
+
+ obj@IterativeSolver();
+ obj.smoother = smoother;
+ end
+
+ function n = get.nLevels(obj)
+ n = obj.smoother.nLevels;
+ end
+
+ %% extend superclass methods
+ function setupSystemMatrix(obj, A, symA, varargin)
+ setupSystemMatrix@IterativeSolver(obj, A, symA);
+ projectedA = obj.smoother.setup(symA, varargin{:});
+ obj.projectedMatrix{obj.nLevels} = projectedA;
+
+ % initialize Hessenberg matrix and Krylov vectors
+ n = size(A, 1);
+ obj.V = zeros(n, obj.maxIter + 1); % sparse matrix for Krylov vectors
+ obj.PV = zeros(n, obj.maxIter + 1); % sparse matrix for preconditioned Krylov vectors
+ obj.H = sparse(obj.maxIter + 1, obj.maxIter); % sparse matrix for Hessenberg matrix
+ end
+
+ function setupRhs(obj, b, varargin)
+ setupRhs@IterativeSolver(obj, b, varargin{:});
+
+ % initialize residual & search direction
+ obj.b = b;
+ obj.residual = b - obj.A * obj.x;
+ [obj.Cresidual, ~] = obj.Vcycle(obj.residual);
+ obj.searchDirection = obj.Cresidual;
+ obj.residualCNorm = sqrt(dot(obj.residual, obj.Cresidual, 1));
+ obj.normb = vecnorm(b, 2, 1);
+
+ % initialize modified residual vector
+ obj.g = zeros(obj.maxIter + 1, 1);
+ obj.g(1) = obj.residualCNorm;
+
+ obj.cosine = zeros(obj.maxIter-1, 1); % % cosine of Givens rotation
+ obj.sine = zeros(obj.maxIter-1, 1); % % sine of Givens rotation
+ end
+ end
+
+ %% implement abstract superclass methods
+ methods (Access=public)
+ function tf = isConverged(obj)
+ tf = ((obj.iterationCount >= obj.maxIter) ...
+ | (obj.residualCNorm < obj.tol) ...
+ | (obj.residualCNorm ./ obj.normb < obj.tol));
+ end
+ end
+
+ methods (Access=protected)
+ function computeUpdate(obj)
+
+ % compute first Krylov vector
+ if obj.iterationCount == 0
+ obj.V(:, 1) = obj.Cresidual / obj.residualCNorm; % normalize search direction
+ end
+
+ % current basis vector
+ k = obj.iterationCount + 1;
+ AVk = obj.A * obj.V(:, k);
+
+ obj.PV(:, 1) = obj.residual/obj.residualCNorm;
+
+ % Arnoldi process
+ for i = 1:k
+ obj.H(i, k) = obj.V(:, i)' * AVk;
+ end
+
+ w = AVk - obj.PV(:, 1:k) * obj.H(1:k, k);
+ %Cw = obj.C.apply(w);
+ [Cw, ~] = obj.Vcycle(w);
+
+ obj.H(k+1, k) = sqrt(Cw'*w); % next off-diagonal element
+
+ % Apply Givens rotation
+ for i = 1:k-1
+ temp = obj.cosine(i) * obj.H(i, k) + obj.sine(i) * obj.H(i + 1, k);
+ obj.H(i + 1, k) = -obj.sine(i) * obj.H(i, k) + obj.cosine(i) * obj.H(i + 1, k);
+ obj.H(i, k) = temp;
+ end
+
+ beta = sqrt( obj.H(k, k)^2 + obj.H(k + 1, k)^2);
+ obj.cosine(k) = obj.H(k, k) / beta;
+ obj.sine(k) = obj.H(k + 1, k) / beta;
+ obj.H(k, k) = beta;
+
+ % Update residual vector
+ obj.g(k + 1) = -obj.sine(k) * obj.g(k);
+ obj.g(k) = obj.cosine(k) * obj.g(k);
+
+ % solve upper triangular system
+ y = zeros(k,1);
+ y(k) = obj.g(k)/obj.H(k,k);
+ for i=k-1:-1:1
+ tmp = obj.H(i,i+1:k) * y(i+1:k);
+ y(i) = (obj.g(i)-tmp) / obj.H(i,i);
+ end
+ obj.x = obj.x0 + obj.V(:, 1:k) * y;
+
+ if(obj.H(k+1,k)==0)
+ obj.residual = 0;
+ return;
+ end
+
+ % update residual
+ obj.residual = obj.b - obj.A * obj.x;
+ [obj.Cresidual, ~] = obj.Vcycle(obj.residual);
+ obj.residualCNorm = sqrt(dot(obj.residual, obj.Cresidual, 1));
+
+ % set new search direction
+ obj.V(:, k+1) = Cw / obj.H(k+1,k);
+ obj.PV(:, k+1) = w / obj.H(k+1,k);
+ 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
\ No newline at end of file
diff --git a/lib/solvers/PgmresSolver.m b/lib/solvers/PgmresSolver.m
index d519fb25e6f6cd85d127a1681bbbc57a137d3b59..9cbb8dcb4b6ccd013dff6ac0a28f5d3a6b64215d 100644
--- a/lib/solvers/PgmresSolver.m
+++ b/lib/solvers/PgmresSolver.m
@@ -101,17 +101,16 @@ classdef PgmresSolver < IterativeSolver
w = AVk - obj.PV(:, 1:k) * obj.H(1:k, k);
Cw = obj.C.apply(w);
+
+ obj.H(k+1, k) = sqrt(Cw'*w); % next off-diagonal element
- obj.H(k+1, k) = sqrt(Cw'*w); % next diagonal element
-
- % Givens-Rotation anwenden
+ % Apply Givens rotation
for i = 1:k-1
temp = obj.cosine(i) * obj.H(i, k) + obj.sine(i) * obj.H(i + 1, k);
obj.H(i + 1, k) = -obj.sine(i) * obj.H(i, k) + obj.cosine(i) * obj.H(i + 1, k);
obj.H(i, k) = temp;
end
- % Apply Givens rotation
beta = sqrt( obj.H(k, k)^2 + obj.H(k + 1, k)^2);
obj.cosine(k) = obj.H(k, k) / beta;
obj.sine(k) = obj.H(k + 1, k) / beta;
@@ -135,7 +134,6 @@ classdef PgmresSolver < IterativeSolver
return;
end
-
% update residual
obj.residual = obj.b - obj.A * obj.x;
obj.Cresidual = obj.C.apply(obj.residual);
@@ -145,84 +143,5 @@ classdef PgmresSolver < IterativeSolver
obj.V(:, k+1) = Cw / obj.H(k+1,k);
obj.PV(:, k+1) = w / obj.H(k+1,k);
end
-
- % function PgmresIteration(obj)
- % %% only iterate on active components
- % idx = find(obj.activeComponents);
- %
- % %%%%%%%%
- % % Begin construction site : rename later
- % %%%%%%%%%%
- % x0 = obj.x;
- %
- % % Initialisierung
- % N = size(obj.residual,1);
- % maxIt = obj.maxIter;
- %
- % %V = zeros(N,maxIt+1);
- % %H = zeros(maxIt+1,maxIt);
- % cs = zeros(maxIt+1,1);
- % sn = zeros(maxIt+1,1);
- % g = zeros(maxIt+1,1);
- % g(1) = obj.residualCNorm;
- %
- %
- % % Start
- % V(:,1) = obj.Cresidual / obj.residualCNorm;
- % PV = zeros(N,maxIt+1);
- % PV(:,1) = obj.residual/obj.residualCNorm;
- %
- % for j=1:maxIt
- % AVj = obj.A * V(:,j);
- % for i=1:j
- % H(i,j) = V(:,i)'*AVj;
- % end
- %
- % w = AVj-PV(:,1:j) * H(:,j);
- % Cw = obj.C.apply(w);
- %
- % H(j+1,j) = sqrt(Cw'*w);
- %
- %
- % for i=1:j-1
- % Hip1jtmp = -sn(i+1)*H(i,j) + cs(i+1)* H(i+1,j);
- % H(i,j) = cs(i+1) * H(i,j) + sn(i+1) * H(i+1,j);
- % H(i+1,j) = Hip1jtmp;
- % end
- %
- % beta = sqrt(H(j,j)^2 + H(j+1,j)^2);
- % sn(j+1) = H(j+1,j)/beta;
- % cs(j+1) = H(j,j)/beta;
- % H(j,j) = beta;
- % g(j+1) = -sn(j+1)*g(j);
- % g(j) = cs(j+1)*g(j);
- %
- % y = zeros(j,1);
- % y(j) = g(j)/H(j,j);
- % for i=j-1:-1:1
- % tmp = H(i,i+1:j)*y(i+1:j);
- % y(i) = (g(i)-tmp)/H(i,i);
- % end
- %
- % obj.x = x0 + V(:,1:j)*y;
- %
- % if(H(j+1,j)==0)
- % res = 0;
- % relres = 0;
- % return;
- % end
- %
- % runprec = obj.b-obj.A * obj.x;
- % r = obj.C.apply(runprec);
- % res = sqrt(r'*runprec);
- %
- % if(res<obj.tol*obj.residualCNorm || res<5e-14)
- % return
- % end
- %
- % V(:,j+1) = Cw/H(j+1,j);
- % PV(:,j+1) = w/H(j+1,j);
- % end
- % end
end
end
\ No newline at end of file