1 {-
    2 The second module of Choleski decomposition.
    3 Contains Choleski factorization function.
    4    
    5 XZ, 24/10/91
    6 -}
    7 
    8 {-
    9         Modified to employ S_array.
   10 
   11         The decomposition is now coded in the column-Choleski
   12         feshion.  It was in row-Choleski feshion.  It might be
   13         benifitial if numerical evaluation is forced (access
   14         locality).
   15 
   16         XZ, 7/2/92
   17 -}
   18    
   19 module Chl_decomp ( chl_factor ) where
   20    
   21 import Defs
   22 import S_Array  -- not needed w/ proper module handling
   23 import Norm     -- ditto
   24 import Asb_routs
   25 import Ix--1.3
   26 infix 1 =:
   27 (=:) a b = (a,b)
   28    
   29 -----------------------------------------------------------
   30 -- Choleski factorization:                               --
   31 -- it adopts the matrix struct derived from Liu's so     --
   32 -- called generalized envelope method.  Called at the    --
   33 -- data setup stage.                                     --
   34 -----------------------------------------------------------
   35 
   36 cal_one_elem (j_off,v_off) (_,(j_s,v_s)) old_v =
   37         case v_s of
   38                 [] -> old_v
   39                 _ ->
   40                         old_v -
   41                         list_inner_prod
   42                         (drop (max_j-j_s) v_s)
   43                         (drop (max_j-j_off) v_off)
   44         where max_j = max j_s j_off
   45 
   46 cal_one_seg segs ((_,v_old),off@(j_off,v_off)) =
   47         zipWith (cal_one_elem off) segs v_old
   48 
   49 match_pairs a@((k1,l1):as) b@(b1@(k2,l2):bs)
   50         | k1<k2 = match_pairs as b
   51         | k2<k1 = match_pairs a bs
   52         | otherwise = (k1=:(l1,l2)):match_pairs as bs
   53 match_pairs _ _ = []
   54 
   55 chl_factor :: S_array (S_array Float,S_array (Int,[Float]))
   56         -> S_array (S_array Float,S_array (Int,[Float]))
   57 
   58 chl_factor init_L = foldl f init_L (range (s_bounds init_L))
   59   where
   60     f old_l j = old_l //^ [j=:step2 step1]
   61       where
   62         (l,u) = block_bnds
   63         block_bnds = s_bounds (fst this_block)
   64         this_block = (old_l!^j)
   65         bindTo x k = k x -- old Haskell 1.0 "let"
   66         step1 = foldl step1_f
   67                       this_block
   68                       -- previous_blocks...
   69                       [ k=: ((snd (old_l!^k)) `bindTo` ( \ b ->
   70                             filter (\ (i,_) -> (l<=i)&&(i<=u)) (s_assocs b) ))
   71                               | k <- range (1,j-1)
   72                       ]
   73           where
   74             step1_f (old_diag,old_off_diag) (k,segs_jk) = (new_diag,new_off_diag)
   75                         where
   76                         subs_segs ((i,pair@((j1,_),_)):rest) =
   77                                 (i=:(j1,cal_one_seg (drop (j1-l) segs_jk) pair)):
   78                                 subs_segs rest
   79                         subs_segs _ = []
   80                         new_diag =
   81                                 s_accum (-) old_diag 
   82                                 [ i =: list_inner_prod vs vs
   83                                         | (i,(_,vs@(_:_))) <- segs_jk
   84                                 ]
   85                         new_off_diag =
   86                                 old_off_diag //^
   87                                 (
   88                                         subs_segs
   89                                         (
   90                                                 match_pairs
   91                                                 (sparse_assocs old_off_diag)
   92                                                 (sparse_assocs (snd (old_l!^k)))
   93                                         )
   94                                 )
   95 
   96         step2 (old_diag,old_off_diag) =
   97                 (
   98                         s_listArray block_bnds (fst ass_res),
   99                         old_off_diag //^ (tail (snd ass_res))
  100                 )
  101                 where
  102                 ass_res =
  103                         gen_assocs (sparse_assocs old_off_diag)
  104                         ([sqrt (old_diag!^l)],[l=:(l+1,[])])
  105                 gen_assocs (old_line@(i,(j,vs)):rest)
  106                                                                 (t_diag,t_off_diag) =
  107                         if i <= u
  108                         then
  109                                 gen_assocs rest
  110                                 (t_diag++[new_diag],t_off_diag++[i=:(j,new_off_diag)])
  111                         else
  112                                 gen_assocs rest
  113                                 (t_diag,t_off_diag++[i=:(j,new_off_diag)])
  114                         where
  115                         new_diag =
  116                                 sqrt
  117                                 ( (old_diag!^i) -
  118                                         list_inner_prod new_off_diag new_off_diag)
  119                         new_off_diag =
  120                                 do_recur vs (drop (j-l) t_diag,drop (j-l) t_off_diag) []
  121                         do_recur (o_v:o_v_r) ((d:d_r),(off:off_r)) res =
  122                                 do_recur o_v_r (d_r,off_r)
  123                                 ( res ++ [ (cal_one_elem (j,res) off o_v)/d ] )
  124                         do_recur _ _ res = res
  125                 gen_assocs _ res = res