1module cholesky;
2
3%**********************************************************************%
4%                                                                      %
5% Computation of the Cholesky decomposition of dense positive definite %
6% matrices containing numeric entries.                                 %
7%                                                                      %
8% Author: Matt Rebbeck, May 1994.                                      %
9%                                                                      %
10% The algorithm was taken from "Linear Algebra" - J.H.Wilkinson        %
11%                                                  & C. Reinsch        %
12%                                                                      %
13%                                                                      %
14% NB: By using the same rounded number techniques as used in svd this  %
15%     could be made a lot faster.                                      %
16%                                                                      %
17%**********************************************************************%
18
19% Redistribution and use in source and binary forms, with or without
20% modification, are permitted provided that the following conditions are met:
21%
22%    * Redistributions of source code must retain the relevant copyright
23%      notice, this list of conditions and the following disclaimer.
24%    * Redistributions in binary form must reproduce the above copyright
25%      notice, this list of conditions and the following disclaimer in the
26%      documentation and/or other materials provided with the distribution.
27%
28% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
29% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
30% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
31% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR
32% CONTRIBUTORS
33% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
34% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
35% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
36% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
37% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
38% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
39% POSSIBILITY OF SUCH DAMAGE.
40%
41
42
43
44
45symbolic procedure cholesky(mat1);
46  %
47  % A must be a positive definite symmetric matrix.
48  %
49  % LU decomposition of matrix A. ie: A=LU, where U is the transpose
50  % of L. The determinant of A is also computed as a side effect but
51  % has been commented out as it is not necessary. The procedure will
52  % fail if A is unsymmetric. It will also fail if A, modified by
53  % rounding errors, is not positive definite.
54  %
55  % The reciprocals of the diagonal elements are stored in p and the
56  % matrix is then 'dragged' out and 'glued' back together in get_l.
57  %
58  %
59  begin
60    scalar x,p,in_mat,l,u,i_turned_rounded_on;  % d1,d2;
61    integer i,j,n;
62    if not !*rounded then << i_turned_rounded_on := t; on rounded; >>;
63    if not matrixp(mat1) then
64     rederr "Error in cholesky:  non matrix input.";
65    if not symmetricp(mat1) then
66     rederr "Error in cholesky: input matrix is not symmetric.";
67    in_mat := copy_mat(mat1);
68    n := first size_of_matrix(in_mat);
69    p := mkvect(n);
70%    d1 := 1;
71%    d2 := 0;
72    for i:=1:n do
73    <<
74      for j:=i:n do
75      <<
76        x := innerprod(1,1,i-1,{'minus,getmat(in_mat,i,j)},
77                       row_vec(in_mat,i,n),row_vec(in_mat,j,n));
78        x := reval{'minus,x};
79        if j=i then
80        <<
81%          d1 := reval{'times,d1,x};
82          if get_num_part(my_reval(x))<=0 then rederr
83           "Error in cholesky: input matrix is not positive definite.";
84%          while abs(get_num_part(d1)) >= 1 do
85%          <<
86%            d1 := reval{'times,d1,0.0625};
87%            d2 := d2+4;
88%          >>;
89%          while abs(get_num_part(d1)) < 0.0625 do
90%          <<
91%            d1 := reval{'times,d1,16};
92%            d2 := d2-4;
93%          >>;
94          putv(p,i,reval{'quotient,1,{'sqrt,x}});
95        >>
96        else
97        <<
98          setmat(in_mat,j,i,reval{'times,x,getv(p,i)});
99        >>;
100      >>;
101    >>;
102    l := get_l(in_mat,p,n);
103    u := algebraic tp(l);
104    if i_turned_rounded_on then off rounded;
105    return {'list,l,u};
106  end;
107
108flag('(cholesky),'opfn);  % So it can be used from algebraic mode.
109
110
111
112symbolic procedure get_l(in_mat,p,sq_size);
113  %
114  % Pulls out L from in_mat and p.
115  %
116  begin
117    scalar l;
118    integer i,j;
119    l := mkmatrix(sq_size,sq_size);
120    for i:=1:sq_size do
121    <<
122      setmat(l,i,i,{'quotient,1,getv(p,i)});
123      for j:=1:i-1 do
124      <<
125        setmat(l,i,j,getmat(in_mat,i,j));
126      >>;
127    >>;
128    return l;
129  end;
130
131
132
133symbolic procedure symmetricp(in_mat);
134  %
135  % Checks input is symmetric. ie: transpose(A) = A. (boolean).
136  %
137  if algebraic (tp(in_mat)) neq in_mat then nil else t;
138
139flag('(symmetricp),'boolean);
140flag('(symmetricp),'opfn);
141
142
143endmodule;  % cholesky.
144
145end;
146