1% Factorize some matrices to increase performance.
2%
3% s = divand_factorize(s)
4%
5% Create the inverse of the posteriori covariance matrix
6% and factorize
7%
8% Input:
9%   s: structure created by divand_background
10%
11% Output:
12%   s.P: factorized P
13
14
15function s = divand_factorize(s)
16
17R = s.R;
18iB = s.iB;
19H = s.H;
20
21if s.primal
22    if strcmp(s.inversion,'chol')
23        iP = iB + H'*(R\H);
24        P = CovarIS(iP);
25
26        % Cholesky factor of the inverse of a posteriori
27        % error covariance iP
28        if s.factorize
29            P = factorize(P);
30        end
31
32        s.P = P;
33    else
34      tic
35        [s.M1,s.M2] = s.compPC(iB,H,R);
36      s.pc_time = toc();
37    end
38else % dual
39    %C = H * (iB \ H') + R;
40    s.B = CovarIS(iB);
41    % pre-conditioning function for conjugate gradient
42    s.funPC = [];
43
44    if s.factorize
45        s.B = factorize(s.B);
46
47        % iM = inv(H * B * H' + sparse_diag(diag(R)))
48        % iM * (H * B * H' + R) is identify matrix if R is diagonal
49
50        M = H * (s.B * H') + sparse_diag(diag(R));
51        %M = H * (s.B * H') + sparse_diag(sum(R,1));
52        iM = CovarIS(M);
53        iM = factorize(iM);
54
55        % pre-conditioning function for conjugate gradient
56        s.funPC = @(x) iM*x;
57
58        if 0
59            C = H * (s.B * H') + sparse_diag(diag(R));
60            [RC, q, QC] = chol (C);
61            s.funPC = @(x) RC' \ (QC*x);
62
63        end
64    end
65end
66
67% LocalWords:  divand
68
69% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
70%
71% This program is free software; you can redistribute it and/or modify it under
72% the terms of the GNU General Public License as published by the Free Software
73% Foundation; either version 2 of the License, or (at your option) any later
74% version.
75%
76% This program is distributed in the hope that it will be useful, but WITHOUT
77% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
78% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
79% details.
80%
81% You should have received a copy of the GNU General Public License along with
82% this program; if not, see <http://www.gnu.org/licenses/>.
83