[project @ 1996-07-25 21:02:03 by partain]
[nofib.git] / real / fluid / Chl_method.hs
1 {-
2 The third part of the Choleski decomposition.
3 Contains forward, backward substitution and their
4 driver.
5
6 XZ, 24/10/91
7 -}
8
9 {-
10 Modified to employ S_array.
11 (Forward and backward substitution functions have been
12 recoded.)
13
14 XZ, 7/2/92
15 -}
16
17 module Chl_method ( chl_method ) where
18
19 import Defs
20 import S_Array -- not needed w/ proper module handling
21 import Norm -- ditto
22 import Asb_routs
23 import Ix--1.3
24 infix 1 =:
25 (=:) a b = (a,b)
26
27 -----------------------------------------------------------
28 -- Forward substitution for the Choleski method. Called --
29 -- in "chl_method". --
30 -----------------------------------------------------------
31
32 lower_part u off_diag =
33 dropWhile (\(i,_)->i<=u) (sparse_assocs off_diag)
34
35 -- forward substitution
36 fwd_sbs chl_fac b_old =
37 s_listArray (s_bounds b_old) (gen_x (1::Int) b_old [])
38 where
39 b_up = snd (s_bounds chl_fac)
40 -- generate solution
41 gen_x blck b x =
42 if blck > b_up
43 then x
44 else gen_x (blck+1) new_b (x++new_x)
45 where
46 diag = fst this_block
47 off_diag = snd this_block
48 this_block = chl_fac!^blck
49 block_bounds = s_bounds diag
50 (l,u) = block_bounds
51 new_x = gen_block_x (l+1) [(b!^l)/(diag!^l)]
52 -- update RHS
53 new_b =
54 s_accum (-) b
55 [
56 k =: list_inner_prod (drop (j-l) new_x) vs
57 | (k,(j,vs)) <- lower_part u off_diag
58 ]
59 -- generate solution for one block
60 gen_block_x i x_res =
61 if i>u
62 then x_res
63 else
64 gen_block_x (i+1)
65 ( x_res ++
66 [ ((b!^i)-list_inner_prod (drop (j-l) x_res) vs)/(diag!^i) ]
67 )
68 where (j,vs) = off_diag!^i
69
70 -----------------------------------------------------------
71 -- Backward substitution for the Choleski method. --
72 -- It works in a column by column manner. --
73 -- Called in "chl_methold". --
74 -----------------------------------------------------------
75
76 -- backward substitution
77 bwd_sbs chl_fac y =
78 gen_x b_up ((s_array (s_bounds y) [])::(My_Array Int Frac_type))
79 where
80 b_up = snd (s_bounds chl_fac)
81 -- generate solution
82 gen_x blck x_res =
83 if blck < (1::Int)
84 then x_res
85 else gen_x (blck-1) new_x
86 where
87 diag = fst this_block
88 off_diag = snd this_block
89 this_block = chl_fac!^blck
90 block_bounds = s_bounds diag
91 (l,u) = block_bounds
92 new_x = gen_block_x u new_b x_res
93 -- update RHS
94 new_b =
95 s_accum (-)
96 (s_listArray block_bounds [y!^i|i<-range block_bounds])
97 ( concat
98 [
99 zipWith (\l v->l=:v) (range (j,u))
100 (map ((*) (x_res!^k)) vs)
101 | (k,(j,vs)) <- lower_part u off_diag
102 ]
103 )
104 -- generate solution for one block
105 gen_block_x i b x_res1 =
106 if i<l
107 then x_res1
108 else
109 gen_block_x (i-1) new_b1 (x_res1//^[i=:new_x])
110 where
111 new_x = (b!^i) / (diag!^i)
112 (j,vs) = off_diag!^i
113 new_b1 =
114 s_accum (-) b
115 (zipWith (\l v->l=:v) (range (j,i-1)) (map ((*) new_x) vs))
116
117 -----------------------------------------------------------
118 -- The driving function for the Choleski method. --
119 -- Because the used generalized envelope mathod reorders --
120 -- the system matrix, the right-hand-side and result are --
121 -- also reordered to match the internal and external --
122 -- forms. --
123 -- Called in the TG iteration. --
124 -- Calls "fwd_sbs" and "bwd_sbs" --
125 -----------------------------------------------------------
126
127 chl_method (chl_fac,o_to_n) b scalor =
128 -- parameters: (Choleski_factor,ordering) right_hand_side
129 -- constant_in_front_of_the_system
130 s_amap ((*) scalor) (s_ixmap bnds ((!^) o_to_n) x)
131 where
132 x = bwd_sbs chl_fac (fwd_sbs chl_fac new_b)
133 new_b =
134 s_array bnds
135 (map (\(i,v)->(o_to_n!^i)=:v) (s_assocs b))
136 bnds = s_bounds b