blas (empty) → 0.4
raw patch · 55 files changed
+5626/−0 lines, 55 filesdep +QuickCheckdep +basedep +ieeebuild-type:Customsetup-changed
Dependencies added: QuickCheck, base, ieee, storable-complex
Files
- BLAS/Access.hs +20/−0
- BLAS/C.hs +20/−0
- BLAS/C/Double.hs +112/−0
- BLAS/C/Level1.hs +125/−0
- BLAS/C/Level2.hs +82/−0
- BLAS/C/Level3.hs +82/−0
- BLAS/C/Types.hs +87/−0
- BLAS/C/Zomplex.hs +122/−0
- BLAS/Elem.hs +20/−0
- BLAS/Elem/Base.hs +51/−0
- BLAS/Internal.hs +198/−0
- BLAS/Matrix.hs +20/−0
- BLAS/Matrix/Base.hs +26/−0
- BLAS/Matrix/Immutable.hs +34/−0
- BLAS/Matrix/ReadOnly.hs +30/−0
- BLAS/Matrix/Solve.hs +16/−0
- BLAS/Matrix/Solve/Immutable.hs +28/−0
- BLAS/Matrix/Solve/ReadOnly.hs +25/−0
- BLAS/Tensor.hs +24/−0
- BLAS/Tensor/Base.hs +25/−0
- BLAS/Tensor/Dense.hs +16/−0
- BLAS/Tensor/Dense/Immutable.hs +23/−0
- BLAS/Tensor/Dense/ReadOnly.hs +25/−0
- BLAS/Tensor/Immutable.hs +70/−0
- BLAS/Tensor/Mutable.hs +76/−0
- BLAS/Tensor/ReadOnly.hs +62/−0
- BLAS/Tensor/Scalable.hs +20/−0
- BLAS/Types.hs +53/−0
- BLAS/Vector.hs +25/−0
- Data/Matrix/Dense.hs +168/−0
- Data/Matrix/Dense/IO.hs +89/−0
- Data/Matrix/Dense/Internal.hs +536/−0
- Data/Matrix/Dense/Operations.hs +347/−0
- Data/Matrix/Herm.hs +64/−0
- Data/Matrix/Herm/Dense.hs +114/−0
- Data/Matrix/Tri.hs +73/−0
- Data/Matrix/Tri/Dense.hs +171/−0
- Data/Vector/Dense.hs +135/−0
- Data/Vector/Dense/IO.hs +64/−0
- Data/Vector/Dense/Internal.hs +422/−0
- Data/Vector/Dense/Operations.hs +373/−0
- LICENSE +30/−0
- Setup.lhs +17/−0
- Test/QuickCheck/Complex.hs +24/−0
- Test/QuickCheck/Matrix.hs +41/−0
- Test/QuickCheck/Matrix/Dense.hs +157/−0
- Test/QuickCheck/Matrix/Herm/Dense.hs +72/−0
- Test/QuickCheck/Matrix/Tri/Dense.hs +132/−0
- Test/QuickCheck/Vector.hs +50/−0
- Test/QuickCheck/Vector/Dense.hs +92/−0
- blas.cabal +95/−0
- tests/HermMatrix.hs +96/−0
- tests/Matrix.hs +334/−0
- tests/TriMatrix.hs +157/−0
- tests/Vector.hs +256/−0
+ BLAS/Access.hs view
@@ -0,0 +1,20 @@+{-# LANGUAGE EmptyDataDecls #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Access+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--+-- Types for access control.+module BLAS.Access (+ Mut,+ Imm+ ) where++-- | Tag for mutable types.+data Mut++-- | Tag for immutable types.+data Imm
+ BLAS/C.hs view
@@ -0,0 +1,20 @@+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.C+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.C (+ module BLAS.C.Types,+ module BLAS.C.Level1,+ module BLAS.C.Level2,+ module BLAS.C.Level3,+ ) where+ +import BLAS.C.Types+import BLAS.C.Level1+import BLAS.C.Level2+import BLAS.C.Level3
+ BLAS/C/Double.hs view
@@ -0,0 +1,112 @@+{-# LANGUAGE ForeignFunctionInterface #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.C.Double+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.C.Double + where+ +import Foreign.Ptr ( Ptr )+import BLAS.C.Types+++---------------------------- Level 1 Routines -------------------------------++foreign import ccall unsafe "cblas.h cblas_ddot"+ ddot :: Int -> Ptr Double -> Int -> Ptr Double -> Int -> IO Double++foreign import ccall unsafe "cblas.h cblas_dnrm2"+ dnrm2 :: Int -> Ptr Double -> Int -> IO Double++foreign import ccall unsafe "cblas.h cblas_dasum"+ dasum :: Int -> Ptr Double -> Int -> IO Double++foreign import ccall unsafe "cblas.h cblas_idamax"+ idamax :: Int -> Ptr Double -> Int -> IO Int++foreign import ccall unsafe "cblas.h cblas_dscal"+ dscal :: Int -> Double -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_dswap"+ dswap :: Int -> Ptr Double -> Int -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_dcopy"+ dcopy :: Int -> Ptr Double -> Int -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_daxpy"+ daxpy :: Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_drotg"+ drotg :: Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> IO ()++foreign import ccall unsafe "cblas.h cblas_drot"+ drot :: Int -> Ptr Double -> Int -> Ptr Double -> Int -> Double -> Double -> IO ()++foreign import ccall unsafe "cblas.h cblas_drotmg"+ drotmg :: Ptr Double -> Ptr Double -> Ptr Double -> Double -> Ptr Double -> IO ()++foreign import ccall unsafe "cblas.h cblas_drotm"+ drotm :: Int -> Ptr Double -> Int -> Ptr Double -> Int -> Ptr Double -> IO ()+++---------------------------- Level 2 Routines -------------------------------++foreign import ccall unsafe "cblas.h cblas_dgemv"+ dgemv :: CBLASOrder -> CBLASTrans -> Int -> Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> Double -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_dgbmv"+ dgbmv :: CBLASOrder -> CBLASTrans -> Int -> Int -> Int -> Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> Double -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_dtrmv"+ dtrmv :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Ptr Double -> Int -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_dtbmv"+ dtbmv :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Int -> Ptr Double -> Int -> Ptr Double -> Int -> IO ()+ +foreign import ccall unsafe "cblas.h cblas_dtrsv"+ dtrsv :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Ptr Double -> Int -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_dtbsv"+ dtbsv :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Int -> Ptr Double -> Int -> Ptr Double -> Int -> IO ()+ +foreign import ccall unsafe "cblas.h cblas_dsymv"+ dsymv :: CBLASOrder -> CBLASUpLo -> Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> Double -> Ptr Double -> Int -> IO ()+ +foreign import ccall unsafe "cblas.h cblas_dsbmv"+ dsbmv :: CBLASOrder -> CBLASUpLo -> Int -> Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> Double -> Ptr Double -> Int -> IO ()+ +foreign import ccall unsafe "cblas.h cblas_dger"+ dger :: CBLASOrder -> Int -> Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> Ptr Double -> Int -> IO ()+ +foreign import ccall unsafe "cblas.h cblas_dsyr"+ dsyr :: CBLASOrder -> CBLASUpLo -> Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_dsyr2"+ dsyr2 :: CBLASOrder -> CBLASUpLo -> Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> Ptr Double -> Int -> IO ()+++---------------------------- Level 3 Routines -------------------------------++foreign import ccall unsafe "cblas.h cblas_dgemm"+ dgemm :: CBLASOrder -> CBLASTrans -> CBLASTrans -> Int -> Int -> Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> Double -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_dsymm"+ dsymm :: CBLASOrder -> CBLASSide -> CBLASUpLo -> Int -> Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> Double -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_dtrmm"+ dtrmm :: CBLASOrder -> CBLASSide -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_dtrsm"+ dtrsm :: CBLASOrder -> CBLASSide -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_dsyrk"+ dsyrk :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> Int -> Int -> Double -> Ptr Double -> Int -> Double -> Ptr Double -> Int -> IO ()+ +foreign import ccall unsafe "cblas.h cblas_dsyr2k" + dsyr2k :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> Int -> Int -> Double -> Ptr Double -> Int -> Ptr Double -> Int -> Double -> Ptr Double -> Int -> IO ()+
+ BLAS/C/Level1.hs view
@@ -0,0 +1,125 @@+{-# LANGUAGE FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.C.Level1+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.C.Level1+ where+ +import Foreign ( Ptr, Storable, advancePtr, castPtr, peek, poke, with )+import Foreign.Storable.Complex ()+import Data.Complex++import BLAS.Elem.Base+import BLAS.C.Double +import BLAS.C.Zomplex+ +class (Elem a) => BLAS1 a where+ dotu :: Int -> Ptr a -> Int -> Ptr a -> Int -> IO a+ dotc :: Int -> Ptr a -> Int -> Ptr a -> Int -> IO a+ nrm2 :: Int -> Ptr a -> Int -> IO Double+ asum :: Int -> Ptr a -> Int -> IO Double+ iamax :: Int -> Ptr a -> Int -> IO Int+ + scal :: Int -> a -> Ptr a -> Int -> IO () ++ swap :: Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()+ copy :: Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()+ axpy :: Int -> a -> Ptr a -> Int -> Ptr a -> Int -> IO ()++ rotg :: Ptr a -> Ptr a -> Ptr a -> Ptr a -> IO ()+ rot :: Int -> Ptr a -> Int -> Ptr a -> Int -> Double -> Double -> IO ()++ -- | conjugate all elements of a vector+ conj :: Int -> Ptr a -> Int -> IO ()++ -- | Replaces @y@ with @alpha (conj x) + y@+ acxpy :: Int -> a -> Ptr a -> Int -> Ptr a -> Int -> IO ()+++instance BLAS1 Double where+ dotu = ddot+ dotc = ddot+ nrm2 = dnrm2+ asum = dasum+ iamax = idamax+ swap = dswap+ copy = dcopy+ axpy = daxpy+ scal = dscal+ rotg = drotg+ rot = drot+ conj _ _ _ = return ()+ acxpy = daxpy++instance BLAS1 (Complex Double) where+ dotu n pX incX pY incY =+ with 0 $ \pDotu -> do+ zdotu_sub n pX incX pY incY pDotu+ peek pDotu++ dotc n pX incX pY incY =+ with 0 $ \pDotc -> do+ zdotc_sub n pX incX pY incY pDotc+ peek pDotc++ nrm2 = znrm2+ asum = zasum+ iamax = izamax+ swap = zswap+ copy = zcopy+ + axpy n alpha pX incX pY incY = + with alpha $ \pAlpha ->+ zaxpy n pAlpha pX incX pY incY+ + scal n alpha pX incX =+ with alpha $ \pAlpha ->+ zscal n pAlpha pX incX+ + rotg = zrotg++ rot = zdrot++ conj n pX incX =+ let pXI = (castPtr pX) `advancePtr` 1+ alpha = -1+ incXI = 2 * incX+ in dscal n alpha pXI incXI+ + acxpy n a pX incX pY incY =+ let pXR = castPtr pX+ pYR = castPtr pY+ pXI = pXR `advancePtr` 1+ pYI = pYR `advancePtr` 1+ incX' = 2 * incX+ incY' = 2 * incY+ in case a of+ (ra :+ 0) -> do+ daxpy n ( ra) pXR incX' pYR incY'+ daxpy n (-ra) pXI incX' pYI incY'+ (0 :+ ia) -> do+ daxpy n ( ia) pXR incX' pYI incY'+ daxpy n ( ia) pXI incX' pYR incY'+ _ -> go n pX pY+ where+ go n' pX' pY'+ | n' `seq` pX' `seq` pY' `seq` False = undefined+ | n' <= 0 =+ return ()+ | otherwise = do+ x <- peek pX'+ y <- peek pY'+ poke pY' (a * (conjugate x) + y)+ + let n'' = n' - 1+ pX'' = pX' `advancePtr` incX+ pY'' = pY' `advancePtr` incY+ + go n'' pX'' pY''+
+ BLAS/C/Level2.hs view
@@ -0,0 +1,82 @@+{-# LANGUAGE FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.C.Level2+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.C.Level2+ where+ +import Data.Complex +import Foreign ( Ptr, with ) ++import BLAS.C.Types+import BLAS.C.Level1+import BLAS.C.Double +import BLAS.C.Zomplex+ +class (BLAS1 a) => BLAS2 a where+ gemv :: CBLASOrder -> CBLASTrans -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+ gbmv :: CBLASOrder -> CBLASTrans -> Int -> Int -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+ trmv :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()+ tbmv :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()+ trsv :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()+ tbsv :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()+ hemv :: CBLASOrder -> CBLASUpLo -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+ hbmv :: CBLASOrder -> CBLASUpLo -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+ geru :: CBLASOrder -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()+ gerc :: CBLASOrder -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()+ her :: CBLASOrder -> CBLASUpLo -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> IO ()+ her2 :: CBLASOrder -> CBLASUpLo -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()++instance BLAS2 Double where+ gemv = dgemv+ gbmv = dgbmv+ trmv = dtrmv+ tbmv = dtbmv+ trsv = dtrsv+ tbsv = dtbsv+ hemv = dsymv+ hbmv = dsbmv+ geru = dger+ gerc = dger+ her = dsyr+ her2 = dsyr2+ +instance BLAS2 (Complex Double) where+ gemv order transA m n alpha pA ldA pX incX beta pY incY =+ with alpha $ \pAlpha -> with beta $ \pBeta ->+ zgemv order transA m n pAlpha pA ldA pX incX pBeta pY incY+ + gbmv order transA m n kl ku alpha pA ldA pX incX beta pY incY =+ with alpha $ \pAlpha -> with beta $ \pBeta ->+ zgbmv order transA m n kl ku pAlpha pA ldA pX incX pBeta pY incY++ trmv = ztrmv+ tbmv = ztbmv+ trsv = ztrsv+ tbsv = ztbsv + + hemv order uplo n alpha pA ldA pX incX beta pY incY =+ with alpha $ \pAlpha -> with beta $ \pBeta -> + zhemv order uplo n pAlpha pA ldA pX incX pBeta pY incY+ + hbmv order uplo n k alpha pA ldA pX incX beta pY incY =+ with alpha $ \pAlpha -> with beta $ \pBeta -> + zhbmv order uplo n k pAlpha pA ldA pX incX pBeta pY incY++ geru order m n alpha pX incX pY incY pA ldA = + with alpha $ \pAlpha -> zgeru order m n pAlpha pX incX pY incY pA ldA++ gerc order m n alpha pX incX pY incY pA ldA = + with alpha $ \pAlpha -> zgerc order m n pAlpha pX incX pY incY pA ldA++ her order uplo n alpha pX incX pA ldA = + with alpha $ \pAlpha -> zher order uplo n pAlpha pX incX pA ldA+ + her2 order uplo n alpha pX incX pY incY pA ldA = + with alpha $ \pAlpha -> zher2 order uplo n pAlpha pX incX pY incY pA ldA
+ BLAS/C/Level3.hs view
@@ -0,0 +1,82 @@+{-# LANGUAGE FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.C.Level3+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.C.Level3+ where+ +import Data.Complex +import Foreign ( Ptr, with ) ++import BLAS.C.Types+import BLAS.C.Level2+import BLAS.C.Double +import BLAS.C.Zomplex + + +class (BLAS2 a) => BLAS3 a where+ gemm :: CBLASOrder -> CBLASTrans -> CBLASTrans -> Int -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+ symm :: CBLASOrder -> CBLASSide -> CBLASUpLo -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+ hemm :: CBLASOrder -> CBLASSide -> CBLASUpLo -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+ trmm :: CBLASOrder -> CBLASSide -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> IO ()+ trsm :: CBLASOrder -> CBLASSide -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> IO ()+ syrk :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> Int -> Int -> a -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+ syr2k :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+ herk :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> Int -> Int -> a -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+ her2k :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+ + +instance BLAS3 Double where+ gemm = dgemm+ symm = dsymm+ hemm = dsymm+ trmm = dtrmm+ trsm = dtrsm+ syrk = dsyrk+ syr2k = dsyr2k+ herk = dsyrk+ her2k = dsyr2k+ + +instance BLAS3 (Complex Double) where+ gemm order transA transB m n k alpha pA ldA pB ldB beta pC ldC =+ with alpha $ \pAlpha -> with beta $ \pBeta ->+ zgemm order transA transB m n k pAlpha pA ldA pB ldB pBeta pC ldC+ + symm order side uplo m n alpha pA ldA pB ldB beta pC ldC =+ with alpha $ \pAlpha -> with beta $ \pBeta ->+ zsymm order side uplo m n pAlpha pA ldA pB ldB pBeta pC ldC++ hemm order side uplo m n alpha pA ldA pB ldB beta pC ldC =+ with alpha $ \pAlpha -> with beta $ \pBeta ->+ zhemm order side uplo m n pAlpha pA ldA pB ldB pBeta pC ldC+ + trmm order side uplo transA diag m n alpha pA ldA pB ldB =+ with alpha $ \pAlpha -> + ztrmm order side uplo transA diag m n pAlpha pA ldA pB ldB+ + trsm order side uplo transA diag m n alpha pA ldA pB ldB =+ with alpha $ \pAlpha -> + ztrsm order side uplo transA diag m n pAlpha pA ldA pB ldB+ + syrk order uplo transA n k alpha pA ldA beta pC ldC =+ with alpha $ \pAlpha -> with beta $ \pBeta ->+ zsyrk order uplo transA n k pAlpha pA ldA pBeta pC ldC+ + syr2k order uplo transA n k alpha pA ldA pB ldB beta pC ldC =+ with alpha $ \pAlpha -> with beta $ \pBeta ->+ zsyr2k order uplo transA n k pAlpha pA ldA pB ldB pBeta pC ldC++ herk order uplo transA n k alpha pA ldA beta pC ldC =+ with alpha $ \pAlpha -> with beta $ \pBeta ->+ zherk order uplo transA n k pAlpha pA ldA pBeta pC ldC+ + her2k order uplo transA n k alpha pA ldA pB ldB beta pC ldC =+ with alpha $ \pAlpha -> with beta $ \pBeta ->+ zher2k order uplo transA n k pAlpha pA ldA pB ldB pBeta pC ldC
+ BLAS/C/Types.hs view
@@ -0,0 +1,87 @@+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.C.Types+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.C.Types (+ CBLASOrder,+ CBLASTrans,+ CBLASUpLo,+ CBLASDiag,+ CBLASSide,++ rowMajor,+ colMajor,+ + noTrans,+ trans,+ conjTrans,+ + upper,+ lower,+ + nonUnit,+ unit,+ + leftSide,+ rightSide,++ cblasOrder,+ cblasTrans,+ cblasUpLo,+ cblasDiag,+ cblasSide,+ ) where+ +import BLAS.Types+ +newtype CBLASOrder = CBLASOrder Int deriving (Eq, Show)+newtype CBLASTrans = CBLASTrans Int deriving (Eq, Show)+newtype CBLASUpLo = CBLASUpLo Int deriving (Eq, Show)+newtype CBLASDiag = CBLASDiag Int deriving (Eq, Show)+newtype CBLASSide = CBLASSide Int deriving (Eq, Show)++rowMajor, colMajor :: CBLASOrder+rowMajor = CBLASOrder 101+colMajor = CBLASOrder 102++cblasOrder :: Order -> CBLASOrder+cblasOrder RowMajor = CBLASOrder 101+cblasOrder ColMajor = CBLASOrder 102++noTrans, trans, conjTrans :: CBLASTrans+noTrans = CBLASTrans 111+trans = CBLASTrans 112+conjTrans = CBLASTrans 113++cblasTrans :: Trans -> CBLASTrans+cblasTrans NoTrans = CBLASTrans 111+cblasTrans ConjTrans = CBLASTrans 113++upper, lower :: CBLASUpLo+upper = CBLASUpLo 121+lower = CBLASUpLo 122++cblasUpLo :: UpLo -> CBLASUpLo+cblasUpLo Upper = CBLASUpLo 121+cblasUpLo Lower = CBLASUpLo 122++nonUnit, unit :: CBLASDiag+nonUnit = CBLASDiag 131+unit = CBLASDiag 132++cblasDiag :: Diag -> CBLASDiag+cblasDiag NonUnit = CBLASDiag 131+cblasDiag Unit = CBLASDiag 132++leftSide, rightSide :: CBLASSide+leftSide = CBLASSide 141+rightSide = CBLASSide 142++cblasSide :: Side -> CBLASSide+cblasSide LeftSide = CBLASSide 141+cblasSide RightSide = CBLASSide 142
+ BLAS/C/Zomplex.hs view
@@ -0,0 +1,122 @@+{-# LANGUAGE ForeignFunctionInterface #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.C.Zomplex+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.C.Zomplex+ where+ +import Data.Complex ( Complex )+import Foreign.Ptr ( Ptr )+import BLAS.C.Types++---------------------------- Level 1 Routines -------------------------------++foreign import ccall unsafe "cblas.h cblas_zdotu_sub"+ zdotu_sub :: Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> IO ()++foreign import ccall unsafe "cblas.h cblas_zdotc_sub"+ zdotc_sub :: Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> IO ()+++foreign import ccall unsafe "cblas.h cblas_dznrm2"+ znrm2 :: Int -> Ptr (Complex Double) -> Int -> IO Double++foreign import ccall unsafe "cblas.h cblas_dzasum"+ zasum :: Int -> Ptr (Complex Double) -> Int -> IO Double++foreign import ccall unsafe "cblas.h cblas_izamax"+ izamax :: Int -> Ptr (Complex Double) -> Int -> IO Int++foreign import ccall unsafe "cblas.h cblas_zscal"+ zscal :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_zswap"+ zswap :: Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_zcopy"+ zcopy :: Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_zaxpy"+ zaxpy :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_zrotg"+ zrotg :: Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()++foreign import ccall unsafe "cblas.h cblas_zdrot"+ zdrot :: Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Double -> Double -> IO ()+++---------------------------- Level 2 Routines -------------------------------++foreign import ccall unsafe "cblas.h cblas_zgemv"+ zgemv :: CBLASOrder -> CBLASTrans -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_zgbmv"+ zgbmv :: CBLASOrder -> CBLASTrans -> Int -> Int -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_ztrmv"+ ztrmv :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_ztbmv"+ ztbmv :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()+ +foreign import ccall unsafe "cblas.h cblas_ztrsv"+ ztrsv :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_ztbsv"+ ztbsv :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()+ +foreign import ccall unsafe "cblas.h cblas_zhemv"+ zhemv :: CBLASOrder -> CBLASUpLo -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_zhbmv"+ zhbmv :: CBLASOrder -> CBLASUpLo -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> IO ()+ +foreign import ccall unsafe "cblas.h cblas_zgeru"+ zgeru :: CBLASOrder -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_zgerc"+ zgerc :: CBLASOrder -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()+ +foreign import ccall unsafe "cblas.h cblas_zher"+ zher :: CBLASOrder -> CBLASUpLo -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_zher2"+ zher2 :: CBLASOrder -> CBLASUpLo -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()+++---------------------------- Level 3 Routines -------------------------------++foreign import ccall unsafe "cblas.h cblas_zgemm"+ zgemm :: CBLASOrder -> CBLASTrans -> CBLASTrans -> Int -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_zsymm"+ zsymm :: CBLASOrder -> CBLASSide -> CBLASUpLo -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_zhemm"+ zhemm :: CBLASOrder -> CBLASSide -> CBLASUpLo -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_ztrmm"+ ztrmm :: CBLASOrder -> CBLASSide -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_ztrsm"+ ztrsm :: CBLASOrder -> CBLASSide -> CBLASUpLo -> CBLASTrans -> CBLASDiag -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_zsyrk"+ zsyrk :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> IO ()+ +foreign import ccall unsafe "cblas.h cblas_zsyr2k" + zsyr2k :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> IO ()++foreign import ccall unsafe "cblas.h cblas_zherk"+ zherk :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> IO ()+ +foreign import ccall unsafe "cblas.h cblas_zher2k" + zher2k :: CBLASOrder -> CBLASUpLo -> CBLASTrans -> Int -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Int -> IO ()+
+ BLAS/Elem.hs view
@@ -0,0 +1,20 @@+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Elem+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Elem (+ module BLAS.Elem.Base,+ module BLAS.C.Level1,+ module BLAS.C.Level2,+ module BLAS.C.Level3+ ) where++import BLAS.Elem.Base+import BLAS.C.Level1 ( BLAS1 )+import BLAS.C.Level2 ( BLAS2 )+import BLAS.C.Level3 ( BLAS3 )
+ BLAS/Elem/Base.hs view
@@ -0,0 +1,51 @@+{-# LANGUAGE FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Elem.Base+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Elem.Base (+ Elem(..)+ ) where++import Data.Complex ( Complex(..), conjugate, magnitude )+import Foreign ( Storable )+import Foreign.Storable.Complex ()++-- | The base class for elements.+class (Storable e, Fractional e) => Elem e where+ -- | Take the complex conjugate of a value. For real values+ -- this is equal to @id@.+ conj :: e -> e+ + -- | Get the magnitude of a value.+ norm :: e -> Double+ + -- | Get the l1 norm of a value.+ norm1 :: e -> Double+ + -- | Convert a double to an element+ fromReal :: Double -> e++ -- | Coerce an element to a double+ toReal :: e -> Double+ ++instance Elem Double where+ conj = id+ norm = abs+ norm1 = abs+ fromReal = id+ toReal = id+ +instance Elem (Complex Double) where+ conj = conjugate+ norm = magnitude+ norm1 (x :+ y) = abs x + abs y+ fromReal x = x :+ 0+ toReal (x :+ _) = x+
+ BLAS/Internal.hs view
@@ -0,0 +1,198 @@+{-# LANGUAGE CPP, ForeignFunctionInterface #-}+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Internal+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--+++module BLAS.Internal (+ clearArray,+ bzero,+ inlinePerformIO,+ checkedSubvector,+ checkedSubvectorWithStride,+ checkVecVecOp,+ checkedRow,+ checkedCol,+ checkedDiag,+ checkedSubmatrix,+ checkMatMatOp,+ checkMatVecMult,+ checkMatMatMult,+ diagStart,+ diagLen, + ) where++import Data.Ix ( inRange )+import Foreign ( Ptr, Storable, castPtr, sizeOf )+import Foreign.C.Types ( CSize )+import Text.Printf ( printf )++#if defined(__GLASGOW_HASKELL__)+import GHC.Base ( realWorld# )+import GHC.IOBase ( IO(IO) )+#else+import System.IO.Unsafe ( unsafePerformIO )+#endif+++clearArray :: Storable e => Ptr e -> Int -> IO ()+clearArray = clearArray' undefined+ where+ clearArray' :: Storable e => e -> Ptr e -> Int -> IO ()+ clearArray' e ptr n =+ let nbytes = (fromInteger . toInteger) (n * sizeOf e)+ in do+ bzero ptr nbytes+{-# INLINE clearArray #-}+++bzero :: Ptr a -> Int -> IO ()+bzero ptr n =+ let ptr' = castPtr ptr+ n' = (fromInteger . toInteger) n+ in bzero_ ptr' n'+ +foreign import ccall "strings.h bzero"+ bzero_ :: Ptr () -> CSize -> IO ()+ ++inlinePerformIO :: IO a -> a+#if defined(__GLASGOW_HASKELL__)+inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r+#else+inlinePerformIO = unsafePerformIO+#endif+{-# INLINE inlinePerformIO #-}++checkedSubvector :: Int -> (Int -> Int -> v) -> Int -> Int -> v+checkedSubvector n sub o n'+ | (o < 0) && (n' /= 0) = + error $ printf + "tried to create a subvector starting at a negative offset: `%d'" o+ | n' < 0 = + error $ printf + "tried to create a subvector with a negative length `%d'" n'+ | n' + o > n = + error $ printf+ ("tried to create a subvector of length `%d' and offset `%d' "+ ++ " from a vector of length `%d'") n' o n+ | otherwise =+ sub o n'+ +checkedSubvectorWithStride :: Int -> Int -> (Int -> Int -> v)+ -> Int -> Int -> v+checkedSubvectorWithStride s n sub o n'+ | (o < 0) && (n' /= 0) =+ error $ printf+ "Tried to create a subvector starting at a negative offset: `%d'" o+ | n' < 0 =+ error $ printf+ "Tried to create a subvector with a negative length `%d'" n'+ | s <= 0 =+ error $ printf+ "Tried to create a subvector with non-positive stride `%d'" s+ | not $ inRange (-1,n) (o + s * n') =+ error $ printf+ ("tried to create a subvector of length `%d', offset `%d',"+ ++ " and stride '%d' from a vector of length `%d'") n' o s n+ | otherwise =+ sub o n'++checkVecVecOp :: String -> Int -> Int -> IO ()+checkVecVecOp name n1 n2+ | n1 /= n2 =+ ioError $ userError $ printf+ ("%s: x and y have different dimensions. x has dimension `%d',"+ ++ " and y has dimension `%d'") name n1 n2+ | otherwise =+ return ()++checkedRow :: (Int,Int) -> (Int -> v) -> Int -> v+checkedRow (m,n) row i + | i < 0 || i >= m =+ error $ printf+ "Error in row index. Tried to get row `%d' in a matrix with shape `(%d,%d)'" i m n+ | otherwise =+ row i++checkedCol :: (Int,Int) -> (Int -> v) -> Int -> v+checkedCol (m,n) col j + | j < 0 || j >= n =+ error $ printf+ "Error in column index. Tried to get column `%d' in a matrix with shape `(%d,%d)'" j m n+ | otherwise =+ col j++checkedDiag :: (Int,Int) -> (Int -> v) -> Int -> v+checkedDiag (m,n) diag i+ | i < 0 && negate i >= m =+ error $ printf+ "Tried to get sub-diagonal `%d' of a matrix with shape `(%d,%d)'" (negate i) m n+ | i > 0 && i >= n =+ error $ printf+ "Tried to get super-diagonal `%d' of a matrix with shape `(%d,%d)'" i m n + | otherwise = + diag i++diagStart :: Int -> (Int,Int)+diagStart i+ | i <= 0 =+ (negate i, 0)+ | otherwise =+ (0, i)+ +diagLen :: (Int,Int) -> Int -> Int+diagLen (m,n) i+ | m <= n =+ if i <= 0 + then max (m + i) 0+ else min (n - i) m+ | otherwise =+ if i > 0+ then max (n - i) 0+ else min (m + i) n++checkedSubmatrix :: (Int,Int) -> ((Int,Int) -> (Int,Int) -> a) -> (Int,Int) -> (Int,Int) -> a+checkedSubmatrix (m,n) sub (i,j) (m',n')+ | or [ i < 0, m' < 0, i + m' > m, + j < 0, n' < 0, j + n' > n ] =+ error $ printf ("tried to create submatrix of a `(%d,%d)' matrix " +++ " using offset `(%d,%d)' and shape (%d,%d)") m n i j m' n'+ | otherwise =+ sub (i,j) (m',n')+++checkMatMatOp :: String -> (Int,Int) -> (Int,Int) -> IO ()+checkMatMatOp name mn1 mn2+ | mn1 /= mn2 =+ ioError $ userError $ printf+ ("%s: x and y have different shapes. x has shape `%s',"+ ++ " and y has shape `%s'") name (show mn1) (show mn2)+ | otherwise =+ return ()+ +checkMatVecMult :: (Int,Int) -> Int -> IO ()+checkMatVecMult mn n+ | snd mn /= n =+ ioError $ userError $ printf+ ("Tried to multiply a matrix with shape `%s' by a vector of dimension `%d'")+ (show mn) n+ | otherwise =+ return ()+ +checkMatMatMult :: (Int,Int) -> (Int,Int) -> IO ()+checkMatMatMult mk kn+ | snd mk /= fst kn =+ ioError $ userError $ printf+ ("Tried to multiply a matrix with shape `%s' by a matrix with shape `%s'")+ (show mk) (show kn)+ | otherwise =+ return ()++
+ BLAS/Matrix.hs view
@@ -0,0 +1,20 @@+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Matrix+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Matrix (+ module BLAS.Matrix.Base,+ module BLAS.Matrix.Immutable,+ module BLAS.Matrix.ReadOnly,+ module BLAS.Matrix.Solve,+ ) where++import BLAS.Matrix.Base+import BLAS.Matrix.Immutable+import BLAS.Matrix.ReadOnly+import BLAS.Matrix.Solve
+ BLAS/Matrix/Base.hs view
@@ -0,0 +1,26 @@+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Matrix.Base+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Matrix.Base (+ Matrix(..)+ ) where++import BLAS.Elem.Base ( Elem )++-- | A base class for matrices.+class Matrix a where+ -- | The number of rows in the matrix.+ numRows :: a (m,n) e -> Int+ + -- | The number of columns in the matrix.+ numCols :: a (m,n) e -> Int+ + -- | Creates a new matrix view that conjugates and transposes the + -- given matrix.+ herm :: Elem e => a (m,n) e -> a (n,m) e
+ BLAS/Matrix/Immutable.hs view
@@ -0,0 +1,34 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Matrix.Immutable+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Matrix.Immutable (+ IMatrix(..)+ ) where++import BLAS.Access+import BLAS.Elem ( BLAS3 )+import qualified BLAS.Matrix.Base as Base+import Data.Vector.Dense+import Data.Matrix.Dense.Internal+import Data.Matrix.Dense.Operations ( apply, applyMat )++infixl 7 <*>, <**>++class Base.Matrix a => IMatrix a e where+ -- | Apply to a vector+ (<*>) :: a (m,n) e -> Vector n e -> Vector m e+ + -- | Apply to a matrix+ (<**>) :: a (m,k) e -> Matrix (k,n) e -> Matrix (m,n) e++instance (BLAS3 e) => IMatrix (DMatrix Imm) e where+ (<*>) = apply+ (<**>) = applyMat+
+ BLAS/Matrix/ReadOnly.hs view
@@ -0,0 +1,30 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Matrix.ReadOnly+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Matrix.ReadOnly (+ RMatrix(..)+ ) where++import BLAS.Elem ( BLAS3 )+import qualified BLAS.Matrix.Base as Base+import Data.Vector.Dense.Internal+import Data.Matrix.Dense.Internal+import qualified Data.Matrix.Dense.Operations as M++class Base.Matrix a => RMatrix a e where+ -- | Apply to a vector+ getApply :: a (m,n) e -> DVector t n e -> IO (DVector r m e)+ + -- | Apply to a matrix+ getApplyMat :: a (m,k) e -> DMatrix t (k,n) e -> IO (DMatrix r (m,n) e)++instance (BLAS3 e) => RMatrix (DMatrix t) e where+ getApply = M.getApply+ getApplyMat = M.getApplyMat
+ BLAS/Matrix/Solve.hs view
@@ -0,0 +1,16 @@+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Matrix.Solve+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Matrix.Solve (+ module BLAS.Matrix.Solve.Immutable,+ module BLAS.Matrix.Solve.ReadOnly,+ ) where++import BLAS.Matrix.Solve.Immutable+import BLAS.Matrix.Solve.ReadOnly
+ BLAS/Matrix/Solve/Immutable.hs view
@@ -0,0 +1,28 @@+{-# LANGUAGE MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Matrix.Solve.Immutable+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Matrix.Solve.Immutable (+ ISolve(..)+ ) where++import BLAS.Matrix.Immutable++import Data.Vector.Dense ( Vector )+import Data.Matrix.Dense ( Matrix )++infixl 7 <\>, <\\>++class IMatrix a e => ISolve a e where+ -- | Solve for a vector+ (<\>) :: a (m,n) e -> Vector m e -> Vector n e+ + -- | Solve for a matrix+ (<\\>) :: a (m,n) e -> Matrix (m,k) e -> Matrix (n,k) e+
+ BLAS/Matrix/Solve/ReadOnly.hs view
@@ -0,0 +1,25 @@+{-# LANGUAGE MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Matrix.Solve.ReadOnly+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Matrix.Solve.ReadOnly (+ RSolve(..)+ ) where++import BLAS.Matrix.ReadOnly++import Data.Vector.Dense.IO ( DVector )+import Data.Matrix.Dense.IO ( DMatrix )++class RMatrix a e => RSolve a e where+ -- | Solve for a vector+ getSolve :: a (m,n) e -> DVector t m e -> IO (DVector r n e)+ + -- | Solve for a matrix+ getSolveMat :: a (m,n) e -> DMatrix t (m,k) e -> IO (DMatrix r (n,k) e)
+ BLAS/Tensor.hs view
@@ -0,0 +1,24 @@+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Tensor+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Tensor (+ module BLAS.Tensor.Base,+ module BLAS.Tensor.Dense,+ module BLAS.Tensor.Scalable,+ module BLAS.Tensor.Immutable,+ module BLAS.Tensor.ReadOnly,+ module BLAS.Tensor.Mutable+ ) where++import BLAS.Tensor.Base+import BLAS.Tensor.Dense+import BLAS.Tensor.Scalable+import BLAS.Tensor.Immutable+import BLAS.Tensor.ReadOnly+import BLAS.Tensor.Mutable
+ BLAS/Tensor/Base.hs view
@@ -0,0 +1,25 @@+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Tensor.Base+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Tensor.Base (+ Tensor(..),+ ) where++import Data.Ix++-- | The base class for tensors (i.e. Vector, Matrix, etc.).+class (Ix i) => Tensor x i e | x -> i where+ -- | Get the shape of the tensor. For vectors this is the dimension.+ -- For matrices, this will be a pair @(m,n)@ of the number of rows+ -- and columns.+ shape :: x e -> i+ + -- | Get the range of valid indices in the tensor.+ bounds :: x e -> (i,i)
+ BLAS/Tensor/Dense.hs view
@@ -0,0 +1,16 @@+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Tensor.Dense+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Tensor.Dense (+ module BLAS.Tensor.Dense.Immutable,+ module BLAS.Tensor.Dense.ReadOnly,+ ) where++import BLAS.Tensor.Dense.Immutable+import BLAS.Tensor.Dense.ReadOnly
+ BLAS/Tensor/Dense/Immutable.hs view
@@ -0,0 +1,23 @@+{-# LANGUAGE MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Tensor.Dense.Immutable+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Tensor.Dense.Immutable (+ IDTensor(..),+ ) where++import BLAS.Tensor.Immutable++-- | Class for immutable dense tensors.+class (ITensor x i e) => (IDTensor x i e) where+ -- | Get a zero tensor of the given shape.+ zero :: i -> x e+ + -- | Get a new constant tensor of the given shape.+ constant :: i -> e -> x e
+ BLAS/Tensor/Dense/ReadOnly.hs view
@@ -0,0 +1,25 @@+{-# LANGUAGE MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Tensor.Dense.ReadOnly+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Tensor.Dense.ReadOnly (+ RDTensor(..),+ ) where++import BLAS.Tensor.ReadOnly++-- | Class for mutable dense read-only tensors.+class (RTensor x i e m) => RDTensor x i e m where+ -- | Creates a new tensor with elements all initialized to zero.+ newZero :: i -> m (x e)+ + -- | Creates a new tensor with elements all initialized to the + -- given value.+ newConstant :: i -> e -> m (x e)+
+ BLAS/Tensor/Immutable.hs view
@@ -0,0 +1,70 @@+{-# LANGUAGE MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Tensor.Immutable+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Tensor.Immutable (+ ITensor(..),+ (!),+ ) where++import BLAS.Tensor.Base+import BLAS.Elem+import Data.Ix++infixl 9 !++-- | A class for immutable tensors.+class (BLAS1 e, Tensor x i e) => ITensor x i e where+ -- | Get the numer of elements stored in the tensor.+ size :: x e -> Int+ + -- | Get a new tensor by replacing the elements at the given indices.+ (//) :: x e -> [(i,e)] -> x e++ -- | Get the value at the given index, without doing any bounds-checking.+ unsafeAt :: x e -> i -> e+ + -- | Same as '(//)' but doesn't do any bounds-checking.+ unsafeReplace :: x e -> [(i,e)] -> x e+ + -- | Get the indices of the elements stored in the tensor.+ indices :: x e -> [i]+ indices = fst . unzip . assocs+ + -- | Get the elements stored in the tensor.+ elems :: x e -> [e]+ elems = snd . unzip . assocs++ -- | Get the list of @(@index@,@ element@)@ pairs stored in the tensor.+ assocs :: x e -> [(i,e)]++ -- accum :: (e -> e' -> e) -> x e -> [(i,e')] -> x e+ + -- | Apply a function elementwise to a tensor.+ amap :: (ITensor x i e') => (e -> e') -> x e -> x e'+ + -- | Apply a function to pairs of elements of tensors that are the + -- same shape.+ azipWith :: (ITensor x i f, ITensor x i g) => (e -> f -> g) -> x e -> x f -> x g+ + -- ixmap :: i -> (i -> i) -> x e -> x e+ -- unsafeIxMap++-- | Get the value at the given index. Range-checks the argument.+(!) :: (ITensor x i e, Show i) => x e -> i -> e+(!) x i =+ case (inRange b i) of+ False -> + error $ "tried to get element at a index `" ++ show i ++ "'"+ ++ " in an object with shape `" ++ show s ++ "'"+ True -> + unsafeAt x i+ where+ b = bounds x+ s = shape x
+ BLAS/Tensor/Mutable.hs view
@@ -0,0 +1,76 @@+{-# LANGUAGE MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Tensor.Mutable+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Tensor.Mutable (+ MTensor(..),+ writeElem,+ modifyElem,+ ) where++import BLAS.Tensor.Base+import BLAS.Tensor.ReadOnly++-- | Class for modifiable mutable tensors.+class (RTensor x i e m) => (MTensor x i e m) where+ -- | Get the maximum number of elements that can be stored in the tensor.+ getMaxSize :: x e -> m Int+ getMaxSize = getSize+ + -- | Sets all stored elements to zero.+ setZero :: x e -> m ()+ + -- | Sets all stored elements to the given value.+ setConstant :: e -> x e -> m ()++ -- | True if the value at a given index can be changed+ canModifyElem :: x e -> i -> m Bool+ + -- | Set the value of the element at the given index, without doing any+ -- range checking.+ unsafeWriteElem :: x e -> i -> e -> m ()+ + -- | Modify the value of the element at the given index, without doing+ -- any range checking.+ unsafeModifyElem :: x e -> i -> (e -> e) -> m ()+ unsafeModifyElem x i f = do+ e <- unsafeReadElem x i+ unsafeWriteElem x i (f e)+ + -- | Replace each element by a function applied to it+ modifyWith :: (e -> e) -> x e -> m ()+ +-- | Set the value of the element at the given index.+writeElem :: (MTensor x i e m, Show i) => x e -> i -> e -> m ()+writeElem x i e = do+ ok <- canModifyElem x i+ case ok of+ False -> + fail $ "tried to set element at index `" ++ show i ++ "'"+ ++ " in an object with shape `" ++ show s ++ "'"+ ++ " but that element cannot be modified"+ True ->+ unsafeWriteElem x i e+ where+ s = shape x++-- | Update the value of the element at the given index.+modifyElem :: (MTensor x i e m, Show i) => x e -> i -> (e -> e) -> m ()+modifyElem x i f = do+ ok <- canModifyElem x i+ case ok of+ False -> + fail $ "tried to modify element at index `" ++ show i ++ "'"+ ++ " in an object with shape `" ++ show s ++ "'"+ ++ " but that element cannot be modified"+ True ->+ unsafeModifyElem x i f+ where+ s = shape x+
+ BLAS/Tensor/ReadOnly.hs view
@@ -0,0 +1,62 @@+{-# LANGUAGE MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Tensor.ReadOnly+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Tensor.ReadOnly (+ RTensor(..),+ readElem,+ ) where++import Data.Ix+import BLAS.Tensor.Base++-- | Class for mutable read-only tensors.+class (Tensor x i e, Monad m) => RTensor x i e m where+ -- | Get the number of elements stored in the tensor.+ getSize :: x e -> m Int+ + -- | Get a copy of the tensor.+ newCopy :: x e -> m (x e)+ + -- | Get the value at the specified index, without doing any + -- range-checking.+ unsafeReadElem :: x e -> i -> m e++ -- | Returns a lazy list of the indices in the tensor. + -- Because of the laziness, this function should be used with care.+ getIndices :: x e -> m [i]+ getIndices x = getAssocs x >>= return . fst . unzip++ -- | Returns a lazy list of the elements in the tensor. + -- Because of the laziness, this function should be used with care.+ getElems :: x e -> m [e]+ getElems x = getAssocs x >>= return . snd . unzip+ + -- | Returns a lazy list of the elements-index pairs in the tensor. + -- Because of the laziness, this function should be used with care.+ getAssocs :: x e -> m [(i,e)]+ getAssocs x = do+ is <- getIndices x+ es <- getElems x+ return $ zip is es+++-- | Gets the value at the specified index after checking that the argument+-- is in bounds.+readElem :: (RTensor x i e m, Show i) => x e -> i -> m e+readElem x i =+ case (inRange b i) of+ False -> + fail $ "tried to get element at a index `" ++ show i ++ "'"+ ++ " in an object with shape `" ++ show s ++ "'"+ True -> + unsafeReadElem x i+ where+ b = bounds x+ s = shape x
+ BLAS/Tensor/Scalable.hs view
@@ -0,0 +1,20 @@+{-# LANGUAGE MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Tensor.Scalable+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Tensor.Scalable (+ Scalable(..)+ ) where++infixl 7 *>++-- | A class for scalable tensors.+class (Num e) => Scalable x e where+ -- | Scale a tensor by the given value.+ (*>) :: e -> x e -> x e
+ BLAS/Types.hs view
@@ -0,0 +1,53 @@+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Types+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Types (+ Order(..),+ Trans(..),+ UpLo(..),+ Diag(..),+ Side(..),+ + flipOrder,+ flipTrans,+ flipUpLo,+ flipSide,+ ) where++-- | Matrix element storage order+data Order = RowMajor | ColMajor deriving (Eq, Show)++-- | Transpose type+data Trans = NoTrans | ConjTrans deriving (Eq,Show)++-- | Lower or upper triangular storage+data UpLo = Upper | Lower deriving (Eq,Show)++-- | Diagonal storage+data Diag = Unit | NonUnit deriving (Eq,Show)++-- | Multiplication side+data Side = LeftSide | RightSide deriving (Eq,Show)+++flipOrder :: Order -> Order+flipOrder RowMajor = ColMajor+flipOrder ColMajor = RowMajor++flipTrans :: Trans -> Trans+flipTrans NoTrans = ConjTrans+flipTrans ConjTrans = NoTrans++flipUpLo :: UpLo -> UpLo+flipUpLo Upper = Lower+flipUpLo Lower = Upper+ +flipSide :: Side -> Side+flipSide LeftSide = RightSide+flipSide RightSide = LeftSide
+ BLAS/Vector.hs view
@@ -0,0 +1,25 @@+-----------------------------------------------------------------------------+-- |+-- Module : BLAS.Vector+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module BLAS.Vector (+ Vector(..)+ ) where++import BLAS.Elem.Base ( Elem )++-- | A class for vectors.+class Vector x where+ -- | Get the dimension of the vector.+ dim :: x n e -> Int+ + -- | Returns a view into an existing vector that takes the complex + -- conjugate of every element. For real vectors, this should compile+ -- to a no-op.+ conj :: (Elem e) => x n e -> x n e+
+ Data/Matrix/Dense.hs view
@@ -0,0 +1,168 @@+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : Data.Matrix.Dense+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Data.Matrix.Dense (+ -- * Dense matrix type+ Matrix,+ + module BLAS.Matrix.Base,+ module BLAS.Matrix.Immutable,+ module BLAS.Tensor.Base,+ module BLAS.Tensor.Dense.Immutable,+ module BLAS.Tensor.Immutable,+ module BLAS.Tensor.Scalable,++ -- * Creating matrices+ matrix, + listMatrix,+ fromCols,+ fromRows,+ + -- * Special matrices+ identity,++ -- * Rows and columns+ row,+ col,+ rows,+ cols,++ -- * Diagonals+ diag,++ -- * Augmenting matrices+ submatrix,++ -- * Matrix multiplication+ apply,+ applyMat,+ sapply,+ sapplyMat,++ -- * Matrix arithmetic+ shift,+ scale,+ invScale,++ -- * Casting matrices+ coerceMatrix,+ + -- * Converting between vectors and matrices+ fromRow,+ fromCol,+ + -- * Unsafe operations+ unsafeMatrix,+ unsafeRow,+ unsafeCol,+ unsafeDiag,+ unsafeSubmatrix,+ + ) where++import Data.Maybe ( fromJust )+import System.IO.Unsafe ( unsafePerformIO )++import BLAS.Access+import BLAS.Elem ( BLAS1, BLAS2 )+import BLAS.Matrix.Base hiding ( Matrix )+import BLAS.Matrix.Immutable+import BLAS.Tensor.Base+import BLAS.Tensor.Dense.Immutable+import BLAS.Tensor.Immutable+import BLAS.Tensor.Scalable++import Data.Matrix.Dense.Internal+import qualified Data.Matrix.Dense.Internal as M+import Data.Matrix.Dense.Operations ( apply, applyMat, sapply, sapplyMat,+ shift, scale, invScale, plus, minus, times, divide )+import Data.Vector.Dense hiding ( scale, invScale, shift )+++-- | Create a new matrix of the given size and initialize the given elements to+-- the given values. All other elements get set to zero.+matrix :: (BLAS1 e) => (Int,Int) -> [((Int,Int), e)] -> Matrix (m,n) e+matrix mn ies = unsafePerformIO $ newMatrix mn ies+{-# NOINLINE matrix #-}++-- | Same as 'matrix' but does not do any bounds checking.+unsafeMatrix :: (BLAS1 e) => (Int,Int) -> [((Int,Int), e)] -> Matrix (m,n) e+unsafeMatrix mn ies = unsafePerformIO $ unsafeNewMatrix mn ies+{-# NOINLINE unsafeMatrix #-}++-- | Create a matrix of the given shape from a list of columns+fromCols :: (BLAS1 e) => (Int,Int) -> [Vector m e] -> Matrix (m,n) e+fromCols mn cs = unsafePerformIO $ newColsMatrix mn cs+{-# NOINLINE fromCols #-}++-- | Create a matrix of the given shape from a list of rows+fromRows :: (BLAS1 e) => (Int,Int) -> [Vector n e] -> Matrix (m,n) e+fromRows mn rs = unsafePerformIO $ newRowsMatrix mn rs+{-# NOINLINE fromRows #-}++-- | Get a new matrix of the given shape with ones along the diagonal and+-- zeroes everywhere else.+identity :: (BLAS1 e) => (Int,Int) -> Matrix (m,n) e+identity mn = unsafePerformIO $ newIdentity mn+{-# NOINLINE identity #-}++-- | Get a matrix from a row vector.+fromRow :: (BLAS1 e) => Vector n e -> Matrix (one,n) e+fromRow x = + case maybeFromRow x of+ Just x' -> x'+ Nothing -> fromJust $ maybeFromRow $ unsafePerformIO $ newCopy x+{-# NOINLINE fromRow #-}++-- | Get a matrix from a column vector.+fromCol :: (BLAS1 e) => Vector m e -> Matrix (m,one) e+fromCol x = + case maybeFromCol x of+ Just x' -> x'+ Nothing -> fromJust $ maybeFromCol $ unsafePerformIO $ newCopy x+{-# NOINLINE fromCol #-}+++instance (BLAS1 e) => Scalable (DMatrix Imm (m,n)) e where+ (*>) = scale+ +instance (BLAS2 e) => Num (DMatrix Imm (m,n) e) where+ (+) = plus+ (-) = minus+ (*) = times+ negate = scale (-1)+ abs = amap abs+ signum = amap signum+ fromInteger n = constant (1,1) (fromInteger n)+ +instance (BLAS2 e) => Fractional (DMatrix Imm (m,n) e) where+ (/) a b = divide a b+ recip = amap recip+ fromRational q = constant (1,1) (fromRational q)+ +instance (BLAS2 e, Floating e) => Floating (DMatrix Imm (m,n) e) where+ pi = constant (1,1) pi+ exp = amap exp + sqrt = amap sqrt+ log = amap log+ (**) = azipWith (**)+ sin = amap sin+ cos = amap cos+ tan = amap tan+ asin = amap asin+ acos = amap acos+ atan = amap atan+ sinh = amap sinh+ cosh = amap cosh+ tanh = amap tanh+ asinh = amap asinh+ acosh = amap acosh+ atanh = amap atanh+
+ Data/Matrix/Dense/IO.hs view
@@ -0,0 +1,89 @@+-----------------------------------------------------------------------------+-- |+-- Module : Data.Matrix.Dense.IO+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--+-- This modules defines a mutable dense matrix and associate operations.+module Data.Matrix.Dense.IO (+ -- * The mutable dense matrix data type+ DMatrix(..),+ IOMatrix,+ + module BLAS.Matrix.Base,+ module BLAS.Matrix.ReadOnly,+ module BLAS.Tensor.Base,+ module BLAS.Tensor.Dense.ReadOnly,+ module BLAS.Tensor.ReadOnly,+ module BLAS.Tensor.Mutable,+ + -- * Creating new matrices+ newMatrix,+ newMatrix_,+ newListMatrix,+ newColsMatrix,+ newRowsMatrix,++ -- * Special matrices+ newIdentity,+ setIdentity,+ + -- * Views+ -- ** Rows and columns+ row,+ col,+ rows,+ cols,+ + -- ** Diagonals+ diag,+ + -- ** Matrix views+ submatrix,++ -- * Operations+ module Data.Matrix.Dense.Operations,+ + -- ** Lifting scalar and vector operations+ liftV,+ liftV2,++ -- * Converting to and from matrices+ -- ** Vectors+ maybeFromRow,+ maybeFromCol,+ maybeToVector,+ + -- ** @ForeignPtr@s+ toForeignPtr,+ fromForeignPtr,+ ldaOf,+ isHerm,+ + -- ** Coercing+ coerceMatrix,+ + -- * Unsafe operations+ unsafeNewMatrix,+ unsafeWithElemPtr,+ unsafeRow,+ unsafeCol,+ unsafeDiag,+ unsafeSubmatrix,+ unsafeFreeze,+ unsafeThaw,+ + ) where++import Data.Matrix.Dense.Internal+import Data.Matrix.Dense.Operations hiding ( gemv, gemm, apply, applyMat,+ sapply, sapplyMat, add, plus, minus, times, divide, getApply, getApplyMat )+ +import BLAS.Matrix.Base hiding ( Matrix )+import BLAS.Matrix.ReadOnly+import BLAS.Tensor.Base+import BLAS.Tensor.Dense.ReadOnly+import BLAS.Tensor.ReadOnly+import BLAS.Tensor.Mutable
+ Data/Matrix/Dense/Internal.hs view
@@ -0,0 +1,536 @@+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : Data.Matrix.Dense.Internal+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Data.Matrix.Dense.Internal (+ -- * Dense matrix data types+ DMatrix(..),+ IOMatrix,+ Matrix,++ module BLAS.Matrix.Base,+ module BLAS.Tensor,++ -- * Converting to and from foreign pointers+ toForeignPtr,+ fromForeignPtr,+ ldaOf,+ isHerm,++ -- * Creating new matrices+ newMatrix,+ newMatrix_,+ newListMatrix,+ newColsMatrix,+ newRowsMatrix,+ + listMatrix,++ -- * Special matrices+ newIdentity,+ setIdentity,+ + -- * Row and column views+ row,+ col,+ rows,+ cols,+ + -- * Diagonal views+ diag,+ + -- * Matrix views+ submatrix,++ -- * Converting to/from vectors+ maybeFromRow,+ maybeFromCol,+ maybeToVector,+ + -- * Lifting scalar and vector operations+ liftV,+ liftV2,+ + -- * Casting matrices+ coerceMatrix,+ + -- * Unsafe operations+ unsafeThaw,+ unsafeFreeze,+ unsafeNewMatrix,+ unsafeWithElemPtr,+ unsafeRow,+ unsafeCol,+ unsafeDiag,+ unsafeSubmatrix,+ ) where++import Control.Monad ( forM_, zipWithM_ )+import Data.Ix ( inRange, range )+import Foreign+import System.IO.Unsafe +import Unsafe.Coerce++import Data.AEq++import Data.Vector.Dense.Internal hiding ( toForeignPtr, fromForeignPtr,+ unsafeFreeze, unsafeThaw, fptr, offset, unsafeWithElemPtr )+import qualified Data.Vector.Dense.Internal as V+import qualified Data.Vector.Dense.Operations as V++import BLAS.Access+import BLAS.Internal ( inlinePerformIO, checkedRow, checkedCol, checkedDiag,+ checkedSubmatrix, diagStart, diagLen )+import BLAS.Elem ( Elem, BLAS1 )+import qualified BLAS.Elem as E+import BLAS.Matrix.Base hiding ( Matrix )+import qualified BLAS.Matrix.Base as C+import BLAS.Tensor+import BLAS.Types+++-- | The mutable dense matrix data type. It can either store elements in +-- column-major order, or provide a view into another matrix. The view +-- transposes and conjugates the underlying matrix.+data DMatrix t mn e =+ DM { fptr :: !(ForeignPtr e) -- ^ a pointer to the storage region+ , offset :: !Int -- ^ an offset (in elements, not bytes) to the first element in the matrix. + , size1 :: !Int -- ^ the number of rows in the matrix+ , size2 :: !Int -- ^ the number of columns in the matrix+ , lda :: !Int -- ^ the leading dimension size of the matrix+ }+ | H !(DMatrix t mn e) -- ^ a transposed and conjugated matrix++type Matrix = DMatrix Imm+type IOMatrix = DMatrix Mut++unsafeFreeze :: DMatrix t mn e -> Matrix mn e+unsafeFreeze = unsafeCoerce++unsafeThaw :: DMatrix t mn e -> IOMatrix mn e+unsafeThaw = unsafeCoerce++-- | Coerce the phantom shape type from one type to another.+coerceMatrix :: DMatrix t mn e -> DMatrix t kl e+coerceMatrix = unsafeCoerce++-- | @fromForeignPtr f o mn l@ creates a matrix view of the data pointed to+-- by @f@ starting at offset @o@ and having shape @mn@ and lda @l@.+fromForeignPtr :: ForeignPtr e -> Int -> (Int,Int) -> Int -> DMatrix t (m,n) e+fromForeignPtr f o (m,n) l = DM f o m n l++-- | Convert a dense matrix to a pointer, offset, size, and lda. Note that this+-- does not give the conjugacy/transpose information. For that, use 'isHerm'.+toForeignPtr :: DMatrix t (m,n) e -> (ForeignPtr e, Int, (Int,Int), Int)+toForeignPtr (H a) = toForeignPtr a+toForeignPtr a@(DM _ _ _ _ _) = (fptr a, offset a, (size1 a, size2 a), lda a)++-- | Get the lda of a matrix, defined as the number of elements in the underlying+-- array that separate two consecutive elements in the same row of the matrix.+ldaOf :: DMatrix t (m,n) e -> Int+ldaOf (H a) = ldaOf a+ldaOf a@(DM _ _ _ _ _) = lda a+++indexOf :: DMatrix t (m,n) e -> (Int,Int) -> Int+indexOf (H a) (i,j) = indexOf a (j,i)+indexOf a@(DM _ _ _ _ _) (i,j) = + let o = offset a+ l = lda a+ in o + i + j*l+ +-- | Get the storage order of the matrix. If 'isTrans' is true, this +-- will be 'RowMajor'. Otherwise, it will be 'ColMajor'.+orderOf :: DMatrix t (m,n) e -> Order+orderOf (H a) = flipOrder (orderOf a)+orderOf (DM _ _ _ _ _) = ColMajor++-- | Get whether or not the matrix is transposed and conjugated.+isHerm :: DMatrix t (m,n) e -> Bool+isHerm (H a) = not (isHerm a)+isHerm (DM _ _ _ _ _) = False++-- | Create a new matrix of the given size and initialize the given elements to+-- the given values. All other elements get initialized to zero.+newMatrix :: (BLAS1 e) => (Int,Int) -> [((Int,Int), e)] -> IO (DMatrix t (m,n) e)+newMatrix = newMatrixHelp writeElem++-- | Same as 'newMatrix' but do not do any bounds-checking.+unsafeNewMatrix :: (BLAS1 e) => (Int,Int) -> [((Int,Int), e)] -> IO (DMatrix t (m,n) e)+unsafeNewMatrix = newMatrixHelp unsafeWriteElem++newMatrixHelp :: (BLAS1 e) => + (IOMatrix (m,n) e -> (Int,Int) -> e -> IO ()) + -> (Int,Int) -> [((Int,Int),e)] -> IO (DMatrix t (m,n) e)+newMatrixHelp set mn ijes = do+ x <- newZero mn+ io <- unsafeInterleaveIO $ mapM_ (uncurry $ set $ unsafeThaw x) ijes+ return $ io `seq` x++-- | Create a new matrix of given shape, but do not initialize the elements.+newMatrix_ :: (Elem e) => (Int,Int) -> IO (DMatrix t (m,n) e)+newMatrix_ (m,n) + | m < 0 || n < 0 =+ ioError $ userError $ + "Tried to create a matrix with shape `" ++ show (m,n) ++ "'"+ | otherwise = do+ f <- mallocForeignPtrArray (m*n)+ return $ fromForeignPtr f 0 (m,n) (max 1 m)++-- | Create a new matrix with the given elements in column-major order.+newListMatrix :: (Elem e) => (Int,Int) -> [e] -> IO (DMatrix t (m,n) e)+newListMatrix (m,n) es = do+ a <- newMatrix_ (m,n)+ unsafeWithElemPtr a (0,0) $ flip pokeArray (take (m*n) es)+ return a++-- | Create a new matrix with the given elements in row-major order.+listMatrix :: (Elem e) => (Int,Int) -> [e] -> Matrix (m,n) e+listMatrix mn es = unsafePerformIO $ newListMatrix mn es+{-# NOINLINE listMatrix #-}++-- | Create a new matrix of the given shape with ones along the diagonal, +-- and zeros everywhere else.+newIdentity :: (BLAS1 e) => (Int,Int) -> IO (DMatrix t (m,n) e)+newIdentity mn = do+ a <- newMatrix_ mn+ setIdentity (unsafeThaw a)+ return a++-- | Set the diagonal to ones, and set everywhere else to zero.+setIdentity :: (BLAS1 e) => IOMatrix (m,n) e -> IO ()+setIdentity a = do+ s <- getSize a+ case s of+ 0 -> return ()+ _ -> setZero a >>+ setConstant 1 (diag a 0)++-- | Form a matrix from a list of column vectors.+newColsMatrix :: (BLAS1 e) => (Int,Int) -> [DVector t m e] -> IO (DMatrix r (m,n) e)+newColsMatrix (m,n) cs = do+ a <- newZero (m,n)+ forM_ (zip [0..(n-1)] cs) $ \(j,c) ->+ V.copyVector (unsafeCol (unsafeThaw a) j) c+ return a++-- | Form a matrix from a list of row vectors.+newRowsMatrix :: (BLAS1 e) => (Int,Int) -> [DVector t n e] -> IO (DMatrix r (m,n) e)+newRowsMatrix (m,n) rs = do+ a <- newZero (m,n)+ forM_ (zip [0..(m-1)] rs) $ \(i,r) ->+ V.copyVector (unsafeRow (unsafeThaw a) i) r+ return a++-- | Evaluate a function with a pointer to the raw storage for the element+-- at the given index. It may be necessary to conjugate or scale values before+-- reading or writing to or from the location.+unsafeWithElemPtr :: (Elem e) => DMatrix t (m,n) e -> (Int,Int) -> (Ptr e -> IO a) -> IO a+unsafeWithElemPtr (H a) (i,j) f = unsafeWithElemPtr a (j,i) f+unsafeWithElemPtr a@(DM _ _ _ _ _) ij f =+ withForeignPtr (fptr a) $ \ptr ->+ let ptr' = ptr `advancePtr` (indexOf a ij)+ in f ptr'+ +-- | Get a vector view of the given row in a matrix.+row :: (Elem e) => DMatrix t (m,n) e -> Int -> DVector t n e+row a = checkedRow (shape a) (unsafeRow a)++-- | Get a list of vector views of the rows of the matrix.+rows :: (Elem e) => DMatrix t (m,n) e -> [DVector t n e]+rows a = [ unsafeRow a i | i <- [0..numRows a - 1] ]++-- | Get a list of vector views of the columns of the matrix.+cols :: (Elem e) => DMatrix t (m,n) e -> [DVector t m e]+cols a = [ unsafeCol a j | j <- [0..numCols a - 1] ]++-- | Get a vector view of the given column in a matrix.+col :: (Elem e) => DMatrix t (m,n) e -> Int -> DVector t m e+col a = checkedCol (shape a) (unsafeCol a)++-- | Same as 'row', but does not do any bounds checking.+unsafeRow :: (Elem e) => DMatrix t (m,n) e -> Int -> DVector t n e+unsafeRow a@(H _) i = conj $ unsafeCol (herm a) i+unsafeRow a@(DM _ _ _ _ _) i =+ let f = fptr a+ o = indexOf a (i,0)+ n = size2 a+ s = lda a+ in V.fromForeignPtr f o n s++-- | Same as 'col', but does not do any bounds checking.+unsafeCol :: (Elem e) => DMatrix t (m,n) e -> Int -> DVector t m e+unsafeCol a@(H _) j = conj $ unsafeRow (herm a) j+unsafeCol a@(DM _ _ _ _ _) j =+ let f = fptr a+ o = indexOf a (0,j)+ m = size1 a+ s = 1+ in V.fromForeignPtr f o m s++-- | @diag a 0@ gets a vector view of the main diagonal of @a@. @diag a k@ for +-- @k@ positive gets a view of the @k@th superdiagonal. For @k@ negative, it+-- gets a view of the @(-k)@th subdiagonal.+diag :: (Elem e) => DMatrix t (m,n) e -> Int -> DVector t k e+diag a = checkedDiag (shape a) (unsafeDiag a)++-- | Same as 'diag', but does not do any bounds checking.+unsafeDiag :: (Elem e) => DMatrix t (m,n) e -> Int -> DVector t k e+unsafeDiag (H a) i = conj $ unsafeDiag a (negate i)+unsafeDiag a@(DM _ _ _ _ _) i =+ let f = fptr a+ o = indexOf a (diagStart i)+ n = diagLen (shape a) i+ s = lda a + 1+ in V.fromForeignPtr f o n s++-- | @submatrix a ij mn@ returns a view of the submatrix of @a@ with element @(0,0)@+-- being element @ij@ in @a@, and having shape @mn@.+submatrix :: (Elem e) => DMatrix t (m,n) e -> (Int,Int) -> (Int,Int) -> DMatrix t (k,l) e+submatrix a = checkedSubmatrix (shape a) (unsafeSubmatrix a)++-- | Same as 'submatrix' but does not do any bounds checking.+unsafeSubmatrix :: (Elem e) => DMatrix t (m,n) e -> (Int,Int) -> (Int,Int) -> DMatrix t (k,l) e+unsafeSubmatrix a@(H _) (i,j) (m',n') = herm $ unsafeSubmatrix (herm a) (j,i) (n',m')+unsafeSubmatrix a@(DM _ _ _ _ _) (i,j) mn' =+ let f = fptr a+ o = indexOf a (i,j)+ l = lda a+ in fromForeignPtr f o mn' l+ ++-- | Create a matrix view of a row vector. This will fail if the+-- stride is not @1@ and the vector is conjugated.+maybeFromRow :: (Elem e) => DVector t m e -> Maybe (DMatrix t (one,m) e)+maybeFromRow (V.C (V.C x)) = maybeFromRow x+maybeFromRow (V.C x@(V.DV _ _ _ _))+ | V.stride x == 1 =+ let f = V.fptr x+ o = V.offset x+ n = V.dim x+ l = max 1 n+ in Just $ herm $ fromForeignPtr f o (n,1) l+ | otherwise =+ Nothing+maybeFromRow x@(V.DV _ _ _ _) =+ let f = V.fptr x+ o = V.offset x+ n = V.dim x+ s = V.stride x+ l = max 1 s+ in Just $ fromForeignPtr f o (1,n) l+++-- | Possibly create a matrix view of a column vector. This will fail+-- if the stride of the vector is not @1@ and the vector is not conjugated.+maybeFromCol :: (Elem e) => DVector t n e -> Maybe (DMatrix t (n,one) e)+maybeFromCol (V.C x) = maybeFromRow x >>= return . herm+maybeFromCol x@(V.DV _ _ _ _)+ | V.stride x == 1 =+ let f = V.fptr x+ o = V.offset x+ m = dim x+ l = max 1 m+ in Just $ fromForeignPtr f o (m,1) l+ | otherwise =+ Nothing++maybeToVector :: (Elem e) => DMatrix t (m,n) e -> Maybe (Order, DVector t k e)+maybeToVector (H a) = maybeToVector a >>= (\(o,x) -> return (flipOrder o, conj x))+maybeToVector (DM f o m n ld)+ | ld == m =+ Just $ (ColMajor, V.fromForeignPtr f o (m*n) 1)+ | m == 1 =+ Just $ (ColMajor, V.fromForeignPtr f o n ld)+ | otherwise =+ Nothing++-- | Modify each element in-place by applying a function to it.+-- modifyWith :: (Elem e) => (e -> e) -> IOMatrix (m,n) e -> IO ()+++-- | Take a unary elementwise vector operation and apply it to the elements of a matrix.+liftV :: (Elem e) => (DVector t k e -> IO ()) -> DMatrix t (m,n) e -> IO ()+liftV f a =+ case maybeToVector a of+ Just (_,x) -> f x+ _ -> + let xs = case orderOf a of+ RowMajor -> rows (coerceMatrix a)+ ColMajor -> cols (coerceMatrix a)+ in mapM_ f xs++-- | Take a binary elementwise vector operation and apply it to the elements of a pair+-- of matrices.+liftV2 :: (Elem e) => (DVector s k e -> DVector t k e -> IO ()) + -> DMatrix s (m,n) e -> DMatrix t (m,n) e -> IO ()+liftV2 f a b =+ case (maybeToVector a, maybeToVector b) of+ (Just (RowMajor,x), Just (RowMajor,y)) -> f x y+ (Just (ColMajor,x), Just (ColMajor,y)) -> f x y+ _ -> + let (xs,ys) = case orderOf a of+ RowMajor -> (rows (coerceMatrix a), rows (coerceMatrix b))+ ColMajor -> (cols (coerceMatrix a), cols (coerceMatrix b))+ in zipWithM_ f xs ys+++instance C.Matrix (DMatrix t) where+ numRows = fst . shape+ numCols = snd . shape++ herm a = case a of+ (H a') -> coerceMatrix a'+ _ -> H (coerceMatrix a)+ +instance Tensor (DMatrix 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 (DMatrix Imm (m,n)) (Int,Int) e where+ size a = (numRows a * numCols a)+ + unsafeAt a = inlinePerformIO . unsafeReadElem a+ {-# INLINE unsafeAt #-}+ + indices a = [ (i,j) | j <- range (0,n-1), i <- range (0,m-1) ]+ where (m,n) = shape a+ {-# INLINE indices #-}+ + elems = inlinePerformIO . getElems+ {-# INLINE elems #-}+ + assocs = inlinePerformIO . getAssocs+ {-# INLINE assocs #-}++ (//) = replaceHelp writeElem+ unsafeReplace = replaceHelp unsafeWriteElem++ amap f a = listMatrix (shape a) (map f $ elems a)++ azipWith f a b+ | shape b /= mn =+ error ("azipWith: matrix shapes differ; first has shape `"+ ++ show mn ++ "' and second has shape `" + ++ show (shape b) ++ "'")+ | otherwise =+ listMatrix mn (zipWith f (elems a) (elems b))+ where+ mn = shape a+++replaceHelp :: (BLAS1 e) => + (IOMatrix (m,n) e -> (Int,Int) -> e -> IO ())+ -> Matrix (m,n) e -> [((Int,Int), e)] -> Matrix (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) => IDTensor (DMatrix Imm (m,n)) (Int,Int) e where+ zero = unsafePerformIO . newZero+ {-# NOINLINE zero #-}+ + constant mn = unsafePerformIO . newConstant mn+ {-# NOINLINE constant #-}+ + +instance (BLAS1 e) => RTensor (DMatrix t (m,n)) (Int,Int) e IO where+ getSize a = return (numRows a * numCols a)+ + newCopy a = case a of+ (H a') -> newCopy a' >>= return . H+ _ -> do+ a' <- newMatrix_ (shape a)+ liftV2 V.copyVector (unsafeThaw a') a+ return a'+ + 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))+ {-# INLINE unsafeReadElem #-}+ + getIndices = return . indices . unsafeFreeze+ {-# INLINE getIndices #-}++ getElems a = return $ go (cols a)+ where go cs | cs `seq` False = undefined+ go [] = []+ go (c:cs) =+ let e = inlinePerformIO $ getElems c+ es = go cs+ in e ++ es+ {-# NOINLINE getElems #-}++ getAssocs a = return $ go (cols a) 0+ where go cs j | cs `seq` j `seq` False = undefined+ go [] _ = []+ go (c:cs) j =+ let ie = inlinePerformIO $ getAssocs c+ ije = map (\(i,e) -> ((i,j),e)) ie+ ijes = go cs (j+1)+ in ije ++ ijes+ {-# NOINLINE getAssocs #-}+++instance (BLAS1 e) => RDTensor (DMatrix t (m,n)) (Int,Int) e IO where+ newZero mn = do+ a <- newMatrix_ mn+ setZero (unsafeThaw a)+ return a+ + newConstant mn e = do+ a <- newMatrix_ mn+ setConstant e (unsafeThaw a)+ return a+ + +instance (BLAS1 e) => MTensor (DMatrix Mut (m,n)) (Int,Int) e IO where + setZero = liftV setZero+ setConstant e = liftV (setConstant e)+ + canModifyElem a ij =+ return $ inRange (bounds a) ij+ {-# INLINE canModifyElem #-}++ 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 = liftV (modifyWith f) +++instance (BLAS1 e, Show e) => Show (DMatrix Imm (m,n) e) where+ show a = case a of+ (H a') -> "herm (" ++ show a' ++ ")"+ _ -> "listMatrix " ++ show (shape a) ++ " " ++ show (elems a)+ +compareHelp :: (BLAS1 e) => + (e -> e -> Bool) -> Matrix (m,n) e -> Matrix (m,n) e -> Bool+compareHelp cmp x y+ | isHerm x && isHerm y =+ compareHelp cmp (herm x) (herm y)+compareHelp cmp x y =+ (shape x == shape y) && (and $ zipWith cmp (elems x) (elems y))++instance (BLAS1 e, Eq e) => Eq (DMatrix Imm (m,n) e) where+ (==) = compareHelp (==)++instance (BLAS1 e, AEq e) => AEq (DMatrix Imm (m,n) e) where+ (===) = compareHelp (===)+ (~==) = compareHelp (~==)+
+ Data/Matrix/Dense/Operations.hs view
@@ -0,0 +1,347 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Data.Matrix.Dense.Operations+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Data.Matrix.Dense.Operations (+ -- * Copy and swap+ copyMatrix,+ swapMatrices,++ -- * Matrix multiplication+ -- ** Pure+ apply,+ applyMat,+ sapply,+ sapplyMat,+ + -- ** Impure+ getApply,+ getApplyMat,+ getSApply,+ getSApplyMat,++ -- * Matrix Arithmetic+ -- ** Pure+ shift,+ scale,+ invScale,+ add,+ plus,+ minus,+ times,+ divide,+ + -- ** Impure+ getShifted,+ getScaled,+ getInvScaled,+ getSum,+ getDiff,+ getProduct,+ getRatio,+ + -- * In-place operations+ doConj,+ shiftBy,+ scaleBy,+ invScaleBy,+ axpy,+ (+=),+ (-=),+ (*=),+ (//=),+ + -- * BLAS operations+ gemv,+ gemm,+ + ) where++import System.IO.Unsafe+import Unsafe.Coerce++import BLAS.Internal ( checkMatMatOp, checkMatVecMult, checkMatMatMult )+import Data.Matrix.Dense.Internal+import Data.Vector.Dense.Internal hiding ( unsafeWithElemPtr, unsafeThaw, + unsafeFreeze )+import qualified Data.Vector.Dense.Operations as V+import qualified Data.Vector.Dense.Internal as V++import BLAS.C ( CBLASTrans, colMajor, noTrans, conjTrans )+import qualified BLAS.C as BLAS+import BLAS.Elem ( BLAS1, BLAS2, BLAS3 )+import qualified BLAS.Elem as E++infixl 7 `apply`, `applyMat`, `scale`, `invScale`+infixl 6 `shift`+infixl 1 +=, -=, *=, //=++-- | @copy dst src@ copies the elements from the second argument to the first.+copyMatrix :: (BLAS1 e) => IOMatrix (m,n) e -> DMatrix t (m,n) e -> IO ()+copyMatrix a b = checkMatMatOp "copyMatrix" (shape a) (shape b) >> unsafeCopyMatrix a b++unsafeCopyMatrix :: (BLAS1 e) => IOMatrix (m,n) e -> DMatrix t (m,n) e -> IO ()+unsafeCopyMatrix = liftV2 (V.unsafeCopyVector)++-- | @swap a b@ exchanges the elements stored in two matrices.+swapMatrices :: (BLAS1 e) => IOMatrix (m,n) e -> IOMatrix (m,n) e -> IO ()+swapMatrices a b = checkMatMatOp "swapMatrices" (shape a) (shape b) >> unsafeSwapMatrices a b++unsafeSwapMatrices :: (BLAS1 e) => IOMatrix (m,n) e -> IOMatrix (m,n) e -> IO ()+unsafeSwapMatrices = liftV2 (V.unsafeSwapVectors)++-- | Multiply a matrix by a vector.+getApply :: (BLAS2 e) => DMatrix s (m,n) e -> DVector t n e -> IO (DVector r m e)+getApply = getSApply 1++-- | Multiply a scaled matrix by a vector.+getSApply :: (BLAS2 e) => e -> DMatrix s (m,n) e -> DVector t n e -> IO (DVector r m e)+getSApply alpha a x = checkMatVecMult (shape a) (V.dim x) >> unsafeGetSApply alpha a x++unsafeGetSApply :: (BLAS2 e) => e -> DMatrix s (m,n) e -> DVector t n e -> IO (DVector r m e)+unsafeGetSApply alpha a x = do+ y <- V.newZero (numRows a)+ gemv alpha a x 0 y+ return (unsafeCoerce y)++-- | Multiply a matrix by a matrix.+getApplyMat :: (BLAS3 e) => DMatrix s (m,k) e -> DMatrix t (k,n) e -> IO (DMatrix r (m,n) e)+getApplyMat = getSApplyMat 1++-- | Multiply a scaled matrix by a matrix.+getSApplyMat :: (BLAS3 e) => e -> DMatrix s (m,k) e -> DMatrix t (k,n) e -> IO (DMatrix r (m,n) e)+getSApplyMat alpha a b = checkMatMatMult (shape a) (shape b) >> unsafeGetSApplyMat alpha a b++unsafeGetSApplyMat :: (BLAS3 e) => e -> DMatrix s (m,k) e -> DMatrix t (k,n) e -> IO (DMatrix r (m,n) e)+unsafeGetSApplyMat alpha a b = do+ c <- newZero (numRows a, numCols b)+ gemm alpha a b 0 c+ return (unsafeCoerce c)++-- | Form a new matrix by adding a value to every element in a matrix.+getShifted :: (BLAS1 e) => e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)+getShifted k = unaryOp (shiftBy k)++-- | Form a new matrix by multiplying every element by a value.+getScaled :: (BLAS1 e) => e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)+getScaled k = unaryOp (scaleBy k)++-- | Form a new matrix by dividing every element by a value.+getInvScaled :: (BLAS1 e) => e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)+getInvScaled k = unaryOp (invScaleBy k)++-- | Create a new matrix by taking the elementwise sum of two matrices.+getSum :: (BLAS1 e) => e -> DMatrix s (m,n) e -> e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)+getSum alpha a beta b = checkMatMatOp "getSum" (shape a) (shape b) + >> unsafeGetSum alpha a beta b++unsafeGetSum :: (BLAS1 e) => e -> DMatrix s (m,n) e -> e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)+unsafeGetSum alpha a@(H _) beta b = do+ s <- unsafeGetSum (E.conj alpha) (herm a) (E.conj beta) (herm b)+ return (herm s)+unsafeGetSum alpha a@(DM _ _ _ _ _) beta b = do+ s <- getScaled alpha a+ axpy beta b (unsafeThaw s)+ return (unsafeCoerce s)++-- | Create a new matrix by taking the elementwise difference of two matrices.+getDiff :: (BLAS1 e) => DMatrix s (m,n) e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)+getDiff a b = checkMatMatOp "minus" (shape a) (shape b) >> unsafeGetSum 1 a (-1) b++-- | Create a new matrix by taking the elementwise product of two matrices.+getProduct :: (BLAS2 e) => DMatrix s (m,n) e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)+getProduct = binaryOp "times" (*=)++-- | Create a new matrix by taking the elementwise ratio of two matrices.+getRatio :: (BLAS2 e) => DMatrix s (m,n) e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)+getRatio = binaryOp "getRatio" (//=)++-- | Conjugate every element in a matrix.+doConj :: (BLAS1 e) => IOMatrix (m,n) e -> IO (IOMatrix (m,n) e)+doConj (H a) = do+ a' <- doConj a+ return (H a')+doConj a@(DM _ _ _ _ _) = do+ liftV (\x -> V.doConj x >> return ()) a+ return a++-- | Scale every element in a matrix by the given value.+scaleBy :: (BLAS1 e) => e -> IOMatrix (m,n) e -> IO ()+scaleBy k = liftV (\x -> V.scaleBy k x >> return ())+ +-- | Scale every element by the given value.+shiftBy :: (BLAS1 e) => e -> IOMatrix (m,n) e -> IO ()+shiftBy k = liftV (V.shiftBy k)++-- | Divide every element by the given value.+invScaleBy :: (BLAS1 e) => e -> IOMatrix (m,n) e -> IO ()+invScaleBy k = liftV (V.invScaleBy k)++axpy :: (BLAS1 e) => e -> DMatrix t (m,n) e -> IOMatrix (m,n) e -> IO ()+axpy alpha = liftV2 (V.axpy alpha)++-- | In-place elementwise add.+(+=) :: (BLAS1 e) => IOMatrix (m,n) e -> DMatrix t (m,n) e -> IO ()+(+=) = liftV2 (V.+=)++-- | In-place elementwise subtract.+(-=) :: (BLAS1 e) => IOMatrix (m,n) e -> DMatrix t (m,n) e -> IO ()+(-=) = liftV2 (V.-=)++-- | In-place elementwise product.+(*=) :: (BLAS2 e) => IOMatrix (m,n) e -> DMatrix t (m,n) e -> IO ()+(*=) = liftV2 (V.*=)++-- | In-place elementwise divide.+(//=) :: (BLAS2 e) => IOMatrix (m,n) e -> DMatrix t (m,n) e -> IO ()+(//=) = liftV2 (V.//=)++blasTransOf :: DMatrix t (m,n) e -> CBLASTrans+blasTransOf a = + case (isHerm a) of+ False -> noTrans+ True -> conjTrans++flipShape :: (Int,Int) -> (Int,Int)+flipShape (m,n) = (n,m)++-- | @gemv alpha a x beta y@ replaces @y := alpha a * x + beta y@+gemv :: (BLAS2 e) => e -> DMatrix s (m,n) e -> DVector t n e -> e -> IOVector m e -> IO ()+gemv alpha a x beta y+ | numRows a == 0 || numCols a == 0 =+ return ()+ | V.isConj y = do+ V.doConj y+ gemv alpha a x beta (V.conj y)+ V.doConj y+ | V.isConj x = do+ x' <- V.newCopy (V.unsafeThaw x)+ V.doConj x'+ gemv alpha a (V.conj x') beta y+ | otherwise =+ let order = colMajor+ transA = blasTransOf a+ (m,n) = case (isHerm a) of+ False -> shape a+ True -> (flipShape . shape) a+ ldA = ldaOf a+ incX = V.strideOf x+ incY = V.strideOf y+ in unsafeWithElemPtr a (0,0) $ \pA ->+ V.unsafeWithElemPtr x 0 $ \pX ->+ V.unsafeWithElemPtr y 0 $ \pY -> do+ BLAS.gemv order transA m n alpha pA ldA pX incX beta pY incY+++-- | @gemm alpha a b beta c@ replaces @c := alpha a * b + beta c@.+gemm :: (BLAS3 e) => e -> DMatrix s (m,k) e -> DMatrix t (k,n) e -> e -> IOMatrix (m,n) e -> IO ()+gemm alpha a b beta c+ | numRows a == 0 || numCols a == 0 || numCols b == 0 = return ()+ | isHerm c = gemm (E.conj alpha) (herm b) (herm a) (E.conj beta) (herm c)+ | otherwise =+ let order = colMajor+ transA = blasTransOf a+ transB = blasTransOf b+ (m,n) = shape c+ k = numCols a+ ldA = ldaOf a+ ldB = ldaOf b+ ldC = ldaOf c+ in unsafeWithElemPtr a (0,0) $ \pA ->+ unsafeWithElemPtr b (0,0) $ \pB ->+ unsafeWithElemPtr c (0,0) $ \pC ->+ BLAS.gemm order transA transB m n k alpha pA ldA pB ldB beta pC ldC+++unaryOp :: (BLAS1 e) => (IOMatrix (m,n) e -> IO ()) + -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)+unaryOp f a = do+ a' <- newCopy a+ f (unsafeThaw a')+ return (unsafeCoerce a')++binaryOp :: (BLAS1 e) => String -> (IOMatrix (m,n) e -> DMatrix t (m,n) e -> IO ()) + -> DMatrix s (m,n) e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)+binaryOp name f a b = + checkMatMatOp name (shape a) (shape b) >> do+ a' <- newCopy a+ f (unsafeThaw a') b+ return (unsafeCoerce a')+ + +-- | Multiply a matrix by a vector.+apply :: (BLAS2 e) => Matrix (m,n) e -> Vector n e -> Vector m e+apply = sapply 1++-- | Multiply a scaled matrix by a vector.+sapply :: (BLAS2 e) => e -> Matrix (m,n) e -> Vector n e -> Vector m e+sapply alpha a x = unsafePerformIO $ getSApply alpha a x+{-# NOINLINE sapply #-}++-- | Multiply a scaled matrix by a matrix.+sapplyMat :: (BLAS3 e) => e -> Matrix (m,k) e -> Matrix (k,n) e -> Matrix (m,n) e+sapplyMat alpha a b = unsafePerformIO $ getSApplyMat alpha a b+{-# NOINLINE sapplyMat #-}+ +-- | Multiply a matrix by a matrix.+applyMat :: (BLAS3 e) => Matrix (m,k) e -> Matrix (k,n) e -> Matrix (m,n) e+applyMat = sapplyMat 1++-- | Create a new matrix by scaling another matrix by the given value.+scale :: (BLAS1 e) => e -> Matrix (m,n) e -> Matrix (m,n) e+scale k a = unsafePerformIO $ getScaled k a+{-# NOINLINE scale #-}++-- | Form a new matrix by adding a value to every element in a matrix.+shift :: (BLAS1 e) => e -> Matrix (m,n) e -> Matrix (m,n) e+shift k a = unsafePerformIO $ getShifted k a+{-# NOINLINE shift #-}++-- | Form a new matrix by dividing every element by a value.+invScale :: (BLAS1 e) => e -> Matrix (m,n) e -> Matrix (m,n) e+invScale k a = unsafePerformIO $ getInvScaled k a+{-# NOINLINE invScale #-}++-- | Create a new matrix by taking the elementwise sum of two matrices.+add :: (BLAS1 e) => e -> Matrix (m,n) e -> e -> Matrix (m,n) e -> Matrix (m,n) e+add alpha a beta b = unsafePerformIO $ getSum alpha a beta b+{-# NOINLINE add #-}++-- | Create a new matrix by taking the elementwise sum of two matrices.+plus :: (BLAS1 e) => Matrix (m,n) e -> Matrix (m,n) e -> Matrix (m,n) e+plus a b = add 1 a 1 b++-- | Create a new matrix by taking the elementwise difference of two matrices.+minus :: (BLAS1 e) => Matrix (m,n) e -> Matrix (m,n) e -> Matrix (m,n) e+minus a b = unsafePerformIO $ getDiff a b+{-# NOINLINE minus #-}++-- | Create a new matrix by taking the elementwise product of two matrices.+times :: (BLAS2 e) => Matrix (m,n) e -> Matrix (m,n) e -> Matrix (m,n) e+times a b = unsafePerformIO $ getProduct a b+{-# NOINLINE times #-}++-- | Create a new matrix by taking the elementwise ratio of two matrices.+divide :: (BLAS2 e) => Matrix (m,n) e -> Matrix (m,n) e -> Matrix (m,n) e+divide a b = unsafePerformIO $ getRatio a b+{-# NOINLINE divide #-} ++{-# RULES+"scale.plus/add" forall k l x y. plus (scale k x) (scale l y) = add k x l y+"scale1.plus/add" forall k x y. plus (scale k x) y = add k x 1 y+"scale2.plus/add" forall k x y. plus x (scale k y) = add 1 x k y++"scale.minus/add" forall k l x y. minus (scale k x) (scale l y) = add k x (-l) y+"scale1.minus/add" forall k x y. minus (scale k x) y = add k x (-1) y+"scale2.minus/add" forall k x y. minus x (scale k y) = add 1 x (-k) y++"scale.apply/sapply" forall k a x. apply (scale k a) x = sapply k a x+"scale.applyMat/sapplyMat" forall k a b. applyMat (scale k a) b = sapplyMat k a b+ #-}+
+ Data/Matrix/Herm.hs view
@@ -0,0 +1,64 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : Data.Matrix.Herm+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Data.Matrix.Herm (+ Herm(..),+ UpLo(..),++ fromBase,+ toBase,+ mapHerm,++ hermL,+ hermU,++ ) where++import Unsafe.Coerce++import qualified BLAS.Elem as E+import BLAS.Matrix+import BLAS.Tensor+import BLAS.Types ( UpLo(..) )++data Herm a nn e = Herm UpLo e (a nn e)++mapHerm :: (a (n,n) e -> b (n,n) e) -> Herm a (n,n) e -> Herm b (n,n) e+mapHerm f (Herm u e a) = Herm u e $ f a++fromBase :: UpLo -> e -> a (n,n) e -> Herm a (n,n) e+fromBase = Herm+ +toBase :: Herm a (n,n) e -> (UpLo, e, a (n,n) e)+toBase (Herm u e a) = (u,e,a)++hermL :: (Num e) => a (n,n) e -> Herm a (n,n) e+hermL = Herm Lower 1++hermU :: (Num e) => a (n,n) e -> Herm a (n,n) e+hermU = Herm Upper 1+ +instance Matrix a => Matrix (Herm a) where+ numRows (Herm _ _ a) = numRows a+ numCols (Herm _ _ a) = numCols a+ herm (Herm u e a) = Herm u (E.conj e) (unsafeCoerce a)+ +instance (Num e) => Scalable (Herm a nn) e where+ (*>) k (Herm u e a) = Herm u (k*e) a++instance (Show (a mn e), Show e, Num e) => Show (Herm a mn e) where+ show (Herm u k a) + | k /= 1 = "(" ++ show k ++ ") *> " ++ show (Herm u 1 a)+ | otherwise =+ constructor ++ " (" ++ show a ++ ")"+ where+ constructor = case u of+ Lower -> "hermL"+ Upper -> "hermU"
+ Data/Matrix/Herm/Dense.hs view
@@ -0,0 +1,114 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : Data.Matrix.Herm.Dense+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Data.Matrix.Herm.Dense (+ module Data.Matrix.Herm,+ module BLAS.Matrix.Immutable,+ module BLAS.Matrix.ReadOnly,++ hemv,+ hemm,+ ) where++import Control.Monad ( zipWithM_ )+import System.IO.Unsafe+import Unsafe.Coerce++import BLAS.Access+import BLAS.Elem ( BLAS2, BLAS3 )+import BLAS.C ( colMajor, rightSide, leftSide, cblasUpLo )+import BLAS.Types ( flipUpLo )+import qualified BLAS.Elem as E+import qualified BLAS.C as BLAS++import Data.Matrix.Dense.Internal+import Data.Vector.Dense.Internal+import qualified Data.Matrix.Dense.Internal as M+import qualified Data.Vector.Dense.Internal as V+import qualified Data.Vector.Dense.Operations as V++import Data.Matrix.Herm+import BLAS.Matrix.Immutable+import BLAS.Matrix.ReadOnly+++instance (BLAS3 e) => IMatrix (Herm (DMatrix Imm)) e where+ (<*>) h x = unsafePerformIO $ getApply h x+ {-# NOINLINE (<*>) #-}++ (<**>) h a = unsafePerformIO $ getApplyMat h a+ {-# NOINLINE (<**>) #-}+++instance (BLAS3 e) => RMatrix (Herm (DMatrix s)) e where+ getApply h x = do+ y <- newZero (dim x)+ hemv 1 (unsafeCoerce h) x 1 y+ return (unsafeCoerce y)+ + getApplyMat h a = do+ b <- newZero (shape a)+ hemm 1 (unsafeCoerce h) a 1 b+ return (unsafeCoerce b)+++hemv :: (BLAS2 e) => e -> Herm (DMatrix t) (n,n) e -> DVector s n e -> e -> IOVector n e -> IO ()+hemv alpha h x beta y+ | numRows h == 0 =+ return ()+ | isConj y = do+ V.doConj y+ hemv alpha h x beta (V.conj y)+ V.doConj y+ | isConj x = do+ x' <- newCopy x+ V.doConj (V.unsafeThaw x')+ hemv alpha h (conj x') beta y+ | otherwise =+ let order = colMajor+ (u,e,a) = toBase h+ alpha' = alpha * e+ n = numCols a+ (u',alpha'') + = case (isHerm a) of+ True -> (flipUpLo u, E.conj alpha')+ False -> (u, alpha')+ uploA = cblasUpLo u'+ ldA = ldaOf a+ incX = strideOf x+ incY = strideOf y+ in M.unsafeWithElemPtr a (0,0) $ \pA ->+ V.unsafeWithElemPtr x 0 $ \pX ->+ V.unsafeWithElemPtr y 0 $ \pY -> do+ BLAS.hemv order uploA n alpha'' pA ldA pX incX beta pY incY++hemm :: (BLAS3 e) => e -> Herm (DMatrix t) (m,m) e -> DMatrix s (m,n) e -> e -> IOMatrix (m,n) e -> IO ()+hemm alpha h b beta c+ | numRows b == 0 || numCols b == 0 || numCols c == 0 = return ()+ | (isHerm a) /= (isHerm c) || (isHerm a) /= (isHerm b) =+ zipWithM_ (\x y -> hemv alpha h x beta y) (cols b) (cols c)+ | otherwise =+ let order = colMajor+ (m,n) = shape c+ alpha' = alpha * e+ (side,u',m',n', alpha'')+ = if isHerm a+ then (rightSide, flipUpLo u, n, m, E.conj alpha')+ else (leftSide, u, m, n, alpha')+ uploA = cblasUpLo u'+ ldA = ldaOf a+ ldB = ldaOf b+ ldC = ldaOf c+ in M.unsafeWithElemPtr a (0,0) $ \pA ->+ M.unsafeWithElemPtr b (0,0) $ \pB ->+ M.unsafeWithElemPtr c (0,0) $ \pC ->+ BLAS.hemm order side uploA m' n' alpha'' pA ldA pB ldB beta pC ldC+ where+ (u,e,a) = toBase h
+ Data/Matrix/Tri.hs view
@@ -0,0 +1,73 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : Data.Matrix.Tri+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Data.Matrix.Tri (+ Tri(..),+ UpLo(..), Diag(..),++ fromBase,+ toBase,+ mapTri,++ lower,+ lowerU,+ upper,+ upperU,++ ) where++import qualified BLAS.Elem as E+import BLAS.Matrix+import BLAS.Tensor+import BLAS.Types ( UpLo(..), Diag(..), flipUpLo )++data Tri a nn e = Tri UpLo Diag e (a nn e)++mapTri :: (a (n,n) e -> b (n,n) e) -> Tri a (n,n) e -> Tri b (n,n) e+mapTri f (Tri u d n a) = Tri u d n $ f a++fromBase :: UpLo -> Diag -> e -> a (n,n) e -> Tri a (n,n) e+fromBase = Tri+ +toBase :: Tri a (n,n) e -> (UpLo, Diag, e, a (n,n) e)+toBase (Tri u d e a) = (u,d,e,a)++lower :: (Num e) => a (n,n) e -> Tri a (n,n) e+lower = Tri Lower NonUnit 1++lowerU :: (Num e) => a (n,n) e -> Tri a (n,n) e+lowerU = Tri Lower Unit 1++upper :: (Num e) => a (n,n) e -> Tri a (n,n) e+upper = Tri Upper NonUnit 1++upperU :: (Num e) => a (n,n) e -> Tri a (n,n) e+upperU = Tri Upper Unit 1+ +instance Matrix a => Matrix (Tri a) where+ numRows (Tri _ _ _ a) = numRows a+ numCols (Tri _ _ _ a) = numCols a+ herm (Tri u d e a) = Tri (flipUpLo u) d (E.conj e) (herm a)+ +instance (Num e) => Scalable (Tri a nn) e where+ (*>) k (Tri u d e a) = Tri u d (k*e) a++instance (Show (a mn e), Show e, Num e) => Show (Tri a mn e) where+ show (Tri u d k a) + | k /= 1 = "(" ++ show k ++ ") *> " ++ show (Tri u d 1 a)+ | otherwise =+ constructor ++ " (" ++ show a ++ ")"+ where+ constructor = case (u,d) of+ (Lower, NonUnit) -> "lower"+ (Lower, Unit ) -> "lowerU"+ (Upper, NonUnit) -> "upper"+ (Upper, Unit ) -> "upperU"+
+ Data/Matrix/Tri/Dense.hs view
@@ -0,0 +1,171 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : Data.Matrix.Tri.Dense+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Data.Matrix.Tri.Dense (+ module Data.Matrix.Tri,+ module BLAS.Matrix.Immutable,+ module BLAS.Matrix.ReadOnly,+ module BLAS.Matrix.Solve,++ trmv,+ trsv,+ trmm, + trsm+ ) where++import Control.Monad ( when )+import Data.Maybe ( fromJust )+import System.IO.Unsafe+import Unsafe.Coerce++import Data.Matrix.Dense.Internal+import qualified Data.Matrix.Dense.Internal as M+import Data.Vector.Dense.Internal+import qualified Data.Vector.Dense.Internal as V+import qualified Data.Vector.Dense.Operations as V++import BLAS.Access+import BLAS.Elem ( BLAS3 )+import qualified BLAS.Elem as E+import BLAS.C.Types ( cblasDiag, cblasUpLo, cblasTrans, colMajor, + noTrans, conjTrans, leftSide, rightSide )+import BLAS.Types ( Trans(..), flipTrans, flipUpLo )+ +import qualified BLAS.C as BLAS+import qualified BLAS.C.Types as BLAS++import BLAS.Matrix.Immutable+import BLAS.Matrix.ReadOnly+import BLAS.Matrix.Solve++import Data.Matrix.Tri++instance (BLAS3 e) => IMatrix (Tri (DMatrix Imm)) e where+ (<*>) t x = unsafePerformIO $ getApply t x+ {-# NOINLINE (<*>) #-}++ (<**>) t a = unsafePerformIO $ getApplyMat t a+ {-# NOINLINE (<**>) #-}++instance (BLAS3 e) => ISolve (Tri (DMatrix Imm)) e where+ (<\>) t x = unsafePerformIO $ getSolve t x+ {-# NOINLINE (<\>) #-}++ (<\\>) t a = unsafePerformIO $ getSolveMat t a+ {-# NOINLINE (<\\>) #-}++instance (BLAS3 e) => RMatrix (Tri (DMatrix s)) e where+ getApply t x = do+ x' <- newCopy x+ trmv (unsafeCoerce t) (V.unsafeThaw x')+ return (unsafeCoerce x')+ + getApplyMat t a = do+ a' <- newCopy a+ trmm (unsafeCoerce t) (M.unsafeThaw a')+ return (unsafeCoerce a')++instance (BLAS3 e) => RSolve (Tri (DMatrix s)) e where+ getSolve t x = do+ x' <- newCopy x+ trsv (unsafeCoerce t) (V.unsafeThaw x')+ return (unsafeCoerce x')+ + getSolveMat t a = do+ a' <- newCopy a+ trsm (unsafeCoerce t) (M.unsafeThaw a')+ return (unsafeCoerce a')+++trmv :: (BLAS3 e) => Tri (DMatrix t) (n,n) e -> IOVector n e -> IO ()+trmv _ x+ | dim x == 0 = return ()+trmv t x+ | isConj x =+ let b = fromJust $ maybeFromCol x+ in trmm t b+trmv t x =+ let (u,d,alpha,a) = toBase t+ order = colMajor+ (transA,u') = if isHerm a then (conjTrans, flipUpLo u) else (noTrans, u)+ uploA = cblasUpLo u'+ diagA = cblasDiag d+ n = numCols a+ ldA = ldaOf a+ incX = strideOf x+ in M.unsafeWithElemPtr a (0,0) $ \pA ->+ V.unsafeWithElemPtr x 0 $ \pX -> do+ BLAS.trmv order uploA transA diagA n pA ldA pX incX+ when (alpha /= 1) $ V.scaleBy alpha x+ + +trmm :: (BLAS3 e) => Tri (DMatrix t) (m,m) e -> IOMatrix (m,n) e -> IO ()+trmm _ b+ | M.numRows b == 0 || M.numCols b == 0 = return ()+trmm t b =+ let (u,d,alpha,a) = toBase t+ order = colMajor+ (h,u') = if isHerm a then (ConjTrans, flipUpLo u) else (NoTrans, u)+ (m,n) = shape b+ (side,h',m',n',alpha')+ = if M.isHerm b+ then (rightSide, flipTrans h, n, m, E.conj alpha)+ else (leftSide , h , m, n, alpha )+ uploA = cblasUpLo u'+ transA = cblasTrans h'+ diagA = cblasDiag d+ ldA = ldaOf a+ ldB = ldaOf b+ in M.unsafeWithElemPtr a (0,0) $ \pA ->+ M.unsafeWithElemPtr b (0,0) $ \pB ->+ BLAS.trmm order side uploA transA diagA m' n' alpha' pA ldA pB ldB++trsv :: (BLAS3 e) =>Tri (DMatrix t) (n,n) e -> IOVector n e -> IO ()+trsv _ x+ | dim x == 0 = return ()+trsv t x+ | isConj x =+ let b = fromJust $ maybeFromCol x+ in trsm t b+trsv t x =+ let (u,d,alpha,a) = toBase t+ order = colMajor+ (transA,u') = if isHerm a then (conjTrans, flipUpLo u) else (noTrans, u)+ uploA = cblasUpLo u'+ diagA = cblasDiag d+ n = numCols a+ ldA = ldaOf a+ incX = strideOf x+ in M.unsafeWithElemPtr a (0,0) $ \pA ->+ V.unsafeWithElemPtr x 0 $ \pX -> do+ BLAS.trsv order uploA transA diagA n pA ldA pX incX+ when (alpha /= 1) $ V.invScaleBy alpha x++trsm :: (BLAS3 e) => Tri (DMatrix t) (m,m) e -> IOMatrix (m,n) e -> IO ()+trsm _ b+ | M.numRows b == 0 || M.numCols b == 0 = return ()+trsm t b =+ let (u,d,alpha,a) = toBase t+ order = colMajor+ (h,u') = if isHerm a then (ConjTrans, flipUpLo u) else (NoTrans, u)+ (m,n) = shape b+ (side,h',m',n',alpha')+ = if isHerm b+ then (rightSide, flipTrans h, n, m, E.conj alpha)+ else (leftSide , h , m, n, alpha )+ uploA = cblasUpLo u'+ transA = cblasTrans h'+ diagA = cblasDiag d+ ldA = ldaOf a+ ldB = ldaOf b+ in M.unsafeWithElemPtr a (0,0) $ \pA ->+ M.unsafeWithElemPtr b (0,0) $ \pB -> do+ BLAS.trsm order side uploA transA diagA m' n' (1/alpha') pA ldA pB ldB+
+ Data/Vector/Dense.hs view
@@ -0,0 +1,135 @@+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : Data.Vector.Dense+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--+module Data.Vector.Dense (+ Vector,+ module BLAS.Vector,+ module BLAS.Tensor.Base,+ module BLAS.Tensor.Dense.Immutable,+ module BLAS.Tensor.Immutable,+ module BLAS.Tensor.Scalable,++ -- * Creating vectors+ vector, + listVector,++ -- * Special vectors+ basis,++ -- * Augmenting vectors+ subvector,+ subvectorWithStride,++ -- * Norms and dot product+ sumAbs,+ norm2,+ whichMaxAbs,+ (<.>),++ -- * Vector arithmetic+ shift,+ scale,+ invScale,++ -- * Converting to and from lists+ toList,+ fromList,+ + -- * Casting vectors+ coerceVector,+ + -- * Unsafe operations+ unsafeVector,+ unsafeSubvector,+ unsafeSubvectorWithStride,+ ) where++import System.IO.Unsafe ++import Data.Vector.Dense.Internal +import Data.Vector.Dense.Operations++import BLAS.Access+import BLAS.Elem ( BLAS1, BLAS2 )+import BLAS.Vector hiding ( Vector )+import BLAS.Tensor.Base+import BLAS.Tensor.Dense.Immutable+import BLAS.Tensor.Immutable+import BLAS.Tensor.Scalable+++-- | Create a vector with the given dimension and elements. The elements+-- given in the association list must all have unique indices, otherwise+-- the result is undefined.+vector :: (BLAS1 e) => Int -> [(Int, e)] -> Vector n e+vector n ies = unsafeFreeze $ unsafePerformIO $ newVector n ies+{-# NOINLINE vector #-}++-- | Same as 'vector', but does not range-check the indices.+unsafeVector :: (BLAS1 e) => Int -> [(Int, e)] -> Vector n e+unsafeVector n ies = unsafeFreeze $ unsafePerformIO $ unsafeNewVector n ies+{-# NOINLINE unsafeVector #-}++-- | Create a vector of the given dimension with elements initialized+-- to the values from the list.+listVector :: (BLAS1 e) => Int -> [e] -> Vector n e+listVector n es = unsafePerformIO $ newListVector n es+{-# NOINLINE listVector #-}++-- | @basis n i@ creates a vector of dimension @n@ with zeros everywhere but+-- position @i@, where there is a one.+basis :: (BLAS1 e) => Int -> Int -> Vector n e+basis n i = unsafeFreeze $ unsafePerformIO $ newBasis n i+{-# NOINLINE basis #-}++-- | Convert a vector to a list. Same as @elems@.+toList :: (BLAS1 e) => Vector n e -> [e]+toList = elems++-- | Convert a list to a vector. @fromList xs = listVector (length xs) xs@.+fromList :: (BLAS1 e) => [e] -> Vector n e+fromList es = listVector (length es) es+{-# INLINE fromList #-}++instance (BLAS1 e) => Scalable (DVector Imm n) e where+ (*>) = scale+ +instance (BLAS2 e) => Num (DVector Imm n e) where+ (+) = plus+ (-) = minus+ (*) = times+ negate = (*>) (-1)+ abs = amap abs+ signum = amap signum+ fromInteger = (constant 1) . fromInteger+ +instance (BLAS2 e) => Fractional (DVector Imm n e) where+ (/) = divide+ recip = amap recip+ fromRational = (constant 1) . fromRational + +instance (BLAS2 e, Floating e) => Floating (DVector Imm n e) where+ pi = constant 1 pi+ exp = amap exp+ sqrt = amap sqrt+ log = amap log+ (**) = azipWith (**)+ sin = amap sin+ cos = amap cos+ tan = amap tan+ asin = amap asin+ acos = amap acos+ atan = amap atan+ sinh = amap sinh+ cosh = amap cosh+ tanh = amap tanh+ asinh = amap asinh+ acosh = amap acosh+ atanh = amap atanh+
+ Data/Vector/Dense/IO.hs view
@@ -0,0 +1,64 @@+-----------------------------------------------------------------------------+-- |+-- Module : Data.Vector.Dense.Mutable+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Data.Vector.Dense.IO (+ -- * Mutable dense vector data type+ IOVector,+ DVector(..),+ + module BLAS.Vector,+ module BLAS.Tensor.Base,+ module BLAS.Tensor.Dense.ReadOnly,+ module BLAS.Tensor.ReadOnly,+ module BLAS.Tensor.Mutable,+ + -- * Creating vectors+ newVector, + newVector_,+ newListVector,+ + -- * Special vectors+ newBasis,+ setBasis,++ -- * Vector views+ subvector,+ subvectorWithStride,++ module Data.Vector.Dense.Operations,+ + -- * Conversion to and from @ForeignPtr@s+ fromForeignPtr,+ toForeignPtr,+ isConj,+ strideOf,+ + -- * Converting between mutable and immutable vectors+ unsafeFreeze,+ unsafeThaw,+ + -- * Casting vectors+ coerceVector,+ + -- * Unsafe operations+ unsafeNewVector,+ unsafeWithElemPtr,+ unsafeSubvector,+ unsafeSubvectorWithStride,+ + ) where++import Data.Vector.Dense.Internal+import Data.Vector.Dense.Operations hiding ( axpy, sumAbs, norm2, whichMaxAbs, + (<.>), shift, scale, invScale, plus, minus, times, divide )+import BLAS.Vector hiding ( Vector )+import BLAS.Tensor.Base+import BLAS.Tensor.Dense.ReadOnly+import BLAS.Tensor.ReadOnly+import BLAS.Tensor.Mutable
+ Data/Vector/Dense/Internal.hs view
@@ -0,0 +1,422 @@+{-# LANGUAGE BangPatterns, FlexibleInstances, MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module : Data.Vector.Dense.Internal+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Data.Vector.Dense.Internal (+ -- * Vector data types+ Vector,+ IOVector,+ DVector(..),++ module BLAS.Vector,+ module BLAS.Tensor,++ -- * Conversion to and from @ForeignPtr@s.+ fromForeignPtr,+ toForeignPtr,+ isConj,+ strideOf,+ + -- * Creating new vectors+ newVector, + newVector_,+ newListVector,+ + -- * Special vectors+ newBasis,+ setBasis,++ -- * Vector views+ subvector,+ subvectorWithStride,+ + -- * Casting vectors+ coerceVector,+ + -- * Unsafe operations+ unsafeNewVector,+ unsafeWithElemPtr,+ unsafeSubvector,+ unsafeSubvectorWithStride,+ + unsafeFreeze,+ unsafeThaw,+ + ) where+ + +import Control.Monad+import Data.Ix+import Foreign +import System.IO.Unsafe +import Unsafe.Coerce++import Data.AEq++import BLAS.Access+import BLAS.Elem.Base ( Elem )+import qualified BLAS.Elem.Base as E+import BLAS.Vector hiding ( Vector )+import qualified BLAS.Vector as C+import BLAS.Tensor++import BLAS.Internal ( clearArray, inlinePerformIO, checkedSubvector,+ checkedSubvectorWithStride )+import BLAS.C.Level1 ( BLAS1, copy )+++-- | A dense vector. @t@ is a type that will usually be @Imm@ or @Mut@. +-- @n@ is a phantom type for the dimension of the vector, and @e@ is the +-- element type. A @DVector@ @x@ stores @dim x@ elements. Indices into+-- the vector are @0@-based.+data DVector t n e = + DV { fptr :: !(ForeignPtr e) -- ^ a pointer to the storage region+ , offset :: !Int -- ^ an offset (in elements, not bytes) to the first element in the vector. + , len :: !Int -- ^ the length of the vector+ , stride :: !Int -- ^ the stride (in elements, not bytes) between elements.+ }+ | C !(DVector t n e) -- ^ a conjugated vector++type Vector n e = DVector Imm n e+type IOVector n e = DVector Mut n e++-- | Cast the phantom length type.+coerceVector :: DVector t n e -> DVector t m e+coerceVector = unsafeCoerce++-- | Gets the pointer to the storage block+storageOf :: DVector t n e -> ForeignPtr e+storageOf (C x) = storageOf x+storageOf x@(DV _ _ _ _) = fptr x++-- | Gets the stride of the vector.+strideOf :: DVector t n e -> Int+strideOf (C x) = strideOf x+strideOf x@(DV _ _ _ _) = stride x+{-# INLINE strideOf #-}++-- | Indicates whether or not the vector has been conjugated. For +-- newly-created vectors, this will be @False@.+isConj :: (Elem e) => DVector t n e -> Bool+isConj (C x) = not (isConj x)+isConj (DV _ _ _ _) = False+{-# INLINE isConj #-}++-- | @fromForeignPtr fptr offset n inc@ creates a vector view of a+-- region in memory starting at the given offset and having dimension @n@,+-- with a stride of @inc@.+fromForeignPtr :: ForeignPtr e -> Int -> Int -> Int -> DVector t n e+fromForeignPtr = DV+{-# INLINE fromForeignPtr #-}++-- | Gets the tuple @(fptr,offset,n,inc)@, where @n@ is the dimension and +-- @inc@ is the stride of the vector. Note that this does not return the+-- conjugacy information of the vector. For that information, use @isConj@.+toForeignPtr :: DVector t n e -> (ForeignPtr e, Int, Int, Int)+toForeignPtr (C x) = toForeignPtr x+toForeignPtr (DV f o n s) = (f, o, n, s)+{-# INLINE toForeignPtr #-}++-- | @subvector x o n@ creates a subvector view of @x@ starting at index @o@ +-- and having length @n@.+subvector :: DVector t n e -> Int -> Int -> DVector t m e+subvector x = checkedSubvector (dim x) (unsafeSubvector x)++-- | Same as 'subvector' but arguments are not range-checked.+unsafeSubvector :: DVector t n e -> Int -> Int -> DVector t m e+unsafeSubvector = unsafeSubvectorWithStride 1++-- | @subvectorWithStride s x o n@ creates a subvector view of @x@ starting +-- at index @o@, having length @n@ and stride @s@.+subvectorWithStride :: Int -> DVector t n e -> Int -> Int -> DVector t m e+subvectorWithStride s x = + checkedSubvectorWithStride s (dim x) (unsafeSubvectorWithStride s x)+ +-- | Same as 'subvectorWithStride' but arguments are not range-checked.+unsafeSubvectorWithStride :: Int -> DVector t n e -> Int -> Int -> DVector t m e+unsafeSubvectorWithStride s (C x) o n = C $ unsafeSubvectorWithStride s x o n+unsafeSubvectorWithStride s x@(DV _ _ _ _) o n =+ let f = fptr x+ o' = indexOf x o+ n' = n+ s' = s * (stride x)+ in + fromForeignPtr f o' n' s'++-- | Creates a new vector of the given length. The elements will be +-- uninitialized.+newVector_ :: (Elem e) => Int -> IO (DVector t n e)+newVector_ n+ | n < 0 = + ioError $ userError $ + "Tried to create a vector with `" ++ show n ++ "' elements."+ | otherwise = do+ arr <- mallocForeignPtrArray n+ return $ fromForeignPtr arr 0 n 1++-- | Creates a new vector of the given dimension with the given elements.+-- If the list has length less than the passed-in dimenson, the tail of+-- the vector will be uninitialized.+newListVector :: (Elem e) => Int -> [e] -> IO (DVector t n e)+newListVector n es = do+ x <- newVector_ n+ withForeignPtr (fptr x) $ flip pokeArray $ take n es+ return x++-- | @listVector n es@ is equivalent to @vector n (zip [0..(n-1)] es)@, except+-- that the result is undefined if @length es@ is less than @n@.+listVector :: (Elem e) => Int -> [e] -> Vector n e+listVector n es = unsafeFreeze $ unsafePerformIO $ newListVector n es+{-# NOINLINE listVector #-}+++-- | Creates a new vector with the given association list. Unspecified+-- indices will get initialized to zero.+newVector :: (BLAS1 e) => Int -> [(Int,e)] -> IO (DVector t n e)+newVector =+ newVectorHelp writeElem++-- | Same as 'newVector' but indices are not range-checked.+unsafeNewVector :: (BLAS1 e) => Int -> [(Int,e)] -> IO (DVector t n e)+unsafeNewVector =+ newVectorHelp unsafeWriteElem++newVectorHelp :: (BLAS1 e) => + (IOVector n e -> Int -> e -> IO ()) + -> Int -> [(Int,e)] -> IO (DVector t n e)+newVectorHelp set n ies = do+ x <- newZero n+ mapM_ (uncurry $ set x) ies+ return (unsafeCoerce x)++-- | @newBasis n i@ creates a vector of length @n@ that is all zero except for+-- at position @i@, where it equal to one.+newBasis :: (BLAS1 e) => Int -> Int -> IO (IOVector n e)+newBasis n i = do+ x <- newVector_ n+ setBasis i x+ return x++-- | @setBasis x i@ sets the @i@th coordinate of @x@ to @1@, and all other+-- coordinates to @0@. If the vector has been scaled, it is possible that+-- @readVector x i@ will not return exactly @1@. See 'setElem'.+setBasis :: (BLAS1 e) => Int -> IOVector n e -> IO ()+setBasis i x+ | i < 0 || i >= dim x =+ ioError $ userError $ + "tried to set a vector of dimension `" ++ show (dim x) ++ "'"+ ++ " to basis vector `" ++ show i ++ "'"+ | otherwise = do+ setZero x+ unsafeWriteElem x i 1 ++indexOf :: DVector t n e -> Int -> Int+indexOf (C x) i = indexOf x i+indexOf x@(DV _ _ _ _) i = offset x + i * stride x+{-# INLINE indexOf #-}++-- | Evaluate a function with a pointer to the value stored at the given+-- index. Note that the value may need to conjugated before using it. See+-- 'isConj'.+unsafeWithElemPtr :: (Elem e) => DVector t n e -> Int -> (Ptr e -> IO a) -> IO a+unsafeWithElemPtr x i f =+ withForeignPtr (storageOf x) $ \ptr ->+ let elemPtr = ptr `advancePtr` (indexOf x i)+ in f elemPtr+{-# INLINE unsafeWithElemPtr #-}++-- | Cast the access type to @Imm@.+unsafeFreeze :: DVector t n e -> Vector n e+unsafeFreeze = unsafeCoerce++-- | Cast the access type to @Mut@.+unsafeThaw :: DVector t n e -> IOVector n e+unsafeThaw = unsafeCoerce++instance C.Vector (DVector t) where+ dim x = case x of+ (C x') -> dim x'+ _ -> len x+ {-# INLINE dim #-}++ conj x = case x of+ (C x') -> x'+ _ -> C x+ {-# INLINE conj #-}++{-# RULES + "conj/Float" conj = conjFloat + "conj/Double" conj = conjDouble+ #-}+conjFloat :: DVector t n Float -> DVector t n Float+conjFloat = id++conjDouble :: DVector t n Double -> DVector t n Double+conjDouble = id++++instance Tensor (DVector t n) Int e where+ shape = dim+ {-# INLINE shape #-}++ bounds x = (0, dim x - 1)+ {-# INLINE bounds #-}++instance (BLAS1 e) => ITensor (DVector Imm n) Int e where+ size = dim+ + indices = range . bounds+ {-# INLINE indices #-}++ elems = inlinePerformIO . getElems . unsafeThaw+ assocs = inlinePerformIO . getAssocs . unsafeThaw++ unsafeAt x = inlinePerformIO . unsafeReadElem (unsafeThaw x)+ {-# INLINE unsafeAt #-}+ + amap f x = listVector (dim x) (map f $ elems x)+ + azipWith f x y+ | dim y /= n =+ error ("amap2: vector lengths differ; first has length `" +++ show n ++ "' and second has length `" +++ show (dim y) ++ "'")+ | otherwise =+ listVector n (zipWith f (elems x) (elems y))+ where+ n = dim x+ + (//) = replaceHelp writeElem+ unsafeReplace = replaceHelp unsafeWriteElem++replaceHelp :: (BLAS1 e) => + (IOVector n e -> Int -> e -> IO ())+ -> Vector n e -> [(Int, e)] -> Vector 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) => IDTensor (DVector Imm n) Int e where+ zero n = unsafeFreeze $ unsafePerformIO $ newZero n+ {-# NOINLINE zero #-}+ + constant n e = unsafeFreeze $ unsafePerformIO $ newConstant n e+ {-# NOINLINE constant #-}+++instance (BLAS1 e) => RTensor (DVector t n) Int e IO where+ getSize = return . dim+ + newCopy x = case x of+ (C x') -> newCopy x' >>= return . C+ _ -> do+ y <- newVector_ (dim x)+ unsafeWithElemPtr x 0 $ \pX ->+ unsafeWithElemPtr y 0 $ \pY ->+ let n = dim x+ incX = strideOf x+ incY = strideOf y+ in copy n pX incX pY incY >> + return y+ + getIndices = return . indices . unsafeFreeze+ {-# INLINE getIndices #-}+ + unsafeReadElem x i = case x of+ (C x') -> unsafeReadElem x' i >>= return . E.conj+ _ -> withForeignPtr (fptr x) $ \ptr ->+ peekElemOff ptr (indexOf x i) ++ getAssocs x = case x of+ (C x') -> getAssocs x' >>= return . map (\(i,e) -> (i,E.conj e))+ _ -> let (f,o,n,incX) = toForeignPtr x+ ptr = (unsafeForeignPtrToPtr f) `advancePtr` o+ in return $ go n f incX ptr 0+ where+ go !n !f !incX !ptr !i + | i >= n = + -- This is very important since we are doing unsafe IO.+ -- Otherwise, the DVector might get discared and the+ -- memory freed before all of the elements are read+ inlinePerformIO $ do+ touchForeignPtr f+ return []+ | otherwise =+ let e = inlinePerformIO $ peek ptr+ ptr' = ptr `advancePtr` incX+ i' = i + 1+ ies = go n f incX ptr' i'+ in e `seq` ((i,e):ies)+ {-# NOINLINE getAssocs #-}+++instance (BLAS1 e) => RDTensor (DVector t n) Int e IO where+ newZero n = newVector_ n >>= (\x -> setZero (unsafeThaw x) >> return x)+ + newConstant n e = newVector_ n >>= (\x -> setConstant e (unsafeThaw x) >> return x)+ + +instance (BLAS1 e) => MTensor (DVector Mut n) Int e IO where+ setZero x + | strideOf x == 1 = unsafeWithElemPtr x 0 $ flip clearArray (dim x)+ | otherwise = setConstant 0 x++ setConstant e x = case x of+ (C x') -> setConstant (E.conj e) x'+ _ -> unsafeWithElemPtr x 0 $ go (dim x)+ where+ go !n !ptr | n <= 0 = return ()+ | otherwise = let ptr' = ptr `advancePtr` (strideOf x)+ n' = n - 1+ in poke ptr e >> + go n' ptr'+ + unsafeWriteElem x i e = case x of+ (C x') -> unsafeWriteElem x' i $ E.conj e+ _ -> withForeignPtr (fptr x) $ \ptr -> + pokeElemOff ptr (indexOf x i) e+ + canModifyElem x i = return $ inRange (bounds x) i+ {-# INLINE canModifyElem #-}+ + modifyWith f x = case x of+ (C x') -> modifyWith (E.conj . f . E.conj) x'+ _ -> withForeignPtr (fptr x) $ go (dim x)+ where+ go !n !ptr | n <= 0 = return ()+ | otherwise = do+ peek ptr >>= poke ptr . f+ go (n-1) (ptr `advancePtr` incX)++ incX = strideOf x+ +compareHelp :: (BLAS1 e) => + (e -> e -> Bool) -> Vector n e -> Vector n e -> Bool+compareHelp cmp (C x) (C y) =+ compareHelp cmp x y+compareHelp cmp x y =+ (dim x == dim y) && (and $ zipWith cmp (elems x) (elems y))++instance (BLAS1 e, Eq e) => Eq (DVector Imm n e) where+ (==) = compareHelp (==)++instance (BLAS1 e, AEq e) => AEq (DVector Imm n e) where+ (===) = compareHelp (===)+ (~==) = compareHelp (~==)++instance (BLAS1 e, Show e) => Show (DVector Imm n e) where+ show x = case x of+ (C x') -> "conj (" ++ show x' ++ ")"+ _ -> "listVector " ++ show (dim x) ++ " " ++ show (elems x)+
+ Data/Vector/Dense/Operations.hs view
@@ -0,0 +1,373 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Data.Vector.Dense.Operations+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Data.Vector.Dense.Operations (+ -- * Copy and swap+ copyVector,+ swapVectors,++ -- * Vector norms and inner products+ -- ** Pure+ sumAbs,+ norm2,+ whichMaxAbs,+ (<.>),+ + getSumAbs,+ getNorm2,+ getWhichMaxAbs,+ getDot,++ -- * Vector arithmetic+ -- ** Pure+ shift,+ scale,+ invScale,+ add,+ plus,+ minus,+ times,+ divide,+ + -- ** Impure+ getShifted,+ getScaled,+ getInvScaled,+ getSum,+ getDiff,+ getProduct,+ getRatio,+ + -- * In-place arithmetic+ doConj,+ scaleBy,+ shiftBy,+ invScaleBy,+ (+=),+ (-=),+ (*=),+ (//=),+ + -- * Unsafe operations+ unsafeCopyVector,+ unsafeSwapVectors,+ + -- * BLAS calls+ axpy,+ + ) where+ +import Control.Monad ( forM_ )+import Data.Vector.Dense.Internal+import BLAS.Tensor+import BLAS.Elem.Base ( Elem )+import qualified BLAS.Elem.Base as E++import Foreign ( Ptr )+import System.IO.Unsafe+import Unsafe.Coerce++import BLAS.Internal ( inlinePerformIO, checkVecVecOp )+import BLAS.C hiding ( copy, swap, iamax, conj, axpy, acxpy )+import qualified BLAS.C as BLAS+import qualified BLAS.C.Types as T++infixl 7 <.>, `times`, `divide`, `scale`, `invScale`+infixl 6 `plus`, `minus`, `shift`+infixl 1 +=, -=, *=, //=+++-- | @copyVector dst src@ replaces the elements in @dst@ with the values from @src@.+-- This may result in a loss of precision.+copyVector :: (BLAS1 e) => IOVector n e -> DVector t n e -> IO ()+copyVector x y = checkVecVecOp "copyVector" (dim x) (dim y) >> unsafeCopyVector x y++-- | Same as 'copyVector' but does not check the dimensions of the arguments.+unsafeCopyVector :: (BLAS1 e) => IOVector n e -> DVector t n e -> IO ()+unsafeCopyVector (C x) (C y) =+ unsafeCopyVector x y+unsafeCopyVector x@(DV _ _ _ _) y@(DV _ _ _ _) =+ call2 BLAS.copy y x+unsafeCopyVector x y = do+ forM_ [0..(dim x - 1)] $ \i -> do+ unsafeReadElem y i >>= unsafeWriteElem x i+++-- | @swapVectors x y@ replaces the elements in @x@ with the values from @y@, and +-- replaces the elements in @y@ with the values from @x@. This may result in +-- a loss of precision.+swapVectors :: (BLAS1 e) => IOVector n e -> IOVector n e -> IO ()+swapVectors x y = checkVecVecOp "swapVectors" (dim x) (dim y) >> unsafeSwapVectors x y++-- | Same as 'swap' but does not check the dimensions of the arguments.+unsafeSwapVectors :: (BLAS1 e) => IOVector n e -> IOVector n e -> IO ()+unsafeSwapVectors (C x) (C y) =+ unsafeSwapVectors x y+unsafeSwapVectors x@(DV _ _ _ _) y@(DV _ _ _ _) =+ call2 BLAS.swap x y+unsafeSwapVectors x y = do+ forM_ [0..(dim x - 1)] $ \i -> do+ tmp <- unsafeReadElem x i+ unsafeReadElem y i >>= unsafeWriteElem x i+ unsafeWriteElem y i tmp++-- | Gets the sum of the absolute values of the vector entries.+getSumAbs :: (BLAS1 e) => DVector t n e -> IO Double+getSumAbs = call asum+ +-- | Gets the 2-norm of a vector.+getNorm2 :: (BLAS1 e) => DVector t n e -> IO Double+getNorm2 = call nrm2++-- | Gets the index and norm of the element with maximum magnitude. This is undefined+-- if any of the elements are @NaN@. It will throw an exception if the +-- dimension of the vector is 0.+getWhichMaxAbs :: (BLAS1 e) => DVector t n e -> IO (Int, e)+getWhichMaxAbs x =+ case (dim x) of+ 0 -> ioError $ userError $ "getWhichMaxAbs of an empty vector"+ _ -> do+ i <- call BLAS.iamax x+ e <- unsafeReadElem x i+ return (i,e)++-- | Computes the dot product of two vectors.+getDot :: (BLAS1 e) => DVector s n e -> DVector t n e -> IO e+getDot x y = checkVecVecOp "dot" (dim x) (dim y) >> unsafeGetDot x y++unsafeGetDot :: (BLAS1 e) => DVector s n e -> DVector t n e -> IO e+unsafeGetDot x@(DV _ _ _ _) y@(DV _ _ _ _) =+ call2 dotc x y+unsafeGetDot (C x@(DV _ _ _ _)) (y@(DV _ _ _ _)) =+ call2 dotu x y+unsafeGetDot (x@(DV _ _ _ _)) (C y@(DV _ _ _ _)) =+ call2 dotu x y >>= return . E.conj+unsafeGetDot x@(DV _ _ _ _) (C (C y)) = + unsafeGetDot x y+unsafeGetDot (C x) y = + unsafeGetDot x (conj y) >>= return . E.conj++-- | Create a new vector by adding a value to every element in a vector.+getShifted :: (BLAS1 e) => e -> DVector t n e -> IO (DVector r n e)+getShifted k x = do+ y <- newCopy x+ shiftBy k (unsafeThaw y)+ return (unsafeCoerce y)++-- | Create a new vector by scaling every element by a value. @scale'k x@+-- is equal to @newCopy' (scale k x)@.+getScaled :: (BLAS1 e) => e -> DVector t n e -> IO (DVector r n e)+getScaled k x = do+ y <- newCopy x+ scaleBy k (unsafeThaw y)+ return (unsafeCoerce y)++-- | Create a new vector by dividing every element by a value.+getInvScaled :: (BLAS1 e) => e -> DVector t n e -> IO (DVector r n e)+getInvScaled k x = do+ y <- newCopy x+ invScaleBy k (unsafeThaw y)+ return (unsafeCoerce y)++-- | Computes the sum of two vectors.+getSum :: (BLAS1 e) => e -> DVector s n e -> e -> DVector t n e -> IO (DVector r n e)+getSum alpha x beta y = checkVecVecOp "getSum" (dim x) (dim y) >> unsafeGetSum alpha x beta y++unsafeGetSum :: (BLAS1 e) => e -> DVector s n e -> e -> DVector t n e -> IO (DVector r n e)+unsafeGetSum 1 x beta y+ | beta /= 1 = unsafeGetSum beta y 1 x+unsafeGetSum alpha (C x) beta y = do+ s <- unsafeGetSum (E.conj alpha) x (E.conj beta) (conj y)+ return (conj s)+unsafeGetSum alpha x@(DV _ _ _ _) beta y = do+ s <- newCopy y+ scaleBy beta (unsafeThaw s)+ axpy alpha x (unsafeThaw s)+ return (unsafeCoerce s)+ +-- | Computes the difference of two vectors.+getDiff :: (BLAS1 e) => DVector s n e -> DVector t n e -> IO (DVector r n e)+getDiff x y = checkVecVecOp "getDiff" (dim x) (dim y) >> unsafeGetSum 1 x (-1) y++-- | Computes the elementwise product of two vectors.+getProduct :: (BLAS2 e) => DVector s n e -> DVector t n e -> IO (DVector r n e)+getProduct = binaryOp "getProduct" (*=)++-- | Computes the elementwise ratio of two vectors.+getRatio :: (BLAS2 e) => DVector s n e -> DVector t n e -> IO (DVector r n e)+getRatio = binaryOp "getRatio" (//=)++-- | Replace every element with its complex conjugate. See also 'conj'.+doConj :: (BLAS1 e) => IOVector n e -> IO ()+doConj = call BLAS.conj+ +-- | Add a value to every element in a vector.+shiftBy :: (BLAS1 e) => e -> IOVector n e -> IO ()+shiftBy alpha (C x) = shiftBy (E.conj alpha) x+shiftBy alpha x = modifyWith (alpha+) x++-- | Scale every element in the vector by the given value, and return a view+-- to the scaled vector. See also 'scale'.+scaleBy :: (BLAS1 e) => e -> IOVector n e -> IO ()+scaleBy 1 _ = return ()+scaleBy k (C x) = do+ scaleBy (E.conj k) x+scaleBy k x@(DV _ _ _ _) =+ call (flip scal k) x++-- | Divide every element by a value.+invScaleBy :: (BLAS1 e) => e -> IOVector n e -> IO ()+invScaleBy k (C x) = invScaleBy (E.conj k) x+invScaleBy k x = modifyWith (/k) x++-- | @y += x@ replaces @y@ by @y + x@.+(+=) :: (BLAS1 e) => IOVector n e -> DVector t n e -> IO ()+(+=) y x = checkVecVecOp "(+=)" (dim y) (dim x) >> axpy 1 x y++axpy :: (BLAS1 e) => e -> DVector t n e -> IOVector n e -> IO ()+axpy alpha x@(DV _ _ _ _) y@(DV _ _ _ _) =+ call2 (flip BLAS.axpy alpha) x y+axpy alpha (C x@(DV _ _ _ _)) y@(DV _ _ _ _) =+ call2 (flip BLAS.acxpy alpha) x y+axpy alpha (C (C x)) y = + axpy alpha x y+axpy alpha x (C y) =+ axpy (E.conj alpha) (conj x) y++-- | @y -= x@ replaces @y@ by @y - x@.+(-=) :: (BLAS1 e) => IOVector n e -> DVector t n e -> IO ()+(-=) y x = checkVecVecOp "(-=)" (dim y) (dim x) >> axpy (-1) x y++-- | @y *= x@ replaces @y@ by @x * y@, the elementwise product.+(*=) :: (BLAS2 e) => IOVector n e -> DVector t n e -> IO ()+(*=) y x = checkVecVecOp "(*=)" (dim y) (dim x) >> timesEquals y x++timesEquals :: (BLAS2 e) => IOVector n e -> DVector t n e -> IO ()+timesEquals y@(DV _ _ _ _) x@(DV _ _ _ _) =+ call2 (flip (tbmv T.colMajor T.upper T.noTrans T.nonUnit) 0) x y+timesEquals y@(DV _ _ _ _) (C x@(DV _ _ _ _)) =+ call2 (flip (tbmv T.colMajor T.upper T.conjTrans T.nonUnit) 0) x y +timesEquals y@(DV _ _ _ _) (C (C x)) = + timesEquals y x+timesEquals (C y) x =+ timesEquals y (conj x)++-- | @y //= x@ replaces @y@ by @y / x@, the elementwise ratio.+(//=) :: (BLAS2 e) => IOVector n e -> DVector t n e -> IO ()+(//=) y x = checkVecVecOp "(//=)" (dim y) (dim x) >> divideEquals y x++divideEquals :: (BLAS2 e) => IOVector n e -> DVector t n e -> IO ()+divideEquals y@(DV _ _ _ _) x@(DV _ _ _ _) =+ call2 (flip (tbsv T.colMajor T.upper T.noTrans T.nonUnit) 0) x y+divideEquals y@(DV _ _ _ _) (C x@(DV _ _ _ _)) =+ call2 (flip (tbsv T.colMajor T.upper T.conjTrans T.nonUnit) 0) x y+divideEquals y@(DV _ _ _ _) (C (C x)) = + divideEquals y x+divideEquals (C y) x =+ divideEquals y (conj x)+ +call :: (Elem e) => (Int -> Ptr e -> Int -> IO a) -> DVector t n e -> IO a+call f x =+ let n = dim x+ incX = strideOf x+ in unsafeWithElemPtr x 0 $ \pX -> f n pX incX++call2 :: (Elem e) => + (Int -> Ptr e -> Int -> Ptr e -> Int -> IO a) + -> DVector s n e -> DVector t m e -> IO a+call2 f x y =+ let n = dim x+ incX = strideOf x+ incY = strideOf y+ in unsafeWithElemPtr x 0 $ \pX ->+ unsafeWithElemPtr y 0 $ \pY ->+ f n pX incX pY incY++binaryOp :: (BLAS1 e) => String -> (IOVector n e -> DVector t n e -> IO ()) + -> DVector s n e -> DVector t n e -> IO (DVector r n e)+binaryOp name f x y =+ checkVecVecOp name (dim x) (dim y) >> do+ x' <- newCopy x >>= return . unsafeThaw+ f x' y+ return $! (unsafeCoerce x')++++-- | Compute the sum of absolute values of entries in the vector.+sumAbs :: (BLAS1 e) => Vector n e -> Double+sumAbs x = inlinePerformIO $ getSumAbs x+{-# NOINLINE sumAbs #-}++-- | Compute the 2-norm of a vector.+norm2 :: (BLAS1 e) => Vector n e -> Double+norm2 x = inlinePerformIO $ getNorm2 x+{-# NOINLINE norm2 #-}++-- | Get the index and norm of the element with absulte value. Not valid +-- if any of the vector entries are @NaN@. Raises an exception if the +-- vector has length @0@.+whichMaxAbs :: (BLAS1 e) => Vector n e -> (Int, e)+whichMaxAbs x = inlinePerformIO $ getWhichMaxAbs x+{-# NOINLINE whichMaxAbs #-}++-- | Compute the dot product of two vectors.+(<.>) :: (BLAS1 e) => Vector n e -> Vector n e -> e+(<.>) x y = inlinePerformIO $ getDot x y+{-# NOINLINE (<.>) #-}++-- | Add a value to every element in a vector.+shift :: (BLAS1 e) => e -> Vector n e -> Vector n e+shift k x = unsafePerformIO $ getShifted k x+{-# NOINLINE shift #-}++-- | Scale every element by the given value.+scale :: (BLAS1 e) => e -> Vector n e -> Vector n e+scale k x = unsafePerformIO $ getScaled k x+{-# NOINLINE scale #-}++-- | Divide every element by a value.+invScale :: (BLAS1 e) => e -> Vector n e -> Vector n e+invScale k x = unsafePerformIO $ getInvScaled k x+{-# NOINLINE invScale #-}++add :: (BLAS1 e) => e -> Vector n e -> e -> Vector n e -> Vector n e+add alpha x beta y = unsafePerformIO $ getSum alpha x beta y+{-# NOINLINE add #-}++-- | Sum of two vectors.+plus :: (BLAS1 e) => Vector n e -> Vector n e -> Vector n e+plus x y = add 1 x 1 y++-- | Difference of vectors.+minus :: (BLAS1 e) => Vector n e -> Vector n e -> Vector n e+minus x y = unsafePerformIO $ getDiff x y+{-# NOINLINE minus #-}++-- | Elementwise vector product.+times :: (BLAS2 e) => Vector n e -> Vector n e -> Vector n e+times x y = unsafePerformIO $ getProduct x y+{-# NOINLINE times #-}++-- | Elementwise vector ratio.+divide :: (BLAS2 e) => Vector n e -> Vector n e -> Vector n e+divide x y = unsafePerformIO $ getRatio x y+{-# NOINLINE divide #-}+++{-# RULES+"scale/plus" forall k l x y. plus (scale k x) (scale l y) = add k x l y+"scale1/plus" forall k x y. plus (scale k x) y = add k x 1 y+"scale2/plus" forall k x y. plus x (scale k y) = add 1 x k y++"scale/minus" forall k l x y. minus (scale k x) (scale l y) = add k x (-l) y+"scale1/minus" forall k x y. minus (scale k x) y = add k x (-1) y+"scale2/minus" forall k x y. minus x (scale k y) = add 1 x (-k) y+ #-}+
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) Patrick Perry <patperry@stanford.edu> 2008++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.++3. Neither the name of the author nor the names of his contributors+ may be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.
+ Setup.lhs view
@@ -0,0 +1,17 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> 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"+> return ()+>+> main = defaultMainWithHooks defaultUserHooks{ runTests=testing }+>
+ Test/QuickCheck/Complex.hs view
@@ -0,0 +1,24 @@+-----------------------------------------------------------------------------+-- |+-- Module : Test.QuickCheck.Complex+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Test.QuickCheck.Complex+ where++import Test.QuickCheck ++import Data.Complex++newtype TestComplex e = TestComplex (Complex e)++instance (Arbitrary e, RealFloat e) => Arbitrary (TestComplex e) where+ arbitrary = do+ (x,y) <- arbitrary+ return $ TestComplex (x :+ y)+ coarbitrary (TestComplex (x :+ y)) = + coarbitrary (x,y)
+ Test/QuickCheck/Matrix.hs view
@@ -0,0 +1,41 @@+-----------------------------------------------------------------------------+-- |+-- Module : Test.QuickCheck.Matrix+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Test.QuickCheck.Matrix+ where+ +import Data.List ( nub )+import Test.QuickCheck.Vector ( Index(..) )++import Test.QuickCheck hiding ( vector )+import qualified Test.QuickCheck as QC++newtype IndexPair = IndexPair (Int,Int) deriving (Eq, Show)+instance Arbitrary IndexPair where+ arbitrary = sized $ \k ->+ let k' = floor (sqrt $ fromInteger $ toInteger k :: Double)+ in resize k' $+ two arbitrary >>= (\((Index i), (Index j)) -> return $ IndexPair (i,j))+ coarbitrary (IndexPair ij) = coarbitrary ij+ +data Assocs e = Assocs (Int,Int) [((Int,Int),e)] deriving (Eq, Show)+instance Arbitrary e => Arbitrary (Assocs e) where+ arbitrary = sized $ \k -> + let k' = floor (sqrt $ fromInteger $ toInteger k :: Double)+ in do+ m <- choose (0,k')+ n <- choose (0,k')+ ijs <- QC.vector (m*n)+ let ijs' = filter (\(i,j) -> i < m && j < n) + $ nub + $ map (\(IndexPair ij) -> ij) ijs+ es <- QC.vector (length ijs')+ return $ Assocs (m,n) (zip ijs' es)+ + coarbitrary (Assocs mn ijes) = coarbitrary (mn,ijes)
+ Test/QuickCheck/Matrix/Dense.hs view
@@ -0,0 +1,157 @@+{-# LANGUAGE FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : Test.QuickCheck.Matrix.Dense+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Test.QuickCheck.Matrix.Dense+ where++import Test.QuickCheck hiding ( vector )+import qualified Test.QuickCheck as QC+import Test.QuickCheck.Vector.Dense ( TestVector(..), dvector )++import Data.Vector.Dense ( Vector, dim )+import Data.Matrix.Dense+import BLAS.Elem ( Elem, BLAS1 )++newtype TestMatrix m n e = TestMatrix (Matrix (m,n) e) deriving (Eq, Show)+data SubMatrix m n e = SubMatrix (Matrix (m,n) e) (Int,Int) (Int,Int) deriving (Eq, Show)++dmatrix :: (Elem e, Arbitrary e) => (Int,Int) -> Gen (Matrix (m,n) e)+dmatrix mn = frequency [ (3, rawMatrix mn) + , (2, hermedMatrix mn)+ , (1, subMatrix mn + >>= \(SubMatrix a ij _) -> + return $ submatrix a ij mn)+ ]++rawMatrix :: (Elem e, Arbitrary e) => (Int,Int) -> Gen (Matrix (m,n) e)+rawMatrix (m,n) = do+ es <- QC.vector (m*n)+ return $ listMatrix (m,n) es++hermedMatrix :: (Elem e, Arbitrary e) => (Int,Int) -> Gen (Matrix (m,n) e)+hermedMatrix (m,n) = do+ x <- dmatrix (n,m)+ return $ (herm x)++subMatrix :: (Elem e, Arbitrary e) => (Int,Int) -> Gen (SubMatrix m n e)+subMatrix (m,n) = do+ i <- choose (0,5)+ j <- choose (0,5)+ e <- choose (0,5)+ f <- choose (0,5)+ x <- dmatrix (i+m+e, j+n+f)++ return $ SubMatrix x (i,j) (m,n)++++instance (Arbitrary e, BLAS1 e) => Arbitrary (TestMatrix m n e) where+ arbitrary = sized $ \k -> + let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)+ in do+ m <- choose (0,k')+ n <- choose (0,k')+ dmatrix (m,n) >>= return . TestMatrix+ + coarbitrary (TestMatrix x) =+ coarbitrary (elems x)++data MatAt m n e = MatAt (Matrix (m,n) e) (Int,Int) deriving (Eq, Show)+instance (Arbitrary e, Elem e) => Arbitrary (MatAt m n e) where+ arbitrary = sized $ \k ->+ let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)+ in do+ m <- choose (1,k'+1)+ n <- choose (1,k'+1)+ i <- choose (0,m-1)+ j <- choose (0,n-1)+ a <- dmatrix (m,n)+ return $ MatAt a (i,j)++ coarbitrary = undefined+ ++instance (Arbitrary e, BLAS1 e) => Arbitrary (SubMatrix m n e) where+ arbitrary = sized $ \k -> + let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double) + in do+ m <- choose (0,k')+ n <- choose (0,k')+ (SubMatrix a ij mn) <- subMatrix (m,n)+ return $ SubMatrix a ij mn+ + coarbitrary (SubMatrix a ij mn) =+ coarbitrary (elems a, ij, mn)+ +data MatrixPair m n e = Pair (Matrix (m,n) e) (Matrix (m,n) e) deriving (Eq, Show)++instance (Arbitrary e, Elem e) => Arbitrary (MatrixPair m n e) where+ arbitrary = sized $ \k -> + let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)+ in do+ m <- choose (0,k')+ n <- choose (0,k')+ a <- dmatrix (m,n)+ b <- dmatrix (m,n)+ return $ Pair a b+ + coarbitrary = undefined+ +data MultMV m n e = MultMV (Matrix (m,n) e) (Vector n e) deriving (Eq, Show)++instance (Arbitrary e, BLAS1 e) => Arbitrary (MultMV m n e) where+ arbitrary = sized $ \k -> + let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)+ in do+ m <- choose (0,k')+ n <- choose (0,k')+ a <- dmatrix (m,n)+ x <- dvector n+ return $ MultMV a x+ + coarbitrary = undefined++data MultMVPair m n e = MultMVPair (Matrix (m,n) e) (Vector n e) (Vector n e) + deriving (Eq, Show)+ +instance (Arbitrary e, BLAS1 e) => Arbitrary (MultMVPair m n e) where+ arbitrary = do+ (MultMV a x) <- arbitrary+ y <- dvector (dim x)+ return $ MultMVPair a x y+ + coarbitrary (MultMVPair a x y) =+ coarbitrary (MultMV a x, TestVector y)+ +data MultMM m n k e = MultMM (Matrix (m,k) e) (Matrix (k,n) e) deriving (Eq, Show)++instance (Arbitrary e, Elem e) => Arbitrary (MultMM m n k e) where+ arbitrary = sized $ \s ->+ let s' = ceiling (sqrt $ fromInteger $ toInteger s :: Double)+ in do+ m <- choose (0,s')+ n <- choose (0,s')+ k <- choose (0,s')+ a <- dmatrix (m,k)+ b <- dmatrix (k,n)+ return $ MultMM a b+ + coarbitrary = undefined+ +data MultMMPair m n k e = MultMMPair (Matrix (m,k) e) (Matrix (k,n) e) (Matrix (k,n) e)+ deriving (Eq, Show)+ +instance (Arbitrary e, Elem e) => Arbitrary (MultMMPair m n k e) where+ arbitrary = do+ (MultMM a b) <- arbitrary+ c <- dmatrix (shape b)+ return $ MultMMPair a b c+ + coarbitrary = undefined
+ Test/QuickCheck/Matrix/Herm/Dense.hs view
@@ -0,0 +1,72 @@+-----------------------------------------------------------------------------+-- |+-- Module : Test.QuickCheck.Matrix.Herm.Dense+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Test.QuickCheck.Matrix.Herm.Dense + where++import Test.QuickCheck hiding ( vector )+import qualified Test.QuickCheck as QC+import Test.QuickCheck.Vector.Dense ( dvector )+import Test.QuickCheck.Matrix.Dense ( dmatrix, rawMatrix )++import BLAS.Elem ( BLAS2 )++import Data.Vector.Dense+import Data.Matrix.Dense+import Data.Matrix.Herm++++hermMatrix :: (BLAS2 e, Arbitrary e) => Int -> Gen (Matrix (n,n) e)+hermMatrix n = do+ a <- rawMatrix (n,n)+ return $ (a + herm a)++data HermMatrixMV n e = + HermMatrixMV UpLo (Matrix (n,n) e) (Matrix (n,n) e) (Vector n e) deriving (Show)++instance (Arbitrary e, BLAS2 e) => Arbitrary (HermMatrixMV n e) where+ arbitrary = sized $ \k ->+ let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)+ in do+ u <- elements [ Upper, Lower ]+ n <- choose (0,k')+ h <- hermMatrix n+ let f = case u of+ Upper -> \(i,j) -> i > j+ Lower -> \(i,j) -> i < j+ zs = zip (filter f (indices h)) (repeat 0)+ a = h // zs+ + x <- dvector n+ return $ HermMatrixMV u a h x+ coarbitrary = undefined+ +data HermMatrixMM m n e = + HermMatrixMM UpLo (Matrix (m,m) e) (Matrix (m,m) e) (Matrix (m,n) e) deriving (Show)+ +instance (Arbitrary e, BLAS2 e) => Arbitrary (HermMatrixMM m n e) where+ arbitrary = sized $ \k ->+ let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)+ in do+ u <- elements [ Upper, Lower ]+ m <- choose (0,k')+ n <- choose (0,k')+ h <- hermMatrix m+ let f = case u of+ Upper -> \(i,j) -> i > j+ Lower -> \(i,j) -> i < j+ zs = zip (filter f (indices h)) (repeat 0)+ a = h // zs++ b <- dmatrix (m,n)+ return $ HermMatrixMM u a h b+ + coarbitrary = undefined+
+ Test/QuickCheck/Matrix/Tri/Dense.hs view
@@ -0,0 +1,132 @@+{-# LANGUAGE FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module : Test.QuickCheck.Matrix.Tri.Dense+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Test.QuickCheck.Matrix.Tri.Dense+ where++import Data.Ix ( range )++import Test.QuickCheck hiding ( vector )+import qualified Test.QuickCheck as QC+import Test.QuickCheck.Vector.Dense ( TestVector(..), dvector )+import Test.QuickCheck.Matrix.Dense ( dmatrix )++import Data.Vector.Dense+import Data.Matrix.Dense+import BLAS.Elem ( BLAS1, BLAS3 )++import Data.Matrix.Tri.Dense ( Tri, UpLo(..), Diag(..), upper, lower, upperU, lowerU )++triMatrix :: (BLAS1 e, Arbitrary e) => UpLo -> Diag -> Int -> Gen (Matrix (n,n) e)+triMatrix u d n = + let nz = case d of+ NonUnit -> n * (n + 1) `div` 2+ Unit -> n * (n - 1) `div` 2+ in do+ h <- arbitrary+ let f = case (h,u,d) of+ (False, Upper, NonUnit) -> \(i,j) -> i <= j+ (False, Upper, Unit) -> \(i,j) -> i < j+ (False, Lower, NonUnit) -> \(i,j) -> i >= j+ (False, Lower, Unit) -> \(i,j) -> i > j+ (True, Upper, NonUnit) -> \(i,j) -> i >= j+ (True, Upper, Unit) -> \(i,j) -> i > j+ (True, Lower, NonUnit) -> \(i,j) -> i <= j+ (True, Lower, Unit) -> \(i,j) -> i < j+ ijs = filter f $ range ((0,0), (n-1,n-1))+ es <- QC.vector nz+ let a = matrix (n,n) $ zip ijs es++ a' <- case h of+ False -> return a+ True -> return (herm a)+ return a'+++data TriMatrix n e = TriMatrix UpLo Diag (Matrix (n,n) e) deriving (Eq, Show)++instance (Arbitrary e, BLAS1 e) => Arbitrary (TriMatrix n e) where+ arbitrary = sized $ \k ->+ let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)+ in do+ u <- elements [ Upper, Lower ]+ d <- elements [ Unit, NonUnit ]+ n <- choose (0,k')+ a <- triMatrix u d n+ return $ TriMatrix u d a+ + coarbitrary = undefined++data TriMatrixMV n e = + TriMatrixMV UpLo Diag (Matrix (n,n) e) (Vector n e) deriving (Eq, Show)++instance (Arbitrary e, BLAS1 e) => Arbitrary (TriMatrixMV n e) where+ arbitrary = do+ (TriMatrix u d a) <- arbitrary+ x <- dvector (numCols a)+ return $ TriMatrixMV u d a x+ + coarbitrary (TriMatrixMV u d a x) =+ coarbitrary (TriMatrix u d a, TestVector x)+ +data TriMatrixMM m n e = + TriMatrixMM UpLo Diag (Matrix (m,m) e) (Matrix (m,n) e) deriving (Eq, Show)++instance (Arbitrary e, BLAS1 e) => Arbitrary (TriMatrixMM m n e) where+ arbitrary = sized $ \k ->+ let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)+ in do+ (TriMatrix u d a) <- arbitrary+ n <- choose (0,k')+ b <- dmatrix (numCols a, n)+ return $ TriMatrixMM u d a b+ + coarbitrary = undefined+ +data TriMatrixSV n e = + TriMatrixSV (Tri Matrix (n,n) e) (Vector n e) deriving (Show)+ +instance (Arbitrary e, BLAS3 e) => Arbitrary (TriMatrixSV n e) where+ arbitrary = do+ (TriMatrix u d a) <- arbitrary+ let t = case (u,d) of+ (Lower,NonUnit) -> lower a+ (Lower,Unit) -> lowerU a+ (Upper,NonUnit) -> upper a+ (Upper,Unit) -> upperU a+ k <- arbitrary+ t' <- elements [ t+ , k *> t+ ]+ x <- dvector (numCols t')+ let y = t' <*> x+ return (TriMatrixSV t' y)+ + coarbitrary = undefined+++data TriMatrixSM m n e = + TriMatrixSM (Tri Matrix (m,m) e) (Matrix (m,n) e) + deriving (Show)+ +instance (Arbitrary e, BLAS3 e) => Arbitrary (TriMatrixSM m n e) where+ arbitrary = sized $ \k ->+ let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)+ in do+ (TriMatrixSV t _) <- arbitrary+ n <- choose (0, k')+ a <- dmatrix (numCols t, n)+ + let b = t <**> a+ return (TriMatrixSM t b)+ + coarbitrary = undefined++
+ Test/QuickCheck/Vector.hs view
@@ -0,0 +1,50 @@++-----------------------------------------------------------------------------+-- |+-- Module : Test.QuickCheck.Vector+-- Copyright : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Test.QuickCheck.Vector+ where++import Data.List ( nub )+import Test.QuickCheck hiding ( vector )+import qualified Test.QuickCheck as QC+++newtype Index = Index Int deriving (Eq, Show)+instance Arbitrary Index where+ arbitrary = sized $ \n -> do+ i <- if n == 0 + then return 0+ else choose (0,n-1)+ return (Index i)+ + coarbitrary (Index i) = coarbitrary i+ +data Basis = Basis Int Int deriving (Eq, Show)+instance Arbitrary Basis where+ arbitrary = do+ n <- arbitrary >>= (\(Index x) -> return (x+1))+ i <- choose (0,n-1)+ return $ Basis n i++ coarbitrary (Basis n i) = coarbitrary (n,i)++data Assocs e = Assocs Int [(Int,e)] deriving (Eq, Show)+instance Arbitrary e => Arbitrary (Assocs e) where+ arbitrary = sized $ \n -> do+ m <- choose (0,n)+ is <- QC.vector m + >>= mapM (\(Index i) -> return i) + >>= return . nub + >>= return . filter (<n)+ es <- QC.vector m+ return $ Assocs n (zip is es)+ + coarbitrary (Assocs n ies) = coarbitrary (n,ies)+
+ Test/QuickCheck/Vector/Dense.hs view
@@ -0,0 +1,92 @@+-----------------------------------------------------------------------------+-- |+-- Module : Test.QuickCheck.Vector.Dense+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++module Test.QuickCheck.Vector.Dense (+ TestVector(..),+ SubVector(..),+ VectorPair(..),+ VectorTriple(..),+ + dvector,+ rawVector,+ conjVector,+ subVector+ ) where++import Test.QuickCheck hiding ( vector )+import qualified Test.QuickCheck as QC++import Data.Vector.Dense+import BLAS.Elem ( BLAS1 )++newtype TestVector n e = TestVector (Vector n e) deriving (Eq, Show)+data SubVector n e = SubVector Int (Vector n e) Int Int deriving (Eq, Show)+data VectorPair n e = Pair (Vector n e) (Vector n e) deriving (Eq, Show)+data VectorTriple n e = Triple (Vector n e) (Vector n e) (Vector n e) deriving (Eq, Show)+++dvector :: (BLAS1 e, Arbitrary e) => Int -> Gen (Vector n e)+dvector n =+ frequency [ (3, rawVector n) + , (2, conjVector n)+ , (1, subVector n >>= \(SubVector s x o _) -> + return $ subvectorWithStride s x o n)+ ] ++rawVector :: (BLAS1 e, Arbitrary e) => Int -> Gen (Vector n e)+rawVector n = do+ es <- QC.vector n+ return $ listVector n es++conjVector :: (BLAS1 e, Arbitrary e) => Int -> Gen (Vector n e)+conjVector n = do+ x <- dvector n+ return $ (conj x)++subVector :: (BLAS1 e, Arbitrary e) => Int -> Gen (SubVector n e)+subVector n = do+ o <- choose (0,5)+ s <- choose (1,5)+ e <- choose (0,5)+ x <- dvector (o + s*n + e)+ return (SubVector s x o n)++instance (Arbitrary e, BLAS1 e) => Arbitrary (TestVector n e) where+ arbitrary = sized $ \m ->+ choose (0,m) >>= dvector >>= return . TestVector+ coarbitrary (TestVector x) =+ coarbitrary (elems x)++instance (Arbitrary e, BLAS1 e) => Arbitrary (SubVector n e) where+ arbitrary = sized $ \m -> + choose (0,m) >>= subVector+ coarbitrary (SubVector s x o n) = + coarbitrary (s,TestVector x,o,n)++instance (Arbitrary e, BLAS1 e) => Arbitrary (VectorPair n e) where+ arbitrary = sized $ \m -> do+ n <- choose (0,m)+ x <- dvector n+ y <- dvector n+ return $ Pair x y+ + coarbitrary (Pair x y) = + coarbitrary (TestVector x, TestVector y)+ +instance (Arbitrary e, BLAS1 e) => Arbitrary (VectorTriple n e) where+ arbitrary = sized $ \m -> do+ n <- choose (0,m)+ x <- dvector n+ y <- dvector n+ z <- dvector n+ return $ Triple x y z+ + coarbitrary (Triple x y z) = + coarbitrary (TestVector x, TestVector y, TestVector z)+
+ blas.cabal view
@@ -0,0 +1,95 @@+name: blas+Version: 0.4+homepage: http://stat.stanford.edu/~patperry/code/blas+synopsis: Bindings to the BLAS library+description:+ The BLAS (Basic Linear Algebra Subprograms) are routines that provide+ standard building blocks for performing basic vector and matrix operations. + The Level 1 BLAS perform scalar, vector and vector-vector operations, the + Level 2 BLAS perform matrix-vector operations, and the Level 3 BLAS perform+ matrix-matrix operations. Because the BLAS are efficient, portable, and+ widely available, they are commonly used in the development of high quality+ linear algebra software, LAPACK for example.+ .+ For more information, see the Netlib BLAS webpage:+ <http://www.netlib.org/blas/> + .+category: Math+license: BSD3+license-file: LICENSE+copyright: (c) 2008. Patrick Perry <patperry@stanford.edu>+author: Patrick Perry+maintainer: Patrick Perry <patperry@stanford.edu>+cabal-version: >= 1.2.0+build-type: Custom+tested-with: GHC == 6.8.2++extra-source-files: tests/Matrix.hs+ tests/HermMatrix.hs+ tests/TriMatrix.hs+ tests/Vector.hs++library+ exposed-modules: BLAS.Access+ BLAS.C+ BLAS.C.Level1+ BLAS.C.Level2+ BLAS.C.Level3+ BLAS.C.Types+ BLAS.Elem+ BLAS.Elem.Base+ BLAS.Internal+ BLAS.Matrix+ BLAS.Matrix.Base+ BLAS.Matrix.Immutable+ BLAS.Matrix.ReadOnly+ BLAS.Matrix.Solve+ BLAS.Matrix.Solve.Immutable+ BLAS.Matrix.Solve.ReadOnly+ BLAS.Tensor+ BLAS.Tensor.Base+ BLAS.Tensor.Dense+ BLAS.Tensor.Dense.Immutable+ BLAS.Tensor.Dense.ReadOnly+ BLAS.Tensor.Immutable+ BLAS.Tensor.ReadOnly+ BLAS.Tensor.Mutable+ BLAS.Tensor.Scalable+ BLAS.Types+ BLAS.Vector+ + Data.Matrix.Dense+ Data.Matrix.Dense.IO+ Data.Matrix.Dense.Internal+ Data.Matrix.Dense.Operations+ + Data.Matrix.Herm+ Data.Matrix.Herm.Dense+ + Data.Matrix.Tri+ Data.Matrix.Tri.Dense+ + Data.Vector.Dense+ Data.Vector.Dense.IO+ Data.Vector.Dense.Internal+ Data.Vector.Dense.Operations+ + Test.QuickCheck.Complex+ Test.QuickCheck.Vector+ Test.QuickCheck.Vector.Dense+ Test.QuickCheck.Matrix+ Test.QuickCheck.Matrix.Dense+ Test.QuickCheck.Matrix.Herm.Dense+ Test.QuickCheck.Matrix.Tri.Dense+ + other-modules: BLAS.C.Double+ BLAS.C.Zomplex + + ghc-options: -Wall + extensions: BangPatterns, EmptyDataDecls, FlexibleContexts, + FlexibleInstances, ForeignFunctionInterface, + FunctionalDependencies, MultiParamTypeClasses+ build-depends: base, ieee, storable-complex, QuickCheck++ extra-libraries: blas+
+ tests/HermMatrix.hs view
@@ -0,0 +1,96 @@+{-# OPTIONS -fglasgow-exts -fno-excess-precision -cpp #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++import Data.AEq+import Numeric.IEEE++import System.Environment ( getArgs )+import Test.QuickCheck.Parallel hiding ( vector )+import qualified Test.QuickCheck as QC+++import qualified BLAS.Elem as E+import Data.Complex ( Complex(..) )++import Data.Matrix.Herm.Dense+import Data.Matrix.Dense+import Data.Vector.Dense++import Test.QuickCheck.Complex+import Test.QuickCheck.Matrix.Herm.Dense++#ifdef COMPLEX+field = "Complex Double"+type E = Complex Double+#else+field = "Double"+type E = Double+#endif++type V = Vector Int E+type M = Matrix (Int,Int) E++instance (Arbitrary e, RealFloat e) => Arbitrary (Complex e) where+ arbitrary = arbitrary >>= \(TestComplex x) -> return x+ coarbitrary = coarbitrary . TestComplex++prop_herm_apply (HermMatrixMV u h (a :: M) x) =+ case u of+ Lower -> hermL h <*> x ~== a <*> x+ Upper -> hermU h <*> x ~== a <*> x++prop_scale_herm_apply k (HermMatrixMV u h (a :: M) x) =+ case u of+ Lower -> (k *> hermL h) <*> x ~== (k *> (a <*> x))+ Upper -> (k *> hermU h) <*> x ~== (k *> (a <*> x))++prop_herm_herm_apply (HermMatrixMV u h (a :: M) x) =+ case u of+ Lower -> hermU (herm h) <*> x ~== a <*> x+ Upper -> hermL (herm h) <*> x ~== a <*> x++prop_herm_compose (HermMatrixMM u h (a :: M) b) =+ case u of+ Lower -> hermL h <**> b ~== a <**> b+ Upper -> hermU h <**> b ~== a <**> b++prop_scale_herm_compose k (HermMatrixMM u h (a :: M) b) =+ case u of+ Lower -> (k *> hermL h) <**> b ~== (k *> (a <**> b))+ Upper -> (k *> hermU h) <**> b ~== (k *> (a <**> b))++prop_herm_herm_compose (HermMatrixMM u h (a :: M) b) =+ case u of+ Lower -> hermU (herm h) <**> b ~== a <**> b+ Upper -> hermL (herm h) <**> b ~== a <**> b+++properties =+ [ + ("herm apply" , pDet prop_herm_apply)+ , ("scale herm apply" , pDet prop_scale_herm_apply)+ , ("herm herm apply" , pDet prop_herm_herm_apply)++ , ("herm compose" , pDet prop_herm_compose)+ , ("scale herm compose" , pDet prop_scale_herm_compose)+ , ("herm herm compose" , pDet prop_herm_herm_compose)+ + ]+++main = do+ args <- getArgs+ n <- case args of+ (a:_) -> readIO a+ _ -> return 1+ main' n++main' n = do+ putStrLn $ "Running tests for " ++ field+ pRun n 400 properties
+ tests/Matrix.hs view
@@ -0,0 +1,334 @@+{-# OPTIONS -fglasgow-exts -fno-excess-precision -cpp #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++import Debug.Trace ( trace )++import qualified Data.Array as Array+import Data.Ix ( inRange, range )+import Data.List ( nub, sortBy )+import Data.Ord ( comparing )+import System.Environment ( getArgs )+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 qualified BLAS.Elem as E+import Data.Complex ( Complex(..) )++import Data.AEq+import Numeric.IEEE++import Test.QuickCheck.Complex+import Test.QuickCheck.Vector hiding ( Assocs )+import Test.QuickCheck.Vector.Dense hiding ( Pair )+import Test.QuickCheck.Matrix+import Test.QuickCheck.Matrix.Dense++import Debug.Trace+ ++#ifdef COMPLEX+field = "Complex Double"+type E = Complex Double+#else+field = "Double"+type E = Double+#endif++type V = Vector Int E+type M = Matrix (Int,Int) E++instance (Arbitrary e, RealFloat e) => Arbitrary (Complex e) where+ arbitrary = arbitrary >>= \(TestComplex x) -> return x+ coarbitrary = coarbitrary . TestComplex++instance (Arbitrary e, BLAS1 e) => Arbitrary (DVector Imm n e) where+ arbitrary = arbitrary >>= \(TestVector x) -> return x+ coarbitrary = coarbitrary . TestVector++instance (Arbitrary e, BLAS1 e) => Arbitrary (DMatrix Imm (m,n) e) where+ arbitrary = arbitrary >>= \(TestMatrix x) -> return x+ coarbitrary = coarbitrary . TestMatrix+++assocsEq :: (BLAS1 e, AEq e) => Matrix (m,n) e -> [((Int,Int), e)] -> Bool+assocsEq x ijes =+ let ijs = fst $ unzip ijes+ in filter (\(ij,e) -> ij `elem` ijs) (sortBy (comparing fst) $ assocs x) === sortBy (comparing fst) ijes+ && (all (==0) $ map snd $ filter (\(ij,e) -> not $ ij `elem` ijs) $ assocs x)++prop_matrix_shape (Assocs mn ijes) =+ shape (matrix mn ijes :: M) == mn+prop_matrix_assocs (Assocs mn ijes) =+ (matrix mn ijes :: M) `assocsEq` ijes++prop_listMatrix_shape (IndexPair mn) es =+ shape (listMatrix mn es :: M) == mn+prop_listMatrix_assocs (IndexPair (m,n)) es =+ let es' = repeat es+ in assocs (listMatrix (m,n) es' :: M) === zip [(i,j) | j <- range (0,n-1), i <- range (0,m-1)] es'++prop_zero_shape (IndexPair mn) =+ shape (zero mn :: M) == mn+prop_zero_elems (IndexPair (m,n)) =+ elems (zero (m,n) :: M) == replicate (m*n) 0++prop_constant_shape (IndexPair mn) (e :: E) =+ shape (constant mn e :: M) == mn+prop_constant_elems (IndexPair (m,n)) (e :: E) =+ elems (constant (m,n) e :: M) == replicate (m*n) e++prop_identity_shape (IndexPair mn) =+ shape (identity mn :: M) == mn+prop_identity_diag (IndexPair (m,n)) =+ 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+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++prop_replace_elems (a :: M) (Assocs _ ijes) =+ let ijes' = filter (\((i,j),_) -> i < numRows a && j < numCols a) ijes+ a' = a // ijes'+ mn = (numRows a - 1, numCols a - 1)+ in and $ zipWith (\(ij1,e1) (ij2,e2) -> (ij1 == ij2) && (e1 === e2))+ (sortBy (comparing fst) $ assocs a')+ (Array.assocs $ (Array.//) (Array.array ((0,0),mn) $ assocs a) ijes')+++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)]+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)]++prop_shape (a :: M) = + shape a == (numRows a, numCols a)+prop_size (a :: M) =+ size a == numRows a * numCols a+prop_bounds (a :: M) =+ bounds a == ((0,0), (numRows a - 1, numCols a - 1))+ +prop_at (MatAt (a :: M) (i,j)) =+ let ij = (i,j)+ k = i + j * numRows a+ in (a!ij) === ((elems a) !! k)+ +prop_row_dim (MatAt (a :: M) (i,_)) =+ V.dim (row a i) == numCols a+prop_col_dim (MatAt (a :: M) (_,j)) =+ V.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)+prop_cols_dims (a :: M) =+ map (V.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)+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))+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)++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)+prop_herm_cols (a :: M) = + cols (herm a) === map (V.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)+prop_diag_herm2 (MatAt (a :: M) (_,k)) =+ diag a k === V.conj (diag (herm a) (-k))++prop_fromRow_shape (x :: V) =+ shape (fromRow x :: M) == (1,V.dim x)+prop_fromRow_elems (x :: V) =+ elems (fromRow x :: M) === V.elems x++prop_fromCol_shape (x :: V) =+ shape (fromCol x :: M) == (V.dim x,1)+prop_fromCol_elems (x :: V) =+ elems (fromCol x :: M) === V.elems x+++prop_apply_basis (MatAt (a :: M) (_,j)) =+ a <*> (V.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)+prop_apply_scale k (MultMV (a :: M) x) =+ a <*> (V.scale k x) ~== V.scale k (a <*> x)+prop_apply_linear (MultMVPair (a :: M) x y) =+ a <*> (x + y) ~== a <*> x + a <*> y++prop_compose_id_right (a :: M) =+ let n = numCols a+ in a <**> (identity (n,n) :: M) ~== a+prop_compose_id_left (a :: M) =+ let m = numRows a+ in (identity (m,m) :: M) <**> a ~== a+prop_compose_scale_left (MultMM (a:: M) b) k =+ a <**> (k *> b) ~== k *> (a <**> b) +prop_compose_scale_right (MultMM (a:: M) b) k =+ (k *> a) <**> b ~== k *> (a <**> b)+prop_compose_linear (MultMMPair (a :: M) b c) =+ a <**> (b + c) ~== a <**> b + a <**> c+prop_compose_herm (MultMM (a :: M) b) =+ herm b <**> herm a ~== herm (a <**> b)+prop_compose_cols (MultMM (a :: M) b) =+ cols (a <**> b) ~== map (a <*> ) (cols b)++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+prop_invScale k (a :: M) =+ invScale k a ~== a / constant (shape a) k++prop_plus (Pair (a :: M) b) =+ elems (a + b) ~== zipWith (+) (elems a) (elems b)+prop_minus (Pair (a :: M) b) =+ elems (a - b) ~== zipWith (-) (elems a) (elems b)+prop_times (Pair (a :: M) b) =+ elems (a * b) ~== zipWith (*) (elems a) (elems b)+prop_divide (Pair (a :: M) b) =+ elems (a / b) ~== zipWith (/) (elems a) (elems b)++prop_negate (a :: M) =+ negate a ~== scale (-1) a+ +prop_abs (a :: M) =+ elems (abs a) ~== map abs (elems a)+prop_signum (a :: M) =+ elems (signum a) === map signum (elems a)+prop_recip (a :: M) =+ elems (recip a) ~== (map recip $ elems a)++properties =+ [ ("shape of matrix" , pDet prop_matrix_shape)+ , ("assocs of matrix" , pDet prop_matrix_assocs)+ , ("shape of listMatrix" , pDet prop_listMatrix_shape)+ , ("assocs of listMatrix" , pDet prop_listMatrix_assocs)+ , ("shape of zero" , pDet prop_zero_shape)+ , ("elems of zero" , pDet prop_zero_elems)+ + , ("shape of constant" , pDet prop_constant_shape)+ , ("elems of constant" , pDet prop_constant_elems)+ + , ("shape of identity" , pDet prop_identity_shape)+ , ("diag of identity" , pDet prop_identity_diag)+ , ("row of identity" , pDet prop_identity_row)+ , ("col of identity" , pDet prop_identity_col)+ + , ("elems of replace" , pDet prop_replace_elems)+ + , ("numRows/numCols" , pDet prop_shape)+ , ("size" , pDet prop_size)+ , ("bounds" , pDet prop_bounds)+ , ("at" , pDet prop_at)+ , ("row dim" , pDet prop_row_dim)+ , ("col dim" , pDet prop_col_dim)+ , ("rows length" , pDet prop_rows_len)+ , ("cols length" , pDet prop_cols_len)+ , ("rows dims" , pDet prop_rows_dims)+ , ("cols dims" , pDet prop_cols_dims)++ , ("indices" , pDet prop_indices)+ , ("elems" , pDet prop_elems)+ , ("assocs" , pDet prop_assocs)+ + , ("shape of submatrix" , pDet prop_submatrix_shape)+ , ("rows of submatrix" , pDet prop_submatrix_rows)+ , ("col of submatrix" , pDet prop_submatrix_cols)+ + , ("elems of scale" , pDet prop_scale_elems)+ , ("elem of herm" , pDet prop_herm_elem)+ , ("herm/scale" , pDet prop_herm_scale)+ + , ("shape . herm" , pDet prop_herm_shape)+ , ("rows . herm" , pDet prop_herm_rows)+ , ("cols . herm" , pDet prop_herm_cols)+ + , ("herm . herm == id" , pDet prop_herm_herm)+ + , ("subdiag . herm" , pDet prop_diag_herm1)+ , ("superdiag . herm" , pDet prop_diag_herm2)+ + , ("shape . fromRow" , pDet prop_fromRow_shape)+ , ("elems . fromRow" , pDet prop_fromRow_elems)+ , ("shape . fromCol" , pDet prop_fromCol_shape)+ , ("elems . fromCol" , pDet prop_fromCol_elems)+ + , ("apply basis" , pDet prop_apply_basis)+ , ("apply herm basis" , pDet prop_apply_herm_basis)+ , ("apply scale" , pDet prop_apply_scale)+ , ("apply linear" , pDet prop_apply_linear)+ + , ("compose id left" , pDet prop_compose_id_left)+ , ("compose id right" , pDet prop_compose_id_right)+ , ("compose scale left" , pDet prop_compose_scale_left)+ , ("compose scale right" , pDet prop_compose_scale_right)+ , ("compose linear" , pDet prop_compose_linear)+ , ("compose herm" , pDet prop_compose_herm)+ , ("compose cols" , pDet prop_compose_cols)+ + , ("shift" , pDet prop_shift)+ , ("scale" , pDet prop_scale)+ , ("invScale" , pDet prop_invScale)+ + , ("plus" , pDet prop_plus)+ , ("minus" , pDet prop_minus)+ , ("times" , pDet prop_times)+ , ("divide" , pDet prop_divide)+ + , ("negate" , pDet prop_negate)+ , ("abs" , pDet prop_abs)+ , ("signum" , pDet prop_signum)+ , ("recip" , pDet prop_recip)++ ]+++main = do+ args <- getArgs+ n <- case args of+ (a:_) -> readIO a+ _ -> return 1+ main' n++main' n = do+ putStrLn $ "Running tests for " ++ field+ pRun n 400 properties
+ tests/TriMatrix.hs view
@@ -0,0 +1,157 @@+{-# OPTIONS -fglasgow-exts -fno-excess-precision -cpp #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--++import System.Environment ( getArgs )+import Test.QuickCheck.Parallel hiding ( vector )+import qualified Test.QuickCheck as QC++import Data.Complex ( Complex(..) )++import qualified BLAS.Elem as E+import Data.Vector.Dense+import Data.Matrix.Dense+import Data.Matrix.Tri.Dense ++import Data.AEq+import Numeric.IEEE++import Test.QuickCheck.Complex+import Test.QuickCheck.Matrix.Tri.Dense++isUndefR x = isNaN x || isInfinite x+isUndefC (x :+ y) = isUndefR x || isUndefR y+ +#ifdef COMPLEX+field = "Complex Double"+type E = Complex Double+isUndef = isUndefC+#else+field = "Double"+type E = Double+isUndef = isUndefR+#endif ++type V = Vector Int E+type M = Matrix (Int,Int) E+type TM = Tri (Matrix) (Int,Int) E++instance (Arbitrary e, RealFloat e) => Arbitrary (Complex e) where+ arbitrary = arbitrary >>= \(TestComplex x) -> return x+ coarbitrary = coarbitrary . TestComplex++prop_tri_apply (TriMatrixMV u d (a :: M) x) =+ case (u,d) of+ (Lower, NonUnit) -> (lower a) <*> x ~== a <*> x+ (Lower, Unit ) -> (lowerU a) <*> x ~== a <*> x + x+ (Upper, NonUnit) -> (upper a) <*> x ~== a <*> x+ (Upper, Unit ) -> (upperU a) <*> x ~== a <*> x + x++prop_herm_tri_apply (TriMatrixMV u d (a :: M) x) =+ case (u,d) of+ (Lower, NonUnit) -> (herm $ lower a) <*> x ~== herm a <*> x+ (Lower, Unit ) -> (herm $ lowerU a) <*> x ~== herm a <*> x + x+ (Upper, NonUnit) -> (herm $ upper a) <*> x ~== herm a <*> x+ (Upper, Unit ) -> (herm $ upperU a) <*> x ~== herm a <*> x + x++prop_scale_tri_apply k (TriMatrixMV u d (a :: M) x) =+ case (u,d) of+ (Lower, NonUnit) -> (k *> lower a) <*> x ~== (k *> a) <*> x+ (Lower, Unit ) -> (k *> lowerU a) <*> x ~== (k *> a) <*> x + (k *> x)+ (Upper, NonUnit) -> (k *> upper a) <*> x ~== (k *> a) <*> x+ (Upper, Unit ) -> (k *> upperU a) <*> x ~== (k *> a) <*> x + (k *> x)++prop_scale_herm_tri_apply k (TriMatrixMV u d (a :: M) x) =+ case (u,d) of+ (Lower, NonUnit) -> (k *> (herm $ lower a)) <*> x ~== (k *> herm a) <*> x+ (Lower, Unit ) -> (k *> (herm $ lowerU a)) <*> x ~== (k *> herm a) <*> x + (k *> x)+ (Upper, NonUnit) -> (k *> (herm $ upper a)) <*> x ~== (k *> herm a) <*> x+ (Upper, Unit ) -> (k *> (herm $ upperU a)) <*> x ~== (k *> herm a) <*> x + (k *> x)++prop_herm_scale_tri_apply k (TriMatrixMV u d (a :: M) x) =+ case (u,d) of+ (Lower, NonUnit) -> (herm $ k *> lower a) <*> x ~== (herm $ k *> a) <*> x+ (Lower, Unit ) -> (herm $ k *> lowerU a) <*> x ~== (herm $ k *> a) <*> x + ((E.conj k) *> x)+ (Upper, NonUnit) -> (herm $ k *> upper a) <*> x ~== (herm $ k *> a) <*> x+ (Upper, Unit ) -> (herm $ k *> upperU a) <*> x ~== (herm $ k *> a) <*> x + ((E.conj k) *> x)+++prop_tri_compose (TriMatrixMM u d (a :: M) b) =+ case (u,d) of+ (Lower, NonUnit) -> (lower a) <**> b ~== a <**> b+ (Lower, Unit ) -> (lowerU a) <**> b ~== a <**> b + b+ (Upper, NonUnit) -> (upper a) <**> b ~== a <**> b+ (Upper, Unit ) -> (upperU a) <**> b ~== a <**> b + b++prop_herm_tri_compose (TriMatrixMM u d (a :: M) b) =+ case (u,d) of+ (Lower, NonUnit) -> (herm $ lower a) <**> b ~== herm a <**> b+ (Lower, Unit ) -> (herm $ lowerU a) <**> b ~== herm a <**> b + b+ (Upper, NonUnit) -> (herm $ upper a) <**> b ~== herm a <**> b+ (Upper, Unit ) -> (herm $ upperU a) <**> b ~== herm a <**> b + b++prop_scale_tri_compose k (TriMatrixMM u d (a :: M) b) =+ case (u,d) of+ (Lower, NonUnit) -> (k *> lower a) <**> b ~== (k *> a) <**> b+ (Lower, Unit ) -> (k *> lowerU a) <**> b ~== (k *> a) <**> b + (k *> b)+ (Upper, NonUnit) -> (k *> upper a) <**> b ~== (k *> a) <**> b+ (Upper, Unit ) -> (k *> upperU a) <**> b ~== (k *> a) <**> b + (k *> b)++prop_scale_herm_tri_compose k (TriMatrixMM u d (a :: M) b) =+ case (u,d) of+ (Lower, NonUnit) -> (k *> (herm $ lower a)) <**> b ~== (k *> herm a) <**> b+ (Lower, Unit ) -> (k *> (herm $ lowerU a)) <**> b ~== (k *> herm a) <**> b + (k *> b)+ (Upper, NonUnit) -> (k *> (herm $ upper a)) <**> b ~== (k *> herm a) <**> b+ (Upper, Unit ) -> (k *> (herm $ upperU a)) <**> b ~== (k *> herm a) <**> b + (k *> b)++prop_herm_scale_tri_compose k (TriMatrixMM u d (a :: M) b) =+ case (u,d) of+ (Lower, NonUnit) -> (herm $ k *> lower a) <**> b ~== (herm $ k *> a) <**> b+ (Lower, Unit ) -> (herm $ k *> lowerU a) <**> b ~== (herm $ k *> a) <**> b + ((E.conj k) *> b)+ (Upper, NonUnit) -> (herm $ k *> upper a) <**> b ~== (herm $ k *> a) <**> b+ (Upper, Unit ) -> (herm $ k *> upperU a) <**> b ~== (herm $ k *> a) <**> b + ((E.conj k) *> b)+++prop_tri_solve (TriMatrixSV (t :: TM) y) =+ let x = t <\> y+ in t <*> x ~== y || (any isUndef $ elems x)++prop_tri_invCompose (TriMatrixSM (t :: TM) b) =+ let a = t <\\> b+ in t <**> a ~== b || (any isUndef $ elems a)+++properties =+ [ ("tri apply" , pDet prop_tri_apply)+ , ("scale tri apply" , pDet prop_scale_tri_apply)+ , ("herm tri apply" , pDet prop_herm_tri_apply)+ , ("scale herm tri apply" , pDet prop_scale_herm_tri_apply)+ , ("herm scale tri apply" , pDet prop_herm_scale_tri_apply)+ + , ("tri compose" , pDet prop_tri_compose)+ , ("scale tri compose" , pDet prop_scale_tri_compose)+ , ("herm tri compose" , pDet prop_herm_tri_compose)+ , ("scale herm tri compose", pDet prop_scale_herm_tri_compose)+ , ("herm scale tri compose", pDet prop_herm_scale_tri_compose)+ + , ("tri solve" , pDet prop_tri_solve)+ , ("tri invCompose" , pDet prop_tri_invCompose)+ + ]+++main = do+ args <- getArgs+ n <- case args of+ (a:_) -> readIO a+ _ -> return 1+ main' n++main' n = do+ putStrLn $ "Running tests for " ++ field+ pRun n 400 properties
+ tests/Vector.hs view
@@ -0,0 +1,256 @@+{-# OPTIONS -fglasgow-exts -cpp #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--+++import qualified Data.Array as Array+import Data.Ix ( inRange )+import Data.List ( nub, sortBy )+import Data.Ord ( comparing )+import System.Environment ( getArgs )+import Test.QuickCheck.Parallel hiding ( vector )+import qualified Test.QuickCheck as QC++import Data.Vector.Dense+import BLAS.Elem ( BLAS1 )+import qualified BLAS.Elem as E+import Data.Complex ( Complex(..) )++import Data.AEq+import Numeric.IEEE++import Test.QuickCheck.Complex+import Test.QuickCheck.Vector+import Test.QuickCheck.Vector.Dense+++import Debug.Trace++#ifdef COMPLEX+field = "Complex Double"+type E = Complex Double+#else+field = "Double"+type E = Double+#endif++type V = Vector Index E++instance (Arbitrary e, RealFloat e) => Arbitrary (Complex e) where+ arbitrary = arbitrary >>= \(TestComplex x) -> return x+ coarbitrary = coarbitrary . TestComplex++instance (Arbitrary e, BLAS1 e) => Arbitrary (Vector n e) where+ arbitrary = arbitrary >>= \(TestVector x) -> return x+ coarbitrary = coarbitrary . TestVector++assocsEq x ies =+ let is = fst $ unzip ies+ in filter (\(i,e) -> i `elem` is) (assocs x) === sortBy (comparing fst) ies+ && (all (==0) $ map snd $ filter (\(i,e) -> not $ i `elem` is) $ assocs x)++prop_vector_dim (Assocs n ies) =+ dim (vector n ies :: V) == n+prop_vector_assocs (Assocs n ies) =+ (vector n ies :: V) `assocsEq` ies++prop_listVector_dim es =+ let n = length es+ in dim (listVector n es :: V) == n+prop_listVector_assocs es =+ let n = length es+ in (listVector n es :: V) `assocsEq` (zip [0..] es)++prop_zero_dim (Index n) =+ dim (zero n :: V) == n+prop_zero_elems (Index n) =+ elems (zero n :: V) == replicate n 0+ +prop_constant_dim (Index n) e =+ dim (constant n e :: V) == n+prop_constant_elems (Index n) (e :: E) =+ elems (constant n e :: V) === replicate n e++prop_basis_dim (Basis n i) =+ dim (basis n i :: V) == n+prop_basis_elems (Basis n i) =+ elems (basis n i :: V) == (replicate i 0) ++ [1] ++ (replicate (n-i-1) 0)++prop_replace_elems (x :: V) (Assocs _ ies) =+ let ies' = filter (\(i,e) -> i < dim x) ies+ x' = x // ies'+ n = dim x+ in and $ zipWith (\(i1,e1) (i2,e2) -> (i1 == i2) && (e1 === e2))+ (assocs x')+ (Array.assocs $ (Array.//) (Array.array (0,n-1) $ assocs x) ies')++prop_subvector_dim (SubVector _ (x :: V) o n) =+ dim (subvector x o n) == n+prop_subvector_elems (SubVector _ (x :: V) o n) =+ elems (subvector x o n) === (take n $ drop o $ elems x)+ +prop_subvectorWithStride_dim (SubVector s (x :: V) o n) =+ dim (subvectorWithStride s x o n) == n+ +prop_subvectorWithStride_elems (SubVector s (x :: V) o n) =+ elems (subvectorWithStride s x o n) === + (map snd $ filter (\(i,_) -> (i - o >= 0) + && ((i - o) `mod` s == 0) + && ((i - o) `div` s <= n))+ (assocs x))++prop_dim (x :: V) = + dim x == length (elems x)+prop_bounds (x :: V) =+ bounds x == (0, dim x - 1)+prop_at (x :: V) (Index i) = + i < dim x ==> (x!i) === ((elems x) !! i)+++prop_indices (x :: V) =+ indices x == [0..(dim x - 1)]+prop_assocs (x :: V) = + assocs x === zip (indices x) (elems x)++prop_scale_elems k (x :: V) =+ (elems $ k *> x) ~== (map (k*) $ elems x)+prop_conj_elems (x :: V) =+ and $ zipWith (===) (elems (conj x)) (map (E.conj) (elems x))+prop_conj_scale k (x :: V) =+ conj (k *> x) === (E.conj k *> (conj x))++prop_to_from_list es =+ toList (fromList es :: V) === es++prop_sumAbs (x :: V) =+ sumAbs x ~== (sum $ map E.norm1 $ elems x)+prop_norm2 (x :: V) =+ norm2 x ~== (sqrt $ sum $ map (^2) $ map E.norm $ elems x)+prop_whichMaxAbs1 (x :: V) =+ (dim x > 0) && all (not . isNaN) (map E.norm1 $ elems x) ==>+ let (i,e) = whichMaxAbs x+ in x ! i === e+prop_whichMaxAbs2 (x :: V) =+ (dim x > 0) && all (not . isNaN) (map E.norm1 $ elems x) ==>+ let a = E.norm1 $ snd $ whichMaxAbs x+ in all (<= a) (map E.norm1 $ elems x)+ +prop_dot_self (x :: V) =+ (sqrt $ x <.> x) ~== (E.fromReal $ norm2 x)+prop_dot_conj (Pair (x :: V) y) =+ (x <.> y) ~== (E.conj $ y <.> x)+prop_dot_scale1 k (Pair (x :: V) y) =+ (x <.> (k *> y)) ~== k * (x <.> y)+prop_dot_scale2 k (Pair (x :: V) y) =+ ((k *> x) <.> y) ~== (E.conj k) * (x <.> y)+prop_dot_linear1 (Triple (x :: V) y z) =+ (x <.> (y + z)) ~== (x <.> y + x <.> z)+prop_dot_linear2 (Triple (x :: V) y z) =+ ((x + y) <.> z) ~== (x <.> z + y <.> z)++prop_shift k (x :: V) =+ shift k x ~== x + constant (dim x) k+prop_scale k (x :: V) =+ scale k x ~== x * constant (dim x) k+prop_invScale k (x :: V) =+ invScale k x ~== x / constant (dim x) k++prop_plus (Pair (x :: V) y) =+ elems (x + y) ~== zipWith (+) (elems x) (elems y)+prop_minus (Pair (x :: V) y) =+ elems (x - y) ~== zipWith (-) (elems x) (elems y)+prop_times (Pair (x :: V) y) =+ elems (x * y) ~== zipWith (*) (elems x) (elems y)+prop_divide (Pair (x :: V) y) =+ elems (x / y) ~== zipWith (/) (elems x) (elems y)++prop_negate (x :: V) =+ negate x ~== (-1) *> x+prop_abs (x :: V) =+ elems (abs x) ~== map abs (elems x)+prop_signum (x :: V) =+ elems (signum x) === map signum (elems x)+prop_recip (x :: V) =+ elems (recip x) ~== (map recip $ elems x)++properties =+ [ ("dim of vector" , pDet prop_vector_dim)+ , ("assocs of vector" , pDet prop_vector_assocs)+ + , ("dim of listVector" , pDet prop_listVector_dim)+ , ("assocs of listVector" , pDet prop_listVector_assocs)+ + , ("dim of zero" , pDet prop_zero_dim)+ , ("elems of zero" , pDet prop_zero_elems)+ + , ("dim of constant" , pDet prop_constant_dim)+ , ("elems of constant" , pDet prop_constant_elems)++ , ("dim of basis" , pDet prop_basis_dim)+ , ("elems of basis" , pDet prop_basis_elems)++ , ("dim of subvector" , pDet prop_subvector_dim)+ , ("elems of subvector" , pDet prop_subvector_elems)+ + , ("dim of subvectorWithStride" + , pDet prop_subvectorWithStride_dim)+ , ("elems of subvectorWithStride" + , pDet prop_subvectorWithStride_elems)+ + , ("elems of replace" , pDet prop_replace_elems)+ , ("elems of scale" , pDet prop_scale_elems)+ , ("elems of conj" , pDet prop_conj_elems)+ , ("conj/scale" , pDet prop_conj_scale)+ + , ("dim" , pDet prop_dim)+ , ("bounds" , pDet prop_bounds)+ , ("at" , pDet prop_at)++ , ("indices" , pDet prop_indices)+ , ("assocs" , pDet prop_assocs)+ + , ("to/from list" , pDet prop_to_from_list)+ + , ("sumAbs" , pDet prop_sumAbs)+ , ("norm2" , pDet prop_norm2)+ , ("whichMaxAbs1" , pDet prop_whichMaxAbs1)+ , ("whichMaxAbs2" , pDet prop_whichMaxAbs2)+ + , ("dot self" , pDet prop_dot_self)+ , ("dot conj" , pDet prop_dot_conj)+ , ("dot scale1" , pDet prop_dot_scale1)+ , ("dot scale2" , pDet prop_dot_scale2)+ , ("dot linear1" , pDet prop_dot_linear1)+ , ("dot linear2" , pDet prop_dot_linear2)+ + , ("shift" , pDet prop_shift)+ , ("scale" , pDet prop_scale)+ , ("invScale" , pDet prop_invScale)+ + , ("plus" , pDet prop_plus)+ , ("minus" , pDet prop_minus)+ , ("times" , pDet prop_times)+ , ("divide" , pDet prop_divide)+ + , ("negate" , pDet prop_negate)+ , ("abs" , pDet prop_abs)+ , ("signum" , pDet prop_signum)+ , ("recip" , pDet prop_recip)+ ]+++main = do+ args <- getArgs+ n <- case args of+ (a:_) -> readIO a+ _ -> return 1+ main' n++main' n = do+ putStrLn $ "Running tests for " ++ field+ pRun n 500 properties