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