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