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

Revmove reference to gpl2.txt and further shorten copyright and license statement

File size: 1.8 KB
Line 
1% -------------------------------------------------------------------
2% ckbs: Constrained Kalman-Bucy Smoother Program: Copyright (C) 2006
3% Authors: Bradlely Bell:        bradbell at washington dot edu
4%          Gianluigi Pillonetto: giapi at dei dot unipd dot it
5% License: GNU General Public License Version 2
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%
25% $head Source Code$$
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.