diff --git a/BLAS/C/Level1.hs b/BLAS/C/Level1.hs --- a/BLAS/C/Level1.hs +++ b/BLAS/C/Level1.hs @@ -1,4 +1,5 @@ {-# LANGUAGE FlexibleInstances #-} +{-# OPTIONS_GHC -fno-excess-precision #-} ----------------------------------------------------------------------------- -- | -- Module : BLAS.C.Level1 diff --git a/Data/Matrix/Banded/Internal.hs b/Data/Matrix/Banded/Internal.hs new file mode 100644 --- /dev/null +++ b/Data/Matrix/Banded/Internal.hs @@ -0,0 +1,419 @@ +{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Data.Matrix.Banded.Internal +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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' diff --git a/INSTALL b/INSTALL new file mode 100644 --- /dev/null +++ b/INSTALL @@ -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. diff --git a/Setup.lhs b/Setup.lhs --- a/Setup.lhs +++ b/Setup.lhs @@ -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 } diff --git a/blas.cabal b/blas.cabal --- a/blas.cabal +++ b/blas.cabal @@ -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 diff --git a/tests/Makefile b/tests/Makefile new file mode 100644 --- /dev/null +++ b/tests/Makefile @@ -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 diff --git a/tests/Matrix.hs b/tests/Matrix.hs --- a/tests/Matrix.hs +++ b/tests/Matrix.hs @@ -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) diff --git a/tests/Vector.hs b/tests/Vector.hs --- a/tests/Vector.hs +++ b/tests/Vector.hs @@ -29,8 +29,6 @@ import Test.QuickCheck.Vector.Dense -import Debug.Trace - #ifdef COMPLEX field = "Complex Double" type E = Complex Double