[project @ 1996-07-25 21:02:03 by partain]
[nofib.git] / real / fluid / Chl_routs.hs
1 {-
2 The first part of Choleski decomposition.
3 Contains a matrix reodering function.
4 The generalized envelope method is implemented here.
5
6 XZ, 24/10/91
7 -}
8
9 {-
10 Modified to adopt S_arrays.
11
12 More efficient algorithms have been adopted.
13 They include:
14 a) minimum degree ordering (in module Min_degree.hs);
15 b) K matrix assembly.
16
17 Also, the output format has been changed.
18
19 XZ, 19/2/92
20 -}
21
22 module Chl_routs ( orded_mat ) where
23
24 import Defs
25 import S_Array -- not needed w/ proper module handling
26 import Norm -- ditto
27 import Min_degree
28 import Ix--1.3
29 infix 1 =:
30 (=:) a b = (a,b)
31 -----------------------------------------------------------
32 -- Liu's generalized envelope method adopted here. --
33 -- Reordering the system matric by firstly applying --
34 -- minimum degree ordering ( to minimize fill-ins ) and --
35 -- secondly applying postordering ( to optimize matrix --
36 -- structure ). The system matrix structure is found --
37 -- using the elimination tree. Used at the data setup --
38 -- stage. --
39 -----------------------------------------------------------
40
41 orded_mat
42 :: Int
43 -> (My_Array Int (Frac_type,((Frac_type,Frac_type,Frac_type),
44 (Frac_type,Frac_type,Frac_type))))
45 -> (My_Array Int [Int])
46 -> [Int]
47 -> (My_Array Int (My_Array Int Frac_type,My_Array Int (Int,[Frac_type])),My_Array Int Int)
48
49 orded_mat p_total el_det_fac p_steer fixed =
50 (init_L,o_to_n)
51 where
52 bindTo x k = k x -- old Haskell 1.0 "let", essentially
53 remove x = filter ((/=) x) -- also old Haskell 1.0 thing
54
55 n_bnds = (1,p_total)
56 n_bnds' = (0,p_total)
57 -- the inverse of an 1-D Int array.
58 inv_map = \a ->
59 s_array n_bnds' (map (\(i,j)->j=:i) (s_assocs a))
60 -- find the column indecies of nonzero entries in a row
61 get_js old_i map_f =
62 filter (\j->j<=i) (map ((!^) map_f) (old_rows!^old_i))
63 where i = map_f!^old_i
64 -- children of individual elimination tree nodes
65 chldrn = \e_tree ->
66 s_accumArray (++) [] n_bnds'
67 (map (\(i,j)->j=:[i]) (s_assocs e_tree))
68 -- the entry map from the input matrix to the output matrix
69 -- ( combination of o_to_min and min_to_n )
70 o_to_n :: (My_Array Int Int)
71 o_to_n = s_amap ((!^) min_to_n) o_to_min
72 n_to_o = inv_map o_to_n
73 -- the entry map of the minimum degree ordering
74 o_to_min = inv_map min_to_o
75 min_to_o = s_listArray n_bnds' (0:min_degree old_rows)
76 -- the entry map of postordering
77 -- switch off ordering
78 min_to_n :: My_Array Int Int
79 -- min_to_n = s_listArray n_bnds' (range n_bnds')
80 -- min_to_o = min_to_n
81 min_to_n =
82 s_array n_bnds' ((0=:0):(fst (recur ([],1) (chn!^0))))
83 where
84 chn = chldrn min_e_tree
85 -- recursive postordering
86 recur =
87 foldl
88 (
89 -- pattern before entering a loop
90 \ res r ->
91 -- current result of post-reordering
92 (recur res (chn!^r)) `bindTo` ( \ (new_reord,label) ->
93 ((r=:label):new_reord,label+1) )
94 )
95 -- the elimination tree of the reordered matrix
96 new_e_tree =
97 s_array n_bnds
98 ( map (\(i,j)-> (min_to_n!^i =: min_to_n!^j))
99 ( s_assocs min_e_tree ))
100 -- elimination tree of the matrix after minimum degree
101 -- ordering
102 min_e_tree =
103 s_def_array n_bnds (0::Int)
104 (all_rs (1::Int) init_arr [])
105 where
106 init_arr = s_def_array n_bnds (0::Int) []
107 -- implementation of an elimination tree construction
108 -- algorithm
109 all_rs i ance pare =
110 if ( i>p_total )
111 then pare
112 else all_rs (i+1) new_ance pare++rss
113 where
114 root old@(k,old_anc) =
115 if ( (new_k==0) || (new_k==i) )
116 then old
117 else root (new_k,old_anc//^[k=:i])
118 where new_k = old_anc!^k
119 -- finding new parents and ancestors
120 (rss,new_ance) =
121 -- looping over connetions of current node in
122 -- the matrix graph
123 foldl
124 (
125 -- pattern before entering a loop
126 \ (rs,anc) k1 ->
127 -- appending a new parent
128 (root (k1,anc)) `bindTo` ( \ (r,new_anc) ->
129 (r=:i) `bindTo` ( \ new_r ->
130 if new_anc!^r /= 0
131 then (rs, new_anc)
132 else (new_r:rs, new_anc //^ [new_r]) ))
133 )
134 ([],ance) (remove i (get_js (min_to_o!^i) o_to_min))
135 -- initial L
136 init_L =
137 s_listArray (1,length block_ends)
138 [
139 (
140 s_listArray bn [get_v i i|i<-range bn],
141 (filter (\ (_,j)->j<=u)
142 [ (i, find_first bn (find_non0 i))
143 | i <- range (l+1,p_total)
144 ]) `bindTo` ( \ non_emp_set ->
145
146 s_def_array (l+1,p_total) (u+1,[])
147 [ i=:(j',[get_v i j | j<- range (j',min u (i-1))])
148 | (i,j') <- non_emp_set
149 ] )
150 )
151 | bn@(l,u) <- block_bnds
152 ]
153 where
154 get_v i j =
155 if ( i'<j' )
156 then (old_mat!^j')!^i'
157 else (old_mat!^i')!^j'
158 where
159 i' = n_to_o!^i
160 j' = n_to_o!^j
161 find_non0 i =
162 foldl ( \ar j -> all_non0s j ar )
163 (s_def_array (1,i) False [])
164 (get_js (n_to_o!^i) o_to_n)
165 where
166 all_non0s j arr =
167 if ( j>i || j==0 || arr!^j )
168 then arr
169 else all_non0s (new_e_tree!^j) (arr//^[j=:True])
170 -- finding the first non-zero entry between l and u of the ith line
171 find_first :: (Int,Int) -> (My_Array Int Bool) -> Int
172 find_first (j1,u) non0_line = f' j1
173 where
174 f' j =
175 if (j>u) || non0_line!^j
176 then j
177 else f' (j+1)
178 -- reordered matrix in a new sparse form
179 block_ends =
180 [ i | (i,j)<-s_assocs new_e_tree, j/=(i+1) ]
181 block_bnds = zip (1:(map ((+) 1) (init block_ends))) block_ends
182 -- descendants of nodes of elimination tree
183 decnd :: My_Array Int [Int]
184 decnd =
185 s_listArray n_bnds
186 [ chn_n ++ concat [ decnd!^i | i <- chn_n ]
187 | chn_n <- s_elems (chldrn new_e_tree)
188 ]
189 -- rows of the K matrix (before ordering)
190 old_rows =
191 s_accumArray (++) [] n_bnds
192 ( concat
193 [
194 [j|(j,_)<-sparse_assocs (old_mat!^i)] `bindTo` ( \ j_set ->
195 (i=:j_set):[j'=:[i]|j'<-j_set,i/=j'] )
196 | i <- range n_bnds
197 ]
198 )
199 -- Value and index pairs of the original matrix.
200 -- This is found by assembling system K.
201 -- Fixed entries are multiplied by a large number
202 old_mat :: My_Array Int (My_Array Int Frac_type)
203 old_mat =
204 arr //^
205 [ (arr!^i) `bindTo` ( \ ar ->
206 i =: ar //^ [i=:(ar!^i)*large_scalor] )
207 | i <- fixed
208 ]
209 where
210 arr =
211 s_listArray n_bnds
212 [
213 s_accumArray (+) (0::Frac_type) (1,i) (temp!^i)
214 | i<-range n_bnds
215 ]
216 temp :: My_Array Int [(Int,Frac_type)]
217 temp =
218 s_accumArray (++) [] n_bnds
219 ( concat
220 [
221 (el_det_fac!^e) `bindTo` ( \ d_f ->
222 (zip (range (1,p_nodel)) (p_steer!^e)) `bindTo` ( \ pairs ->
223 concat
224 [
225 (dd_mat!^ii) `bindTo` ( \ dd_m ->
226 [ i =: [j =: (dd_m!^jj) d_f]
227 | (jj,j) <- pairs, j<=i
228 ] )
229 | (ii,i) <- pairs
230 ] ))
231 | e <- s_indices el_det_fac
232 ]
233 )
234 -- element contribution matrix
235 dd_mat =
236 s_listArray (1,p_nodel) [
237 s_listArray (1,p_nodel) [f11,f12,f13],
238 s_listArray (1,p_nodel) [f12,f22,f23],
239 s_listArray (1,p_nodel) [f13,f23,f33]
240 ]
241 where
242 f = \x y u v d -> (x*y+u*v)*d
243 s1 = \(x,_,_) -> x
244 s2 = \(_,y,_) -> y
245 s3 = \(_,_,z) -> z
246 f11 (det,(x,y)) = f c1 c1 c2 c2 det
247 where
248 c1 = s1 x
249 c2 = s1 y
250 f12 = \(det,(x,y)) -> f (s1 x) (s2 x) (s1 y) (s2 y) det
251 f13 = \(det,(x,y)) -> f (s1 x) (s3 x) (s1 y) (s3 y) det
252 f22 (det,(x,y)) = f c1 c1 c2 c2 det
253 where
254 c1 = s2 x
255 c2 = s2 y
256 f23 = \(det,(x,y)) -> f (s2 x) (s3 x) (s2 y) (s3 y) det
257 f33 (det,(x,y)) = f c1 c1 c2 c2 det
258 where
259 c1 = s3 x
260 c2 = s3 y