1% Solve the variational problem. 2% 3% fi = divand_solve(s,yo) 4% 5% Derive the analysis based on all contraints included in s and using the 6% observations yo 7% 8% Input: 9% s: structure created by divand_factorize 10% yo: value of the observations 11% 12% Output: 13% fi: analyzed field 14 15function [fi,s] = divand_solve(s,yo) 16 17H = s.H; 18sv = s.sv; 19R = s.R; 20% fix me: ignore argument 21yo = s.yo; 22 23 24if s.primal 25 if strcmp(s.inversion,'chol') 26 P = s.P; 27 fpi = P * (H'* (R \ yo(:))); 28 else 29 %H = double(H); 30 HiRyo = H'* (R \ yo(:)); 31 32 fun = @(x) s.iB*x + H'*(R\ (H * x)); 33 34 if ~s.keepLanczosVectors 35 [fpi,s.flag,s.relres,s.iter] = pcg(fun,HiRyo,s.tol,s.maxit,s.M1,s.M2); 36 37 if s.flag ~= 0 38 warning('divand-noconvergence','pcg did not converge'); 39 end 40 41 else 42 %pc = @(x) s.M2 \ (s.M1 \ x); 43 pc = []; 44 x0 = zeros(size(s.iB,1),1); 45 [fpi,Q,T,diag] = conjugategradient(fun,HiRyo,'tol',s.tol,... 46 'maxit',s.maxit,... 47 'minit',s.minit,... 48 'x0',x0,'renorm',1,'pc',pc); 49 s.iter = diag.iter; 50 s.relres = diag.relres; 51 s.P = CovarLanczos(Q,T); 52 end 53 end 54else % dual 55 B = s.B; 56 % fun(x) computes (H B H' + R)*x 57 fun = @(x) H * (B * (H'*x)) + R*x; 58 59 [tmp,flag,relres,iter] = pcg(fun, yo,s.tol,s.maxit,s.funPC); 60 61 if (flag ~= 0) 62 error('divand:pcg', ['Preconditioned conjugate gradients method'... 63 ' did not converge %d %g %g'],flag,relres,iter); 64 end 65 66 fpi = B * (H'*tmp); 67end 68 69 70[fi] = statevector_unpack(sv,fpi); 71 72fi(~s.mask) = NaN; 73 74% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be> 75% 76% This program is free software; you can redistribute it and/or modify it under 77% the terms of the GNU General Public License as published by the Free Software 78% Foundation; either version 2 of the License, or (at your option) any later 79% version. 80% 81% This program is distributed in the hope that it will be useful, but WITHOUT 82% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 83% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 84% details. 85% 86% You should have received a copy of the GNU General Public License along with 87% this program; if not, see <http://www.gnu.org/licenses/>. 88