diff --git a/BLAS/Access.hs b/BLAS/Access.hs new file mode 100644 --- /dev/null +++ b/BLAS/Access.hs @@ -0,0 +1,20 @@ +{-# LANGUAGE EmptyDataDecls #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Access +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- +-- Types for access control. +module BLAS.Access ( + Mut, + Imm + ) where + +-- | Tag for mutable types. +data Mut + +-- | Tag for immutable types. +data Imm diff --git a/BLAS/C.hs b/BLAS/C.hs new file mode 100644 --- /dev/null +++ b/BLAS/C.hs @@ -0,0 +1,20 @@ +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.C +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/C/Double.hs b/BLAS/C/Double.hs new file mode 100644 --- /dev/null +++ b/BLAS/C/Double.hs @@ -0,0 +1,112 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.C.Double +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 () + diff --git a/BLAS/C/Level1.hs b/BLAS/C/Level1.hs new file mode 100644 --- /dev/null +++ b/BLAS/C/Level1.hs @@ -0,0 +1,125 @@ +{-# LANGUAGE FlexibleInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.C.Level1 +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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'' + diff --git a/BLAS/C/Level2.hs b/BLAS/C/Level2.hs new file mode 100644 --- /dev/null +++ b/BLAS/C/Level2.hs @@ -0,0 +1,82 @@ +{-# LANGUAGE FlexibleInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.C.Level2 +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/C/Level3.hs b/BLAS/C/Level3.hs new file mode 100644 --- /dev/null +++ b/BLAS/C/Level3.hs @@ -0,0 +1,82 @@ +{-# LANGUAGE FlexibleInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.C.Level3 +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/C/Types.hs b/BLAS/C/Types.hs new file mode 100644 --- /dev/null +++ b/BLAS/C/Types.hs @@ -0,0 +1,87 @@ +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.C.Types +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/C/Zomplex.hs b/BLAS/C/Zomplex.hs new file mode 100644 --- /dev/null +++ b/BLAS/C/Zomplex.hs @@ -0,0 +1,122 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.C.Zomplex +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 () + diff --git a/BLAS/Elem.hs b/BLAS/Elem.hs new file mode 100644 --- /dev/null +++ b/BLAS/Elem.hs @@ -0,0 +1,20 @@ +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Elem +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 ) diff --git a/BLAS/Elem/Base.hs b/BLAS/Elem/Base.hs new file mode 100644 --- /dev/null +++ b/BLAS/Elem/Base.hs @@ -0,0 +1,51 @@ +{-# LANGUAGE FlexibleInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Elem.Base +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 + diff --git a/BLAS/Internal.hs b/BLAS/Internal.hs new file mode 100644 --- /dev/null +++ b/BLAS/Internal.hs @@ -0,0 +1,198 @@ +{-# LANGUAGE CPP, ForeignFunctionInterface #-} +{-# OPTIONS_GHC -fglasgow-exts #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Internal +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 () + + diff --git a/BLAS/Matrix.hs b/BLAS/Matrix.hs new file mode 100644 --- /dev/null +++ b/BLAS/Matrix.hs @@ -0,0 +1,20 @@ +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Matrix +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/Matrix/Base.hs b/BLAS/Matrix/Base.hs new file mode 100644 --- /dev/null +++ b/BLAS/Matrix/Base.hs @@ -0,0 +1,26 @@ +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Matrix.Base +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/Matrix/Immutable.hs b/BLAS/Matrix/Immutable.hs new file mode 100644 --- /dev/null +++ b/BLAS/Matrix/Immutable.hs @@ -0,0 +1,34 @@ +{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Matrix.Immutable +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 + diff --git a/BLAS/Matrix/ReadOnly.hs b/BLAS/Matrix/ReadOnly.hs new file mode 100644 --- /dev/null +++ b/BLAS/Matrix/ReadOnly.hs @@ -0,0 +1,30 @@ +{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Matrix.ReadOnly +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/Matrix/Solve.hs b/BLAS/Matrix/Solve.hs new file mode 100644 --- /dev/null +++ b/BLAS/Matrix/Solve.hs @@ -0,0 +1,16 @@ +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Matrix.Solve +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/Matrix/Solve/Immutable.hs b/BLAS/Matrix/Solve/Immutable.hs new file mode 100644 --- /dev/null +++ b/BLAS/Matrix/Solve/Immutable.hs @@ -0,0 +1,28 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Matrix.Solve.Immutable +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 + diff --git a/BLAS/Matrix/Solve/ReadOnly.hs b/BLAS/Matrix/Solve/ReadOnly.hs new file mode 100644 --- /dev/null +++ b/BLAS/Matrix/Solve/ReadOnly.hs @@ -0,0 +1,25 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Matrix.Solve.ReadOnly +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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) diff --git a/BLAS/Tensor.hs b/BLAS/Tensor.hs new file mode 100644 --- /dev/null +++ b/BLAS/Tensor.hs @@ -0,0 +1,24 @@ +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Tensor +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/Tensor/Base.hs b/BLAS/Tensor/Base.hs new file mode 100644 --- /dev/null +++ b/BLAS/Tensor/Base.hs @@ -0,0 +1,25 @@ +{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Tensor.Base +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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) diff --git a/BLAS/Tensor/Dense.hs b/BLAS/Tensor/Dense.hs new file mode 100644 --- /dev/null +++ b/BLAS/Tensor/Dense.hs @@ -0,0 +1,16 @@ +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Tensor.Dense +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/Tensor/Dense/Immutable.hs b/BLAS/Tensor/Dense/Immutable.hs new file mode 100644 --- /dev/null +++ b/BLAS/Tensor/Dense/Immutable.hs @@ -0,0 +1,23 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Tensor.Dense.Immutable +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/Tensor/Dense/ReadOnly.hs b/BLAS/Tensor/Dense/ReadOnly.hs new file mode 100644 --- /dev/null +++ b/BLAS/Tensor/Dense/ReadOnly.hs @@ -0,0 +1,25 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Tensor.Dense.ReadOnly +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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) + diff --git a/BLAS/Tensor/Immutable.hs b/BLAS/Tensor/Immutable.hs new file mode 100644 --- /dev/null +++ b/BLAS/Tensor/Immutable.hs @@ -0,0 +1,70 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Tensor.Immutable +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/Tensor/Mutable.hs b/BLAS/Tensor/Mutable.hs new file mode 100644 --- /dev/null +++ b/BLAS/Tensor/Mutable.hs @@ -0,0 +1,76 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Tensor.Mutable +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 + diff --git a/BLAS/Tensor/ReadOnly.hs b/BLAS/Tensor/ReadOnly.hs new file mode 100644 --- /dev/null +++ b/BLAS/Tensor/ReadOnly.hs @@ -0,0 +1,62 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Tensor.ReadOnly +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/Tensor/Scalable.hs b/BLAS/Tensor/Scalable.hs new file mode 100644 --- /dev/null +++ b/BLAS/Tensor/Scalable.hs @@ -0,0 +1,20 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Tensor.Scalable +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/Types.hs b/BLAS/Types.hs new file mode 100644 --- /dev/null +++ b/BLAS/Types.hs @@ -0,0 +1,53 @@ +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Types +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/BLAS/Vector.hs b/BLAS/Vector.hs new file mode 100644 --- /dev/null +++ b/BLAS/Vector.hs @@ -0,0 +1,25 @@ +----------------------------------------------------------------------------- +-- | +-- Module : BLAS.Vector +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 + diff --git a/Data/Matrix/Dense.hs b/Data/Matrix/Dense.hs new file mode 100644 --- /dev/null +++ b/Data/Matrix/Dense.hs @@ -0,0 +1,168 @@ +{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Data.Matrix.Dense +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 + diff --git a/Data/Matrix/Dense/IO.hs b/Data/Matrix/Dense/IO.hs new file mode 100644 --- /dev/null +++ b/Data/Matrix/Dense/IO.hs @@ -0,0 +1,89 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Data.Matrix.Dense.IO +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/Data/Matrix/Dense/Internal.hs b/Data/Matrix/Dense/Internal.hs new file mode 100644 --- /dev/null +++ b/Data/Matrix/Dense/Internal.hs @@ -0,0 +1,536 @@ +{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Data.Matrix.Dense.Internal +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 (~==) + diff --git a/Data/Matrix/Dense/Operations.hs b/Data/Matrix/Dense/Operations.hs new file mode 100644 --- /dev/null +++ b/Data/Matrix/Dense/Operations.hs @@ -0,0 +1,347 @@ +{-# OPTIONS_GHC -fglasgow-exts #-} +----------------------------------------------------------------------------- +-- | +-- Module : Data.Matrix.Dense.Operations +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 + #-} + diff --git a/Data/Matrix/Herm.hs b/Data/Matrix/Herm.hs new file mode 100644 --- /dev/null +++ b/Data/Matrix/Herm.hs @@ -0,0 +1,64 @@ +{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, FlexibleInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : Data.Matrix.Herm +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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" diff --git a/Data/Matrix/Herm/Dense.hs b/Data/Matrix/Herm/Dense.hs new file mode 100644 --- /dev/null +++ b/Data/Matrix/Herm/Dense.hs @@ -0,0 +1,114 @@ +{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : Data.Matrix.Herm.Dense +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/Data/Matrix/Tri.hs b/Data/Matrix/Tri.hs new file mode 100644 --- /dev/null +++ b/Data/Matrix/Tri.hs @@ -0,0 +1,73 @@ +{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, FlexibleInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : Data.Matrix.Tri +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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" + diff --git a/Data/Matrix/Tri/Dense.hs b/Data/Matrix/Tri/Dense.hs new file mode 100644 --- /dev/null +++ b/Data/Matrix/Tri/Dense.hs @@ -0,0 +1,171 @@ +{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : Data.Matrix.Tri.Dense +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 + diff --git a/Data/Vector/Dense.hs b/Data/Vector/Dense.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Dense.hs @@ -0,0 +1,135 @@ +{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Dense +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 + diff --git a/Data/Vector/Dense/IO.hs b/Data/Vector/Dense/IO.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Dense/IO.hs @@ -0,0 +1,64 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Dense.Mutable +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/Data/Vector/Dense/Internal.hs b/Data/Vector/Dense/Internal.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Dense/Internal.hs @@ -0,0 +1,422 @@ +{-# LANGUAGE BangPatterns, FlexibleInstances, MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Dense.Internal +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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) + diff --git a/Data/Vector/Dense/Operations.hs b/Data/Vector/Dense/Operations.hs new file mode 100644 --- /dev/null +++ b/Data/Vector/Dense/Operations.hs @@ -0,0 +1,373 @@ +{-# OPTIONS_GHC -fglasgow-exts #-} +----------------------------------------------------------------------------- +-- | +-- Module : Data.Vector.Dense.Operations +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 + #-} + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) Patrick Perry 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. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -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 } +> diff --git a/Test/QuickCheck/Complex.hs b/Test/QuickCheck/Complex.hs new file mode 100644 --- /dev/null +++ b/Test/QuickCheck/Complex.hs @@ -0,0 +1,24 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Test.QuickCheck.Complex +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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) diff --git a/Test/QuickCheck/Matrix.hs b/Test/QuickCheck/Matrix.hs new file mode 100644 --- /dev/null +++ b/Test/QuickCheck/Matrix.hs @@ -0,0 +1,41 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Test.QuickCheck.Matrix +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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) diff --git a/Test/QuickCheck/Matrix/Dense.hs b/Test/QuickCheck/Matrix/Dense.hs new file mode 100644 --- /dev/null +++ b/Test/QuickCheck/Matrix/Dense.hs @@ -0,0 +1,157 @@ +{-# LANGUAGE FlexibleInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : Test.QuickCheck.Matrix.Dense +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/Test/QuickCheck/Matrix/Herm/Dense.hs b/Test/QuickCheck/Matrix/Herm/Dense.hs new file mode 100644 --- /dev/null +++ b/Test/QuickCheck/Matrix/Herm/Dense.hs @@ -0,0 +1,72 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Test.QuickCheck.Matrix.Herm.Dense +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 + diff --git a/Test/QuickCheck/Matrix/Tri/Dense.hs b/Test/QuickCheck/Matrix/Tri/Dense.hs new file mode 100644 --- /dev/null +++ b/Test/QuickCheck/Matrix/Tri/Dense.hs @@ -0,0 +1,132 @@ +{-# LANGUAGE FlexibleInstances #-} +----------------------------------------------------------------------------- +-- | +-- Module : Test.QuickCheck.Matrix.Tri.Dense +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 + + diff --git a/Test/QuickCheck/Vector.hs b/Test/QuickCheck/Vector.hs new file mode 100644 --- /dev/null +++ b/Test/QuickCheck/Vector.hs @@ -0,0 +1,50 @@ + +----------------------------------------------------------------------------- +-- | +-- Module : Test.QuickCheck.Vector +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 ( +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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) + diff --git a/blas.cabal b/blas.cabal new file mode 100644 --- /dev/null +++ b/blas.cabal @@ -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: + + . +category: Math +license: BSD3 +license-file: LICENSE +copyright: (c) 2008. Patrick Perry +author: Patrick Perry +maintainer: Patrick Perry +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 + diff --git a/tests/HermMatrix.hs b/tests/HermMatrix.hs new file mode 100644 --- /dev/null +++ b/tests/HermMatrix.hs @@ -0,0 +1,96 @@ +{-# OPTIONS -fglasgow-exts -fno-excess-precision -cpp #-} +----------------------------------------------------------------------------- +-- | +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/tests/Matrix.hs b/tests/Matrix.hs new file mode 100644 --- /dev/null +++ b/tests/Matrix.hs @@ -0,0 +1,334 @@ +{-# OPTIONS -fglasgow-exts -fno-excess-precision -cpp #-} +----------------------------------------------------------------------------- +-- | +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/tests/TriMatrix.hs b/tests/TriMatrix.hs new file mode 100644 --- /dev/null +++ b/tests/TriMatrix.hs @@ -0,0 +1,157 @@ +{-# OPTIONS -fglasgow-exts -fno-excess-precision -cpp #-} +----------------------------------------------------------------------------- +-- | +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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 diff --git a/tests/Vector.hs b/tests/Vector.hs new file mode 100644 --- /dev/null +++ b/tests/Vector.hs @@ -0,0 +1,256 @@ +{-# OPTIONS -fglasgow-exts -cpp #-} +----------------------------------------------------------------------------- +-- | +-- Copyright : Copyright (c) 2008, Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- 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