# source:projects/ckbs/trunk/test/tridiag_solve_ok.m@95

Last change on this file since 95 was 35, checked in by bradbell, 11 years ago

File size: 1.8 KB
Line
1% -------------------------------------------------------------------
2% ckbs: Constrained Kalman-Bucy Smoother Program: Copyright (C) 2006
4%          Gianluigi Pillonetto: giapi at dei dot unipd dot it
6% -------------------------------------------------------------------
7% \$begin tridiag_solve_ok.m\$\$ \$newlinech %\$\$
8% \$spell
9%       ckbs
10%       ak
11%       blk
12%       lhs
13%       qk
14%       tridiag
15%       uk
16% \$\$
17%
18% \$section ckbs_tridiag_solve Example and Test\$\$
19%
20% \$index ckbs_tridiag_solve, example and test\$\$
21% \$index tridiag_solve, example and test\$\$
22% \$index example, tridiag_solve\$\$
23% \$index test, tridiag_solve\$\$
24%
26% \$newlinech \$\$ \$codep
27function [ok] = tridiag_solve_ok()
28ok = true;
29m  = 2;
30n  = 3;
31N  = 4;
32% case where uk = 0, qk = I, and ak is random
33rand('seed', 123);
34a = rand(n, n, N);
35r = rand(n * N, m);
36c = zeros(n, n, N);
37for k = 1 : N
38        ak         = a(:, :, k);
39        b(:, :, k) = 2 * eye(n) + ak * ak';
40        c(:, :, k) = ak';
41end
42% ----------------------------------------
43[e, lambda] = ckbs_tridiag_solve(b, c, r);
44% ----------------------------------------
45%
46% check equation corresponding to first row
47blk_n = 1 : n;
48lhs   = b(:, :, 1) * e(blk_n, :) + c(:, :, 2)' * e(blk_n + n, :);
49ok    = ok & ( max( max( abs( lhs - r(blk_n, :) ) ) ) < 1e-10 );
50%
51% check equation corresponding to last row
52blk_n = (N-1) * n + (1 : n);
53lhs   = c(:, :, N) * e(blk_n - n, :) + b(:, :, N) * e(blk_n, :);
54ok    = ok & ( max( max( abs( lhs - r(blk_n, :) ) ) ) < 1e-10 );
55%
56% check other rows
57blk_n = 1 : n;
58for k = 2 : N-1
59        blk_n = blk_n + n;
60        lhs   = c(:, :, k)    * e(blk_n - n, :)  ...
61              + b(:, :, k)    * e(blk_n, :)    ...
62              + c(:, :, k+1)' * e(blk_n + n, :);
63        ok    = ok & ( max( max( abs( lhs - r(blk_n, :) ) ) ) < 1e-10 );
64end
65return
66end
67% \$\$ \$newlinech %\$\$
68% \$end
Note: See TracBrowser for help on using the repository browser.