1 {- 2 Implementation of the Jacobi iteration 3 4 XZ, 24/10/91 5 -} 6 7 {- 8 Modified to adopt S_array 9 10 XZ, 19/2/92 11 -} 12 13 module Jcb_method ( jcb_method ) where 14 15 import Defs 16 import S_Array -- not needed w/ proper module handling 17 import Norm -- ditto 18 import Asb_routs 19 import Tol_cal 20 import Ix--1.3 21 22 ----------------------------------------------------------- 23 -- Jacobi iteration method: -- 24 -- solves: scalor*MX = B -- 25 -- with iteration counter (r) as -- 26 -- X(r+1) = X(r) + relax*(B/scalor - MX(r))/D -- 27 -- The system matrix M is never really assembled -- 28 ----------------------------------------------------------- 29 30 jcb_method 31 :: Int 32 -> My_Array Int (Frac_type,((Frac_type,Frac_type,Frac_type), 33 (Frac_type,Frac_type,Frac_type))) 34 -> My_Array Int [(Int,Int)] 35 -> My_Array Int [Int] 36 -> (My_Array Int Bool,(My_Array Int Bool, My_Array Int Bool)) 37 -> (My_Array Int Frac_type, My_Array Int Frac_type) 38 -> Frac_type 39 -> Int 40 -> Frac_type 41 -> Frac_type 42 -> (My_Array Int Frac_type, My_Array Int Frac_type) 43 44 sparse_elems = \y -> map (\(_,x)->x) (sparse_assocs y) 45 46 jcb_method f_step el_det_fac asb_table v_steer 47 (all_bry,(x_fixed,y_fixed)) (b1,b2) scalor 48 max_iter m_iter_toler relax = 49 sub_jacobi (init_u,init_u) max_iter 50 where 51 n_bnds = s_bounds asb_table 52 init_u = s_array n_bnds [] 53 -- the recursive function of the Jacobi iteration 54 sub_jacobi old_x n = -- X(r) iteration_counter 55 if 56 -- checking the iteration limit 57 ( n <= 1 ) || 58 -- checking the tolerance 59 ( 60 (n /= max_iter) && 61 ( tol_cal 62 ((s_elems (fst new_x)) ++ (s_elems (snd new_x))) 63 ((sparse_elems (fst diff)) ++ (sparse_elems (snd diff))) 64 True 65 ) < m_iter_toler 66 ) 67 {- 68 ( 69 (sum ( map abs ( s_elems (fst diff) ) )) + 70 (sum ( map abs ( s_elems (snd diff) ) )) 71 ) 72 < m_iter_toler 73 -} 74 then new_x -- X(r+1) 75 else sub_jacobi new_x (n-1) 76 where 77 -- new adaptive x 78 new_x = 79 if ( n == max_iter ) 80 -- first iteration 81 then diff 82 -- otherwise 83 else add_u old_x diff 84 diff = 85 if ( f_step == 1 ) 86 then 87 -- For step 1. Only fixes some boundry nodes. 88 ( 89 find_diff x_fixed (fst old_x) b1, 90 find_diff y_fixed (snd old_x) b2 91 ) 92 else 93 -- For step 3. Fixes all boundry nodes. 94 ( 95 find_diff all_bry (fst old_x) b1, 96 find_diff all_bry (snd old_x) b2 97 ) 98 bindTo x k = k x -- old haskell 1.0 "let" 99 100 -- function which does one adaptation, 101 -- ie, calculates: 102 -- relax*(B - MX(r))/D 103 find_diff fixed x b = -- list_of_unfixed_nodes X(r) B 104 s_def_listArray n_bnds (0::Frac_type) 105 [ 106 if fixed!^i 107 then 0 108 else 109 ((b!^i) / scalor) `bindTo` ( \ b' -> 110 (( 111 if ( n == max_iter ) 112 then b' 113 else 114 b' - 115 sum [ 116 (list_inner_prod 117 (get_val x (v_steer!^e)) 118 (map (mult (fst (el_det_fac!^e))) (m_mat!^id))) 119 | (e,id)<-asb_table!^i 120 ] 121 ) * relax) / (mat_diag!^i) ) 122 | i <- range n_bnds 123 ] 124 -- row sums of system M 125 mat_diag = 126 s_listArray n_bnds 127 [ 128 sum 129 [ (fst ( el_det_fac!^e)) * (m_r_sum!^id) 130 | (e,id)<-asb_table!^i 131 ] 132 | i <- range n_bnds 133 ] 134 135 ----------------------------------------------------------- 136 -- row sums of element matrix m_mat -- 137 ----------------------------------------------------------- 138 139 m_r_sum = 140 s_def_listArray (1,v_nodel) def_v 141 [ def_v,def_v,def_v,v',v',v' ] 142 where 143 def_v = (12 / 180)::Frac_type 144 v' = (68 / 180)::Frac_type 145 146 ----------------------------------------------------------- 147 -- element matrix -- 148 ----------------------------------------------------------- 149 150 m_mat = 151 s_listArray (1,v_nodel) 152 (map (map (mult inv_180)) 153 [ 154 [6,(-1),(-1),(-4),0,0], 155 [(-1),6,(-1),0,(-4),0], 156 [(-1),(-1),6,0,0,(-4)], 157 [(-4),0,0,32,16,16], 158 [0,(-4),0,16,32,16], 159 [0,0,(-4),16,16,32] 160 ]) 161 where 162 inv_180 = 1 / 180