From 8595bbf391549f11b1a64d908ccd218c02c4c2cb Mon Sep 17 00:00:00 2001 From: Angela Date: Tue, 28 Oct 2025 11:22:44 -0400 Subject: [PATCH 1/2] adding new version of GMRES for operators of the form A = zI+K. --- chunkie/+chnk/gmresF2K.m | 88 +++++++++++++++++++++++++ devtools/test/gmresF2KconvergenceTest.m | 32 +++++++++ 2 files changed, 120 insertions(+) create mode 100644 chunkie/+chnk/gmresF2K.m create mode 100644 devtools/test/gmresF2KconvergenceTest.m diff --git a/chunkie/+chnk/gmresF2K.m b/chunkie/+chnk/gmresF2K.m new file mode 100644 index 0000000..06ee80d --- /dev/null +++ b/chunkie/+chnk/gmresF2K.m @@ -0,0 +1,88 @@ +function [x,k,g] = gmresF2K(z, K, b, x, opts) + % Solve (z*I + K) x = b + % z is a complex number + % K is a matrix or function handle + m = opts.max_iterations; + tol = opts.threshold; + n = numel(b); + if isa(K,'function_handle') + Kop = @(v) K(v); + else + Kop = @(v) K*v; + end + + r = b - z*x - Kop(x); + beta = norm(r); + bnorm = norm(b); + if bnorm==0, bnorm=1; end + if beta/bnorm <= tol, return; end + + Q = zeros(n,m+1); % Krylov basis + H = complex(zeros(m+1,m)); % Hessenberg + c = complex(zeros(m,1)); + s = complex(zeros(m,1)); + g = complex(zeros(m+1,1)); + + Q(:,1) = r/beta; g(1) = beta; + kend = m; + for k = 1:m + % Arnoldi / MGS + w = Kop(Q(:,k)); + for i = 1:k + H(i,k) = Q(:,i)'*w; + w = w - H(i,k)*Q(:,i); + end + H(k+1,k) = norm(w); + if H(k+1,k) ~= 0 + Q(:,k+1) = w/H(k+1,k); + else + Q(:,k+1) = Q(:,k); + end + + H(k,k) = H(k,k)+z; + + % Apply previous rotations + for i = 1:k-1 + Hik = H(i,k); + Hik1 = H(i+1,k); + H(i,k) = conj(c(i))*Hik - s(i)*Hik1; + H(i+1,k) = conj(s(i))*Hik + c(i)*Hik1; + end + + % New rotation to zero H(k+1,k) + [c(k),s(k),rlen] = cgivens(H(k,k), H(k+1,k)); + H(k,k) = rlen; H(k+1,k)=0; + + % Update g (||r_k|| = |g(k+1)|) + gk = g(k); + g(k) = conj(c(k))*gk - s(k)*g(k+1); + g(k+1) = conj(s(k))*gk + c(k)*g(k+1); + + % Stop if reached tolerance + if abs(g(k+1))/bnorm <= tol, kend = k; break; end + end + + if kend == m + kend = m; + end + + % Solve using least squares + y = zeros(kend,1); + for i = kend:-1:1 + ssum = g(i); + for j = i+1:kend + ssum = ssum - H(i,j)*y(j); + end + y(i) = ssum / H(i,i); + end + + x = x + Q(:,1:kend)*y; +end + +function [c,s,r] = cgivens(a,b) + % Complex Givens + % [conj(c) -s; conj(s) c]*[a;b] = [r;0] + aa = abs(a); bb = abs(b); + r = hypot(aa,bb); + if r==0, c=1; s=0; else, c=a/r; s=b/r; end +end \ No newline at end of file diff --git a/devtools/test/gmresF2KconvergenceTest.m b/devtools/test/gmresF2KconvergenceTest.m new file mode 100644 index 0000000..b34791f --- /dev/null +++ b/devtools/test/gmresF2KconvergenceTest.m @@ -0,0 +1,32 @@ +clearvars; close all; +iseed = 8675309; +rng(iseed); +npts = 1e3; +% A = zI+K +matK = (rand(npts)-0.5)*1e-8; % K is a small random matrix +z = 1i; % z is any complex number +matA = z*eye(npts) + matK; +rhs = matA * rand(npts,1); +rhs = rhs(:); + +start = tic; [sol1,flag,relres,iter1,resvec1] = gmres(matA,rhs,[],1e-15,100); t1 = toc(start); +% test new gmres +% set up parameters +opts = []; +opts.max_iterations = 1e4; +opts.threshold = 1e-100; +x = zeros(size(rhs)); +% test gmres for +[sol3,iter3,resvec3] = gmresF2K(z,matK,rhs,x, opts); +resvec3 = abs(resvec3); +resvec3 = resvec3(resvec3~=0); + +clf; +plot(log10(resvec1),Color='b',LineWidth=2); +hold on; +plot(log10(resvec3),Color='r',LineWidth=2); +legend('Native GMRES', 'GMRES for F2K operators'); + +fprintf("Native GMRES exited with flag %i (3=stalled), converged to %f digits. \n",flag,log10(resvec1(end))); +fprintf("New GMRES for F2K operators converged to %f digits. \n",log10(resvec3(end))); +assert(resvec3(end) Date: Wed, 29 Oct 2025 16:28:28 -0400 Subject: [PATCH 2/2] edit returned residue --- chunkie/+chnk/gmresF2K.m | 200 ++++++++++++++++-------- devtools/test/gmresF2KconvergenceTest.m | 91 +++++++---- 2 files changed, 199 insertions(+), 92 deletions(-) diff --git a/chunkie/+chnk/gmresF2K.m b/chunkie/+chnk/gmresF2K.m index 06ee80d..a4c3880 100644 --- a/chunkie/+chnk/gmresF2K.m +++ b/chunkie/+chnk/gmresF2K.m @@ -1,87 +1,159 @@ -function [x,k,g] = gmresF2K(z, K, b, x, opts) - % Solve (z*I + K) x = b - % z is a complex number - % K is a matrix or function handle - m = opts.max_iterations; - tol = opts.threshold; +function [x,iter,resvec,resfin] = gmresF2K(K, b, z, x, opts) +%GMRESF2K Stall-free Generalized Minimum Residual Method +% for operators of the form A*x = (z*I+K)*x = b, +% where z is a complex number, and +% K is a matrix or function handle. +% returns X: the solution, +% ITER: number of iterations, +% RESVEC: residue for all iterations, +% RESFIN: final value for residue. +% +% X = GMRES(K,B) attempts to solve the system of linear equations (I+K)*X = B +% for X. The N-by-N matrix K must be square and the right hand side +% column vector B must have length N. This uses the unrestarted +% method with MIN(N,100) total iterations. +% +% X = GMRES(KFUN,B) accepts a function handle KFUN instead of the matrix +% K. KFUN(X) accepts a vector input X and returns the matrix-vector +% product K*X. In all of the following syntaxes, you can replace K by +% KFUN. +% +% X = GMRES(KFUN,B,Z) attempts to solve the system of linear equations +% (Z*I+K)*X = B for X, where Z is a complex number. +% +% X = GMRES(KFUN,B,Z,X0) specifies the first initial guess. +% If X0 is [] then GMRESF2K uses the default, an all zero vector. +% +% X = GMRES(KFUN,B,Z,X0,OPTS) specifies the tolerance and maximum +% number of iterations of the method. +% If OPTS.TOL is [] then GMRESF2K uses the default, 1e-16. +% If OPTS.MAXIT is [] then GMRESF2K uses the default, min(N,100). +% +% Example: +% n = 21; A = gallery('wilk',n); b = sum(A,2); +% tol = 1e-12; maxit = 15; M = diag([10:-1:1 1 1:10]); +% opts = []; opts.max_iterations = 1e4; opts.threshold = 1e-100; +% x = gmresF2K(A,b,1,[],opts); +% Or, use this matrix-vector product function +% %-----------------------------------------------------------------% +% function y = kfun(x,n) +% y = [0; x(1:n-1)] + [((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x+[x(2:n); 0]; +% %-----------------------------------------------------------------% +% as inputs to GMRESF2K: +% x1 = gmresF2K(@(x)kfun(x,n),b,1,[],opts); + + if nargin < 2, error('gmresF2K: need at least K and b'); end n = numel(b); + + if nargin < 5 || isempty(opts), opts = struct; end + if ~isfield(opts,'maxit') && ~isfield(opts,'max_iterations') + opts.maxit = min(n,100); + end + if ~isfield(opts,'tol') && ~isfield(opts,'threshold') + opts.tol = 1e-16; + end + + if nargin < 4 || isempty(x), x = zeros(size(b)); end + if nargin < 3 || isempty(z), z = 1; end + if isa(K,'function_handle') Kop = @(v) K(v); else Kop = @(v) K*v; end - r = b - z*x - Kop(x); - beta = norm(r); - bnorm = norm(b); - if bnorm==0, bnorm=1; end - if beta/bnorm <= tol, return; end + if isreal(z) + usecomplex = false; + else + usecomplex = true; + end + + n = length(b); + if isempty(x), x = zeros(size(b)); end + m = opts.maxit; + tol = opts.tol; + + r = b-z*x-Kop(x); - Q = zeros(n,m+1); % Krylov basis - H = complex(zeros(m+1,m)); % Hessenberg - c = complex(zeros(m,1)); - s = complex(zeros(m,1)); - g = complex(zeros(m+1,1)); + beta = norm(r); + b_norm = norm(b); + + % Preallocate + Q = zeros(n, m+1); + H = zeros(m+1, m); + cs = zeros(m,1); + sn = zeros(m,1); + g = zeros(m+1,1); + resvec = zeros(m,1); + % Init + Q(:,1) = r / beta; + g(1) = beta; + + if b_norm < tol || abs(g(1)) / b_norm < tol + return + end - Q(:,1) = r/beta; g(1) = beta; - kend = m; for k = 1:m - % Arnoldi / MGS - w = Kop(Q(:,k)); - for i = 1:k - H(i,k) = Q(:,i)'*w; - w = w - H(i,k)*Q(:,i); - end - H(k+1,k) = norm(w); - if H(k+1,k) ~= 0 - Q(:,k+1) = w/H(k+1,k); - else - Q(:,k+1) = Q(:,k); - end - - H(k,k) = H(k,k)+z; - % Apply previous rotations - for i = 1:k-1 - Hik = H(i,k); - Hik1 = H(i+1,k); - H(i,k) = conj(c(i))*Hik - s(i)*Hik1; - H(i+1,k) = conj(s(i))*Hik + c(i)*Hik1; + % Arnoldi + [H(1:k+1, k), Q(:, k+1)] = arnoldi(K, Q, k); + H(k, k) = H(k, k) + z; + [H(1:k+1, k), cs(k), sn(k)] = apply_givens_rotation(H(1:k+1,k), cs, sn, k, usecomplex); + + % Apply Givens rotation + g(k + 1) = -sn(k) * g(k); + g(k) = cs(k) * g(k); + + % Find residue and error + %err(k) = norm(Q(:, 1:k+1)*g(k + 1)); + + yk = H(1:k, 1:k) \ g(1:k); + xk = x + Q(:, 1:k) * yk; + resvec(k) = norm(b-z*xk-Kop(xk)); + + current_relative_tol = abs(g(k + 1)) / b_norm; + if (current_relative_tol <= tol) + break; end - % New rotation to zero H(k+1,k) - [c(k),s(k),rlen] = cgivens(H(k,k), H(k+1,k)); - H(k,k) = rlen; H(k+1,k)=0; + end - % Update g (||r_k|| = |g(k+1)|) - gk = g(k); - g(k) = conj(c(k))*gk - s(k)*g(k+1); - g(k+1) = conj(s(k))*gk + c(k)*g(k+1); + y = H(1:k, 1:k) \ g(1:k); + x = x + Q(:, 1:k) * y; + iter = k; + resfin = norm(b-z*x-Kop(x)); + resvec = resvec(1:k); +end - % Stop if reached tolerance - if abs(g(k+1))/bnorm <= tol, kend = k; break; end +function [h, cs_k, sn_k] = apply_givens_rotation(h, cs, sn, k, usecomplex) + for i = 1:k-1 + temp = cs(i)*h(i) + sn(i)*h(i+1); + h(i+1) = -sn(i)*h(i) + cs(i)*h(i+1); + h(i) = temp; end - - if kend == m - kend = m; + if usecomplex + rotate = @cgivens_rotation; + else + rotate = @givens_rotation; end - - % Solve using least squares - y = zeros(kend,1); - for i = kend:-1:1 - ssum = g(i); - for j = i+1:kend - ssum = ssum - H(i,j)*y(j); - end - y(i) = ssum / H(i,i); + [cs_k, sn_k] = rotate(h(k), h(k+1)); + h(k) = cs_k*h(k) + sn_k*h(k+1); + h(k+1) = 0.0; +end + +function [cs, sn] = givens_rotation(v1, v2) + if (v1 == 0) + cs = 0; + sn = 1; + else + t = sqrt(v1^2 + v2^2); + cs = v1 / t; + sn = v2 / t; end - - x = x + Q(:,1:kend)*y; end -function [c,s,r] = cgivens(a,b) - % Complex Givens - % [conj(c) -s; conj(s) c]*[a;b] = [r;0] +function [c,s,r] = cgivens_rotation(a,b) + % Complex Givens: [conj(c) -s; conj(s) c]*[a;b] = [r;0], aa = abs(a); bb = abs(b); r = hypot(aa,bb); if r==0, c=1; s=0; else, c=a/r; s=b/r; end diff --git a/devtools/test/gmresF2KconvergenceTest.m b/devtools/test/gmresF2KconvergenceTest.m index b34791f..fa599d6 100644 --- a/devtools/test/gmresF2KconvergenceTest.m +++ b/devtools/test/gmresF2KconvergenceTest.m @@ -1,32 +1,67 @@ -clearvars; close all; +gmresF2KconvergenceTest0(); + + +function gmresF2KconvergenceTest0() + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% . . . builds a simple chunker +% and tests new gmres using chunkermat +% from the Helmholtz transmission problem +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + iseed = 8675309; rng(iseed); -npts = 1e3; -% A = zI+K -matK = (rand(npts)-0.5)*1e-8; % K is a small random matrix -z = 1i; % z is any complex number -matA = z*eye(npts) + matK; -rhs = matA * rand(npts,1); -rhs = rhs(:); - -start = tic; [sol1,flag,relres,iter1,resvec1] = gmres(matA,rhs,[],1e-15,100); t1 = toc(start); -% test new gmres -% set up parameters -opts = []; -opts.max_iterations = 1e4; -opts.threshold = 1e-100; -x = zeros(size(rhs)); -% test gmres for -[sol3,iter3,resvec3] = gmresF2K(z,matK,rhs,x, opts); -resvec3 = abs(resvec3); -resvec3 = resvec3(resvec3~=0); + +kvec = 20*[1;-1.5]; +zk = norm(kvec); + +% discretize domain +cparams = []; +cparams.eps = 1.0e-6; +cparams.nover = 0; +cparams.maxchunklen = 4.0/zk; % setting a chunk length helps when the + % frequency is known + +pref = []; +pref.k = 16; +narms = 5; +amp = 0.25; +chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); +chnkr = refine(chnkr); + +fkern = kernel('helm','c',zk,[1,-zk*1i]); +sysmatK = chunkermat(chnkr,fkern)*1e-6; +sysmatA = 0.5*eye(chnkr.k*chnkr.nch) + sysmatK; +rhs = -planewave(kvec(:),chnkr.r(:,:)); rhs = rhs(:); + +tol = 1e-40; +maxit = 200; +opts.tol = tol; +opts.maxit = maxit; +dens0 = sysmatA\rhs; +[dens1,~,~,iter1,resvec1] = gmres(sysmatA,rhs,[],tol,maxit); +[dens2,iter2,resvec2,~] = gmresF2K(sysmatK,rhs,0.5,[],opts); clf; -plot(log10(resvec1),Color='b',LineWidth=2); -hold on; -plot(log10(resvec3),Color='r',LineWidth=2); -legend('Native GMRES', 'GMRES for F2K operators'); - -fprintf("Native GMRES exited with flag %i (3=stalled), converged to %f digits. \n",flag,log10(resvec1(end))); -fprintf("New GMRES for F2K operators converged to %f digits. \n",log10(resvec3(end))); -assert(resvec3(end)