Adapt `nofib` code to Foldable-generalised Prelude
[nofib.git] / real / fluid / Chl_decomp.hs
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 Data.Ix
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