From f796508520b652d9af412817fe6ffb3aa1c8379a Mon Sep 17 00:00:00 2001 From: Travis Askham Date: Thu, 17 Jun 2021 22:24:09 -0400 Subject: [PATCH] add chunkerflam utility, passes some limited tests --- chunkie/+chnk/+flam/nproxy_square.m | 15 ++- chunkie/FLAM | 2 +- chunkie/chunkerflam.m | 196 ++++++++++++++++++++++++++++ chunkie/chunkermat.m | 5 - devtools/test/flamutilitiesTest.m | 15 ++- 5 files changed, 220 insertions(+), 13 deletions(-) create mode 100644 chunkie/chunkerflam.m diff --git a/chunkie/+chnk/+flam/nproxy_square.m b/chunkie/+chnk/+flam/nproxy_square.m index ea5ca4a..64c4321 100644 --- a/chunkie/+chnk/+flam/nproxy_square.m +++ b/chunkie/+chnk/+flam/nproxy_square.m @@ -9,23 +9,26 @@ % width - box/cube width (of center box, proxy surface at % 1.5*width) % opts - options structure -% opts.eps - tolerance (default 1e-13) +% opts.rank_or_tol - tolerance or maximum rank (default 1e-13) +% opts.eps - alias for rank_or_tol for backward compat (default 1e-13) % opts.nsrc - number of sources to use in test (default 200) % % Output: -% npxy - number of points on side of box to be passed -% to proxy surface routine +% npxy - number of points on perimeter of box to be sent to proxy routine % nsrc = 200; - eps = 1e-13; + rank_or_tol = 1e-13; if nargin < 3 opts = []; end if isfield(opts,'eps') - eps = opts.eps; + rank_or_tol = opts.eps; + end + if isfield(opts,'rank_or_tol') + rank_or_tol = opts.rank_or_tol; end if isfield(opts,'nsrc') nsrc = opts.nsrc; @@ -50,7 +53,7 @@ mat = kern(srcinfo,targinfo); - [sk,~] = id(mat,eps); + [sk,~] = id(mat,rank_or_tol); if length(sk) < min(nsrc,npxy) npxy = (floor((length(sk)-1)/4+0.1)+1)*4; diff --git a/chunkie/FLAM b/chunkie/FLAM index f278b10..8a67c78 160000 --- a/chunkie/FLAM +++ b/chunkie/FLAM @@ -1 +1 @@ -Subproject commit f278b10a385fd9430ce20bf0810c3cfe7065f284 +Subproject commit 8a67c788b053e409f1655627d15f1e5f14f9f489 diff --git a/chunkie/chunkerflam.m b/chunkie/chunkerflam.m new file mode 100644 index 0000000..c7d9bc8 --- /dev/null +++ b/chunkie/chunkerflam.m @@ -0,0 +1,196 @@ +function [F] = chunkerflam(chnkr,kern,dval,opts) +%CHUNKERFLAM build the requested FLAM compressed representation +% (e.g. a recursive skeletonization factorization) of the system matrix +% for given kernel and chunker description of boundary. This routine +% has the same quadrature options as CHUNKERMAT. +% +% Syntax: sysmat = chunkerflam(chnkr,kern,opts) +% +% Input: +% chnkr - chunker object describing boundary +% kern - kernel function. By default, this should be a function handle +% accepting input of the form kern(srcinfo,targinfo), where +% srcinfo and targinfo are in the ptinfo struct format, i.e. +% ptinfo.r - positions (2,:) array +% ptinfo.d - first derivative in underlying +% parameterization (2,:) +% ptinfo.d2 - second derivative in underlying +% parameterization (2,:) +% dval - (default 0.0) float or float array. Let A be the matrix +% corresponding to on-curve convolution with the provided kernel. +% If a scalar is provided, the system matrix is +% A + dval*eye(size(A)) +% If a vector is provided, it should be length size(A,1). The +% system matrix is then +% A + diag(dval) +% +% Optional input: +% opts - options structure. available options (default settings) +% opts.flamtype = string ('rskelf'), type of compressed +% representation to compute. Available: +% +% The recursive skeletonization routines +% should be sufficient curves which are not +% approximately space filling. +% +% - 'rskelf', recursive skeletonization +% factorization. Can be immediately used +% for system solves (via RSKELF_SV) and +% determinants (via RSKELF_LOGDET) +% - 'rskel', recursive skeletonization +% compression. Can be immediately used for +% matvecs or embedded in a sparse matrix. +% Not recommended unless you have a +% compelling reason. +% opts.useproxy = boolean (true), use a proxy function to speed +% up the FLAM compression. It may be desirable to +% turn this off if there appear to be precision +% issues. +% +% opts.quad = string ('ggqlog'), specify quadrature routine to +% use. +% +% - 'ggqlog' uses a generalized Gaussian quadrature +% designed for logarithmically singular kernels and +% smooth kernels with removable singularities +% - 'native' selects standard scaled Gauss-Legendre +% quadrature for native functions +% opts.l2scale = boolean (false), if true scale rows by +% sqrt(whts) and columns by 1/sqrt(whts) +% opts.occ = integer "occupancy" parameter (200) determines +% how many sources are in any leaf of the tree used to +% sort points. Determines the smaller dimensions of the +% first set of compressions. Sent to FLAM +% opts.rank_or_tol = integer or float "rank_or_tol" +% parameter (1e-14). Lower precision increases speed. +% Sent to FLAM +% +% Output: +% F - the requested FLAM compressed representation of the +% system matrix +% +% Examples: +% F = chunkermat(chnkr,kern,dval); % standard options +% sysmat = chunkermat(chnkr,kern,dval,opts); +% + +F = []; + +if length(chnkr) > 1 + chnkr = merge(chnkr); +end + +if nargin < 3 + dval = 0.0; +end + +if nargin < 4 + opts = []; +end + +quad = 'ggqlog'; +flamtype = 'rskelf'; +useproxy = true; +occ = 200; +rank_or_tol = 1e-14; + + +if or(chnkr.nch < 1,chnkr.k < 1) + warning('empty chunker, doing nothing') + sp = []; + return +end + +% determine operator dimensions using first two points + +srcinfo = []; targinfo = []; +srcinfo.r = chnkr.r(:,1); srcinfo.d = chnkr.d(:,1); +srcinfo.d2 = chnkr.d2(:,1); +i2 = min(2,chnkr.npt); +targinfo.r = chnkr.r(:,i2); targinfo.d = chnkr.d(:,i2); +targinfo.d2 = chnkr.d2(:,i2); + +ftemp = kern(srcinfo,targinfo); +opdims = size(ftemp); + +if (length(dval) == 1) + dval = dval*ones(opdims(1)*chnkr.npt,1); +end + +if (length(dval) ~= opdims(1)*chnkr.npt) + warning('provided dval array is length %d. must be scalar or length %d',... + length(dval),opdims(1)*chnkr.npt); + return +end + +% get opts from struct if available + +if isfield(opts,'quad') + quad = opts.quad; +end +if isfield(opts,'l2scale') + l2scale = opts.l2scale; +end + +if isfield(opts,'occ') + occ = opts.occ; +end +if isfield(opts,'rank_or_tol') + rank_or_tol = opts.rank_or_tol; +end + +% check if chosen FLAM routine implemented before doing real work + +if ~ (strcmpi(flamtype,'rskelf') || strcmpi(flamtype,'rskel')) + warning('selected flamtype %s not available, doing nothing',flamtype) + return +end + +% get nonsmooth quadrature + +if strcmpi(quad,'ggqlog') + + type = 'log'; + sp = chnk.quadggq.buildmattd(chnkr,kern,opdims,type); + +elseif strcmpi(quad,'native') + + sp = sparse(chnkr.npt,chnkr.npt); + +else + warning('specified quadrature method not available'); + return; +end + +sp = sp + spdiags(dval,0,chnkr.npt,chnkr.npt); + +% prep and call flam + +wts = weights(chnkr); +xflam = chnkr.r(:,:); + +width = max(max(chnkr)-min(chnkr)); +optsnpxy = []; optsnpxy.rank_or_tol = rank_or_tol; +optsnpxy.nsrc = occ; + +npxy = chnk.flam.nproxy_square(kern,width,optsnpxy); + +matfun = @(i,j) chnk.flam.kernbyindex(i,j,chnkr,wts,kern,opdims,sp); +[pr,ptau,pw,pin] = chnk.flam.proxy_square_pts(npxy); + +if strcmpi(flamtype,'rskelf') + ifaddtrans = true; + pxyfun = @(x,slf,nbr,l,ctr) chnk.flam.proxyfun(slf,nbr,l,ctr,chnkr,wts, ... + kern,opdims,pr,ptau,pw,pin,ifaddtrans); + F = rskelf(matfun,xflam,occ,rank_or_tol,pxyfun); +end + +if strcmpi(flamtype,'rskel') + pxyfunr = @(rc,rx,cx,slf,nbr,l,ctr) chnk.flam.proxyfunr(rc,rx,slf,nbr,l, ... + ctr,chnkr,wts,kern,opdims,pr,ptau,pw,pin); + F = rskel(matfun,xflam,xflam,occ,rank_or_tol,pxyfunr); +end + + + +end diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m index 7f06684..4b80bdc 100644 --- a/chunkie/chunkermat.m +++ b/chunkie/chunkermat.m @@ -28,7 +28,6 @@ % - 'native' selects standard scaled Gauss-Legendre % quadrature for native functions % -% opts.quadorder = integer (chnkr.k), desired quadrature order. % opts.nonsmoothonly = boolean (false), if true, only compute the % entries for which a special quadrature is used % (e.g. self and neighbor interactoins) and return @@ -53,7 +52,6 @@ opts = []; end -quadorder = chnkr.k; quad = 'ggqlog'; nonsmoothonly = false; l2scale = false; @@ -76,9 +74,6 @@ % get opts from struct if available -if isfield(opts,'quadorder') - quadorder = opts.quadorder; -end if isfield(opts,'quad') quad = opts.quad; end diff --git a/devtools/test/flamutilitiesTest.m b/devtools/test/flamutilitiesTest.m index 3b799d6..44439d1 100644 --- a/devtools/test/flamutilitiesTest.m +++ b/devtools/test/flamutilitiesTest.m @@ -103,11 +103,24 @@ ifaddtrans = true; pxyfun = @(x,slf,nbr,l,ctr) chnk.flam.proxyfun(slf,nbr,l,ctr,chnkr,wts, ... fkern,opdims,pr,ptau,pw,pin,ifaddtrans); + + start = tic; F = rskelf(matfun,xflam,200,1e-14,pxyfun); t1 = toc(start); +F = chunkerflam(chnkr,fkern,-0.5); + +fprintf('%5.2e s : time for flam rskelf compress\n',t1) + +pxyfunr = @(rc,rx,cx,slf,nbr,l,ctr) chnk.flam.proxyfunr(rc,rx,slf,nbr,l, ... + ctr,chnkr,wts,fkern,opdims,pr,ptau,pw,pin); + +opts = []; +start = tic; F2 = rskel(matfun,xflam,xflam,200,1e-14,pxyfunr,opts); t1 = toc(start); + +fprintf('%5.2e s : time for flam rskel compress\n',t1) afun = @(x) rskelf_mv(F,x); -fprintf('%5.2e s : time for flam rskelf compress\n',t1) + err2 = norm(sys2-sys,'fro')/norm(sys,'fro'); fprintf('%5.2e : fro error of build \n',err2)