1% Testing divand in 2 dimensions with pcg solver. 2 3% grid of background field 4[xi,yi] = ndgrid(linspace(0,1,20)); 5 6mask = ones(size(xi)); 7pm = ones(size(xi)) / (xi(2,1)-xi(1,1)); 8pn = ones(size(xi)) / (yi(1,2)-yi(1,1)); 9 10% grid of observations 11[x,y] = ndgrid(linspace(eps,1-eps,10)); 12x = x(:); 13y = y(:); 14v = sin( x*6 ) .* cos( y*6); 15 16 17len = 0.1; 18lambda = 1; 19n = sum(mask(:)); 20 21% constrain 22fun = @(x) mean(x,1); 23funt = @(x) repmat(x,[n 1])/n; 24 25H = ones(1,n)/n; 26c.H = MatFun([1 n],fun,funt); 27%c.H = H; 28c.R = 0.001; 29c.yo = 1; 30 31% bug 32%[va,err,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'diagnostics',1,'primal',1,'factorize',0,'constraint',c,'inversion','pcg'); 33 34[va0,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,... 35 'diagnostics',1,'primal',1); 36 37[va,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,... 38 'diagnostics',1,'primal',1,'inversion','pcg',... 39 'tol',1e-6,'maxit',10000); 40 41 42rms_diff = []; 43rms_diff(end+1) = sqrt(mean((va(:) - va0(:)).^2)); 44 45 46if any(rms_diff > 1e-4) 47 error('unexpected large difference'); 48end 49 50fprintf('(max rms=%g) ',max(rms_diff)); 51 52% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be> 53% 54% This program is free software; you can redistribute it and/or modify it under 55% the terms of the GNU General Public License as published by the Free Software 56% Foundation; either version 2 of the License, or (at your option) any later 57% version. 58% 59% This program is distributed in the hope that it will be useful, but WITHOUT 60% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 61% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 62% details. 63% 64% You should have received a copy of the GNU General Public License along with 65% this program; if not, see <http://www.gnu.org/licenses/>. 66