1 {- 2 The second module of Choleski decomposition. 3 Contains Choleski factorization function. 4 5 XZ, 24/10/91 6 -} 7 8 {- 9 Modified to employ S_array. 10 11 The decomposition is now coded in the column-Choleski 12 feshion. It was in row-Choleski feshion. It might be 13 benifitial if numerical evaluation is forced (access 14 locality). 15 16 XZ, 7/2/92 17 -} 18 19 module Chl_decomp ( chl_factor ) where 20 21 import Defs 22 import S_Array -- not needed w/ proper module handling 23 import Norm -- ditto 24 import Asb_routs 25 import Ix--1.3 26 infix 1 =: 27 (=:) a b = (a,b) 28 29 ----------------------------------------------------------- 30 -- Choleski factorization: -- 31 -- it adopts the matrix struct derived from Liu's so -- 32 -- called generalized envelope method. Called at the -- 33 -- data setup stage. -- 34 ----------------------------------------------------------- 35 36 cal_one_elem (j_off,v_off) (_,(j_s,v_s)) old_v = 37 case v_s of 38 [] -> old_v 39 _ -> 40 old_v - 41 list_inner_prod 42 (drop (max_j-j_s) v_s) 43 (drop (max_j-j_off) v_off) 44 where max_j = max j_s j_off 45 46 cal_one_seg segs ((_,v_old),off@(j_off,v_off)) = 47 zipWith (cal_one_elem off) segs v_old 48 49 match_pairs a@((k1,l1):as) b@(b1@(k2,l2):bs) 50 | k1<k2 = match_pairs as b 51 | k2<k1 = match_pairs a bs 52 | otherwise = (k1=:(l1,l2)):match_pairs as bs 53 match_pairs _ _ = [] 54 55 chl_factor :: S_array (S_array Float,S_array (Int,[Float])) 56 -> S_array (S_array Float,S_array (Int,[Float])) 57 58 chl_factor init_L = foldl f init_L (range (s_bounds init_L)) 59 where 60 f old_l j = old_l //^ [j=:step2 step1] 61 where 62 (l,u) = block_bnds 63 block_bnds = s_bounds (fst this_block) 64 this_block = (old_l!^j) 65 bindTo x k = k x -- old Haskell 1.0 "let" 66 step1 = foldl step1_f 67 this_block 68 -- previous_blocks... 69 [ k=: ((snd (old_l!^k)) `bindTo` ( \ b -> 70 filter (\ (i,_) -> (l<=i)&&(i<=u)) (s_assocs b) )) 71 | k <- range (1,j-1) 72 ] 73 where 74 step1_f (old_diag,old_off_diag) (k,segs_jk) = (new_diag,new_off_diag) 75 where 76 subs_segs ((i,pair@((j1,_),_)):rest) = 77 (i=:(j1,cal_one_seg (drop (j1-l) segs_jk) pair)): 78 subs_segs rest 79 subs_segs _ = [] 80 new_diag = 81 s_accum (-) old_diag 82 [ i =: list_inner_prod vs vs 83 | (i,(_,vs@(_:_))) <- segs_jk 84 ] 85 new_off_diag = 86 old_off_diag //^ 87 ( 88 subs_segs 89 ( 90 match_pairs 91 (sparse_assocs old_off_diag) 92 (sparse_assocs (snd (old_l!^k))) 93 ) 94 ) 95 96 step2 (old_diag,old_off_diag) = 97 ( 98 s_listArray block_bnds (fst ass_res), 99 old_off_diag //^ (tail (snd ass_res)) 100 ) 101 where 102 ass_res = 103 gen_assocs (sparse_assocs old_off_diag) 104 ([sqrt (old_diag!^l)],[l=:(l+1,[])]) 105 gen_assocs (old_line@(i,(j,vs)):rest) 106 (t_diag,t_off_diag) = 107 if i <= u 108 then 109 gen_assocs rest 110 (t_diag++[new_diag],t_off_diag++[i=:(j,new_off_diag)]) 111 else 112 gen_assocs rest 113 (t_diag,t_off_diag++[i=:(j,new_off_diag)]) 114 where 115 new_diag = 116 sqrt 117 ( (old_diag!^i) - 118 list_inner_prod new_off_diag new_off_diag) 119 new_off_diag = 120 do_recur vs (drop (j-l) t_diag,drop (j-l) t_off_diag) [] 121 do_recur (o_v:o_v_r) ((d:d_r),(off:off_r)) res = 122 do_recur o_v_r (d_r,off_r) 123 ( res ++ [ (cal_one_elem (j,res) off o_v)/d ] ) 124 do_recur _ _ res = res 125 gen_assocs _ res = res