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