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 |
---|

27 | function [ok] = tridiag_solve_ok() |
---|

28 | ok = true; |
---|

29 | m = 2; |
---|

30 | n = 3; |
---|

31 | N = 4; |
---|

32 | % case where uk = 0, qk = I, and ak is random |
---|

33 | rand('seed', 123); |
---|

34 | a = rand(n, n, N); |
---|

35 | r = rand(n * N, m); |
---|

36 | c = zeros(n, n, N); |
---|

37 | for k = 1 : N |
---|

38 | ak = a(:, :, k); |
---|

39 | b(:, :, k) = 2 * eye(n) + ak * ak'; |
---|

40 | c(:, :, k) = ak'; |
---|

41 | end |
---|

42 | % ---------------------------------------- |
---|

43 | [e, lambda] = ckbs_tridiag_solve(b, c, r); |
---|

44 | % ---------------------------------------- |
---|

45 | % |
---|

46 | % check equation corresponding to first row |
---|

47 | blk_n = 1 : n; |
---|

48 | lhs = b(:, :, 1) * e(blk_n, :) + c(:, :, 2)' * e(blk_n + n, :); |
---|

49 | ok = ok & ( max( max( abs( lhs - r(blk_n, :) ) ) ) < 1e-10 ); |
---|

50 | % |
---|

51 | % check equation corresponding to last row |
---|

52 | blk_n = (N-1) * n + (1 : n); |
---|

53 | lhs = c(:, :, N) * e(blk_n - n, :) + b(:, :, N) * e(blk_n, :); |
---|

54 | ok = ok & ( max( max( abs( lhs - r(blk_n, :) ) ) ) < 1e-10 ); |
---|

55 | % |
---|

56 | % check other rows |
---|

57 | blk_n = 1 : n; |
---|

58 | for 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 ); |
---|

64 | end |
---|

65 | return |
---|

66 | end |
---|

67 | % $$ $newlinech %$$ |
---|

68 | % $end |
---|