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