From e97361c46180f6793fc3ee035e32df46e3971ee8 Mon Sep 17 00:00:00 2001
From: Julian Streitberger <julian.streitberger@tuwien.ac.at>
Date: Mon, 20 Jan 2025 16:30:18 +0100
Subject: [PATCH] Add new Pgmres with multigrid

---
 examples/afemWithGmres.m                      |  15 +-
 .../@IterativeSolver/IterativeSolver.m        |   1 -
 lib/solvers/@IterativeSolver/choose.m         |  14 +-
 .../@PgmresSolverSetup/PgmresSolverSetup.m    | 111 ---------
 lib/solvers/@PgmresSolverSetup/choose.m       |  54 -----
 lib/solvers/PgmresMultigrid.m                 | 218 ++++++++++++++++++
 lib/solvers/PgmresSolver.m                    |  87 +------
 7 files changed, 241 insertions(+), 259 deletions(-)
 delete mode 100644 lib/solvers/@PgmresSolverSetup/PgmresSolverSetup.m
 delete mode 100644 lib/solvers/@PgmresSolverSetup/choose.m
 create mode 100644 lib/solvers/PgmresMultigrid.m

diff --git a/examples/afemWithGmres.m b/examples/afemWithGmres.m
index 0a19331..52eeb67 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 496a7d6..c7f4889 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 c5aa8e5..5f29b28 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 b61e09f..0000000
--- 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 77900ed..0000000
--- 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 0000000..c8343d6
--- /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 d519fb2..9cbb8dc 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
-- 
GitLab