From 01641c05f162e4fef8c54013281a8907c70729a6 Mon Sep 17 00:00:00 2001 From: Ken Ho Date: Tue, 1 May 2018 08:12:10 -0700 Subject: [PATCH] Minor fixes in solve error estimation. --- hifde/test/fd_cube1.m | 4 ++-- hifde/test/fd_cube1_diag.m | 2 +- hifde/test/fd_cube1x.m | 4 ++-- hifde/test/fd_cube1x_diag.m | 2 +- hifde/test/fd_cube2.m | 4 ++-- hifde/test/fd_cube2x.m | 4 ++-- hifde/test/fd_cube3.m | 4 ++-- hifde/test/fd_cube3x.m | 4 ++-- hifde/test/fd_cube4x.m | 4 ++-- hifde/test/fd_square1.m | 4 ++-- hifde/test/fd_square1_diag.m | 2 +- hifde/test/fd_square1x.m | 4 ++-- hifde/test/fd_square1x_diag.m | 2 +- hifde/test/fd_square2.m | 4 ++-- hifde/test/fd_square2x.m | 4 ++-- hifde/test/fd_square3.m | 4 ++-- hifde/test/fd_square3x.m | 4 ++-- hifde/test/fd_square4x.m | 4 ++-- hifde/test/fd_square5x.m | 4 ++-- hifie/test/cov_cube1.m | 3 ++- hifie/test/cov_cube2.m | 3 ++- hifie/test/cov_square1.m | 3 ++- hifie/test/cov_square2.m | 3 ++- hifie/test/fd_cube.m | 5 +++-- hifie/test/fd_square.m | 5 +++-- hifie/test/ie_cube1.m | 5 +++-- hifie/test/ie_cube2.m | 3 ++- hifie/test/ie_cube2x.m | 3 ++- hifie/test/ie_square1.m | 5 +++-- hifie/test/ie_square2.m | 3 ++- hifie/test/ie_square2x.m | 3 ++- hifie/test/ie_square3.m | 2 +- hifie/test/ie_square3x.m | 2 +- ifmm/test/ie_circle.m | 2 +- ifmm/test/ie_cube.m | 2 +- ifmm/test/ie_sphere.m | 2 +- ifmm/test/ie_square.m | 2 +- mf/test/fd_cube1.m | 4 ++-- mf/test/fd_cube1_diag.m | 2 +- mf/test/fd_cube1x.m | 4 ++-- mf/test/fd_cube1x_diag.m | 2 +- mf/test/fd_cube2.m | 4 ++-- mf/test/fd_cube2x.m | 4 ++-- mf/test/fd_cube3.m | 4 ++-- mf/test/fd_cube3x.m | 4 ++-- mf/test/fd_line1x.m | 4 ++-- mf/test/fd_line1x_diag.m | 2 +- mf/test/fd_line2x.m | 4 ++-- mf/test/fd_line2x_diag.m | 2 +- mf/test/fd_square1.m | 4 ++-- mf/test/fd_square1_diag.m | 2 +- mf/test/fd_square1x.m | 4 ++-- mf/test/fd_square1x_diag.m | 2 +- mf/test/fd_square2.m | 4 ++-- mf/test/fd_square2x.m | 4 ++-- mf/test/fd_square3.m | 4 ++-- mf/test/fd_square3x.m | 4 ++-- mf/test/fd_square4x.m | 4 ++-- mf/test/fd_square5x.m | 4 ++-- rskel/test/cov_line1.m | 18 +++++++++++------- rskel/test/cov_line2.m | 18 +++++++++++------- rskel/test/fd_square.m | 16 ++++++++++------ rskel/test/ie_circle.m | 16 ++++++++++------ rskel/test/ie_square1.m | 16 ++++++++++------ rskel/test/ie_square2.m | 14 +++++++++----- rskel/test/ie_square3.m | 2 +- rskelf/test/cov_circle1.m | 3 ++- rskelf/test/cov_circle2.m | 3 ++- rskelf/test/cov_line1.m | 3 ++- rskelf/test/cov_line2.m | 3 ++- rskelf/test/cov_square1.m | 3 ++- rskelf/test/cov_square2.m | 3 ++- rskelf/test/fd_cube.m | 5 +++-- rskelf/test/fd_square.m | 5 +++-- rskelf/test/ie_circle.m | 3 ++- rskelf/test/ie_cube1.m | 5 +++-- rskelf/test/ie_cube2.m | 3 ++- rskelf/test/ie_square1.m | 5 +++-- rskelf/test/ie_square2.m | 3 ++- rskelf/test/ie_square3.m | 2 +- 80 files changed, 199 insertions(+), 150 deletions(-) diff --git a/hifde/test/fd_cube1.m b/hifde/test/fd_cube1.m index 42202fd..e958407 100644 --- a/hifde/test/fd_cube1.m +++ b/hifde/test/fd_cube1.m @@ -100,7 +100,7 @@ function fd_cube1(n,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -128,7 +128,7 @@ function fd_cube1(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_cube1_diag.m b/hifde/test/fd_cube1_diag.m index 5e99809..87606ba 100644 --- a/hifde/test/fd_cube1_diag.m +++ b/hifde/test/fd_cube1_diag.m @@ -103,7 +103,7 @@ function fd_cube1_diag(n,occ,rank_or_tol,skip,symm,spdiag) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % prepare for diagonal extracation diff --git a/hifde/test/fd_cube1x.m b/hifde/test/fd_cube1x.m index 2b414d3..2e9ffb0 100644 --- a/hifde/test/fd_cube1x.m +++ b/hifde/test/fd_cube1x.m @@ -103,7 +103,7 @@ function fd_cube1x(n,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -131,7 +131,7 @@ function fd_cube1x(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_cube1x_diag.m b/hifde/test/fd_cube1x_diag.m index 6a57991..b74057f 100644 --- a/hifde/test/fd_cube1x_diag.m +++ b/hifde/test/fd_cube1x_diag.m @@ -106,7 +106,7 @@ function fd_cube1x_diag(n,occ,rank_or_tol,skip,symm,spdiag) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % prepare for diagonal extracation diff --git a/hifde/test/fd_cube2.m b/hifde/test/fd_cube2.m index 4c61686..55b7f97 100644 --- a/hifde/test/fd_cube2.m +++ b/hifde/test/fd_cube2.m @@ -127,7 +127,7 @@ function fd_cube2(n,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -155,7 +155,7 @@ function fd_cube2(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_cube2x.m b/hifde/test/fd_cube2x.m index 0c09647..2c9c0c3 100644 --- a/hifde/test/fd_cube2x.m +++ b/hifde/test/fd_cube2x.m @@ -130,7 +130,7 @@ function fd_cube2x(n,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -158,7 +158,7 @@ function fd_cube2x(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_cube3.m b/hifde/test/fd_cube3.m index 48b3a5d..602ea95 100644 --- a/hifde/test/fd_cube3.m +++ b/hifde/test/fd_cube3.m @@ -103,7 +103,7 @@ function fd_cube3(n,k,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -111,7 +111,7 @@ function fd_cube3(n,k,occ,rank_or_tol,skip,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_cube3x.m b/hifde/test/fd_cube3x.m index 146340b..3e1c232 100644 --- a/hifde/test/fd_cube3x.m +++ b/hifde/test/fd_cube3x.m @@ -106,7 +106,7 @@ function fd_cube3x(n,k,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -114,7 +114,7 @@ function fd_cube3x(n,k,occ,rank_or_tol,skip,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_cube4x.m b/hifde/test/fd_cube4x.m index 633dbe9..bf4d918 100644 --- a/hifde/test/fd_cube4x.m +++ b/hifde/test/fd_cube4x.m @@ -140,7 +140,7 @@ function fd_cube4x(n,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -168,7 +168,7 @@ function fd_cube4x(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_square1.m b/hifde/test/fd_square1.m index 87bf12e..2f705fd 100644 --- a/hifde/test/fd_square1.m +++ b/hifde/test/fd_square1.m @@ -90,7 +90,7 @@ function fd_square1(n,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -118,7 +118,7 @@ function fd_square1(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_square1_diag.m b/hifde/test/fd_square1_diag.m index 026767c..573beee 100644 --- a/hifde/test/fd_square1_diag.m +++ b/hifde/test/fd_square1_diag.m @@ -93,7 +93,7 @@ function fd_square1_diag(n,occ,rank_or_tol,skip,symm,spdiag) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % prepare for diagonal extracation diff --git a/hifde/test/fd_square1x.m b/hifde/test/fd_square1x.m index 289ef48..46dd6ed 100644 --- a/hifde/test/fd_square1x.m +++ b/hifde/test/fd_square1x.m @@ -93,7 +93,7 @@ function fd_square1x(n,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -121,7 +121,7 @@ function fd_square1x(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_square1x_diag.m b/hifde/test/fd_square1x_diag.m index 2b8325c..cebc65c 100644 --- a/hifde/test/fd_square1x_diag.m +++ b/hifde/test/fd_square1x_diag.m @@ -96,7 +96,7 @@ function fd_square1x_diag(n,occ,rank_or_tol,skip,symm,spdiag) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % prepare for diagonal extracation diff --git a/hifde/test/fd_square2.m b/hifde/test/fd_square2.m index 01948cb..bf3b6ee 100644 --- a/hifde/test/fd_square2.m +++ b/hifde/test/fd_square2.m @@ -112,7 +112,7 @@ function fd_square2(n,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -140,7 +140,7 @@ function fd_square2(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_square2x.m b/hifde/test/fd_square2x.m index 9c43a8d..370ddda 100644 --- a/hifde/test/fd_square2x.m +++ b/hifde/test/fd_square2x.m @@ -115,7 +115,7 @@ function fd_square2x(n,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -143,7 +143,7 @@ function fd_square2x(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_square3.m b/hifde/test/fd_square3.m index 93dcf4e..31c3210 100644 --- a/hifde/test/fd_square3.m +++ b/hifde/test/fd_square3.m @@ -93,7 +93,7 @@ function fd_square3(n,k,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -101,7 +101,7 @@ function fd_square3(n,k,occ,rank_or_tol,skip,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_square3x.m b/hifde/test/fd_square3x.m index bd71782..944f09e 100644 --- a/hifde/test/fd_square3x.m +++ b/hifde/test/fd_square3x.m @@ -96,7 +96,7 @@ function fd_square3x(n,k,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -104,7 +104,7 @@ function fd_square3x(n,k,occ,rank_or_tol,skip,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_square4x.m b/hifde/test/fd_square4x.m index a4eaf31..0521fd2 100644 --- a/hifde/test/fd_square4x.m +++ b/hifde/test/fd_square4x.m @@ -116,7 +116,7 @@ function fd_square4x(n,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -144,7 +144,7 @@ function fd_square4x(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifde/test/fd_square5x.m b/hifde/test/fd_square5x.m index 9ab36bb..d85fad0 100644 --- a/hifde/test/fd_square5x.m +++ b/hifde/test/fd_square5x.m @@ -143,7 +143,7 @@ function fd_square5x(n,occ,rank_or_tol,skip,symm) tic Y = hifde_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifde_sv(F,x)),@(x)(x - hifde_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -171,7 +171,7 @@ function fd_square5x(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifde_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifde_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifie/test/cov_cube1.m b/hifie/test/cov_cube1.m index f8f99f6..fe19f0d 100644 --- a/hifie/test/cov_cube1.m +++ b/hifie/test/cov_cube1.m @@ -85,7 +85,8 @@ function cov_cube1(n,occ,p,rank_or_tol,skip,symm,noise,scale,spdiag) tic hifie_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))), ... + @(x)(x - hifie_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') diff --git a/hifie/test/cov_cube2.m b/hifie/test/cov_cube2.m index efdbd1c..9dac7a9 100644 --- a/hifie/test/cov_cube2.m +++ b/hifie/test/cov_cube2.m @@ -85,7 +85,8 @@ function cov_cube2(n,occ,p,rank_or_tol,skip,symm,noise,scale,spdiag) tic hifie_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))), ... + @(x)(x - hifie_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') diff --git a/hifie/test/cov_square1.m b/hifie/test/cov_square1.m index 0f180b6..4b58add 100644 --- a/hifie/test/cov_square1.m +++ b/hifie/test/cov_square1.m @@ -80,7 +80,8 @@ function cov_square1(n,occ,p,rank_or_tol,skip,symm,noise,scale,spdiag) tic hifie_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))), ... + @(x)(x - hifie_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') diff --git a/hifie/test/cov_square2.m b/hifie/test/cov_square2.m index 87d6632..c64fc8b 100644 --- a/hifie/test/cov_square2.m +++ b/hifie/test/cov_square2.m @@ -80,7 +80,8 @@ function cov_square2(n,occ,p,rank_or_tol,skip,symm,noise,scale,spdiag) tic hifie_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))), ... + @(x)(x - hifie_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') diff --git a/hifie/test/fd_cube.m b/hifie/test/fd_cube.m index 14cec8b..be95e84 100644 --- a/hifie/test/fd_cube.m +++ b/hifie/test/fd_cube.m @@ -81,7 +81,8 @@ function fd_cube(n,occ,rank_or_tol,skip,symm) tic Y = hifie_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifie_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifie_sv(F,x)), ... + @(x)(x - hifie_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run CG @@ -89,7 +90,7 @@ function fd_cube(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifie_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifie_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifie/test/fd_square.m b/hifie/test/fd_square.m index 39dacad..561566c 100644 --- a/hifie/test/fd_square.m +++ b/hifie/test/fd_square.m @@ -75,7 +75,8 @@ function fd_square(n,occ,rank_or_tol,skip,symm) tic Y = hifie_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*hifie_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*hifie_sv(F,x)), ... + @(x)(x - hifie_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run CG @@ -83,7 +84,7 @@ function fd_square(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(hifie_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)hifie_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/hifie/test/ie_cube1.m b/hifie/test/ie_cube1.m index d7e0ed5..3c89e66 100644 --- a/hifie/test/ie_cube1.m +++ b/hifie/test/ie_cube1.m @@ -77,7 +77,8 @@ function ie_cube1(n,occ,p,rank_or_tol,skip,symm) tic Y = hifie_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))), ... + @(x)(x - hifie_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -85,7 +86,7 @@ function ie_cube1(n,occ,p,rank_or_tol,skip,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)(hifie_sv(F,x))); + [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)hifie_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - mv(Z))/norm(X); diff --git a/hifie/test/ie_cube2.m b/hifie/test/ie_cube2.m index 3684e0c..2b099e3 100644 --- a/hifie/test/ie_cube2.m +++ b/hifie/test/ie_cube2.m @@ -77,7 +77,8 @@ function ie_cube2(n,occ,p,rank_or_tol,skip,symm) tic hifie_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))), ... + @(x)(x - hifie_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) end diff --git a/hifie/test/ie_cube2x.m b/hifie/test/ie_cube2x.m index e149519..9d39457 100644 --- a/hifie/test/ie_cube2x.m +++ b/hifie/test/ie_cube2x.m @@ -77,7 +77,8 @@ function ie_cube2x(n,occ,p,rank_or_tol,skip,symm) tic hifie_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))), ... + @(x)(x - hifie_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) end diff --git a/hifie/test/ie_square1.m b/hifie/test/ie_square1.m index d5cded3..fc2f634 100644 --- a/hifie/test/ie_square1.m +++ b/hifie/test/ie_square1.m @@ -71,7 +71,8 @@ function ie_square1(n,occ,p,rank_or_tol,skip,symm) tic Y = hifie_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))), ... + @(x)(x - hifie_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -79,7 +80,7 @@ function ie_square1(n,occ,p,rank_or_tol,skip,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)(hifie_sv(F,x))); + [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)hifie_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - mv(Z))/norm(X); diff --git a/hifie/test/ie_square2.m b/hifie/test/ie_square2.m index 686759d..a6abefa 100644 --- a/hifie/test/ie_square2.m +++ b/hifie/test/ie_square2.m @@ -71,7 +71,8 @@ function ie_square2(n,occ,p,rank_or_tol,skip,symm) tic hifie_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))), ... + @(x)(x - hifie_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) end diff --git a/hifie/test/ie_square2x.m b/hifie/test/ie_square2x.m index 95bc32c..ed839aa 100644 --- a/hifie/test/ie_square2x.m +++ b/hifie/test/ie_square2x.m @@ -71,7 +71,8 @@ function ie_square2x(n,occ,p,rank_or_tol,skip,symm) tic hifie_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(hifie_sv(F,x))), ... + @(x)(x - hifie_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) end diff --git a/hifie/test/ie_square3.m b/hifie/test/ie_square3.m index a0cfd86..8c3ddb4 100644 --- a/hifie/test/ie_square3.m +++ b/hifie/test/ie_square3.m @@ -96,7 +96,7 @@ function ie_square3(n,k,occ,p,rank_or_tol,skip,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)(hifie_sv(F,x))); + [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)hifie_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - mv(Z))/norm(X); diff --git a/hifie/test/ie_square3x.m b/hifie/test/ie_square3x.m index 85c8415..e17ca6c 100644 --- a/hifie/test/ie_square3x.m +++ b/hifie/test/ie_square3x.m @@ -96,7 +96,7 @@ function ie_square3x(n,k,occ,p,rank_or_tol,skip,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)(hifie_sv(F,x))); + [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)hifie_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - mv(Z))/norm(X); diff --git a/ifmm/test/ie_circle.m b/ifmm/test/ie_circle.m index 78817be..33c5914 100644 --- a/ifmm/test/ie_circle.m +++ b/ifmm/test/ie_circle.m @@ -66,7 +66,7 @@ function ie_circle(n,occ,p,rank_or_tol,near,store,symm) % solve for surface density tic - [X,~,~,iter] = gmres(@(x)(ifmm_mv(F,x,Afun)),B,[],1e-12,32); + [X,~,~,iter] = gmres(@(x)ifmm_mv(F,x,Afun),B,[],1e-12,32); t = toc; e = norm(B - mv(X))/norm(B); fprintf('gmres: %10.4e / %4d / %10.4e (s)\n',e,iter(2),t) diff --git a/ifmm/test/ie_cube.m b/ifmm/test/ie_cube.m index b212bb7..558a3d8 100644 --- a/ifmm/test/ie_cube.m +++ b/ifmm/test/ie_cube.m @@ -78,7 +78,7 @@ function ie_cube(n,occ,p,rank_or_tol,near,store,symm) % run GMRES tic - [Y,~,~,iter] = gmres(@(x)(ifmm_mv(F,x,Afun)),X,[],1e-12,32); + [Y,~,~,iter] = gmres(@(x)ifmm_mv(F,x,Afun),X,[],1e-12,32); t = toc; e = norm(X - mv(Y))/norm(X); fprintf('gmres: %10.4e / %4d / %10.4e (s)\n',e,iter(2),t) diff --git a/ifmm/test/ie_sphere.m b/ifmm/test/ie_sphere.m index edbac20..83b2d56 100644 --- a/ifmm/test/ie_sphere.m +++ b/ifmm/test/ie_sphere.m @@ -152,7 +152,7 @@ function ie_sphere(n,nquad,occ,p,rank_or_tol,near,store) % solve for surface density tic - [X,~,~,iter] = gmres(@(x)(ifmm_mv(F,x,Afun)),B,[],1e-6,32); + [X,~,~,iter] = gmres(@(x)ifmm_mv(F,x,Afun),B,[],1e-6,32); t = toc; r = randperm(N); r = r(1:min(N,128)); diff --git a/ifmm/test/ie_square.m b/ifmm/test/ie_square.m index b6e130a..8d2bd92 100644 --- a/ifmm/test/ie_square.m +++ b/ifmm/test/ie_square.m @@ -72,7 +72,7 @@ function ie_square(n,occ,p,rank_or_tol,near,store,symm) % run GMRES tic - [Y,~,~,iter] = gmres(@(x)(ifmm_mv(F,x,Afun)),X,[],1e-12,32); + [Y,~,~,iter] = gmres(@(x)ifmm_mv(F,x,Afun),X,[],1e-12,32); t = toc; e = norm(X - mv(Y))/norm(X); fprintf('gmres: %10.4e / %4d / %10.4e (s)\n',e,iter(2),t) diff --git a/mf/test/fd_cube1.m b/mf/test/fd_cube1.m index e909cc0..b6b2ecc 100644 --- a/mf/test/fd_cube1.m +++ b/mf/test/fd_cube1.m @@ -94,7 +94,7 @@ function fd_cube1(n,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -122,7 +122,7 @@ function fd_cube1(n,occ,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_cube1_diag.m b/mf/test/fd_cube1_diag.m index 8473f28..e0a29ec 100644 --- a/mf/test/fd_cube1_diag.m +++ b/mf/test/fd_cube1_diag.m @@ -97,7 +97,7 @@ function fd_cube1_diag(n,occ,symm,spdiag) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % prepare for diagonal extracation diff --git a/mf/test/fd_cube1x.m b/mf/test/fd_cube1x.m index 1965eb6..5e0d227 100644 --- a/mf/test/fd_cube1x.m +++ b/mf/test/fd_cube1x.m @@ -97,7 +97,7 @@ function fd_cube1x(n,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -125,7 +125,7 @@ function fd_cube1x(n,occ,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_cube1x_diag.m b/mf/test/fd_cube1x_diag.m index d75864d..6f21b5d 100644 --- a/mf/test/fd_cube1x_diag.m +++ b/mf/test/fd_cube1x_diag.m @@ -100,7 +100,7 @@ function fd_cube1x_diag(n,occ,symm,spdiag) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % prepare for diagonal extracation diff --git a/mf/test/fd_cube2.m b/mf/test/fd_cube2.m index 25631a5..e171b54 100644 --- a/mf/test/fd_cube2.m +++ b/mf/test/fd_cube2.m @@ -121,7 +121,7 @@ function fd_cube2(n,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -149,7 +149,7 @@ function fd_cube2(n,occ,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_cube2x.m b/mf/test/fd_cube2x.m index ce5235c..62886e7 100644 --- a/mf/test/fd_cube2x.m +++ b/mf/test/fd_cube2x.m @@ -124,7 +124,7 @@ function fd_cube2x(n,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -152,7 +152,7 @@ function fd_cube2x(n,occ,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_cube3.m b/mf/test/fd_cube3.m index 6ef0b98..68fd050 100644 --- a/mf/test/fd_cube3.m +++ b/mf/test/fd_cube3.m @@ -97,7 +97,7 @@ function fd_cube3(n,k,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -105,7 +105,7 @@ function fd_cube3(n,k,occ,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_cube3x.m b/mf/test/fd_cube3x.m index 2e10e25..1908005 100644 --- a/mf/test/fd_cube3x.m +++ b/mf/test/fd_cube3x.m @@ -100,7 +100,7 @@ function fd_cube3x(n,k,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -108,7 +108,7 @@ function fd_cube3x(n,k,occ,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_line1x.m b/mf/test/fd_line1x.m index f9e16b2..a82cfde 100644 --- a/mf/test/fd_line1x.m +++ b/mf/test/fd_line1x.m @@ -75,7 +75,7 @@ function fd_line1x(n,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -103,7 +103,7 @@ function fd_line1x(n,occ,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_line1x_diag.m b/mf/test/fd_line1x_diag.m index e0e73e6..d991ac4 100644 --- a/mf/test/fd_line1x_diag.m +++ b/mf/test/fd_line1x_diag.m @@ -78,7 +78,7 @@ function fd_line1x_diag(n,occ,symm,spdiag) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % prepare for diagonal extracation diff --git a/mf/test/fd_line2x.m b/mf/test/fd_line2x.m index 7fd6829..114bb5b 100644 --- a/mf/test/fd_line2x.m +++ b/mf/test/fd_line2x.m @@ -87,7 +87,7 @@ function fd_line2x(n,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -115,7 +115,7 @@ function fd_line2x(n,occ,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_line2x_diag.m b/mf/test/fd_line2x_diag.m index ffb808c..40c90bd 100644 --- a/mf/test/fd_line2x_diag.m +++ b/mf/test/fd_line2x_diag.m @@ -90,7 +90,7 @@ function fd_line2x_diag(n,occ,symm,spdiag) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % prepare for diagonal extracation diff --git a/mf/test/fd_square1.m b/mf/test/fd_square1.m index 046a7d7..3060190 100644 --- a/mf/test/fd_square1.m +++ b/mf/test/fd_square1.m @@ -84,7 +84,7 @@ function fd_square1(n,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -112,7 +112,7 @@ function fd_square1(n,occ,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_square1_diag.m b/mf/test/fd_square1_diag.m index a23312c..05e0104 100644 --- a/mf/test/fd_square1_diag.m +++ b/mf/test/fd_square1_diag.m @@ -87,7 +87,7 @@ function fd_square1_diag(n,occ,symm,spdiag) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % prepare for diagonal extracation diff --git a/mf/test/fd_square1x.m b/mf/test/fd_square1x.m index ea7f825..2867644 100644 --- a/mf/test/fd_square1x.m +++ b/mf/test/fd_square1x.m @@ -87,7 +87,7 @@ function fd_square1x(n,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -115,7 +115,7 @@ function fd_square1x(n,occ,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_square1x_diag.m b/mf/test/fd_square1x_diag.m index a10b308..ea7b674 100644 --- a/mf/test/fd_square1x_diag.m +++ b/mf/test/fd_square1x_diag.m @@ -90,7 +90,7 @@ function fd_square1x_diag(n,occ,symm,spdiag) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % prepare for diagonal extracation diff --git a/mf/test/fd_square2.m b/mf/test/fd_square2.m index 4eb53ea..5d3454f 100644 --- a/mf/test/fd_square2.m +++ b/mf/test/fd_square2.m @@ -106,7 +106,7 @@ function fd_square2(n,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -134,7 +134,7 @@ function fd_square2(n,occ,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_square2x.m b/mf/test/fd_square2x.m index 28e68fe..049fc71 100644 --- a/mf/test/fd_square2x.m +++ b/mf/test/fd_square2x.m @@ -109,7 +109,7 @@ function fd_square2x(n,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -137,7 +137,7 @@ function fd_square2x(n,occ,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_square3.m b/mf/test/fd_square3.m index 8744df1..e7018cf 100644 --- a/mf/test/fd_square3.m +++ b/mf/test/fd_square3.m @@ -87,7 +87,7 @@ function fd_square3(n,k,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -95,7 +95,7 @@ function fd_square3(n,k,occ,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_square3x.m b/mf/test/fd_square3x.m index 5bf3224..b1e6abf 100644 --- a/mf/test/fd_square3x.m +++ b/mf/test/fd_square3x.m @@ -90,7 +90,7 @@ function fd_square3x(n,k,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -98,7 +98,7 @@ function fd_square3x(n,k,occ,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = gmres(@(x)(A*x),X,[],1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_square4x.m b/mf/test/fd_square4x.m index 31ab794..553085f 100644 --- a/mf/test/fd_square4x.m +++ b/mf/test/fd_square4x.m @@ -110,7 +110,7 @@ function fd_square4x(n,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -138,7 +138,7 @@ function fd_square4x(n,occ,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/mf/test/fd_square5x.m b/mf/test/fd_square5x.m index a4b7b0c..c0ee14f 100644 --- a/mf/test/fd_square5x.m +++ b/mf/test/fd_square5x.m @@ -137,7 +137,7 @@ function fd_square5x(n,occ,symm) tic Y = mf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*mf_sv(F,x)),@(x)(x - mf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') @@ -165,7 +165,7 @@ function fd_square5x(n,occ,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(mf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)mf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/rskel/test/cov_line1.m b/rskel/test/cov_line1.m index 0fd2722..b283757 100644 --- a/rskel/test/cov_line1.m +++ b/rskel/test/cov_line1.m @@ -71,7 +71,7 @@ function cov_line1(n,occ,p,rank_or_tol,symm,noise,scale) spmem = (spmem + w.bytes)/1e6; end fprintf('lu/ldl: %10.4e (s) / %6.2f (MB)\n',t,spmem) - sv = @(x)sv2(FA,x); + sv = @(x,trans)sv2(FA,x,trans); % set up FFT multiplication a = Afun(1:n,1); @@ -95,9 +95,9 @@ function cov_line1(n,occ,p,rank_or_tol,symm,noise,scale) % NORM(INV(A) - INV(F))/NORM(INV(A)) <= NORM(I - A*INV(F)) tic - sv(X); + sv(X,'n'); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(sv(x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(sv(x,'n'))),@(x)(x - sv(mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % prepare for selected inversion @@ -150,7 +150,7 @@ function cov_line1(n,occ,p,rank_or_tol,symm,noise,scale) end end t = toc/m; - Y = sv(X); + Y = sv(X,'n'); for i = 1:m T(i) = Y(r(i,1),i); end @@ -171,7 +171,7 @@ function cov_line1(n,occ,p,rank_or_tol,symm,noise,scale) end end t = toc/m; - Y = sv(X); + Y = sv(X,'n'); for i = 1:m T(i) = Y(r(i,2),i); end @@ -211,11 +211,15 @@ function cov_line1(n,occ,p,rank_or_tol,symm,noise,scale) end % sparse LU solve -function Y = sv2(F,X) +function Y = sv2(F,X,trans) N = size(X,1); X = [X; zeros(size(F.L,1)-N,size(X,2))]; if F.lu - Y = F.U\(F.L\X); + if strcmpi(trans,'n') + Y = F.U\(F.L\X); + else + Y = F.L'\(F.U'\X); + end else Y = F.P*(F.L'\(F.D\(F.L\(F.P'*X)))); end diff --git a/rskel/test/cov_line2.m b/rskel/test/cov_line2.m index 0689ca4..36fc71c 100644 --- a/rskel/test/cov_line2.m +++ b/rskel/test/cov_line2.m @@ -71,7 +71,7 @@ function cov_line2(n,occ,p,rank_or_tol,symm,noise,scale) spmem = (spmem + w.bytes)/1e6; end fprintf('lu/ldl: %10.4e (s) / %6.2f (MB)\n',t,spmem) - sv = @(x)sv2(FA,x); + sv = @(x,trans)sv2(FA,x,trans); % set up FFT multiplication a = Afun(1:n,1); @@ -95,9 +95,9 @@ function cov_line2(n,occ,p,rank_or_tol,symm,noise,scale) % NORM(INV(A) - INV(F))/NORM(INV(A)) <= NORM(I - A*INV(F)) tic - sv(X); + sv(X,'n'); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(sv(x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(sv(x,'n'))),@(x)(x - sv(mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % prepare for selected inversion @@ -150,7 +150,7 @@ function cov_line2(n,occ,p,rank_or_tol,symm,noise,scale) end end t = toc/m; - Y = sv(X); + Y = sv(X,'n'); for i = 1:m T(i) = Y(r(i,1),i); end @@ -171,7 +171,7 @@ function cov_line2(n,occ,p,rank_or_tol,symm,noise,scale) end end t = toc/m; - Y = sv(X); + Y = sv(X,'n'); for i = 1:m T(i) = Y(r(i,2),i); end @@ -211,11 +211,15 @@ function cov_line2(n,occ,p,rank_or_tol,symm,noise,scale) end % sparse LU solve -function Y = sv2(F,X) +function Y = sv2(F,X,trans) N = size(X,1); X = [X; zeros(size(F.L,1)-N,size(X,2))]; if F.lu - Y = F.U\(F.L\X); + if strcmpi(trans,'n') + Y = F.U\(F.L\X); + else + Y = F.L'\(F.U'\X); + end else Y = F.P*(F.L'\(F.D\(F.L\(F.P'*X)))); end diff --git a/rskel/test/fd_square.m b/rskel/test/fd_square.m index a343481..806489a 100644 --- a/rskel/test/fd_square.m +++ b/rskel/test/fd_square.m @@ -87,7 +87,7 @@ function fd_square(n,occ,rank_or_tol,symm) spmem = (spmem + w.bytes)/1e6; end fprintf('lu/ldl: %10.4e (s) / %6.2f (MB)\n',t,spmem) - sv = @(x)sv2(FS,x); + sv = @(x,trans)sv2(FS,x,trans); % test accuracy using randomized power method X = rand(N,1); @@ -103,9 +103,9 @@ function fd_square(n,occ,rank_or_tol,symm) % NORM(INV(A) - INV(F))/NORM(INV(A)) <= NORM(I - A*INV(F)) tic - Y = sv(X); + Y = sv(X,'n'); t = toc; - [e,niter] = snorm(N,@(x)(x - A*sv(x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*sv(x,'n')),@(x)(x - sv(A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run CG @@ -113,7 +113,7 @@ function fd_square(n,occ,rank_or_tol,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,sv); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)sv(x,'n')); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); @@ -152,11 +152,15 @@ function fd_square(n,occ,rank_or_tol,symm) end % sparse LU solve -function Y = sv2(F,X) +function Y = sv2(F,X,trans) N = size(X,1); X = [X; zeros(size(F.L,1)-N,size(X,2))]; if F.lu - Y = F.U\(F.L\X); + if strcmpi(trans,'n') + Y = F.U\(F.L\X); + else + Y = F.L'\(F.U'\X); + end else Y = F.P*(F.L'\(F.D\(F.L\(F.P'*X)))); end diff --git a/rskel/test/ie_circle.m b/rskel/test/ie_circle.m index 782a0ff..cf625fd 100644 --- a/rskel/test/ie_circle.m +++ b/rskel/test/ie_circle.m @@ -66,7 +66,7 @@ function ie_circle(n,occ,p,rank_or_tol,symm) spmem = (spmem + w.bytes)/1e6; end fprintf('lu/ldl: %10.4e (s) / %6.2f (MB)\n',t,spmem) - sv = @(x)sv2(FA,x); + sv = @(x,trans)sv2(FA,x,trans); % set up FFT multiplication G = fft(Afun(1:N,1)); @@ -86,9 +86,9 @@ function ie_circle(n,occ,p,rank_or_tol,symm) % NORM(INV(A) - INV(F))/NORM(INV(A)) <= NORM(I - A*INV(F)) tic - sv(X); + sv(X,'n'); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(sv(x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(sv(x,'n'))),@(x)(x - sv(mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % generate field due to exterior sources @@ -99,7 +99,7 @@ function ie_circle(n,occ,p,rank_or_tol,symm) B = Kfun(x,src,'s')*q; % solve for surface density - X = sv(B); + X = sv(B,'n'); % evaluate field at interior targets trg = 0.5*[cos(theta); sin(theta)]; @@ -155,11 +155,15 @@ function ie_circle(n,occ,p,rank_or_tol,symm) end % sparse LU solve -function Y = sv2(F,X) +function Y = sv2(F,X,trans) N = size(X,1); X = [X; zeros(size(F.L,1)-N,size(X,2))]; if F.lu - Y = F.U\(F.L\X); + if strcmpi(trans,'n') + Y = F.U\(F.L\X); + else + Y = F.L'\(F.U'\X); + end else Y = F.P*(F.L'\(F.D\(F.L\(F.P'*X)))); end diff --git a/rskel/test/ie_square1.m b/rskel/test/ie_square1.m index 9492aa7..ad1f9f2 100644 --- a/rskel/test/ie_square1.m +++ b/rskel/test/ie_square1.m @@ -71,7 +71,7 @@ function ie_square1(n,occ,p,rank_or_tol,symm) spmem = (spmem + w.bytes)/1e6; end fprintf('lu/ldl: %10.4e (s) / %6.2f (MB)\n',t,spmem) - sv = @(x)sv2(FA,x); + sv = @(x,trans)sv2(FA,x,trans); % set up FFT multiplication a = reshape(Afun(1:N,1),n,n); @@ -99,9 +99,9 @@ function ie_square1(n,occ,p,rank_or_tol,symm) % NORM(INV(A) - INV(F))/NORM(INV(A)) <= NORM(I - A*INV(F)) tic - Y = sv(X); + Y = sv(X,'n'); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(sv(x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(sv(x,'n'))),@(x)(x - sv(mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -109,7 +109,7 @@ function ie_square1(n,occ,p,rank_or_tol,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,sv); + [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)sv(x,'n')); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - mv(Z))/norm(X); @@ -158,11 +158,15 @@ function ie_square1(n,occ,p,rank_or_tol,symm) end % sparse LU solve -function Y = sv2(F,X) +function Y = sv2(F,X,trans) N = size(X,1); X = [X; zeros(size(F.L,1)-N,size(X,2))]; if F.lu - Y = F.U\(F.L\X); + if strcmpi(trans,'n') + Y = F.U\(F.L\X); + else + Y = F.L'\(F.U'\X); + end else Y = F.P*(F.L'\(F.D\(F.L\(F.P'*X)))); end diff --git a/rskel/test/ie_square2.m b/rskel/test/ie_square2.m index c7c7379..df8e557 100644 --- a/rskel/test/ie_square2.m +++ b/rskel/test/ie_square2.m @@ -71,7 +71,7 @@ function ie_square2(n,occ,p,rank_or_tol,symm) spmem = (spmem + w.bytes)/1e6; end fprintf('lu/ldl: %10.4e (s) / %6.2f (MB)\n',t,spmem) - sv = @(x)sv2(FA,x); + sv = @(x,trans)sv2(FA,x,trans); % set up FFT multiplication a = reshape(Afun(1:N,1),n,n); @@ -99,9 +99,9 @@ function ie_square2(n,occ,p,rank_or_tol,symm) % NORM(INV(A) - INV(F))/NORM(INV(A)) <= NORM(I - A*INV(F)) tic - sv(X); + sv(X,'n'); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(sv(x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(sv(x,'n'))),@(x)(x - sv(mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) end @@ -146,11 +146,15 @@ function ie_square2(n,occ,p,rank_or_tol,symm) end % sparse LU solve -function Y = sv2(F,X) +function Y = sv2(F,X,trans) N = size(X,1); X = [X; zeros(size(F.L,1)-N,size(X,2))]; if F.lu - Y = F.U\(F.L\X); + if strcmpi(trans,'n') + Y = F.U\(F.L\X); + else + Y = F.L'\(F.U'\X); + end else Y = F.P*(F.L'\(F.D\(F.L\(F.P'*X)))); end diff --git a/rskel/test/ie_square3.m b/rskel/test/ie_square3.m index b6f0aef..c2ea8c6 100644 --- a/rskel/test/ie_square3.m +++ b/rskel/test/ie_square3.m @@ -115,7 +115,7 @@ function ie_square3(n,k,occ,p,rank_or_tol,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)(sv(x,'n'))); + [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)sv(x,'n')); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - mv(Z))/norm(X); diff --git a/rskelf/test/cov_circle1.m b/rskelf/test/cov_circle1.m index b638a39..43c1d04 100644 --- a/rskelf/test/cov_circle1.m +++ b/rskelf/test/cov_circle1.m @@ -68,7 +68,8 @@ function cov_circle1(n,occ,p,rank_or_tol,symm,noise,scale,spdiag) tic rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))), ... + @(x)(x - rskelf_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') diff --git a/rskelf/test/cov_circle2.m b/rskelf/test/cov_circle2.m index 2fcfeed..82d6816 100644 --- a/rskelf/test/cov_circle2.m +++ b/rskelf/test/cov_circle2.m @@ -68,7 +68,8 @@ function cov_circle2(n,occ,p,rank_or_tol,symm,noise,scale,spdiag) tic rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))), ... + @(x)(x - rskelf_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') diff --git a/rskelf/test/cov_line1.m b/rskelf/test/cov_line1.m index ba0ab37..7673d6c 100644 --- a/rskelf/test/cov_line1.m +++ b/rskelf/test/cov_line1.m @@ -67,7 +67,8 @@ function cov_line1(n,occ,p,rank_or_tol,symm,noise,scale,spdiag) tic rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))), ... + @(x)(x - rskelf_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') diff --git a/rskelf/test/cov_line2.m b/rskelf/test/cov_line2.m index e5f26df..32abcad 100644 --- a/rskelf/test/cov_line2.m +++ b/rskelf/test/cov_line2.m @@ -67,7 +67,8 @@ function cov_line2(n,occ,p,rank_or_tol,symm,noise,scale,spdiag) tic rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))), ... + @(x)(x - rskelf_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') diff --git a/rskelf/test/cov_square1.m b/rskelf/test/cov_square1.m index 34db8cf..1bb86d3 100644 --- a/rskelf/test/cov_square1.m +++ b/rskelf/test/cov_square1.m @@ -77,7 +77,8 @@ function cov_square1(n,occ,p,rank_or_tol,symm,noise,scale,spdiag) tic rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))), ... + @(x)(x - rskelf_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') diff --git a/rskelf/test/cov_square2.m b/rskelf/test/cov_square2.m index f2d44d8..98b27b9 100644 --- a/rskelf/test/cov_square2.m +++ b/rskelf/test/cov_square2.m @@ -77,7 +77,8 @@ function cov_square2(n,occ,p,rank_or_tol,symm,noise,scale,spdiag) tic rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))), ... + @(x)(x - rskelf_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) if strcmpi(symm,'p') diff --git a/rskelf/test/fd_cube.m b/rskelf/test/fd_cube.m index d11f90e..9d0c311 100644 --- a/rskelf/test/fd_cube.m +++ b/rskelf/test/fd_cube.m @@ -80,7 +80,8 @@ function fd_cube(n,occ,rank_or_tol,skip,symm) tic Y = rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*rskelf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*rskelf_sv(F,x)), ... + @(x)(x - rskelf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run CG @@ -88,7 +89,7 @@ function fd_cube(n,occ,rank_or_tol,skip,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(rskelf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)rskelf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/rskelf/test/fd_square.m b/rskelf/test/fd_square.m index c1a92ec..11cbf82 100644 --- a/rskelf/test/fd_square.m +++ b/rskelf/test/fd_square.m @@ -71,7 +71,8 @@ function fd_square(n,occ,rank_or_tol,symm) tic Y = rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - A*rskelf_sv(F,x)),[],[],1); + [e,niter] = snorm(N,@(x)(x - A*rskelf_sv(F,x)), ... + @(x)(x - rskelf_sv(F,A*x,'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run CG @@ -79,7 +80,7 @@ function fd_square(n,occ,rank_or_tol,symm) % run PCG tic - [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)(rskelf_sv(F,x))); + [Z,~,~,piter] = pcg(@(x)(A*x),X,1e-12,32,@(x)rskelf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - A*Z)/norm(X); diff --git a/rskelf/test/ie_circle.m b/rskelf/test/ie_circle.m index 525841a..efc6740 100644 --- a/rskelf/test/ie_circle.m +++ b/rskelf/test/ie_circle.m @@ -55,7 +55,8 @@ function ie_circle(n,occ,p,rank_or_tol,symm) tic rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))), ... + @(x)(x - rskelf_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % generate field due to exterior sources diff --git a/rskelf/test/ie_cube1.m b/rskelf/test/ie_cube1.m index d6e3b6c..19cb2c4 100644 --- a/rskelf/test/ie_cube1.m +++ b/rskelf/test/ie_cube1.m @@ -74,7 +74,8 @@ function ie_cube1(n,occ,p,rank_or_tol,symm) tic Y = rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))), ... + @(x)(x - rskelf_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -82,7 +83,7 @@ function ie_cube1(n,occ,p,rank_or_tol,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)(rskelf_sv(F,x))); + [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)rskelf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - mv(Z))/norm(X); diff --git a/rskelf/test/ie_cube2.m b/rskelf/test/ie_cube2.m index 1e1d027..f79f4c7 100644 --- a/rskelf/test/ie_cube2.m +++ b/rskelf/test/ie_cube2.m @@ -74,7 +74,8 @@ function ie_cube2(n,occ,p,rank_or_tol,symm) tic rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))), ... + @(x)(x - rskelf_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) end diff --git a/rskelf/test/ie_square1.m b/rskelf/test/ie_square1.m index cc1f115..ef43ce9 100644 --- a/rskelf/test/ie_square1.m +++ b/rskelf/test/ie_square1.m @@ -68,7 +68,8 @@ function ie_square1(n,occ,p,rank_or_tol,symm) tic Y = rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))), ... + @(x)(x - rskelf_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) % run unpreconditioned GMRES @@ -76,7 +77,7 @@ function ie_square1(n,occ,p,rank_or_tol,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)(rskelf_sv(F,x))); + [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)rskelf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - mv(Z))/norm(X); diff --git a/rskelf/test/ie_square2.m b/rskelf/test/ie_square2.m index d941246..c6caf9a 100644 --- a/rskelf/test/ie_square2.m +++ b/rskelf/test/ie_square2.m @@ -68,7 +68,8 @@ function ie_square2(n,occ,p,rank_or_tol,symm) tic rskelf_sv(F,X); t = toc; - [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))),[],[],1); + [e,niter] = snorm(N,@(x)(x - mv(rskelf_sv(F,x))), ... + @(x)(x - rskelf_sv(F,mv(x),'c'))); fprintf('sv: %10.4e / %4d / %10.4e (s)\n',e,niter,t) end diff --git a/rskelf/test/ie_square3.m b/rskelf/test/ie_square3.m index da5bf66..42a1595 100644 --- a/rskelf/test/ie_square3.m +++ b/rskelf/test/ie_square3.m @@ -93,7 +93,7 @@ function ie_square3(n,k,occ,p,rank_or_tol,symm) % run preconditioned GMRES tic - [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)(rskelf_sv(F,x))); + [Z,~,~,piter] = gmres(mv,X,[],1e-12,32,@(x)rskelf_sv(F,x)); t = toc; e1 = norm(Z - Y)/norm(Z); e2 = norm(X - mv(Z))/norm(X);