{-
The second module of Choleski decomposition.
Contains Choleski factorization function.
XZ, 24/10/91
-}
{-
Modified to employ S_array.
The decomposition is now coded in the column-Choleski
feshion. It was in row-Choleski feshion. It might be
benifitial if numerical evaluation is forced (access
locality).
XZ, 7/2/92
-}
module Chl_decomp ( chl_factor ) where
import Defs
import S_Array -- not needed w/ proper module handling
import Norm -- ditto
import Asb_routs
import Ix--1.3
infix 1 =:
(=:) a b = (a,b)
-----------------------------------------------------------
-- Choleski factorization: --
-- it adopts the matrix struct derived from Liu's so --
-- called generalized envelope method. Called at the --
-- data setup stage. --
-----------------------------------------------------------
cal_one_elem (j_off,v_off) (_,(j_s,v_s)) old_v =
case v_s of
[] -> old_v
_ ->
old_v -
list_inner_prod
(drop (max_j-j_s) v_s)
(drop (max_j-j_off) v_off)
where max_j = max j_s j_off
cal_one_seg segs ((_,v_old),off@(j_off,v_off)) =
zipWith (cal_one_elem off) segs v_old
match_pairs a@((k1,l1):as) b@(b1@(k2,l2):bs)
| k1<k2 = match_pairs as b
| k2<k1 = match_pairs a bs
| otherwise = (k1=:(l1,l2)):match_pairs as bs
match_pairs _ _ = []
chl_factor :: S_array (S_array Float,S_array (Int,[Float]))
-> S_array (S_array Float,S_array (Int,[Float]))
chl_factor init_L = foldl f init_L (range (s_bounds init_L))
where
f old_l j = old_l //^ [j=:step2 step1]
where
(l,u) = block_bnds
block_bnds = s_bounds (fst this_block)
this_block = (old_l!^j)
bindTo x k = k x -- old Haskell 1.0 "let"
step1 = foldl step1_f
this_block
-- previous_blocks...
[ k=: ((snd (old_l!^k)) `bindTo` ( \ b ->
filter (\ (i,_) -> (l<=i)&&(i<=u)) (s_assocs b) ))
| k <- range (1,j-1)
]
where
step1_f (old_diag,old_off_diag) (k,segs_jk) = (new_diag,new_off_diag)
where
subs_segs ((i,pair@((j1,_),_)):rest) =
(i=:(j1,cal_one_seg (drop (j1-l) segs_jk) pair)):
subs_segs rest
subs_segs _ = []
new_diag =
s_accum (-) old_diag
[ i =: list_inner_prod vs vs
| (i,(_,vs@(_:_))) <- segs_jk
]
new_off_diag =
old_off_diag //^
(
subs_segs
(
match_pairs
(sparse_assocs old_off_diag)
(sparse_assocs (snd (old_l!^k)))
)
)
step2 (old_diag,old_off_diag) =
(
s_listArray block_bnds (fst ass_res),
old_off_diag //^ (tail (snd ass_res))
)
where
ass_res =
gen_assocs (sparse_assocs old_off_diag)
([sqrt (old_diag!^l)],[l=:(l+1,[])])
gen_assocs (old_line@(i,(j,vs)):rest)
(t_diag,t_off_diag) =
if i <= u
then
gen_assocs rest
(t_diag++[new_diag],t_off_diag++[i=:(j,new_off_diag)])
else
gen_assocs rest
(t_diag,t_off_diag++[i=:(j,new_off_diag)])
where
new_diag =
sqrt
( (old_diag!^i) -
list_inner_prod new_off_diag new_off_diag)
new_off_diag =
do_recur vs (drop (j-l) t_diag,drop (j-l) t_off_diag) []
do_recur (o_v:o_v_r) ((d:d_r),(off:off_r)) res =
do_recur o_v_r (d_r,off_r)
( res ++ [ (cal_one_elem (j,res) off o_v)/d ] )
do_recur _ _ res = res
gen_assocs _ res = res
|