1 -- Glasgow Haskell 0.403 : BLAS ( Basic Linear Algebra System ) 2 -- ********************************************************************** 3 -- * * 4 -- * FILE NAME : vblldecomp.hs DATE : 5-3-1991 * 5 -- * * 6 -- * CONTENTS : LL decomposition method for variable bandwidth matrix * 7 -- * * 8 -- ********************************************************************** 9 10 module VBlldecomp(vblldecomp, vbllsolution, vbllsolution') where 11 12 import Basics 13 14 import Vector 15 16 import VBmatrix 17 18 vblldecomp :: Vbm Float -> Vbm Float 19 20 -- L[i,i] = sqrt(A[i,i] - ***) 21 -- where 22 -- *** = SUM [ L[i,k]*L[i,k] | k <- [1..i-1] ] 23 -- L[i,j] = (A[i,j] - ***) / L[j,j] 24 -- where 25 -- *** = SUM [ L[i,k]*L[j,k] | k <- [1..j-1] ] 26 27 vbllsolution :: Vbm Float -> Vec Float -> Vec Float 28 -- Solve equation system, the variable bandwidth matrix is 29 -- undecomposed, ie it is the original matrix. 30 31 vbllsolution' :: Vbm Float -> Vec Float -> Vec Float 32 -- Solve equation system, the variable bandwidth matrix is decomposed. 33 34 vblldecomp mA = m 35 where 36 m = makevbmat (boundvbmat mA) (diagadrvbm mA) f 37 f (i,j) = 38 if (i == j) then 39 sqrt ( vbmatsub mA (i,i) - 40 sum ( map ( \ k -> vbmatsub m (i,k) * 41 vbmatsub m (i,k) ) 42 [(fstclvbmat mA i) .. (i-1)] ) ) 43 44 else 45 ( vbmatsub mA (i,j) - 46 sum ( map ( \ k -> (vbmatsub m (i,k) * 47 vbmatsub m (j,k) ) ) 48 [ (k0 i j) .. (j - 1)] ) 49 ) / ( vbmatsub m (j,j) ) 50 k0 i j = if ( (fstclvbmat mA i) >= (fstclvbmat mA j) ) then 51 (fstclvbmat mA i) 52 else 53 (fstclvbmat mA j) 54 55 vbllsolution a b = 56 fst (backwarding ( forwarding (b, vblldecomp a) ) ) 57 58 vbllsolution' a' b = 59 fst (backwarding ( forwarding (b, a') ) ) 60 61 forwarding ( b, mVB ) = 62 (b',mVB) 63 where 64 b' = makevec n f 65 n = boundvec b 66 f i = ( vecsub b i - 67 sum [(vbmatsub mVB (i,j)) * (vecsub b' j) | j<-[l i..i-1] ] 68 ) / (vbmatsub mVB (i,i)) 69 l i = fstclvbmat mVB i 70 71 backwarding ( b, mVB ) = 72 (b',mVB) 73 where 74 n = boundvec b 75 b' = makevec n f 76 f i = ( vecsub b i - 77 sum [( vbmatsub mVB (j,i) ) * ( vecsub b' j ) 78 | j <- ( validj i [i+1..n] ) ] 79 ) / (vbmatsub mVB (i,i)) 80 validj i (j:js) = if ( i >= fstclvbmat mVB j ) then 81 j : (validj i js) 82 else 83 validj i js 84 validj i [] = [] 85