blas 0.4 → 0.4.1
raw patch · 8 files changed
+655/−49 lines, 8 filessetup-changed
Files
- BLAS/C/Level1.hs +1/−0
- Data/Matrix/Banded/Internal.hs +419/−0
- INSTALL +100/−0
- Setup.lhs +8/−8
- blas.cabal +59/−3
- tests/Makefile +37/−0
- tests/Matrix.hs +31/−36
- tests/Vector.hs +0/−2
BLAS/C/Level1.hs view
@@ -1,4 +1,5 @@ {-# LANGUAGE FlexibleInstances #-}+{-# OPTIONS_GHC -fno-excess-precision #-} ----------------------------------------------------------------------------- -- | -- Module : BLAS.C.Level1
+ Data/Matrix/Banded/Internal.hs view
@@ -0,0 +1,419 @@+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : Data.Matrix.Banded.Internal+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Data.Matrix.Banded.Internal (+ -- * Banded matrix data types+ BMatrix(..),+ Banded,+ IOBanded,++ module BLAS.Matrix.Base,+ module BLAS.Tensor,++ -- * Converting to and from foreign pointers+ toForeignPtr,+ fromForeignPtr,+ ldaOf,+ isHerm,+ + -- * To and from the underlying storage matrix+ toMatrix,+ fromMatrix,+ + -- * Bandwith properties+ bandwidth,+ numLower,+ numUpper,+ + -- * Creating new matrices+ -- ** Pure+ banded,+ listsBanded,+ -- ** Impure+ newBanded_,+ newBanded,+ newListsBanded,+ + -- * Getting rows and columns+ row,+ col,+ getRow,+ getCol,+ + -- * Vector views+ diag,+ rowView,+ colView,+ + -- * Casting matrices+ coerceMatrix,+ + -- * Unsafe operations+ unsafeBanded,+ unsafeNewBanded,+ unsafeFreeze,+ unsafeThaw,+ unsafeWithElemPtr,+ unsafeDiag,+ unsafeGetRow,+ unsafeGetCol,+ unsafeRow,+ unsafeCol,+ unsafeRowView,+ unsafeColView,+ ) where++import Control.Arrow ( second )+import Control.Monad ( zipWithM_ )+import Data.Ix ( inRange, range )+import Data.List ( foldl' )+import Foreign+import System.IO.Unsafe+import Unsafe.Coerce ++import BLAS.Access+import BLAS.Elem ( Elem, BLAS1 )+import qualified BLAS.Elem as E+import BLAS.Internal ( checkedRow, checkedCol, checkedDiag, diagStart, + diagLen, clearArray, inlinePerformIO )+ +import BLAS.Matrix.Base hiding ( Matrix )+import qualified BLAS.Matrix.Base as C+import BLAS.Tensor++import Data.Matrix.Dense.Internal ( DMatrix )+import qualified Data.Matrix.Dense.Internal as M+ +import Data.Vector.Dense.Internal ( DVector, Vector, conj, dim, newListVector )+import qualified Data.Vector.Dense.Internal as V++ +data BMatrix t mn e + = BM { fptr :: !(ForeignPtr e)+ , offset :: !Int+ , size1 :: !Int+ , size2 :: !Int+ , lowBW :: !Int+ , upBW :: !Int+ , lda :: !Int+ }+ | H !(BMatrix t mn e) ++type Banded = BMatrix Imm+type IOBanded = BMatrix Mut++fromForeignPtr :: ForeignPtr e -> Int -> (Int,Int) -> (Int,Int) -> Int + -> BMatrix t (m,n) e+fromForeignPtr f o (m,n) (kl,ku) l = BM f o m n kl ku l++toForeignPtr :: BMatrix t (m,n) e -> (ForeignPtr e, Int, (Int,Int), (Int,Int), Int)+toForeignPtr (H a) = toForeignPtr a+toForeignPtr (BM f o m n kl ku l) = (f, o, (m,n), (kl,ku), l)++ldaOf :: BMatrix t (m,n) e -> Int+ldaOf (H a) = ldaOf a+ldaOf a = lda a++isHerm :: BMatrix t (m,n) e -> Bool+isHerm (H a) = not (isHerm a)+isHerm _ = False++unsafeFreeze :: BMatrix t mn e -> Banded mn e+unsafeFreeze = unsafeCoerce++unsafeThaw :: BMatrix t mn e -> IOBanded mn e+unsafeThaw = unsafeCoerce++-- | Coerce the phantom shape type from one type to another.+coerceMatrix :: BMatrix t mn e -> BMatrix t kl e+coerceMatrix = unsafeCoerce++toMatrix :: (Elem e) => BMatrix t (m,n) e -> (DMatrix t (m',n') e, (Int,Int), (Int,Int))+toMatrix (H a) = + case toMatrix (herm a) of (b, (m,n), (kl,ku)) -> (herm b, (n,m), (ku,kl))+toMatrix (BM f o m n kl ku ld) = + (M.fromForeignPtr f o (kl+1+ku,n) ld, (m,n), (kl,ku))++fromMatrix :: (Elem e) => DMatrix t (m,n) e -> (Int,Int) -> (Int,Int) -> BMatrix t (m',n') e+fromMatrix a (m,n) (kl,ku) = case a of+ (M.H a') -> herm (fromMatrix a' (n,m) (ku,kl))+ _ -> + let (f,o,(m',n'),ld) = M.toForeignPtr a+ in case undefined of+ _ | m' /= kl+1+ku -> + error $ "fromMatrix: number of rows must be equal to number of diagonals"+ _ | n' /= n ->+ error $ "fromMatrix: numbers of columns must be equal"+ _ ->+ BM f o m n kl ku ld+ ++bandwidth :: BMatrix t (m,n) e -> (Int,Int)+bandwidth a = + let (kl,ku) = (numLower a, numUpper a)+ in (negate kl, ku)++numLower :: BMatrix t (m,n) e -> Int+numLower (H a) = numUpper a+numLower a = lowBW a++numUpper :: BMatrix t (m,n) e -> Int+numUpper (H a) = numLower a+numUpper a = upBW a+++newBanded_ :: (Elem e) => (Int,Int) -> (Int,Int) -> IO (BMatrix t (m,n) e)+newBanded_ (m,n) (kl,ku)+ | m < 0 || n < 0 =+ err "dimensions must be non-negative."+ | kl < 0 =+ err "lower bandwdth must be non-negative."+ | m /= 0 && kl >= m =+ err "lower bandwidth must be less than m."+ | ku < 0 =+ err "upper bandwidth must be non-negative."+ | n /= 0 && ku >= n =+ err "upper bandwidth must be less than n."+ | otherwise =+ let off = 0+ m' = kl + 1 + ku+ l = m'+ in do + ptr <- mallocForeignPtrArray (m' * n)+ return $ fromForeignPtr ptr off (m,n) (kl,ku) l+ where+ err s = ioError $ userError $ + "newBanded_ " ++ show (m,n) ++ " " ++ show (kl,ku) ++ ": " ++ s++banded :: (BLAS1 e) => (Int,Int) -> (Int,Int) -> [((Int,Int), e)] -> Banded (m,n) e+banded mn kl ijes = unsafePerformIO $ newBanded mn kl ijes+{-# NOINLINE banded #-}++unsafeBanded :: (BLAS1 e) => (Int,Int) -> (Int,Int) -> [((Int,Int), e)] -> Banded (m,n) e+unsafeBanded mn kl ijes = unsafePerformIO $ unsafeNewBanded mn kl ijes+{-# NOINLINE unsafeBanded #-}++newBanded :: (BLAS1 e) => (Int,Int) -> (Int,Int) -> [((Int,Int), e)] -> IO (BMatrix t (m,n) e)+newBanded = newBandedHelp writeElem++unsafeNewBanded :: (BLAS1 e) => (Int,Int) -> (Int,Int) -> [((Int,Int), e)] -> IO (BMatrix t (m,n) e)+unsafeNewBanded = newBandedHelp unsafeWriteElem++newBandedHelp :: (BLAS1 e) => + (IOBanded (m,n) e -> (Int,Int) -> e -> IO ()) + -> (Int,Int) -> (Int,Int) -> [((Int,Int),e)] -> IO (BMatrix t (m,n) e)+newBandedHelp set (m,n) (kl,ku) ijes = do+ x <- newBanded_ (m,n) (kl,ku)+ withForeignPtr (fptr x) $ flip clearArray ((kl+1+ku)*n)+ mapM_ (uncurry $ set $ unsafeThaw x) ijes+ return x++listsBanded :: (BLAS1 e) => (Int,Int) -> (Int,Int) -> [[e]] -> Banded (m,n) e+listsBanded mn kl xs = unsafePerformIO $ newListsBanded mn kl xs+{-# NOINLINE listsBanded #-}++newListsBanded :: (BLAS1 e) => (Int,Int) -> (Int,Int) -> [[e]] -> IO (BMatrix t (m,n) e)+newListsBanded (m,n) (kl,ku) xs = do+ a <- newBanded_ (m,n) (kl,ku)+ zipWithM_ (writeDiagElems (unsafeThaw a)) [(negate kl)..ku] xs+ return a+ where+ writeDiagElems a i es =+ let d = diag a i+ nb = max 0 (negate i)+ es' = drop nb es+ in zipWithM_ (unsafeWriteElem d) [0..(dim d - 1)] es'++unsafeDiag :: (Elem e) => BMatrix t (m,n) e -> Int -> DVector t k e+unsafeDiag (H a) d = + conj $ unsafeDiag a (negate d)+unsafeDiag a d =+ let f = fptr a+ off = indexOf a (diagStart d)+ len = diagLen (shape a) d+ stride = lda a+ in V.fromForeignPtr f off len stride+ +diag :: (Elem e) => BMatrix t (m,n) e -> Int -> DVector t k e+diag a = checkedDiag (shape a) (unsafeDiag a) +++indexOf :: BMatrix t mn e -> (Int,Int) -> Int+indexOf (H a) (i,j) = indexOf a (j,i)+indexOf (BM _ off _ _ _ ku ld) (i,j) =+ off + ku + (i - j) + j * ld+ --off + i * tda + (j - i + kl)+ ++ +unsafeWithElemPtr :: (Elem e) => BMatrix t (m,n) e -> (Int,Int) -> (Ptr e -> IO a) -> IO a+unsafeWithElemPtr a (i,j) f = case a of+ (H a') -> unsafeWithElemPtr a' (j,i) f+ _ -> withForeignPtr (fptr a) $ \ptr ->+ f $ ptr `advancePtr` (indexOf a (i,j))++row :: (BLAS1 e) => Banded (m,n) e -> Int -> Vector n e+row a = checkedRow (shape a) (unsafeRow a)++unsafeRow :: (BLAS1 e) => Banded (m,n) e -> Int -> Vector n e+unsafeRow a i = unsafePerformIO $ getRow a i+{-# NOINLINE unsafeRow #-}++getRow :: (BLAS1 e) => BMatrix t (m,n) e -> Int -> IO (DVector r n e)+getRow a = checkedRow (shape a) (unsafeGetRow a)++unsafeGetRow :: (BLAS1 e) => BMatrix t (m,n) e -> Int -> IO (DVector r n e)+unsafeGetRow a i = + let (nb,x,na) = unsafeRowView a i+ n = numCols a+ in do+ es <- getElems x+ newListVector n $ (replicate nb 0) ++ es ++ (replicate na 0)++col :: (BLAS1 e) => Banded (m,n) e -> Int -> Vector m e+col a = checkedCol (shape a) (unsafeCol a)++unsafeCol :: (BLAS1 e) => Banded (m,n) e -> Int -> Vector m e+unsafeCol a i = unsafePerformIO $ getCol a i+{-# NOINLINE unsafeCol #-}++getCol :: (BLAS1 e) => BMatrix t (m,n) e -> Int -> IO (DVector r m e)+getCol a = checkedCol (shape a) (unsafeGetCol a)++unsafeGetCol :: (BLAS1 e) => BMatrix t (m,n) e -> Int -> IO (DVector r m e)+unsafeGetCol a j = unsafeGetRow (herm a) j >>= return . conj+++unsafeColView :: (Elem e) => BMatrix t (m,n) e -> Int -> (Int, DVector t k e, Int)+unsafeColView (BM f off m _ kl ku ld) j =+ let nb = max (j - ku) 0+ na = max (m - 1 - j - kl) 0+ r = max (ku - j) 0 + c = j + off' = off + r + c * ld+ stride = 1+ len = m - (nb + na)+ in if len >= 0+ then (nb, V.fromForeignPtr f off' len stride, na)+ else (m , V.fromForeignPtr f off' 0 stride, 0)+ +unsafeColView a j = + case unsafeRowView (herm a) j of (nb, v, na) -> (nb, conj v, na)++unsafeRowView :: (Elem e) => BMatrix t (m,n) e -> Int -> (Int, DVector t k e, Int)+unsafeRowView (BM f off _ n kl ku ld) i =+ let nb = max (i - kl) 0+ na = max (n - 1 - i - ku) 0+ r = min (ku + i) (kl + ku)+ c = max (i - kl) 0 + off' = off + r + c * ld+ stride = ld - 1+ len = n - (nb + na)+ in if len >= 0 + then (nb, V.fromForeignPtr f off' len stride, na)+ else (n , V.fromForeignPtr f off' 0 stride, 0)++unsafeRowView a i =+ case unsafeColView (herm a) i of (nb, v, na) -> (nb, conj v, na)+++rowView :: (Elem e) => BMatrix t (m,n) e -> Int -> (Int, DVector t k e, Int)+rowView a = checkedRow (shape a) (unsafeRowView a) ++colView :: (Elem e) => BMatrix t (m,n) e -> Int -> (Int, DVector t k e, Int)+colView a = checkedCol (shape a) (unsafeColView a)+++instance C.Matrix (BMatrix t) where+ numRows = fst . shape+ numCols = snd . shape++ herm a = case a of+ (H a') -> coerceMatrix a'+ _ -> H (coerceMatrix a)+ +instance Tensor (BMatrix t (m,n)) (Int,Int) e where+ shape a = case a of+ (H a') -> case shape a' of (m,n) -> (n,m)+ _ -> (size1 a, size2 a)++ bounds a = let (m,n) = shape a in ((0,0), (m-1,n-1))+ +instance (BLAS1 e) => ITensor (BMatrix Imm (m,n)) (Int,Int) e where+ size = inlinePerformIO . getSize+ + unsafeAt a = inlinePerformIO . (unsafeReadElem a)+ + indices = inlinePerformIO . getIndices+ elems = inlinePerformIO . getElems+ assocs = inlinePerformIO . getAssocs++ (//) = replaceHelp writeElem+ unsafeReplace = replaceHelp unsafeWriteElem++ amap f a = banded (shape a) (bandwidth a) ies+ where+ ies = map (second f) (assocs a)++replaceHelp :: (BLAS1 e) => + (IOBanded (m,n) e -> (Int,Int) -> e -> IO ())+ -> Banded (m,n) e -> [((Int,Int), e)] -> Banded (m,n) e+replaceHelp set x ies =+ unsafeFreeze $ unsafePerformIO $ do+ y <- newCopy (unsafeThaw x)+ mapM_ (uncurry $ set y) ies+ return y+{-# NOINLINE replaceHelp #-}+ + + +instance (BLAS1 e) => RTensor (BMatrix t (m,n)) (Int,Int) e IO where+ newCopy b = + let (a,mn,kl) = toMatrix b+ in do+ a' <- newCopy a+ return $ fromMatrix a' mn kl+ + getSize a = case a of+ (H a') -> getSize a'+ (BM _ _ m n kl ku _) -> return $ foldl' (+) 0 $ + map (diagLen (m,n)) [(-kl)..ku]+ + unsafeReadElem a (i,j) = case a of+ (H a') -> unsafeReadElem a' (j,i) >>= return . E.conj+ _ -> withForeignPtr (fptr a) $ \ptr ->+ peekElemOff ptr (indexOf a (i,j))++ getIndices a =+ return $ filter (\ij -> inlinePerformIO $ canModifyElem (unsafeThaw a) ij)+ (range $ bounds a)++ getElems a = getAssocs a >>= return . (map snd)++ getAssocs a = do+ is <- unsafeInterleaveIO $ getIndices a+ mapM (\i -> unsafeReadElem a i >>= \e -> return (i,e)) is++instance (BLAS1 e) => MTensor (BMatrix Mut (m,n)) (Int,Int) e IO where+ setZero a = case toMatrix a of (a',_,_) -> setZero a'+ setConstant e a = case toMatrix a of (a',_,_) -> setConstant e a'+ + canModifyElem a (i,j) = case a of+ (H a') -> canModifyElem a' (j,i)+ (BM _ _ m n kl ku _) -> return $ inRange (0,m-1) i && + inRange (0,n-1) j &&+ inRange (max 0 (j-ku), min (m-1) (j+kl)) i+ + unsafeWriteElem a (i,j) e = case a of+ (H a') -> unsafeWriteElem a' (j,i) (E.conj e)+ _ -> withForeignPtr (fptr a) $ \ptr ->+ pokeElemOff ptr (indexOf a (i,j)) e++ modifyWith f a = case toMatrix a of (a',_,_) -> modifyWith f a'
+ INSTALL view
@@ -0,0 +1,100 @@++ INSTALLING THE blas LIBRARY++This document tells you how to build and install the blas library, and how+to build and run the tests for the library. ++TABLE OF CONTENTS+I. Prerequisites+II. Configuring+III. Building and Installing+IV. Verifying the Installation+V. Building and Running the Tests+VI. Custom CBLAS Configurations+++I. PREREQUISITES++Installing blas requires a working CBLAS to be installed on your system. If you+do not have a CBLAS installed and you are not too concerned about performance,+probably the easiest one to install is the one that comes with the GNU +Scientific Library (GSL), available at http://www.gnu.org/software/gsl/. If+you care about performance, a better option is to use ATLAS, available at+http://math-atlas.sourceforge.net/. Other options include the Goto BLAS and+Intel's MKL library. If you are running Mac OS X, you can use vecLib, which+is part of the Accelerate framework.+++II. CONFIGURING++To configure blas, first decide which CBLAS you want to link to and figure out+the appropriate flag. The default flags are as follows:++ CBLAS Vendor Flag+ ------------------------------------------------------+ ATLAS -fatlas+ GSL -fgsl+ Intel MKL -fmkl+ vecLib -fveclib+ Other vendor or non-standard installation -fcustom+ +Now, run the command "runhaskell Setup.lhs configure" with the appropriate +flag. For example, to link with ATLAS, you would run the command:++# runhaskell Setup.lhs configure -fatlas++You can also pass additional arguments to the configure script. To see a full+list, do "runhaskell Setup.lhs configure --help". If you do not pass a flag+to the configure script, ATLAS is used by default. Do not try passing more than+one CBLAS flag to the script.+++III. BUILDING AND INSTALLING++To build the library, execute the command++# runhaskell Setup.lhs build++To build the haddock documentation, run++# runhaskell Setup.lhs haddock++Lastly, to install the library, do++# runhaskell Setup.lhs install+++IV. VERIFYING THE INSTALLATION++Once you have installed blas, try running the command "ghci -package blas". You+should not get any "unknown symbol" errors.+++V. BUILDING AND RUNNING THE TESTS++If you have verified that blas is installed correctly, you can try building and+running the tests. These are in the "tests" directory, and are in the form of+QuickCheck properties. To build and run them, just run "make" in the tests+directory.+++VI. CUSTOM CBLAS CONFIGURATIONS++If your CBLAS is not one of the defaults listed above or it is not installed in+a standard location, you need to edit the blas.cabal file to tell the +installation script where it lives. At the end of the file, you will see the+following lines:++ if flag(custom)+ -- CUSTOM CBLAS LIBS GO HERE+ extra-libraries: + + -- PATH TO CUSTOM LIB DIR GOES HERE + extra-lib-dirs: + +You need to fill in the names and locations of your CBLAS libraries. If your+CBLAS is in a library with dependencies, make sure you list the dependencies+*after* the libraries that depend on them. ++Once you have edited blas.cabal, build and install as normal, using the +"-fcustom" flag when you configure the library.
Setup.lhs view
@@ -3,14 +3,14 @@ > import System.Cmd > > testing _ _ _ _ = do-> system "runhaskell -lblas -DREAL tests/Vector.hs"-> system "runhaskell -lblas -DCOMPLEX tests/Vector.hs"-> system "runhaskell -lblas -DREAL tests/Matrix.hs"-> system "runhaskell -lblas -DCOMPLEX tests/Matrix.hs"-> system "runhaskell -lblas -DREAL tests/HermMatrix.hs"-> system "runhaskell -lblas -DCOMPLEX tests/HermMatrix.hs"-> system "runhaskell -lblas -DREAL tests/TriMatrix.hs"-> system "runhaskell -lblas -DCOMPLEX tests/TriMatrix.hs"+> system "runhaskell -lcblas -DREAL tests/Vector.hs"+> system "runhaskell -lcblas -DCOMPLEX tests/Vector.hs"+> system "runhaskell -lcblas -DREAL tests/Matrix.hs"+> system "runhaskell -lcblas -DCOMPLEX tests/Matrix.hs"+> system "runhaskell -lcblas -DREAL tests/HermMatrix.hs"+> system "runhaskell -lcblas -DCOMPLEX tests/HermMatrix.hs"+> system "runhaskell -lcblas -DREAL tests/TriMatrix.hs"+> system "runhaskell -lcblas -DCOMPLEX tests/TriMatrix.hs" > return () > > main = defaultMainWithHooks defaultUserHooks{ runTests=testing }
blas.cabal view
@@ -1,5 +1,5 @@ name: blas-Version: 0.4+version: 0.4.1 homepage: http://stat.stanford.edu/~patperry/code/blas synopsis: Bindings to the BLAS library description:@@ -24,11 +24,40 @@ build-type: Custom tested-with: GHC == 6.8.2 -extra-source-files: tests/Matrix.hs+extra-source-files: INSTALL+ tests/Makefile+ tests/Matrix.hs tests/HermMatrix.hs tests/TriMatrix.hs tests/Vector.hs ++-- Below are the flags for specifying which CBLAS to link with. If no flag+-- is specified, the default is to use ATLAS. To use the "custom" flag, you+-- must edit the section at the end of the file with the name and location+-- of the CBLAS library you want to use.++flag atlas+ description: Link with ATLAS.+ default: False++flag gsl+ description: Link with GSL unoptimized CBLAS.+ default: False++flag mkl+ description: Link with Intel MKL.+ default: False++flag veclib+ description: Link with Mac OS X vecLib.+ default: False+ +flag custom+ description: Link with a custom CBLAS. + default: False+ -- You must edit this file below to use this option. + library exposed-modules: BLAS.Access BLAS.C@@ -58,6 +87,8 @@ BLAS.Types BLAS.Vector + Data.Matrix.Banded.Internal+ Data.Matrix.Dense Data.Matrix.Dense.IO Data.Matrix.Dense.Internal@@ -91,5 +122,30 @@ FunctionalDependencies, MultiParamTypeClasses build-depends: base, ieee, storable-complex, QuickCheck - extra-libraries: blas+ + if flag(atlas)+ extra-libraries: cblas atlas+ + if flag(gsl)+ extra-libraries: gslcblas++ if flag(mkl)+ if arch(x86_64)+ extra-libraries: mkl_lapack mkl_intel_lp64 mkl_sequential mkl_core+ else+ extra-libraries: mkl_lapack mkl_intel mkl_sequential mkl_core++ if flag(veclib)+ extra-libraries: cblas++ if flag(custom)+ -- CUSTOM CBLAS LIBS GO HERE+ extra-libraries: + + -- PATH TO CUSTOM LIB DIR GOES HERE + extra-lib-dirs: + + -- fall back to ATLAS if no flag is specified+ if !flag(atlas) && !flag(gsl) && !flag(mkl) && !flag(veclib) && !flag(custom)+ extra-libraries: cblas atlas
+ tests/Makefile view
@@ -0,0 +1,37 @@+# To build and run the tests, make sure you have blas installed and that you+# also have the "pqc" package installed.+++# Change NTHREADS to the number of cores on your machine to run the tests+# in parallel.+NTHREADS = 1+THREADS = $(NTHREADS) +RTS -N$(NTHREADS) -RTS++all:+ ghc --make -fforce-recomp -O -DREAL -threaded Vector.hs -o vector-real+ ./vector-real $(THREADS)+ + ghc --make -fforce-recomp -O -DCOMPLEX -threaded Vector.hs -o vector-cplx+ ./vector-cplx $(THREADS)+ + ghc --make -fforce-recomp -O -DREAL -threaded Matrix.hs -o matrix-real+ ./matrix-real $(THREADS)+ + ghc --make -fforce-recomp -O -DCOMPLEX -threaded Matrix.hs -o matrix-cplx+ ./matrix-cplx $(THREADS)+ + ghc --make -fforce-recomp -O -DREAL -threaded HermMatrix.hs -o herm-matrix-real+ ./herm-matrix-real $(THREADS)+ + ghc --make -fforce-recomp -O -DCOMPLEX -threaded HermMatrix.hs -o herm-matrix-cplx+ ./herm-matrix-cplx $(THREADS)+ + ghc --make -fforce-recomp -O -DREAL -threaded TriMatrix.hs -o tri-matrix-real+ ./tri-matrix-real $(THREADS)+ + ghc --make -fforce-recomp -O -DCOMPLEX -threaded TriMatrix.hs -o tri-matrix-cplx+ ./tri-matrix-cplx $(THREADS)++clean:+ rm -f *~ *.hi *.o vector-real vector-cplx matrix-real matrix-cplx \+ herm-matrix-real herm-matrix-cplx tri-matrix-real tri-matrix-cplx
tests/Matrix.hs view
@@ -7,7 +7,6 @@ -- Stability : experimental -- -import Debug.Trace ( trace ) import qualified Data.Array as Array import Data.Ix ( inRange, range )@@ -17,14 +16,12 @@ import Test.QuickCheck.Parallel hiding ( vector ) import qualified Test.QuickCheck as QC -import BLAS.Access import Data.Matrix.Dense-import Data.Vector.Dense.Internal ( DVector )-import Data.Matrix.Dense.Internal ( DMatrix ) import Data.Vector.Dense hiding ( shift, scale, invScale )-import qualified Data.Vector.Dense as V-import BLAS.Elem ( Elem, BLAS1 )++import BLAS.Elem ( BLAS1 ) import qualified BLAS.Elem as E+ import Data.Complex ( Complex(..) ) import Data.AEq@@ -35,8 +32,6 @@ import Test.QuickCheck.Vector.Dense hiding ( Pair ) import Test.QuickCheck.Matrix import Test.QuickCheck.Matrix.Dense--import Debug.Trace #ifdef COMPLEX@@ -54,11 +49,11 @@ arbitrary = arbitrary >>= \(TestComplex x) -> return x coarbitrary = coarbitrary . TestComplex -instance (Arbitrary e, BLAS1 e) => Arbitrary (DVector Imm n e) where+instance (Arbitrary e, BLAS1 e) => Arbitrary (Vector n e) where arbitrary = arbitrary >>= \(TestVector x) -> return x coarbitrary = coarbitrary . TestVector -instance (Arbitrary e, BLAS1 e) => Arbitrary (DMatrix Imm (m,n) e) where+instance (Arbitrary e, BLAS1 e) => Arbitrary (Matrix (m,n) e) where arbitrary = arbitrary >>= \(TestMatrix x) -> return x coarbitrary = coarbitrary . TestMatrix @@ -96,12 +91,12 @@ diag (identity (m,n) :: M) 0 === (constant (min m n) 1) prop_identity_row (Basis m i) (Basis n _) = if i < min m n - then row (identity (m,n) :: M) i === V.basis n i- else row (identity (m,n) :: M) i === V.zero n+ then row (identity (m,n) :: M) i === basis n i+ else row (identity (m,n) :: M) i === zero n prop_identity_col (Basis m _) (Basis n j) = if j < min m n- then col (identity (m,n) :: M) j === V.basis m j- else col (identity (m,n) :: M) j === V.zero m+ then col (identity (m,n) :: M) j === basis m j+ else col (identity (m,n) :: M) j === zero m prop_replace_elems (a :: M) (Assocs _ ijes) = let ijes' = filter (\((i,j),_) -> i < numRows a && j < numCols a) ijes@@ -115,9 +110,9 @@ prop_submatrix_shape (SubMatrix a ij mn) = shape (submatrix a ij mn :: M) == mn prop_submatrix_rows (SubMatrix a (i,j) (m,n)) =- rows (submatrix a (i,j) (m,n) :: M) === map (\k -> V.subvector (row a (i+k)) j n) [0..(m-1)]+ rows (submatrix a (i,j) (m,n) :: M) === map (\k -> subvector (row a (i+k)) j n) [0..(m-1)] prop_submatrix_cols (SubMatrix a (i,j) (m,n)) (Index l) =- cols (submatrix a (i,j) (m,n) :: M) === map (\l -> V.subvector (col a (j+l)) i m) [0..(n-1)]+ cols (submatrix a (i,j) (m,n) :: M) === map (\l -> subvector (col a (j+l)) i m) [0..(n-1)] prop_shape (a :: M) = shape a == (numRows a, numCols a)@@ -132,65 +127,65 @@ in (a!ij) === ((elems a) !! k) prop_row_dim (MatAt (a :: M) (i,_)) =- V.dim (row a i) == numCols a+ dim (row a i) == numCols a prop_col_dim (MatAt (a :: M) (_,j)) =- V.dim (col a j) == numRows a+ dim (col a j) == numRows a prop_rows_len (a :: M) = length (rows a) == numRows a prop_cols_len (a :: M) = length (cols a) == numCols a prop_rows_dims (a :: M) =- map (V.dim) (rows a) == replicate (numRows a) (numCols a)+ map dim (rows a) == replicate (numRows a) (numCols a) prop_cols_dims (a :: M) =- map (V.dim) (cols a) == replicate (numCols a) (numRows a)+ map dim (cols a) == replicate (numCols a) (numRows a) prop_indices (a :: M) = let (m,n) = shape a in indices a == [(i,j) | j <- range (0,n-1), i <- range(0,m-1)] prop_elems (a :: M) =- and $ zipWith (===) (elems a) $ concatMap V.elems (cols a)+ and $ zipWith (===) (elems a) $ concatMap elems (cols a) prop_assocs (a :: M) = assocs a === zip (indices a) (elems a) prop_scale_elems (a :: M) k =- and $ zipWith (===) (elems (scale k a)) (map (k*) (elems a))+ and $ zipWith (~==) (elems (k *> a)) (map (k*) (elems a)) prop_herm_elem (MatAt (a :: M) (i,j)) = (herm a) ! (j,i) == E.conj (a!(i,j)) prop_herm_scale (a :: M) k =- herm (scale k a) === scale (E.conj k) (herm a)+ herm (k *> a) === (E.conj k) *> (herm a) prop_herm_shape (a :: M) = shape (herm a) == (numCols a, numRows a) prop_herm_rows (a :: M) =- rows (herm a) === map (V.conj) (cols a)+ rows (herm a) === map conj (cols a) prop_herm_cols (a :: M) = - cols (herm a) === map (V.conj) (rows a)+ cols (herm a) === map conj (rows a) prop_herm_herm (a :: M) = herm (herm a) === a prop_diag_herm1 (MatAt (a :: M) (k,_)) =- diag a (-k) === V.conj (diag (herm a) k)+ diag a (-k) === conj (diag (herm a) k) prop_diag_herm2 (MatAt (a :: M) (_,k)) =- diag a k === V.conj (diag (herm a) (-k))+ diag a k === conj (diag (herm a) (-k)) prop_fromRow_shape (x :: V) =- shape (fromRow x :: M) == (1,V.dim x)+ shape (fromRow x :: M) == (1,dim x) prop_fromRow_elems (x :: V) =- elems (fromRow x :: M) === V.elems x+ elems (fromRow x :: M) === elems x prop_fromCol_shape (x :: V) =- shape (fromCol x :: M) == (V.dim x,1)+ shape (fromCol x :: M) == (dim x,1) prop_fromCol_elems (x :: V) =- elems (fromCol x :: M) === V.elems x+ elems (fromCol x :: M) === elems x prop_apply_basis (MatAt (a :: M) (_,j)) =- a <*> (V.basis (numCols a) j :: V) ~== col a j+ a <*> (basis (numCols a) j :: V) ~== col a j prop_apply_herm_basis (MatAt (a :: M) (i,_)) =- (herm a) <*> (V.basis (numRows a) i :: V) ~== V.conj (row a i)+ (herm a) <*> (basis (numRows a) i :: V) ~== conj (row a i) prop_apply_scale k (MultMV (a :: M) x) =- a <*> (V.scale k x) ~== V.scale k (a <*> x)+ a <*> (k *> x) ~== k *> (a <*> x) prop_apply_linear (MultMVPair (a :: M) x y) = a <*> (x + y) ~== a <*> x + a <*> y @@ -214,7 +209,7 @@ prop_shift k (a :: M) = shift k a ~== a + constant (shape a) k prop_scale k (a :: M) =- scale k a ~== a * constant (shape a) k+ k *> a ~== a * constant (shape a) k prop_invScale k (a :: M) = invScale k a ~== a / constant (shape a) k @@ -228,7 +223,7 @@ elems (a / b) ~== zipWith (/) (elems a) (elems b) prop_negate (a :: M) =- negate a ~== scale (-1) a+ negate a ~== (-1) *> a prop_abs (a :: M) = elems (abs a) ~== map abs (elems a)
tests/Vector.hs view
@@ -29,8 +29,6 @@ import Test.QuickCheck.Vector.Dense -import Debug.Trace- #ifdef COMPLEX field = "Complex Double" type E = Complex Double