diff --git a/hblas.cabal b/hblas.cabal --- a/hblas.cabal +++ b/hblas.cabal @@ -1,31 +1,31 @@ --- Initial hblas.cabal generated by cabal init. For further +-- Initial hblas.cabal generated by cabal init. For further -- documentation, see http://haskell.org/cabal/users-guide/ -- The name of the package. name: hblas --- The package version. See the Haskell package versioning policy (PVP) +-- The package version. See the Haskell package versioning policy (PVP) -- for standards guiding when and how versions should be incremented. -- http://www.haskell.org/haskellwiki/Package_versioning_policy -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change -version: 0.3.0.1 +version: 0.3.1.0 -- A short (one-line) description of the package. synopsis: Human friendly BLAS and Lapack bindings for Haskell. -- A longer description of the package. -description: - User friendly, simple bindings to BLAS and Lapack. Easy to extend and use. +description: + User friendly, simple bindings to BLAS and Lapack. Easy to extend and use. . HBLAS assumes you have BLAS and LAPACK installed. Users on OSX already - have a BLAS and LAPACK installed for them. Users on LINUX and similar + have a BLAS and LAPACK installed for them. Users on LINUX and similar platforms will need to do the system equivalent of @sudo apt-get install libblas liblapack@. . - HBLAS is not intended to be a library for end user array operations, - but rather an easy adaptor for calling BLAS and LAPACK routines in their - full generality from your end-user haskell linear algebra / array library + HBLAS is not intended to be a library for end user array operations, + but rather an easy adaptor for calling BLAS and LAPACK routines in their + full generality from your end-user haskell linear algebra / array library of choice. . HBLAS roughly mirrors the naming conventions and datatype names of the standard @@ -33,11 +33,11 @@ and moderately comfortable with haskell will feel right at home. . If you are not familiar with the the operations and conventions of BLAS and - Lapack, please read the - , + Lapack, please read the + , you'll learn a lot about what BLAS and LAPACK can do. . - HBLAS doesn't provide every BLAS and LAPACK operation currently, but + HBLAS doesn't provide every BLAS and LAPACK operation currently, but is designed so that adding support for new operations is quite easy! @@ -54,41 +54,45 @@ -- The package author(s). author: Carter Tazio Schonwald --- An email address to which users can send suggestions, bug reports, and +-- An email address to which users can send suggestions, bug reports, and -- patches. maintainer: carter at wellposed dot com -- A copyright notice. --- copyright: +-- copyright: cabal-version: >=1.10 category: Math build-type: Custom -extra-source-files: +extra-source-files: readme.md, - changelog.md + changelog.md source-repository head type: git location: http://github.com/wellposed/hblas.git - --- Extra files to be distributed with the package, such as examples or a + +-- Extra files to be distributed with the package, such as examples or a -- README. flag OpenBLAS - default: False - manual: True + default: False + manual: True -- Constraint on the version of Cabal needed to build this package. -- flag lib-Werror -- default: False -- manual: True --- flag lib-Wall +-- flag lib-Wall -- default: False -- manual: True +flag CBLAS + default: False + manual: True + description: Use on platforms when libcblas isn't part of libblas library @@ -97,30 +101,33 @@ Numerical.HBLAS.MatrixTypes Numerical.HBLAS.UtilsFFI Numerical.HBLAS.BLAS - Numerical.HBLAS.BLAS.Internal + Numerical.HBLAS.BLAS.Internal Numerical.HBLAS.Lapack - Numerical.HBLAS.Lapack.FFI + Numerical.HBLAS.Lapack.FFI Numerical.HBLAS.Constants -- Modules included in this library but not exported. - -- other-modules: - + -- other-modules: + -- LANGUAGE extensions used by modules in this package. - -- other-extensions: + -- other-extensions: --- if flag(lib-Werror) --- ghc-options: -Werror --- if flag(lib-Wall) - --- ghc-options: -Wall + --- ghc-options: -Wall -- -Wall no wall for now -- Other library packages from which modules are imported. - build-depends: base >=4.6 && <4.8, storable-complex >= 0.2.0 && < 0.3.0, - vector , primitive >= 0.5 && < 0.6 + build-depends: + base >= 4.6 && < 4.8, + storable-complex >= 0.2.0 && < 0.3.0, + primitive >= 0.5 && < 0.6, + vector if flag(OpenBLAS) extra-libraries: openblas pthread if flag(OpenBLAS)&& os(OSX) - extra-lib-dirs: /usr/local/lib + extra-lib-dirs: /usr/local/lib if os(OSX)&&!flag(OpenBLAS) frameworks: Accelerate @@ -132,35 +139,32 @@ if !os(windows)&& !os(OSX) && !flag(OpenBLAS) extra-libraries: blas lapack - + if flag(CBLAS) + extra-libraries: cblas + + -- Directories containing source files. hs-source-dirs: src - + -- Base language which the package is written in. default-language: Haskell2010 ghc-options: -Wall -fno-warn-name-shadowing -Test-suite unit-testsuite +Test-suite unit-testsuite default-language: Haskell2010 type: exitcode-stdio-1.0 - build-depends: base >=4.6 && < 5 - --,hspec >= 1.9 && < 1.10 - ,tasty >= 0.8 && < 0.9 - ,tasty-hunit >= 0.8 && < 0.9 - ,HUnit >= 1.2.5 && < 1.3 - ,vector >= 0.5 && < 0.12 - --,QuickCheck >= 2.7 && < 2.8 - ,hblas - hs-source-dirs: tests + build-depends: + hblas + , base >= 4.6 && < 5 + , tasty >= 0.8 && < 0.9 + , tasty-hunit >= 0.8 && < 0.9 + , HUnit >= 1.2.5 && < 1.3 + , vector >= 0.5 && < 0.12 + --,QuickCheck >= 2.7 && < 2.8 + --,hspec >= 1.9 && < 1.10 + hs-source-dirs: tests main-is: MainUnit.hs - - - - - - - diff --git a/readme.md b/readme.md --- a/readme.md +++ b/readme.md @@ -1,46 +1,97 @@ -[![Wellposed](http://www.wellposed.com/mini.png)](http://www.wellposed.com)™ +[![Wellposed](http://www.wellposed.com/mini.png)](http://www.wellposed.com)™ # About hblas -hblas is an open source component of the [Wellposed](http://www.wellposed.com)® mathematical software suite. - -Members of the numerical haskell open source community can be found on irc at `#numerical-haskell` on freenode, -and via the [numericalhaskell mailing list](https://groups.google.com/forum/#!forum/numericalhaskell). +hblas is an open source component of the [Wellposed](http://www.wellposed.com)® mathematical software suite. +Members of the numerical haskell open source community can be found on irc at `#numerical-haskell` on freenode, +and via the [numericalhaskell mailing list](https://groups.google.com/forum/#!forum/numericalhaskell). [![Build Status](https://secure.travis-ci.org/wellposed/hblas.png?branch=master)](http://travis-ci.org/wellposed/hblas) - -hblas is a self contained full (well, not quite yet) BLAS and LAPACK binding that provides the -full BLAS and LAPACKE APIs in a simple, unopinionated, Haskell wrapper. +hblas is a self contained full (well, not quite yet) BLAS and LAPACK binding that provides the +full BLAS and LAPACKE APIs in a simple, unopinionated, Haskell wrapper. -This library is *NOT* meant to be used by end users, it is designed to be +This library is *NOT* meant to be used by end users, it is designed to be an unopinionated, simple, portable, easy to install BLAS/LAPACK substrate for higher level numerical computing libraries to build upon. Morever, this library is strictly a wrapper, and simply makes using the functionality of BLAS and LAPACK more accessible. This library is *NOT* meant to be used a standalone array library (except in desperation), -but rather should be used by a higher level numerical array library to provide -high performance linear algebra routines. +but rather should be used by a higher level numerical array library to provide +high performance linear algebra routines. +## Install -## how to install -By default, hblas will assume you have BLAS and LAPACK built and installed, +By default, hblas will assume you have BLAS and LAPACK built and installed. -* On OS X systems, things will just work. -* On linux and bsd systems, the equivalent of +### OSX + +On OS X systems, things will just work. + +```bash +$ cabal install ``` -sudo apt-get install libblas liblapack + +### Linux + +On linux and bsd systems, you will need to manually install the BLAS and LAPACK libraries beforehand. + +```bash +$ sudo apt-get install libblas liblapack +$ cabal install ``` -is all you should have to do before hand -## getting involved -patches, bug reports, tests, and other contributions welcome. +## Testing -Want to add a new routine, check out the ones listed in the [lapack section](http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm) of the Intel MKL manual to get some human +To run the test suite execute: + +```bash +$ cabal test +``` + +## Linking +If you get an error like `undefined reference to `cblas_sdsdot'` when building or running an HBLAS program, +you might be on a system that builds BLAS and CBLAS separately, such as Arch Linux. + +In which case, be sure to install CBLAS, and also `cabal configure`/ `cabal install` HBLAS +with the `-fCBLAS` flag. + +## Usage + +API is subject to change. + +```haskell +import Foreign.Storable +import Numerical.HBLAS.BLAS +import Numerical.HBLAS.MatrixTypes + +-- Generate the constant mutable square matrix of the given type and dimensions. +constMatrix :: Storable a => Int -> a -> IO (IODenseMatrix Row a) +constMatrix n k = generateMutableDenseMatrix SRow (n,n) (const k) + +example_dgemm :: IO () +example_dgemm = do + left <- constMatrix 2 (2 :: Double) + right <- constMatrix 2 (3 :: Double) + out <- constMatrix 2 (0 :: Double) + + dgemm NoTranspose NoTranspose 1.0 1.0 left right res + + resulting <- mutableVectorToList $ _bufferDenMutMat out + print resulting +``` + +## Getting Involved + +Patches, bug reports, tests, and other contributions welcome. + +If you want to add a new routine, check out the ones listed in the [lapack section](http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm) of the Intel MKL manual to get some human readable documentation. +## Commercial Support -# I have > 32bit size arrays, help! +*I have > 32bit size arrays, help!* + Congrats, you have ``big compute on big data'', contact [Carter](http://www.wellposed.com) -and we'll try to help you out. +and we'll try to help you out. diff --git a/src/Numerical/HBLAS/BLAS.hs b/src/Numerical/HBLAS/BLAS.hs --- a/src/Numerical/HBLAS/BLAS.hs +++ b/src/Numerical/HBLAS/BLAS.hs @@ -1,7 +1,7 @@ {-# LANGUAGE BangPatterns , RankNTypes, GADTs, DataKinds #-} {- | The 'Numerical.HBLAS.BLAS' module provides a fully general -yet type safe BLAS API. +yet type safe BLAS API. When in doubt about the semantics of an operation, consult your system's BLAS api documentation, or just read the documentation @@ -10,15 +10,15 @@ A few basic notes about how to invoke BLAS routines. -Many BLAS operations take one or more arguments of type 'Transpose'. +Many BLAS operations take one or more arguments of type 'Transpose'. 'Tranpose' has the following different constructors, which tell BLAS routines what transformation to implicitly apply to an input matrix @mat@ with dimension @n x m@. -* 'NoTranspose' leaves the matrix @mat@ as is. +* 'NoTranspose' leaves the matrix @mat@ as is. * 'Transpose' treats the @mat@ as being implicitly transposed, with dimension @m x n@. Entry @mat(i,j)@ being treated as actually being the entry - @mat(j,i)@. For Real matrices this is also the matrix adjoint operation. + @mat(j,i)@. For Real matrices this is also the matrix adjoint operation. ie @Tranpose(mat)(i,j)=mat(j,i)@ * 'ConjNoTranspose' will implicitly conjugate @mat@, which is a no op for Real ('Float' or 'Double') matrices, but for @@ -35,7 +35,7 @@ * @'sgemm trLeft trRight alpha beta left right result'@, where @trLeft@ and @trRight@ are values of type 'Transpose' that respectively act on the matrices @left@ and @right@. -* the generalized matrix computation thusly formed can be viewed as being +* the generalized matrix computation thusly formed can be viewed as being @result = alpha * trLeft(left) * trRight(right) + beta * result@ @@ -43,9 +43,9 @@ being vectors rather than matrices. -the *trsv operations solve for @x@ in the equation @A x = y@ given @A@ and @y@. -The 'MatUpLo' argument determines if the matrix should be treated as upper or -lower triangular and 'MatDiag' determines if the triangular solver should treat +the *trsv operations solve for @x@ in the equation @A x = y@ given @A@ and @y@. +The 'MatUpLo' argument determines if the matrix should be treated as upper or +lower triangular and 'MatDiag' determines if the triangular solver should treat the diagonal of the matrix as being all 1's or not. A general pattern of invocation would be @'strsv' matuplo tranposeMatA matdiag matrixA xVector@. A key detail to note is that the input vector is ALSO the result vector, @@ -55,32 +55,32 @@ module Numerical.HBLAS.BLAS( GemvFun - ,GemmFun - ,TrsvFun - + ,GemmFun + ,TrsvFun + ,dgemm ,sgemm ,cgemm ,zgemm - ,sgemv + ,sgemv ,dgemv - ,cgemv - ,zgemv + ,cgemv + ,zgemv ,strsv ,dtrsv ,ctrsv - ,ztrsv - ) where + ,ztrsv + ) where -import Numerical.HBLAS.UtilsFFI -import Numerical.HBLAS.BLAS.FFI -import Numerical.HBLAS.BLAS.Internal +import Numerical.HBLAS.UtilsFFI +import Numerical.HBLAS.BLAS.FFI +import Numerical.HBLAS.BLAS.Internal import Numerical.HBLAS.MatrixTypes import Control.Monad.Primitive -import Data.Complex +import Data.Complex @@ -93,24 +93,24 @@ type TrsvFun el orient s m = MatUpLo -> Transpose -> MatDiag - -> MDenseMatrix s orient el -> MDenseVector s Direct el -> m () + -> MDenseMatrix s orient el -> MDenseVector s Direct el -> m () -sgemm :: PrimMonad m=> GemmFun Float orient (PrimState m) m -sgemm = gemmAbstraction "sgemm" cblas_sgemm_safe cblas_sgemm_unsafe (\x f -> f x ) - +sgemm :: PrimMonad m=> GemmFun Float orient (PrimState m) m +sgemm = gemmAbstraction "sgemm" cblas_sgemm_safe cblas_sgemm_unsafe (\x f -> f x ) -dgemm :: PrimMonad m=> GemmFun Double orient (PrimState m) m -dgemm = gemmAbstraction "dgemm" cblas_dgemm_safe cblas_dgemm_unsafe (\x f -> f x ) - -cgemm :: PrimMonad m=> GemmFun (Complex Float) orient (PrimState m) m -cgemm = gemmAbstraction "cgemm" cblas_cgemm_safe cblas_cgemm_unsafe withRStorable_ +dgemm :: PrimMonad m=> GemmFun Double orient (PrimState m) m +dgemm = gemmAbstraction "dgemm" cblas_dgemm_safe cblas_dgemm_unsafe (\x f -> f x ) -zgemm :: PrimMonad m=> GemmFun (Complex Double) orient (PrimState m) m -zgemm = gemmAbstraction "zgemm" cblas_zgemm_safe cblas_zgemm_unsafe withRStorable_ + +cgemm :: PrimMonad m=> GemmFun (Complex Float) orient (PrimState m) m +cgemm = gemmAbstraction "cgemm" cblas_cgemm_safe cblas_cgemm_unsafe withRStorable_ + +zgemm :: PrimMonad m=> GemmFun (Complex Double) orient (PrimState m) m +zgemm = gemmAbstraction "zgemm" cblas_zgemm_safe cblas_zgemm_unsafe withRStorable_ sgemv :: PrimMonad m => GemvFun Float orient (PrimState m) m diff --git a/src/Numerical/HBLAS/BLAS/FFI.hs b/src/Numerical/HBLAS/BLAS/FFI.hs --- a/src/Numerical/HBLAS/BLAS/FFI.hs +++ b/src/Numerical/HBLAS/BLAS/FFI.hs @@ -2,19 +2,19 @@ {-# LANGUAGE ForeignFunctionInterface, GeneralizedNewtypeDeriving #-} -module Numerical.HBLAS.BLAS.FFI where +module Numerical.HBLAS.BLAS.FFI where import Foreign.Ptr import Foreign() import Foreign.C.Types -import Data.Complex +import Data.Complex -- /*Set the number of threads on runtime.*/ --foreign import ccall unsafe "openblas_set_num_threads" openblas_set_num_threads_unsafe :: CInt -> IO () --foreign import ccall unsafe "goto_set_num_threads" goto_set_num_threads_unsafe :: CInt -> IO () -{- | For All of the BlAS FFI operations, +{- | For All of the BlAS FFI operations, -} @@ -29,20 +29,20 @@ -} -newtype CBLAS_INDEX = CBIndex CSize +newtype CBLAS_INDEX = CBIndex CSize deriving (Eq,Show) newtype CBLAS_ORDERT = CBOInt CInt deriving (Eq,Show) -data BLASOrder = BLASRowMajor | BLASColMajor +data BLASOrder = BLASRowMajor | BLASColMajor deriving (Eq,Show) encodeOrder :: BLASOrder -> CBLAS_ORDERT encodeOrder BLASRowMajor = CBOInt 101 -encodeOrder BLASColMajor = CBOInt 102 +encodeOrder BLASColMajor = CBOInt 102 newtype CBLAS_TRANSPOSET = CBLAS_TransposeT{ unCBLAS_TransposeT :: CInt } deriving (Eq, Show) -data BLAS_Transpose = BlasNoTranspose | BlasTranspose | BlasConjTranspose | BlasConjNoTranspose +data BLAS_Transpose = BlasNoTranspose | BlasTranspose | BlasConjTranspose | BlasConjNoTranspose encodeTranspose :: BLAS_Transpose -> CBLAS_TRANSPOSET encodeTranspose BlasNoTranspose = CBLAS_TransposeT 111 @@ -50,31 +50,31 @@ encodeTranspose BlasConjTranspose = CBLAS_TransposeT 113 encodeTranspose BlasConjNoTranspose = CBLAS_TransposeT 114 -newtype CBLAS_UPLOT = CBlasUPLO CInt +newtype CBLAS_UPLOT = CBlasUPLO CInt deriving (Eq,Show) data BLASUplo = BUpper | BLower deriving (Eq,Show) -encodeUPLO :: BLASUplo -> CBLAS_UPLOT -encodeUPLO BUpper = CBlasUPLO 121 +encodeUPLO :: BLASUplo -> CBLAS_UPLOT +encodeUPLO BUpper = CBlasUPLO 121 encodeUPLO BLower = CBlasUPLO 122 -newtype CBLAS_DIAGT = CBLAS_DiagT CUChar +newtype CBLAS_DIAGT = CBLAS_DiagT CUChar deriving (Show,Eq) -data BlasDiag = BlasNonUnit | BlasUnit +data BlasDiag = BlasNonUnit | BlasUnit deriving (Eq,Show ) - -encodeDiag :: BlasDiag -> CBLAS_DIAGT + +encodeDiag :: BlasDiag -> CBLAS_DIAGT encodeDiag BlasNonUnit = CBLAS_DiagT 131 encodeDiag BlasUnit = CBLAS_DiagT 132 -newtype CBLAS_SIDET = CBLAS_SideT { unCBLAS_SideT :: CUChar } +newtype CBLAS_SIDET = CBLAS_SideT { unCBLAS_SideT :: CUChar } deriving (Eq, Show) -data BlasSide = BlasLeft | BlasRight +data BlasSide = BlasLeft | BlasRight deriving (Eq,Show) -encodeSide :: BlasSide -> CBLAS_SIDET +encodeSide :: BlasSide -> CBLAS_SIDET encodeSide BlasLeft = CBLAS_SideT 141 encodeSide BlasRight = CBLAS_SideT 142 @@ -86,7 +86,7 @@ --typedef enum CBLAS_SIDE {CblasLeft=141, CblasRight=142} CBLAS_SIDE; --dot products -foreign import ccall unsafe "cblas_sdsdot" cblas_sdsdot_unsafe :: CInt -> Float -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO Float +foreign import ccall unsafe "cblas_sdsdot" cblas_sdsdot_unsafe :: CInt -> Float -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO Float foreign import ccall unsafe "cblas_dsdot" cblas_dsdot_unsafe :: CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO Double foreign import ccall unsafe "cblas_sdot" cblas_sdot_unsafe :: CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO Float foreign import ccall unsafe "cblas_ddot" cblas_ddot_unsafe :: CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> IO Double @@ -109,20 +109,20 @@ ---- not sure what to do for these complex +--- not sure what to do for these complex --void cblas_cdotu_sub( CInt n, Float *x, CInt incx, Float *y, CInt incy, openblas_complex_Float *ret); --void cblas_cdotc_sub( CInt n, Float *x, CInt incx, Float *y, CInt incy, openblas_complex_Float *ret); --void cblas_zdotu_sub( CInt n, Double *x, CInt incx, Double *y, CInt incy, openblas_complex_Double *ret); --void cblas_zdotc_sub( CInt n, Double *x, CInt incx, Double *y, CInt incy, openblas_complex_Double *ret); ---- absolute value -foreign import ccall unsafe "cblas_sasum" cblas_sasum_unsafe:: +foreign import ccall unsafe "cblas_sasum" cblas_sasum_unsafe:: CInt -> Ptr Float -> CInt -> IO Float -foreign import ccall unsafe "cblas_dasum" cblas_dasum_unsafe :: +foreign import ccall unsafe "cblas_dasum" cblas_dasum_unsafe :: CInt -> Ptr Double -> CInt -> IO Double -foreign import ccall unsafe "cblas_scasum" cblas_casum_unsafe :: +foreign import ccall unsafe "cblas_scasum" cblas_casum_unsafe :: CInt -> Ptr (Complex Float)-> CInt -> IO Float -foreign import ccall unsafe "cblas_dzasum" cblas_zasum_unsafe :: +foreign import ccall unsafe "cblas_dzasum" cblas_zasum_unsafe :: CInt -> Ptr (Complex Double) -> CInt -> IO Double --Float cblas_sasum ( CInt n, Float *x, CInt incx); --Double cblas_dasum ( CInt n, Double *x, CInt incx); @@ -130,13 +130,13 @@ --Double cblas_dzasum( CInt n, Double *x, CInt incx); -foreign import ccall unsafe "cblas_snrm2" cblas_snrm2_unsafe :: +foreign import ccall unsafe "cblas_snrm2" cblas_snrm2_unsafe :: CInt -> Ptr Float -> CInt -> IO Float -foreign import ccall unsafe "cblas_dnrm2" cblas_dnrm2_unsafe :: +foreign import ccall unsafe "cblas_dnrm2" cblas_dnrm2_unsafe :: CInt -> Ptr Double -> CInt -> IO Double -foreign import ccall unsafe "cblas_scnrm2" cblas_scnrm2_unsafe :: +foreign import ccall unsafe "cblas_scnrm2" cblas_scnrm2_unsafe :: CInt -> Ptr (Complex Float)-> CInt -> IO Float -foreign import ccall unsafe "cblas_dznrm2" cblas_dznrm2_unsafe :: +foreign import ccall unsafe "cblas_dznrm2" cblas_dznrm2_unsafe :: CInt -> Ptr (Complex Double) -> CInt -> IO Double @@ -147,14 +147,14 @@ -foreign import ccall unsafe "cblas_isamax" cblas_isamax_unsafe :: +foreign import ccall unsafe "cblas_isamax" cblas_isamax_unsafe :: CInt -> Ptr Float -> CInt -> IO CInt -foreign import ccall unsafe "cblas_idamax" cblas_idamax_unsafe :: +foreign import ccall unsafe "cblas_idamax" cblas_idamax_unsafe :: CInt -> Ptr Float -> CInt -> IO CInt -foreign import ccall unsafe "cblas_icamax" cblas_icamax_unsafe :: +foreign import ccall unsafe "cblas_icamax" cblas_icamax_unsafe :: CInt -> Ptr (Complex Float) -> CInt -> IO CInt -foreign import ccall unsafe "cblas_izamax" cblas_izamax_unsafe :: - CInt -> Ptr (Complex Double) -> CInt -> IO CInt +foreign import ccall unsafe "cblas_izamax" cblas_izamax_unsafe :: + CInt -> Ptr (Complex Double) -> CInt -> IO CInt --CBLAS_INDEX cblas_isamax( CInt n, Float *x, CInt incx); --CBLAS_INDEX cblas_idamax( CInt n, Double *x, CInt incx); --CBLAS_INDEX cblas_icamax( CInt n, Float *x, CInt incx); @@ -170,14 +170,14 @@ -foreign import ccall unsafe "cblas_scopy" cblas_scopy_unsafe :: - CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO () -foreign import ccall unsafe "cblas_dcopy" cblas_dcopy_unsafe :: +foreign import ccall unsafe "cblas_scopy" cblas_scopy_unsafe :: + CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO () +foreign import ccall unsafe "cblas_dcopy" cblas_dcopy_unsafe :: CInt -> Ptr Double-> CInt -> Ptr Double -> CInt -> IO () -foreign import ccall unsafe "cblas_ccopy" cblas_ccopy_unsafe :: +foreign import ccall unsafe "cblas_ccopy" cblas_ccopy_unsafe :: CInt -> Ptr (Complex Float) -> CInt -> Ptr (Complex Float) -> CInt -> IO () -foreign import ccall unsafe "cblas_zcopy" cblas_zcopy_unsafe :: - CInt -> Ptr (Complex Double)-> CInt -> Ptr (Complex Double) -> CInt -> IO () +foreign import ccall unsafe "cblas_zcopy" cblas_zcopy_unsafe :: + CInt -> Ptr (Complex Double)-> CInt -> Ptr (Complex Double) -> CInt -> IO () --void cblas_scopy( CInt n, Float *x, CInt incx, Float *y, CInt incy); --void cblas_dcopy( CInt n, Double *x, CInt incx, Double *y, CInt incy); @@ -186,14 +186,14 @@ -foreign import ccall unsafe "cblas_sswap" cblas_sswap_unsafe :: - CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO () -foreign import ccall unsafe "cblas_dswap" cblas_dswap_unsafe :: +foreign import ccall unsafe "cblas_sswap" cblas_sswap_unsafe :: + CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO () +foreign import ccall unsafe "cblas_dswap" cblas_dswap_unsafe :: CInt -> Ptr Double-> CInt -> Ptr Double -> CInt -> IO () -foreign import ccall unsafe "cblas_cswap" cblas_cswap_unsafe :: +foreign import ccall unsafe "cblas_cswap" cblas_cswap_unsafe :: CInt -> Ptr (Complex Float) -> CInt -> Ptr (Complex Float) -> CInt -> IO () -foreign import ccall unsafe "cblas_zswap" cblas_zswap_unsafe :: - CInt -> Ptr (Complex Double)-> CInt -> Ptr (Complex Double) -> CInt -> IO () +foreign import ccall unsafe "cblas_zswap" cblas_zswap_unsafe :: + CInt -> Ptr (Complex Double)-> CInt -> Ptr (Complex Double) -> CInt -> IO () --void cblas_sswap( CInt n, Float *x, CInt incx, Float *y, CInt incy); --void cblas_dswap( CInt n, Double *x, CInt incx, Double *y, CInt incy); @@ -367,7 +367,7 @@ -- CInt N, CInt K, Double *A, CInt lda, Double *X, CInt incX); ---------------- ---- | solves Ax=v where A is k+1 banded triangular matrix, and x and +--- | solves Ax=v where A is k+1 banded triangular matrix, and x and ---------------- --void cblas_stbsv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag, -- CInt N, CInt K, Float *A, CInt lda, Float *X, CInt incX); @@ -391,7 +391,7 @@ -- CInt N, Double *Ap, Double *X, CInt incX); -------------------------------------------------- ---- | solve Ax=v where A is a nxn packed triangular matrix, v vector input, writes the solution into x. +--- | solve Ax=v where A is a nxn packed triangular matrix, v vector input, writes the solution into x. -------------------------------------------------- --void cblas_stpsv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag, -- CInt N, Float *Ap, Float *X, CInt incX); @@ -408,7 +408,7 @@ --------------------------------- type SymvFunFFI el = CBLAS_ORDERT -> CBLAS_UPLOT -> CInt -> el -> Ptr el -> CInt -> - Ptr el -> CInt -> el -> Ptr el -> CInt -> IO () + Ptr el -> CInt -> el -> Ptr el -> CInt -> IO () --void cblas_ssymv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Float alpha, Float *A, -- CInt lda, Float *X, CInt incX, Float beta, Float *Y, CInt incY); @@ -470,9 +470,9 @@ -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- ------------------------------- | BLAS LEVEL 3 ROUTINES +------------------------------ | BLAS LEVEL 3 ROUTINES -------------------------------------------------------------------------------- ------------------------ | Level 3 ops are faster than Levels 1 or 2 +----------------------- | Level 3 ops are faster than Levels 1 or 2 -------------------------------------------------------------------------------- @@ -488,35 +488,35 @@ -- | Matrix mult for general dense matrices type GemmFunFFI scale el = CBLAS_ORDERT -> CBLAS_TRANSPOSET -> CBLAS_TRANSPOSET-> - CInt -> CInt -> CInt -> {- scal A * B -} scale -> {- Matrix A-} Ptr el -> CInt -> {- B -} Ptr el -> CInt-> + CInt -> CInt -> CInt -> {- scal A * B -} scale -> {- Matrix A-} Ptr el -> CInt -> {- B -} Ptr el -> CInt-> scale -> {- C -} Ptr el -> CInt -> IO () {- C := alpha*op( A )*op( B ) + beta*C , -} -- matrix mult! -foreign import ccall unsafe "cblas_sgemm" +foreign import ccall unsafe "cblas_sgemm" cblas_sgemm_unsafe :: GemmFunFFI Float Float -foreign import ccall unsafe "cblas_dgemm" +foreign import ccall unsafe "cblas_dgemm" cblas_dgemm_unsafe :: GemmFunFFI Double Double -foreign import ccall unsafe "cblas_cgemm" +foreign import ccall unsafe "cblas_cgemm" cblas_cgemm_unsafe :: GemmFunFFI (Ptr(Complex Float)) (Complex Float) -foreign import ccall unsafe "cblas_zgemm" +foreign import ccall unsafe "cblas_zgemm" cblas_zgemm_unsafe :: GemmFunFFI (Ptr (Complex Double)) (Complex Double) -- safe ffi variant for large inputs -foreign import ccall "cblas_sgemm" +foreign import ccall "cblas_sgemm" cblas_sgemm_safe :: GemmFunFFI Float Float -foreign import ccall "cblas_dgemm" +foreign import ccall "cblas_dgemm" cblas_dgemm_safe :: GemmFunFFI Double Double -foreign import ccall "cblas_cgemm" +foreign import ccall "cblas_cgemm" cblas_cgemm_safe :: GemmFunFFI (Ptr(Complex Float)) (Complex Float) -foreign import ccall "cblas_zgemm" +foreign import ccall "cblas_zgemm" cblas_zgemm_safe :: GemmFunFFI (Ptr (Complex Double)) (Complex Double) ----------------------------------------- @@ -536,35 +536,35 @@ type SymmFunFFI scale el = CBLAS_ORDERT -> CBLAS_SIDET -> CBLAS_UPLOT -> CInt->CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt -> scale ->Ptr el -> CInt -> IO () -foreign import ccall unsafe "cblas_ssymm" - cblas_ssymm_unsafe :: SymmFunFFI Float Float +foreign import ccall unsafe "cblas_ssymm" + cblas_ssymm_unsafe :: SymmFunFFI Float Float -foreign import ccall unsafe "cblas_dsymm" - cblas_dsymm_unsafe :: SymmFunFFI Double Double +foreign import ccall unsafe "cblas_dsymm" + cblas_dsymm_unsafe :: SymmFunFFI Double Double -foreign import ccall unsafe "cblas_csymm" +foreign import ccall unsafe "cblas_csymm" cblas_csymm_unsafe :: SymmFunFFI (Ptr (Complex Float )) (Complex Float) -foreign import ccall unsafe "cblas_zsymm" +foreign import ccall unsafe "cblas_zsymm" cblas_zsymm_unsafe :: SymmFunFFI (Ptr (Complex Double)) (Complex Double) --- safe ffi variant, -foreign import ccall "cblas_ssymm" - cblas_ssymm_safe :: SymmFunFFI Float Float +-- safe ffi variant, +foreign import ccall "cblas_ssymm" + cblas_ssymm_safe :: SymmFunFFI Float Float -foreign import ccall "cblas_dsymm" - cblas_dsymm_safe :: SymmFunFFI Double Double +foreign import ccall "cblas_dsymm" + cblas_dsymm_safe :: SymmFunFFI Double Double -foreign import ccall "cblas_csymm" +foreign import ccall "cblas_csymm" cblas_csymm_safe :: SymmFunFFI (Ptr (Complex Float )) (Complex Float) -foreign import ccall "cblas_zsymm" +foreign import ccall "cblas_zsymm" cblas_zsymm_safe :: SymmFunFFI (Ptr (Complex Double)) (Complex Double) - + ----------------------------------- --- | symmetric rank k matrix update, C := alpha*A*A' + beta*C ---- or C = alpha*A'*A + beta*C +--- or C = alpha*A'*A + beta*C ------------------------------------ @@ -579,13 +579,13 @@ type SyrkFunFFI scale el = CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CInt->CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt ->scale -> Ptr el -> CInt -> IO () -foreign import ccall unsafe "cblas_ssyrk" - cblas_ssyrk_unsafe :: SyrkFunFFI Float Float -foreign import ccall unsafe "cblas_dsyrk" +foreign import ccall unsafe "cblas_ssyrk" + cblas_ssyrk_unsafe :: SyrkFunFFI Float Float +foreign import ccall unsafe "cblas_dsyrk" cblas_dsyrk_unsafe :: SyrkFunFFI Double Double -foreign import ccall unsafe "cblas_csyrk" +foreign import ccall unsafe "cblas_csyrk" cblas_csyrk_unsafe :: SyrkFunFFI (Ptr(Complex Float)) (Complex Float) -foreign import ccall unsafe "cblas_zsyrk" +foreign import ccall unsafe "cblas_zsyrk" cblas_zsyrk_unsafe :: SyrkFunFFI (Ptr(Complex Double)) (Complex Double) ---------------------- @@ -604,17 +604,17 @@ --const blasint N, const blasint K, const double *alpha, const double *A, const blasint lda, const double *B, const blasint ldb, const double *beta, double *C, const blasint ldc); type Syr2kFunFFI scale el = CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> - CInt->CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt -> + CInt->CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt -> scale ->Ptr el -> CInt -> IO () -foreign import ccall unsafe "cblas_ssyr2k" - cblas_ssyr2k_unsafe :: Syr2kFunFFI Float Float -foreign import ccall unsafe "cblas_dsyr2k" +foreign import ccall unsafe "cblas_ssyr2k" + cblas_ssyr2k_unsafe :: Syr2kFunFFI Float Float +foreign import ccall unsafe "cblas_dsyr2k" cblas_dsyr2k_unsafe :: Syr2kFunFFI Double Double -foreign import ccall unsafe "cblas_csyr2k" - cblas_csyr2k_unsafe :: Syr2kFunFFI (Ptr (Complex Float)) Float -foreign import ccall unsafe "cblas_zsyr2k" - cblas_zsyr2k_unsafe :: Syr2kFunFFI (Ptr (Complex Double)) Double +foreign import ccall unsafe "cblas_csyr2k" + cblas_csyr2k_unsafe :: Syr2kFunFFI (Ptr (Complex Float)) Float +foreign import ccall unsafe "cblas_zsyr2k" + cblas_zsyr2k_unsafe :: Syr2kFunFFI (Ptr (Complex Double)) Double @@ -626,16 +626,16 @@ -type TrmmFunFFI scale el = CBLAS_ORDERT -> CBLAS_SIDET -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT -> +type TrmmFunFFI scale el = CBLAS_ORDERT -> CBLAS_SIDET -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT -> CInt->CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO () -foreign import ccall unsafe "cblas_strmm" - cblas_strmm_unsafe :: TrmmFunFFI Float Float -foreign import ccall unsafe "cblas_dtrmm" - cblas_dtrmm_unsafe :: TrmmFunFFI Double Double -foreign import ccall unsafe "cblas_ctrmm" - cblas_ctrmm_unsafe :: TrmmFunFFI (Ptr (Complex Float )) (Complex Float) -foreign import ccall unsafe "cblas_ztrmm" - cblas_ztrmm_unsafe :: TrmmFunFFI (Ptr (Complex Double )) (Complex Double) +foreign import ccall unsafe "cblas_strmm" + cblas_strmm_unsafe :: TrmmFunFFI Float Float +foreign import ccall unsafe "cblas_dtrmm" + cblas_dtrmm_unsafe :: TrmmFunFFI Double Double +foreign import ccall unsafe "cblas_ctrmm" + cblas_ctrmm_unsafe :: TrmmFunFFI (Ptr (Complex Float )) (Complex Float) +foreign import ccall unsafe "cblas_ztrmm" + cblas_ztrmm_unsafe :: TrmmFunFFI (Ptr (Complex Double )) (Complex Double) --void cblas_strmm( enum CBLAS_ORDER Order, enum CBLAS_SIDE Side, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, -- enum CBLAS_DIAG Diag, CInt M, CInt N, Float alpha, Float *A, CInt lda, Float *B, CInt ldb); @@ -647,25 +647,25 @@ -- enum CBLAS_DIAG Diag, CInt M, CInt N, Double *alpha, Double *A, CInt lda, Double *B, CInt ldb); ------------------------ --- | triangular solvers +-- | triangular solvers ----------------------- --- ---TRSM solves op(A)*X = alpha*B or X*op(A) = alpha*B +-- +--TRSM solves op(A)*X = alpha*B or X*op(A) = alpha*B --op(A) is one of op(A) = A, or op(A) = A', or op(A) = conjg(A'). --- A is a unit, or non-unit, upper or lower triangular matrix ----- -type TrsmFunFFI scale el = CBLAS_ORDERT -> CBLAS_SIDET -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT -> +-- A is a unit, or non-unit, upper or lower triangular matrix +---- +type TrsmFunFFI scale el = CBLAS_ORDERT -> CBLAS_SIDET -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT -> CInt->CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO () -foreign import ccall unsafe "cblas_strsm" - cblas_strsm_unsafe :: TrmmFunFFI Float Float -foreign import ccall unsafe "cblas_dtrsm" - cblas_dtrsm_unsafe :: TrmmFunFFI Double Double -foreign import ccall unsafe "cblas_ctrsm" - cblas_ctrsm_unsafe :: TrmmFunFFI (Ptr (Complex Float )) (Complex Float) -foreign import ccall unsafe "cblas_ztrsm" - cblas_ztrsm_unsafe :: TrmmFunFFI (Ptr (Complex Double )) (Complex Double) +foreign import ccall unsafe "cblas_strsm" + cblas_strsm_unsafe :: TrmmFunFFI Float Float +foreign import ccall unsafe "cblas_dtrsm" + cblas_dtrsm_unsafe :: TrmmFunFFI Double Double +foreign import ccall unsafe "cblas_ctrsm" + cblas_ctrsm_unsafe :: TrmmFunFFI (Ptr (Complex Float )) (Complex Float) +foreign import ccall unsafe "cblas_ztrsm" + cblas_ztrsm_unsafe :: TrmmFunFFI (Ptr (Complex Double )) (Complex Double) --void cblas_strsm( enum CBLAS_ORDER Order, enum CBLAS_SIDE Side, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, -- enum CBLAS_DIAG Diag, CInt M, CInt N, Float alpha, Float *A, CInt lda, Float *B, CInt ldb); @@ -683,10 +683,10 @@ type HemmFunFFI el = CBLAS_ORDERT -> CBLAS_SIDET -> CBLAS_UPLOT -> CInt->CInt -> Ptr el -> Ptr el -> CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO () -foreign import ccall unsafe "cblas_chemm" - cblas_chemm_unsafe :: HemmFunFFI (Complex Float) -foreign import ccall unsafe "cblas_zhemm" - cblas_zhemm_unsafe :: HemmFunFFI (Complex Double) +foreign import ccall unsafe "cblas_chemm" + cblas_chemm_unsafe :: HemmFunFFI (Complex Float) +foreign import ccall unsafe "cblas_zhemm" + cblas_zhemm_unsafe :: HemmFunFFI (Complex Double) --void cblas_chemm( enum CBLAS_ORDER Order, enum CBLAS_SIDE Side, enum CBLAS_UPLO Uplo, CInt M, CInt N, -- Float *alpha, Float *A, CInt lda, Float *B, CInt ldb, Float *beta, Float *C, CInt ldc); @@ -696,10 +696,10 @@ type HerkFun scale el = CBLAS_ORDERT -> CBLAS_SIDET-> CBLAS_TRANSPOSET -> CInt->CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt ->scale ->Ptr el -> CInt -> IO () -foreign import ccall unsafe "cblas_cherk" - cblas_cherk_unsafe :: HerkFun Float (Complex Float) -foreign import ccall unsafe "cblas_zherk" - cblas_zherk_unsafe :: HerkFun Double (Complex Double) +foreign import ccall unsafe "cblas_cherk" + cblas_cherk_unsafe :: HerkFun Float (Complex Float) +foreign import ccall unsafe "cblas_zherk" + cblas_zherk_unsafe :: HerkFun Double (Complex Double) --void cblas_cherk( enum CBLAS_ORDER Order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE Trans, CInt N, CInt K, -- Float alpha, Float *A, CInt lda, Float beta, Float *C, CInt ldc); --void cblas_zherk( enum CBLAS_ORDER Order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE Trans, CInt N, CInt K, @@ -708,9 +708,9 @@ type Her2kFunFFI scale el = CBLAS_ORDERT -> CBLAS_SIDET -> CBLAS_TRANSPOSET -> CInt->CInt -> Ptr el -> Ptr el -> CInt -> Ptr el -> CInt ->scale ->Ptr el -> CInt -> IO () -foreign import ccall unsafe "cblas_cher2k" - cblas_cher2k_unsafe :: Her2kFunFFI Float (Complex Float) -foreign import ccall unsafe "cblas_zher2k" +foreign import ccall unsafe "cblas_cher2k" + cblas_cher2k_unsafe :: Her2kFunFFI Float (Complex Float) +foreign import ccall unsafe "cblas_zher2k" cblas_zher2k_unsafe :: Her2kFunFFI Double (Complex Double) diff --git a/src/Numerical/HBLAS/BLAS/Internal.hs b/src/Numerical/HBLAS/BLAS/Internal.hs --- a/src/Numerical/HBLAS/BLAS/Internal.hs +++ b/src/Numerical/HBLAS/BLAS/Internal.hs @@ -4,15 +4,15 @@ gemmAbstraction ,gemvAbstraction ,trsvAbstraction - ) where + ) where import Numerical.HBLAS.Constants -import Numerical.HBLAS.UtilsFFI -import Numerical.HBLAS.BLAS.FFI +import Numerical.HBLAS.UtilsFFI +import Numerical.HBLAS.BLAS.FFI import Numerical.HBLAS.MatrixTypes import Control.Monad.Primitive import qualified Data.Vector.Storable.Mutable as SM -import Data.Int +import Data.Int type GemmFun el orient s m = Transpose ->Transpose -> el -> el -> MDenseMatrix s orient el -> MDenseMatrix s orient el -> MDenseMatrix s orient el -> m () @@ -23,37 +23,37 @@ type TrsvFun el orient s m = MatUpLo -> Transpose -> MatDiag - -> MDenseMatrix s orient el -> MDenseVector s Direct el -> m () + -> MDenseMatrix s orient el -> MDenseVector s Direct el -> m () gemmComplexity :: Integral a => a -> a -> a -> Int64 gemmComplexity a b c = fromIntegral a * fromIntegral b *fromIntegral c -- this will be wrong by some constant factor, albeit a small one gemvComplexity :: Integral a => a -> a -> Int64 -gemvComplexity a b = fromIntegral a * fromIntegral b +gemvComplexity a b = fromIntegral a * fromIntegral b -- this covers the ~6 cases for checking the dimensions for GEMM quite nicely isBadGemm :: (Ord a, Num a) => Transpose -> Transpose -> a -> a -> a -> a -> a -> a -> Bool isBadGemm tra trb ax ay bx by cx cy = isBadGemmHelper (cds tra (ax,ay)) (cds trb (bx,by) ) (cx,cy) - where - cds = coordSwapper - isBadGemmHelper !(ax,ay) !(bx,by) !(cx,cy) = (minimum [ax, ay, bx, by, cx ,cy] <= 0) + where + cds = coordSwapper + isBadGemmHelper !(ax,ay) !(bx,by) !(cx,cy) = (minimum [ax, ay, bx, by, cx ,cy] <= 0) || not ( cy == ay && cx == bx && ax == by) coordSwapper :: Transpose -> (a,a)-> (a,a) coordSwapper NoTranspose (a,b) = (a,b) -coordSwapper ConjNoTranspose (a,b) = (a,b) +coordSwapper ConjNoTranspose (a,b) = (a,b) coordSwapper Transpose (a,b) = (b,a) coordSwapper ConjTranspose (a,b) = (b,a) --- / checks if the size of a matrices rows matches input vector size +-- / checks if the size of a matrices rows matches input vector size -- and the column count matchesresult vector size -isBadGemv :: Transpose -> Int -> Int -> Int -> Int -> Bool +isBadGemv :: Transpose -> Int -> Int -> Int -> Int -> Bool isBadGemv tr ax ay bdim cdim = isBadGemvHelper (cds tr (ax,ay)) - where + where cds = coordSwapper - isBadGemvHelper (realX,realY) = + isBadGemvHelper (realX,realY) = minimum [realY,realX,bdim,cdim] <= 0 || not (realX == bdim && realY == cdim ) @@ -63,10 +63,10 @@ encodeFFITranspose :: Transpose -> CBLAS_TRANSPOSET -encodeFFITranspose x= encodeTranspose $ encodeNiceTranspose x +encodeFFITranspose x= encodeTranspose $ encodeNiceTranspose x encodeNiceTranspose :: Transpose -> BLAS_Transpose -encodeNiceTranspose x = case x of +encodeNiceTranspose x = case x of NoTranspose -> BlasNoTranspose Transpose -> BlasTranspose ConjTranspose -> BlasConjTranspose @@ -93,12 +93,12 @@ {- -A key design goal of this ffi is to provide *safe* throughput guarantees +A key design goal of this ffi is to provide *safe* throughput guarantees for a concurrent application built on top of these apis, while evading any overheads for providing such safety. Accordingly, on inputs sizes where the estimated flops count will be more then 1-10 microseconds, safe ffi calls are used. For inputs whose runtime is under that -unsafe ffi calls are used. +unsafe ffi calls are used. -} @@ -106,41 +106,41 @@ ---- | Matrix mult for general dense matrices --type GemmFunFFI scale el = CBLAS_ORDERT -> CBLAS_TRANSPOSET -> CBLAS_TRANSPOSET-> - --CInt -> CInt -> CInt -> {- scal A * B -} scale -> {- Matrix A-} Ptr el -> CInt -> {- B -} Ptr el -> CInt-> + --CInt -> CInt -> CInt -> {- scal A * B -} scale -> {- Matrix A-} Ptr el -> CInt -> {- B -} Ptr el -> CInt-> --scale -> {- C -} Ptr el -> CInt -> IO () --type GemmFun = MutDenseMatrix or el -> MutDenseMatrix or el -> MutDenseMatrix or el -> m () {-# NOINLINE gemmAbstraction #-} -gemmAbstraction:: (SM.Storable el, PrimMonad m) => String -> - GemmFunFFI scale el -> GemmFunFFI scale el -> (el -> (scale -> m ())->m ()) -> forall orient . GemmFun el orient (PrimState m) m -gemmAbstraction gemmName gemmSafeFFI gemmUnsafeFFI constHandler = go - where - shouldCallFast :: Int -> Int -> Int -> Bool +gemmAbstraction:: (SM.Storable el, PrimMonad m) => String -> + GemmFunFFI scale el -> GemmFunFFI scale el -> (el -> (scale -> m ())->m ()) -> forall orient . GemmFun el orient (PrimState m) m +gemmAbstraction gemmName gemmSafeFFI gemmUnsafeFFI constHandler = go + where + shouldCallFast :: Int -> Int -> Int -> Bool shouldCallFast cy cx ax = flopsThreshold >= gemmComplexity cy cx ax - go tra trb alpha beta - (MutableDenseMatrix ornta ax ay astride abuff) - (MutableDenseMatrix _ bx by bstride bbuff) - (MutableDenseMatrix _ cx cy cstride cbuff) + go tra trb alpha beta + (MutableDenseMatrix ornta ax ay astride abuff) + (MutableDenseMatrix _ bx by bstride bbuff) + (MutableDenseMatrix _ cx cy cstride cbuff) | isBadGemm tra trb ax ay bx by cx cy = error $! "bad dimension args to GEMM: ax ay bx by cx cy: " ++ show [ax, ay, bx, by, cx ,cy] - | SM.overlaps abuff cbuff || SM.overlaps bbuff cbuff = - error $ "the read and write inputs for: " ++ gemmName ++ " overlap. This is a programmer error. Please fix." - | otherwise = + | SM.overlaps abuff cbuff || SM.overlaps bbuff cbuff = + error $ "the read and write inputs for: " ++ gemmName ++ " overlap. This is a programmer error. Please fix." + | otherwise = {- FIXME : Add Sharing check that also errors out for now-} - unsafeWithPrim abuff $ \ap -> - unsafeWithPrim bbuff $ \bp -> - unsafeWithPrim cbuff $ \cp -> - constHandler alpha $ \alphaPtr -> - constHandler beta $ \betaPtr -> + unsafeWithPrim abuff $ \ap -> + unsafeWithPrim bbuff $ \bp -> + unsafeWithPrim cbuff $ \cp -> + constHandler alpha $ \alphaPtr -> + constHandler beta $ \betaPtr -> do (axNew,_) <- return $ coordSwapper tra (ax,ay) --- dont need to swap b, info is in a and c --- c doesn't get implicitly transposed blasOrder <- return $ encodeNiceOrder ornta -- all three are the same orientation - rawTra <- return $ encodeFFITranspose tra + rawTra <- return $ encodeFFITranspose tra rawTrb <- return $ encodeFFITranspose trb -- example of why i want to switch to singletones - unsafePrimToPrim $! (if shouldCallFast cy cx axNew then gemmUnsafeFFI else gemmSafeFFI ) - blasOrder rawTra rawTrb (fromIntegral cy) (fromIntegral cx) (fromIntegral ax) + unsafePrimToPrim $! (if shouldCallFast cy cx axNew then gemmUnsafeFFI else gemmSafeFFI ) + blasOrder rawTra rawTrb (fromIntegral cy) (fromIntegral cx) (fromIntegral ax) alphaPtr ap (fromIntegral astride) bp (fromIntegral bstride) betaPtr cp (fromIntegral cstride) @@ -156,14 +156,14 @@ gemvAbstraction gemvName gemvSafeFFI gemvUnsafeFFI constHandler = gemv where shouldCallFast :: Int -> Int -> Bool - shouldCallFast a b = flopsThreshold >= gemvComplexity a b + shouldCallFast a b = flopsThreshold >= gemvComplexity a b gemv tr alpha beta (MutableDenseMatrix ornta ax ay astride abuff) (MutableDenseVector _ bdim bstride bbuff) (MutableDenseVector _ cdim cstride cbuff) | isBadGemv tr ax ay bdim cdim = error $! "Bad dimension args to GEMV: ax ay xdim ydim: " ++ show [ax, ay, bdim, cdim] | SM.overlaps abuff cbuff || SM.overlaps bbuff cbuff = - error $! "The read and write inputs for: " ++ gemvName ++ " overlap. This is a programmer error. Please fix." + error $! "The read and write inputs for: " ++ gemvName ++ " overlap. This is a programmer error. Please fix." | otherwise = call where (newx,newy) = coordSwapper tr (ax,ay) @@ -174,7 +174,7 @@ constHandler beta $ \betaPtr -> unsafePrimToPrim $! (if shouldCallFast newx newy then gemvUnsafeFFI else gemvSafeFFI) (encodeNiceOrder ornta) (encodeFFITranspose tr) - (fromIntegral newx) (fromIntegral newy) alphaPtr ap (fromIntegral astride) bp + (fromIntegral newx) (fromIntegral newy) alphaPtr ap (fromIntegral astride) bp (fromIntegral bstride) betaPtr cp (fromIntegral cstride) diff --git a/src/Numerical/HBLAS/Constants.lhs b/src/Numerical/HBLAS/Constants.lhs --- a/src/Numerical/HBLAS/Constants.lhs +++ b/src/Numerical/HBLAS/Constants.lhs @@ -1,6 +1,6 @@ \begin{code} -module Numerical.HBLAS.Constants(flopsThreshold) where -import Data.Int +module Numerical.HBLAS.Constants(flopsThreshold) where +import Data.Int flopsThreshold :: Int64 flopsThreshold = 10000 diff --git a/src/Numerical/HBLAS/Lapack.lhs b/src/Numerical/HBLAS/Lapack.lhs --- a/src/Numerical/HBLAS/Lapack.lhs +++ b/src/Numerical/HBLAS/Lapack.lhs @@ -1,12 +1,12 @@ \begin{code} -module Numerical.HBLAS.Lapack where - +module Numerical.HBLAS.Lapack where -import Numerical.HBLAS.UtilsFFI -import Numerical.HBLAS.Lapack.FFI + +import Numerical.HBLAS.UtilsFFI +import Numerical.HBLAS.Lapack.FFI import Numerical.HBLAS.MatrixTypes import Control.Monad.Primitive -import Data.Complex +import Data.Complex import qualified Data.Vector.Storable.Mutable as SM diff --git a/src/Numerical/HBLAS/Lapack/FFI.hs b/src/Numerical/HBLAS/Lapack/FFI.hs --- a/src/Numerical/HBLAS/Lapack/FFI.hs +++ b/src/Numerical/HBLAS/Lapack/FFI.hs @@ -4,8 +4,8 @@ import Foreign.Ptr import Foreign() import Foreign.C.Types -import Data.Complex -import Data.Int +import Data.Complex +import Data.Int @@ -14,7 +14,7 @@ stylenote: we will not use the LAPACKE_* operations, only the LAPACKE_*_work variants that require an explicitly provided work buffer. -This is to ensure that solver routine allocation behavior is transparent +This is to ensure that solver routine allocation behavior is transparent -} @@ -45,13 +45,13 @@ newtype Stride_C = Stride_C Int32 newtype Equilib_C = Equilib_C CChar -type Fun_FFI_GESVX el = Ptr Fact_C {- fact -}-> Ptr Trans_C {- trans -} - -> Ptr Int32 {-n -}-> Ptr Int32 {- NRHS -}-> +type Fun_FFI_GESVX el = Ptr Fact_C {- fact -}-> Ptr Trans_C {- trans -} + -> Ptr Int32 {-n -}-> Ptr Int32 {- NRHS -}-> Ptr el {- a -} -> Ptr Stride_C {- lda -} -> Ptr Double {- af -} -> Ptr Stride_C {- ldf-}-> Ptr Int32 -> Ptr Equilib_C {- equed -} -> Ptr el {- r -} -> Ptr el -> - Ptr el {- b -} -> Ptr Stride_C {- ld b -} -> Ptr el {- x -} -> Ptr Stride_C {- ldx -}-> + Ptr el {- b -} -> Ptr Stride_C {- ld b -} -> Ptr el {- x -} -> Ptr Stride_C {- ldx -}-> Ptr el {-rcond -}-> Ptr el {- ferr-} -> Ptr el {-berr-} -> Ptr el {-work-}-> - Ptr Int32 {-iwork -}-> Ptr Int32 {-info -} -> IO () + Ptr Int32 {-iwork -}-> Ptr Int32 {-info -} -> IO () @@ -61,12 +61,12 @@ -fact will be a 1 character C string -either - "F", then the inputs af and ipiv already contain the permuted LU factorization +fact will be a 1 character C string +either + "F", then the inputs af and ipiv already contain the permuted LU factorization (act as input rather than result params) "E", Matrix input A will be equilibriated if needed, then copied to AF and Factored - "N", matrix input A will be copied to AF + "N", matrix input A will be copied to AF -} @@ -75,13 +75,13 @@ Xgesvx is the s -sing -im assuming for now that any real use of *gesvx routines, or any other +im assuming for now that any real use of *gesvx routines, or any other n^3 complexity algs from LAPACK, are on inputs typically n>=15, which means > 1000 flops, -which is > 1µs, and thus ok to +which is > 1µs, and thus ok to -} --need to get around to wrapping these, but thats for another day -foreign import ccall "sgesvx_" sgesvx :: Fun_FFI_GESVX Float +foreign import ccall "sgesvx_" sgesvx :: Fun_FFI_GESVX Float foreign import ccall "dgesvx_" dgesvx :: Fun_FFI_GESVX Double foreign import ccall "cgesvx_" cgesvx :: Fun_FFI_GESVX (Complex Float) foreign import ccall "zgesvx_" zgesvx :: Fun_FFI_GESVX (Complex Double) @@ -89,20 +89,20 @@ --lapack_int ?syev_( char *jobz, char *uplo, lapack_int *n, ?* a, lapack_int * lda, ?* w ); --- ? is Double or Float +-- ? is Double or Float -newtype JobTy = JBT CChar +newtype JobTy = JBT CChar newtype UploTy = UPLT CChar -newtype Info = Info Int32 +newtype Info = Info Int32 --basic symmetric eigen value solvers type SYEV_FUN_FFI elem = Ptr JobTy -> Ptr UploTy -> Ptr Int32 -> Ptr elem -> Ptr Int32 -> Ptr elem -> Ptr Info-> IO () foreign import ccall "ssyev_" ssyev_ffi :: SYEV_FUN_FFI Float -foreign import ccall "dsyev_" dsyev_ffi :: SYEV_FUN_FFI Double +foreign import ccall "dsyev_" dsyev_ffi :: SYEV_FUN_FFI Double {-unsafe versions of lapack routines are meant to ONLY be used for workspace queries-} foreign import ccall unsafe "ssyev_" ssyev_ffi_unsafe :: SYEV_FUN_FFI Float -foreign import ccall unsafe "dsyev_" dsyev_ffi_unsafe :: SYEV_FUN_FFI Double +foreign import ccall unsafe "dsyev_" dsyev_ffi_unsafe :: SYEV_FUN_FFI Double --lapack_int LAPACKE_gesv( int matrix_order, lapack_int n, lapack_int nrhs, * a, lapack_int lda, lapack_int* ipiv, * b, lapack_int ldb ); --call sgesv( n, nrhs, a, lda, ipiv, b, ldb, info ) @@ -110,14 +110,14 @@ type GESV_FUN_FFI elem = Ptr Int32 {- n -} -> Ptr Int32 {- nrhs -} -> Ptr elem {- a -} -> Ptr Int32 {- lda -} -> Ptr Int32 {- permutation vector -} -> Ptr elem {- b -} -> Ptr Int32 {- ldb -} -> Ptr Info -> IO () --- | basic Linear system solvers. they act inplace -foreign import ccall "sgesv_" sgesv_ffi ::GESV_FUN_FFI Float +-- | basic Linear system solvers. they act inplace +foreign import ccall "sgesv_" sgesv_ffi ::GESV_FUN_FFI Float foreign import ccall "dgesv_" dgesv_ffi :: GESV_FUN_FFI Double foreign import ccall "cgesv_" cgesv_ffi :: GESV_FUN_FFI (Complex Float) foreign import ccall "zgesv_" zgesv_ffi :: GESV_FUN_FFI (Complex Double) -- / not sure if linear solvers ever are run in < 1 microsecond size instances in practice -foreign import ccall unsafe "sgesv_" sgesv_ffi_unsafe ::GESV_FUN_FFI Float +foreign import ccall unsafe "sgesv_" sgesv_ffi_unsafe ::GESV_FUN_FFI Float foreign import ccall unsafe "dgesv_" dgesv_ffi_unsafe :: GESV_FUN_FFI Double foreign import ccall unsafe "cgesv_" cgesv_ffi_unsafe :: GESV_FUN_FFI (Complex Float) foreign import ccall unsafe "zgesv_" zgesv_ffi_unsafe :: GESV_FUN_FFI (Complex Double) diff --git a/src/Numerical/HBLAS/MatrixTypes.hs b/src/Numerical/HBLAS/MatrixTypes.hs --- a/src/Numerical/HBLAS/MatrixTypes.hs +++ b/src/Numerical/HBLAS/MatrixTypes.hs @@ -10,8 +10,8 @@ {-| PSA, the matrix data types used in the hBLAS binding should not be regarded as being general purpose matrices. -They are designed to exactly express only the matrices which are -valid inputs for BLAS. When applicable, such matrices should be easily mapped +They are designed to exactly express only the matrices which are +valid inputs for BLAS. When applicable, such matrices should be easily mapped to and from other matrix libraries. That said, the BLAS and LAPACK matrix formats capture a rich and very expressive subset of Dense Matrix formats. @@ -28,73 +28,76 @@ There is a work in progress binding to help this in the numerical-hblas package (which may not be public yet at the time of this writing) --} +-} module Numerical.HBLAS.MatrixTypes where -import qualified Data.Vector.Storable as S +import qualified Data.Vector.Storable as S import qualified Data.Vector.Storable.Mutable as SM -import Control.Monad.Primitive +import Control.Monad.Primitive -import Data.Typeable +import Data.Typeable - + {- what I really want is this, but its not possible till datakinds works on types that aren't kind *, data Eff :: * -> * where - Pure :: Eff () - Mut :: s -> Eff s - -data EVector :: * -> * -> * where - PureVector :: S.Vector el -> EVector Pure el - MutVector :: SM.MVector s el -> EVector (Mut s) el + Pure :: Eff () + Mut :: s -> Eff s +data EVector :: * -> * -> * where + PureVector :: S.Vector el -> EVector Pure el + MutVector :: SM.MVector s el -> EVector (Mut s) el + -} --data Eff s where --- Pure :: Eff s --- Mut :: s -> Eff s +-- Pure :: Eff s +-- Mut :: s -> Eff s ---data EVector s el where --- PureVector :: S.Vector el -> EVector Pure el --- MutVector :: SM.MVector s e -> EVector (Mut s) el +--data EVector s el where +-- PureVector :: S.Vector el -> EVector Pure el +-- MutVector :: SM.MVector s e -> EVector (Mut s) el ---data Eff = Pure | Mut +--data Eff = Pure | Mut ---data EVector :: Eff -> * -> * -> * where --- PureVector :: S.Vector el -> EVector Pure () el --- MutVector :: SM.MVector s el -> EVector Mut s el +--data EVector :: Eff -> * -> * -> * where +-- PureVector :: S.Vector el -> EVector Pure () el +-- MutVector :: SM.MVector s el -> EVector Mut s el -data Orientation = Row | Column +data Orientation = Row | Column deriving (Eq,Show,Typeable) +type Row = 'Row +type Column = 'Column + data SOrientation :: Orientation -> * where - SRow :: SOrientation Row - SColumn :: SOrientation Column + SRow :: SOrientation Row + SColumn :: SOrientation Column #if defined(__GLASGOW_HASKELL_) && (__GLASGOW_HASKELL__ >= 707) deriving (Typeable) -#endif +#endif instance Show (SOrientation Row) where show _ = "SRow" instance Show (SOrientation Column) where - show _ = "SColumn" + show _ = "SColumn" instance Eq (SOrientation Row) where - (==) _ _ = True + (==) _ _ = True instance Eq (SOrientation Column) where - (==) _ _ = True + (==) _ _ = True -sTranpose :: (x~ TransposeF y, y~TransposeF x ) =>SOrientation x -> SOrientation y +sTranpose :: (x~ TransposeF y, y~TransposeF x ) =>SOrientation x -> SOrientation y sTranpose SColumn = SRow sTranpose SRow = SColumn @@ -105,19 +108,19 @@ -- | For Symmetric, Hermetian or Triangular matrices, which part is modeled. --- Applies to both Padded and Packed variants -data MatUpLo = MatUpper | MatLower +data MatUpLo = MatUpper | MatLower deriving(Typeable,Eq,Show) --- | Many triangular matrix routines expect to know if the matrix is +-- | Many triangular matrix routines expect to know if the matrix is -- all 1 (unit ) on the diagonal or not. Likewise, Many Factorizations routines -- can be assumed to return unit triangular matrices data MatDiag= MatUnit | MatNonUnit deriving(Typeable,Eq,Show) -- | For certain Square matrix product, do you want to Compute A*B or B*A --- only used as an argument -data EquationSide = LeftSide | RightSide +-- only used as an argument +data EquationSide = LeftSide | RightSide deriving(Typeable,Eq,Show) {- should think long and hard before adding implicit transposition to the internal data model @@ -126,7 +129,7 @@ type family TransposeF (x :: Orientation) :: Orientation type instance TransposeF Row = Column -type instance TransposeF Column = Row +type instance TransposeF Column = Row @@ -135,74 +138,74 @@ data Variant = Direct | Implicit deriving(Typeable,Eq,Show) -- | 'Variant' and 'SVariant' are a bit odd looking, --- They crop up when needing to talk about eg the row vectors of a +-- They crop up when needing to talk about eg the row vectors of a -- packed triangular row major matrix wrt both their logical size and manifest sizes --- this notion only makes sense in the 1dim case. +-- this notion only makes sense in the 1dim case. -- If you don't understand this parameter, just use 'SDirect' and 'Direct' --- as they will generally be the correct choice for most users. -data SVariant :: Variant -> * where - SImplicit :: {_frontPadding ::{-UNPACK-} !Int, _endPadding:: {-#UNPACK#-} !Int } -> SVariant Implicit - SDirect :: SVariant Direct +-- as they will generally be the correct choice for most users. +data SVariant :: Variant -> * where + SImplicit :: {_frontPadding ::{-UNPACK-} !Int, _endPadding:: {-#UNPACK#-} !Int } -> SVariant Implicit + SDirect :: SVariant Direct data DenseVector :: Variant -> * -> * where DenseVector :: { _VariantDenseVect :: !(SVariant varnt) - ,_LogicalDimDenseVector :: {-#UNPACK#-}!Int - ,_StrideDenseVector :: {-#UNPACK#-} ! Int - ,_bufferDenseVector :: !(S.Vector elem) - } -> DenseVector varnt elem + ,_LogicalDimDenseVector :: {-#UNPACK#-}!Int + ,_StrideDenseVector :: {-#UNPACK#-} ! Int + ,_bufferDenseVector :: !(S.Vector elem) + } -> DenseVector varnt elem #if defined(__GLASGOW_HASKELL_) && (__GLASGOW_HASKELL__ >= 707) deriving (Typeable) -#endif +#endif data MDenseVector :: * -> Variant -> * -> * where MutableDenseVector :: { _VariantMutDenseVect :: !(SVariant varnt) - ,_LogicalDimMutDenseVector :: {-#UNPACK#-}!Int - ,_StrideMutDenseVector :: {-#UNPACK#-} ! Int - ,_bufferMutDenseVector :: !(S.MVector s elem) - } -> MDenseVector s varnt elem + ,_LogicalDimMutDenseVector :: {-#UNPACK#-}!Int + ,_StrideMutDenseVector :: {-#UNPACK#-} ! Int + ,_bufferMutDenseVector :: !(S.MVector s elem) + } -> MDenseVector s varnt elem #if defined(__GLASGOW_HASKELL_) && (__GLASGOW_HASKELL__ >= 707) deriving (Typeable) -#endif +#endif -- | 'DenseMatrix' is for dense row or column major matrices -data DenseMatrix :: Orientation -> * -> * where +data DenseMatrix :: Orientation -> * -> * where DenseMatrix ::{ _OrientationMat :: SOrientation ornt , - _XdimDenMat :: {-# UNPACK #-}!Int, + _XdimDenMat :: {-# UNPACK #-}!Int, _YdimDenMat :: {-# UNPACK #-}!Int , - _StrideDenMat :: {-# UNPACK #-} !Int , - _bufferDenMat :: {-#UNPACK#-}!(S.Vector elem) }-> DenseMatrix ornt elem + _StrideDenMat :: {-# UNPACK #-} !Int , + _bufferDenMat :: {-#UNPACK#-}!(S.Vector elem) }-> DenseMatrix ornt elem #if defined(__GLASGOW_HASKELL_) && (__GLASGOW_HASKELL__ >= 707) deriving (Typeable) #endif -- | this should never be used in real code, ever ever, but its handy for testing -- but seriously never use this in real code, it doesn't do what you think --- because in the case of a matrix slice, the underlying buffer will have --- additional elements aside from the ones you expect! --- never use this in real code please. :) +-- because in the case of a matrix slice, the underlying buffer will have +-- additional elements aside from the ones you expect! +-- never use this in real code please. :) mutableVectorToList :: (PrimMonad m, S.Storable a) => S.MVector (PrimState m) a -> m [a] mutableVectorToList mv = do - v <- S.unsafeFreeze mv + v <- S.unsafeFreeze mv return (S.toList v ) {-# NOINLINE mutableVectorToList #-} {- -need to handle rendering a slice differently than a direct matrix +need to handle rendering a slice differently than a direct matrix -} instance (Show el,SM.Storable el )=> Show (DenseMatrix Row el) where show mat@(DenseMatrix SRow xdim ydim stride buffer) | stride == xdim = "DenseMatrix SRow " ++ " " ++show xdim ++ " " ++ show ydim ++ " " ++ show stride ++ "(" ++ show buffer ++ ")" - | otherwise = show $ mapDenseMatrix id mat + | otherwise = show $ mapDenseMatrix id mat instance (Show el,SM.Storable el )=> Show (DenseMatrix Column el) where show mat@(DenseMatrix SColumn xdim ydim stride buffer) | stride == xdim = "DenseMatrix SColumn " ++ " " ++show xdim ++ " " ++ show ydim ++ " " ++ show stride ++ "(" ++ show buffer ++ ")" - | otherwise = show $ mapDenseMatrix id mat + | otherwise = show $ mapDenseMatrix id mat --- | 'MDenseMatrix' -data MDenseMatrix :: * ->Orientation -> * -> * where +-- | 'MDenseMatrix' +data MDenseMatrix :: * ->Orientation -> * -> * where MutableDenseMatrix :: { _OrientationMutMat :: SOrientation ornt , - _XdimDenMutMat :: {-# UNPACK #-}!Int , + _XdimDenMutMat :: {-# UNPACK #-}!Int , _YdimDenMutMat :: {-# UNPACK #-}!Int, _StrideDenMutMat :: {-# UNPACK #-} !Int, _bufferDenMutMat :: {-# UNPACK #-} !(SM.MVector s elem) } -> MDenseMatrix s ornt elem @@ -210,20 +213,20 @@ --instance (Show el,SM.Storable el, PrimMonad m )=> Show (DenseMatrix (PrimState m) Row el) where -- show mat@(DenseMatrix SRow xdim ydim stride buffer) -- | stride == xdim = "MutableDenseMatrix SRow " ++ " " ++show xdim ++ " " ++ show ydim ++ " " ++ show stride ++ "(" ++ show buffer ++ ")" --- | otherwise = show $ mapDenseMatrix id mat +-- | otherwise = show $ mapDenseMatrix id mat --instance (Show el,SM.Storable el,PrimMonad m )=> Show (DenseMatrix (PrimState m ) Column el) where -- show mat@(DenseMatrix SColumn xdim ydim stride buffer) -- | stride == xdim = "DenseMatrix SColumn " ++ " " ++show xdim ++ " " ++ show ydim ++ " " ++ show stride ++ "(" ++ show buffer ++ ")" --- | otherwise = show $ mapDenseMatrix id mat +-- | otherwise = show $ mapDenseMatrix id mat -type IODenseMatrix = MDenseMatrix RealWorld ---type MutDenseMatrixIO or elem = +type IODenseMatrix = MDenseMatrix RealWorld +--type MutDenseMatrixIO or elem = -- data PaddedSymmetricMatrix -- data PaddedHermetianMatrix ---data PaddedTriangularMatrix ---- these three may just be wrappers for general dense matrices +--data PaddedTriangularMatrix +--- these three may just be wrappers for general dense matrices --data SymmetricMatrix --data HermitianMatrix -- may just be a wrapper for symmetric? @@ -232,20 +235,20 @@ {-#NOINLINE unsafeFreezeDenseMatrix #-} unsafeFreezeDenseMatrix :: (SM.Storable elem, PrimMonad m)=> MDenseMatrix (PrimState m) or elem -> m (DenseMatrix or elem) unsafeFreezeDenseMatrix (MutableDenseMatrix ornt a b c mv) = do - v <- S.unsafeFreeze mv - return $! DenseMatrix ornt a b c v + v <- S.unsafeFreeze mv + return $! DenseMatrix ornt a b c v {-# NOINLINE unsafeThawDenseMatrix #-} -unsafeThawDenseMatrix :: (SM.Storable elem, PrimMonad m)=> DenseMatrix or elem-> m (MDenseMatrix (PrimState m) or elem) -unsafeThawDenseMatrix (DenseMatrix ornt a b c v) = do +unsafeThawDenseMatrix :: (SM.Storable elem, PrimMonad m)=> DenseMatrix or elem-> m (MDenseMatrix (PrimState m) or elem) +unsafeThawDenseMatrix (DenseMatrix ornt a b c v) = do mv <- S.unsafeThaw v - return $! MutableDenseMatrix ornt a b c mv + return $! MutableDenseMatrix ornt a b c mv ---freezeDenseMatrix - +--freezeDenseMatrix + getDenseMatrixRow :: DenseMatrix or elem -> Int getDenseMatrixRow (DenseMatrix _ _ ydim _ _)= ydim @@ -254,27 +257,27 @@ getDenseMatrixColumn (DenseMatrix _ xdim _ _ _)= xdim -getDenseMatrixLeadingDimStride :: DenseMatrix or elem -> Int +getDenseMatrixLeadingDimStride :: DenseMatrix or elem -> Int getDenseMatrixLeadingDimStride (DenseMatrix _ _ _ stride _ ) = stride -getDenseMatrixArray :: DenseMatrix or elem -> S.Vector elem +getDenseMatrixArray :: DenseMatrix or elem -> S.Vector elem getDenseMatrixArray (DenseMatrix _ _ _ _ arr) = arr -getDenseMatrixOrientation :: DenseMatrix or elem -> SOrientation or -getDenseMatrixOrientation m = _OrientationMat m +getDenseMatrixOrientation :: DenseMatrix or elem -> SOrientation or +getDenseMatrixOrientation m = _OrientationMat m -uncheckedDenseMatrixIndex :: (S.Storable elem )=> DenseMatrix or elem -> (Int,Int) -> elem +uncheckedDenseMatrixIndex :: (S.Storable elem )=> DenseMatrix or elem -> (Int,Int) -> elem uncheckedDenseMatrixIndex (DenseMatrix SRow _ _ ystride arr) = \ (x,y)-> arr `S.unsafeIndex` (x + y * ystride) uncheckedDenseMatrixIndex (DenseMatrix SColumn _ _ xstride arr) = \ (x,y)-> arr `S.unsafeIndex` (y + x* xstride) -uncheckedDenseMatrixIndexM :: (Monad m ,S.Storable elem )=> DenseMatrix or elem -> (Int,Int) -> m elem +uncheckedDenseMatrixIndexM :: (Monad m ,S.Storable elem )=> DenseMatrix or elem -> (Int,Int) -> m elem uncheckedDenseMatrixIndexM (DenseMatrix SRow _ _ ystride arr) = \ (x,y)-> return $! arr `S.unsafeIndex` (x + y * ystride) uncheckedDenseMatrixIndexM (DenseMatrix SColumn _ _ xstride arr) = \ (x,y)-> return $! arr `S.unsafeIndex` (y + x* xstride) -uncheckedMutableDenseMatrixIndexM :: (PrimMonad m ,S.Storable elem )=> MDenseMatrix (PrimState m) or elem -> (Int,Int) -> m elem +uncheckedMutableDenseMatrixIndexM :: (PrimMonad m ,S.Storable elem )=> MDenseMatrix (PrimState m) or elem -> (Int,Int) -> m elem uncheckedMutableDenseMatrixIndexM (MutableDenseMatrix SRow _ _ ystride arr) = \ (x,y)-> arr `SM.unsafeRead` (x + y * ystride) uncheckedMutableDenseMatrixIndexM (MutableDenseMatrix SColumn _ _ xstride arr) = \ (x,y)-> arr `SM.unsafeRead` (y + x* xstride) @@ -283,82 +286,82 @@ {-# INLINE swap #-} -- | `map f matrix` -mapDenseMatrix :: (S.Storable a, S.Storable b) => (a->b) -> DenseMatrix or a -> DenseMatrix or b +mapDenseMatrix :: (S.Storable a, S.Storable b) => (a->b) -> DenseMatrix or a -> DenseMatrix or b mapDenseMatrix f rm@(DenseMatrix SRow xdim ydim _ _) = DenseMatrix SRow xdim ydim xdim $! - S.generate (xdim * ydim) (\ix -> f $! uncheckedDenseMatrixIndex rm (swap $ quotRem ix xdim ) ) -mapDenseMatrix f rm@(DenseMatrix SColumn xdim ydim _ _) = + S.generate (xdim * ydim) (\ix -> f $! uncheckedDenseMatrixIndex rm (swap $ quotRem ix xdim ) ) +mapDenseMatrix f rm@(DenseMatrix SColumn xdim ydim _ _) = DenseMatrix SColumn xdim ydim ydim $! S.generate (xdim * ydim ) (\ix -> f $! uncheckedDenseMatrixIndex rm ( quotRem ix ydim ) ) -imapDenseMatrix :: (S.Storable a, S.Storable b) => ((Int,Int)->a->b) -> DenseMatrix or a -> DenseMatrix or b -imapDenseMatrix f rm@(DenseMatrix sornt xdim ydim _ _) = +imapDenseMatrix :: (S.Storable a, S.Storable b) => ((Int,Int)->a->b) -> DenseMatrix or a -> DenseMatrix or b +imapDenseMatrix f rm@(DenseMatrix sornt xdim ydim _ _) = generateDenseMatrix sornt (xdim,ydim) (\ix -> f ix $! uncheckedDenseMatrixIndex rm ix ) -- | In Matrix format memory order enumeration of the index tuples, for good locality 2dim map uncheckedDenseMatrixNextTuple :: DenseMatrix or elem -> (Int,Int) -> Maybe (Int,Int) -uncheckedDenseMatrixNextTuple (DenseMatrix SRow xdim ydim _ _) = - \(!x,!y)-> if (x >= xdim && y >= ydim) then Nothing else Just $! swap $! quotRem (x+ xdim * y + 1) xdim -uncheckedDenseMatrixNextTuple (DenseMatrix SColumn xdim ydim _ _ ) = - \(!x,!y) -> if (x >= xdim && y >= ydim) then Nothing else Just $! quotRem (y + ydim * x + 1) ydim +uncheckedDenseMatrixNextTuple (DenseMatrix SRow xdim ydim _ _) = + \(!x,!y)-> if (x >= xdim && y >= ydim) then Nothing else Just $! swap $! quotRem (x+ xdim * y + 1) xdim +uncheckedDenseMatrixNextTuple (DenseMatrix SColumn xdim ydim _ _ ) = + \(!x,!y) -> if (x >= xdim && y >= ydim) then Nothing else Just $! quotRem (y + ydim * x + 1) ydim --- dont need the swap for column major -- | generateDenseMatrix Row (k,k) \(i,j)-> if i == j then 1.0 else 0.0 would generate a KxK identity matrix -generateDenseMatrix :: (S.Storable a)=> SOrientation x -> (Int,Int)->((Int,Int)-> a) -> DenseMatrix x a +generateDenseMatrix :: (S.Storable a)=> SOrientation x -> (Int,Int)->((Int,Int)-> a) -> DenseMatrix x a generateDenseMatrix SRow (xdim,ydim) f = DenseMatrix SRow xdim ydim xdim $! - S.generate (xdim * ydim) (\ix -> let !ixtup@(!_,!_) = swap $ quotRem ix xdim in - f ixtup ) + S.generate (xdim * ydim) (\ix -> let !ixtup@(!_,!_) = swap $ quotRem ix xdim in + f ixtup ) generateDenseMatrix SColumn (xdim,ydim) f = DenseMatrix SColumn xdim ydim ydim $! - S.generate (xdim * ydim ) (\ix -> let ixtup@(!_,!_) = ( quotRem ix ydim ) in - f ixtup ) + S.generate (xdim * ydim ) (\ix -> let ixtup@(!_,!_) = ( quotRem ix ydim ) in + f ixtup ) -- | mutable version of generateDenseMatrix {-# NOINLINE generateMutableDenseMatrix #-} -generateMutableDenseMatrix :: (S.Storable a,PrimMonad m)=> - SOrientation x -> (Int,Int)->((Int,Int)-> a) -> m (MDenseMatrix (PrimState m) x a) +generateMutableDenseMatrix :: (S.Storable a,PrimMonad m)=> + SOrientation x -> (Int,Int)->((Int,Int)-> a) -> m (MDenseMatrix (PrimState m) x a) generateMutableDenseMatrix sor dims fun = do - x <- unsafeThawDenseMatrix $! generateDenseMatrix sor dims fun - return x + x <- unsafeThawDenseMatrix $! generateDenseMatrix sor dims fun + return x {-#NOINLINE generateMutableDenseVector#-} generateMutableDenseVector :: (S.Storable a,PrimMonad m) => Int -> (Int -> a) -> m (MDenseVector (PrimState m ) Direct a) generateMutableDenseVector size init = do - mv <- S.unsafeThaw $ S.generate size init - return $! MutableDenseVector SDirect size 1 mv + mv <- S.unsafeThaw $ S.generate size init + return $! MutableDenseVector SDirect size 1 mv --- this (uncheckedMatrixSlice) will need to have its inlining quality checked --- | slice over matrix element in the range (inclusive) [xstart..xend] X [ystart .. yend] --- call as @'uncheckedMatrixSlice' matrix (xstart,ystart) (xend,yend) @ -uncheckedDenseMatrixSlice :: (S.Storable elem)=> DenseMatrix or elem -> (Int,Int)-> (Int,Int)-> DenseMatrix or elem +uncheckedDenseMatrixSlice :: (S.Storable elem)=> DenseMatrix or elem -> (Int,Int)-> (Int,Int)-> DenseMatrix or elem uncheckedDenseMatrixSlice (DenseMatrix SRow xdim _ ystride arr) (xstart,ystart) (xend,yend) = res where !res = DenseMatrix SRow (xend - xstart + 1) -- X : n - 0 + 1, because zero indexed (yend - ystart+1) -- Y : m - 0 + 1, because zero indexed (ystride + xstart + (xdim - xend)) -- how much start and end padding per row (S.slice ixStart (ixEnd - ixStart) arr ) !ixStart = (xstart+ystart*ystride) - !ixEnd = (xend+yend*ystride) + !ixEnd = (xend+yend*ystride) uncheckedDenseMatrixSlice (DenseMatrix SColumn _ ydim xstride arr) (xstart,ystart) (xend,yend) = res - where !res = DenseMatrix SColumn (xend - xstart + 1) - (yend - ystart+1) + where !res = DenseMatrix SColumn (xend - xstart + 1) + (yend - ystart+1) (xstride + ystart + (ydim - yend)) (S.slice ixStart (ixEnd - ixStart) arr ) !ixStart = (ystart+xstart*xstride) !ixEnd = (yend+xend*xstride) -- | tranposeMatrix does a shallow transpose that swaps the format and the x y params, but changes nothing --- in the memory layout. +-- in the memory layout. -- Most applications where transpose is used in a computation need a deep, copying, tranpose operation -transposeDenseMatrix :: (inor ~ (TransposeF outor) , outor ~ (TransposeF inor) ) => DenseMatrix inor elem -> DenseMatrix outor elem +transposeDenseMatrix :: (inor ~ (TransposeF outor) , outor ~ (TransposeF inor) ) => DenseMatrix inor elem -> DenseMatrix outor elem transposeDenseMatrix (DenseMatrix orient x y stride arr)= (DenseMatrix (sTranpose orient) y x stride arr) diff --git a/src/Numerical/HBLAS/UtilsFFI.hs b/src/Numerical/HBLAS/UtilsFFI.hs --- a/src/Numerical/HBLAS/UtilsFFI.hs +++ b/src/Numerical/HBLAS/UtilsFFI.hs @@ -4,48 +4,48 @@ -import Data.Vector.Storable.Mutable as M +import Data.Vector.Storable.Mutable as M import Control.Monad.Primitive import Foreign.ForeignPtr.Safe import Foreign.ForeignPtr.Unsafe import Foreign.Storable.Complex() -import Data.Vector.Storable as S +import Data.Vector.Storable as S import Foreign.Ptr {- the IO version of these various utils is in Base. -but would like to have the +but would like to have the -} -withRWStorable:: (Storable a, PrimMonad m)=> a -> (Ptr a -> m b) -> m a -withRWStorable val fun = do - valVect <- M.replicate 1 val - _ <- unsafeWithPrim valVect fun - M.unsafeRead valVect 0 -{-# INLINE withRWStorable #-} +withRWStorable:: (Storable a, PrimMonad m)=> a -> (Ptr a -> m b) -> m a +withRWStorable val fun = do + valVect <- M.replicate 1 val + _ <- unsafeWithPrim valVect fun + M.unsafeRead valVect 0 +{-# INLINE withRWStorable #-} -withRStorable :: (Storable a, PrimMonad m)=> a -> (Ptr a -> m b) -> m b -withRStorable val fun = do - valVect <- M.replicate 1 val - unsafeWithPrim valVect fun -{-# INLINE withRStorable #-} +withRStorable :: (Storable a, PrimMonad m)=> a -> (Ptr a -> m b) -> m b +withRStorable val fun = do + valVect <- M.replicate 1 val + unsafeWithPrim valVect fun +{-# INLINE withRStorable #-} withRStorable_ :: (Storable a, PrimMonad m)=> a -> (Ptr a -> m ()) -> m () -withRStorable_ val fun = do - valVect <- M.replicate 1 val - unsafeWithPrim valVect fun +withRStorable_ val fun = do + valVect <- M.replicate 1 val + unsafeWithPrim valVect fun - return () -{-# INLINE withRStorable_ #-} + return () +{-# INLINE withRStorable_ #-} withForeignPtrPrim :: PrimMonad m => ForeignPtr a -> (Ptr a -> m b) -> m b withForeignPtrPrim fo act = do r <- act (unsafeForeignPtrToPtr fo) touchForeignPtrPrim fo return r -{-# INLINE withForeignPtrPrim #-} +{-# INLINE withForeignPtrPrim #-} touchForeignPtrPrim ::PrimMonad m => ForeignPtr a -> m () touchForeignPtrPrim fp = unsafePrimToPrim $! touchForeignPtr fp @@ -64,17 +64,17 @@ unsafeWithPurePrim ::( Storable a, PrimMonad m )=> Vector a -> ((Ptr a)-> m b) -> m b {-# INLINE unsafeWithPurePrim #-} -unsafeWithPurePrim v fun = case S.unsafeToForeignPtr0 v of - (fp,_) -> do +unsafeWithPurePrim v fun = case S.unsafeToForeignPtr0 v of + (fp,_) -> do res <- withForeignPtrPrim fp (\x -> fun x) - touchForeignPtrPrim fp - return res + touchForeignPtrPrim fp + return res unsafeWithPurePrimLen ::( Storable a, PrimMonad m )=> Vector a -> ((Ptr a, Int )-> m b) -> m b {-# INLINE unsafeWithPurePrimLen #-} -unsafeWithPurePrimLen v fun = case S.unsafeToForeignPtr0 v of - (fp,n) -> do +unsafeWithPurePrimLen v fun = case S.unsafeToForeignPtr0 v of + (fp,n) -> do res <- withForeignPtrPrim fp (\x -> fun (x,n)) - touchForeignPtrPrim fp - return res + touchForeignPtrPrim fp + return res