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