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