packages feed

hblas 0.1.0.0 → 0.2.0.0

raw patch · 11 files changed

+571/−137 lines, 11 filesdep +HUnitdep +tastydep +tasty-hunitdep ~basedep ~vector

Dependencies added: HUnit, tasty, tasty-hunit

Dependency ranges changed: base, vector

Files

+ changelog.md view
@@ -0,0 +1,9 @@+++# PRE-RELEASE version 0.2.0.0+* Support for all BLAS operations defined on General Dense matrices,+both row and column major. +*  general dense linear and least squares solvers from LAPACK. Simple Drivers.++# version 0.1.0.0+* support basic manipulations of BLAS style dense matrices
hblas.cabal view
@@ -10,13 +10,13 @@ -- PVP summary:      +-+------- breaking API changes --                   | | +----- non-breaking API additions --                   | | | +--- code changes with no API change-version:             0.1.0.0+version:             0.2.0.0  -- A short (one-line) description of the package.-synopsis:            BLAS and Lapack bindings for OpenBLAS+synopsis:            Human friendly BLAS and Lapack bindings for Haskell.  -- A longer description of the package.-description:   User friendly, simple bindings to BLAS and Lapack, favoring OpenBLAS as the substrate.+description:   User friendly, simple bindings to BLAS and Lapack. Easy to extend and use.   -- The license under which the package is released. license:             BSD3@@ -40,17 +40,20 @@   extra-source-files:  -  readme.md +  readme.md,+  changelog.md    source-repository head   type: git-  location: git://github.com/wellposed/hOpenBLAS.git+  location: http://github.com/wellposed/hblas.git    -- Extra files to be distributed with the package, such as examples or a  -- README. -+flag OpenBLAS+  default: False +  manual: True  -- Constraint on the version of Cabal needed to build this package. -- flag lib-Werror --   default: False@@ -63,10 +66,13 @@ library    -- Modules exported by the library.-  exposed-modules:      Numerical.HBLAS.BLAS.FFI,-      Numerical.HBLAS.MatrixTypes,-      Numerical.HBLAS.UtilsFFI-      Numerical.HBLAS.BLAS+  exposed-modules:      Numerical.HBLAS.BLAS.FFI+    Numerical.HBLAS.MatrixTypes+    Numerical.HBLAS.UtilsFFI+    Numerical.HBLAS.BLAS+    Numerical.HBLAS.Lapack+    Numerical.HBLAS.Lapack.FFI +    Numerical.HBLAS.Constants   -- Modules included in this library but not exported.   -- other-modules:          @@ -80,48 +86,49 @@   -- Other library packages from which modules are imported.   build-depends:       base >=4.5 && <4.8, storable-complex >= 0.2.0 && < 0.3.0,           vector , primitive >= 0.5 && < 0.6-  -- look, no external deps for openblas! -  -- extra-libraries: openblas pthread -  if os(OSX) +  if flag(OpenBLAS)+    extra-libraries: openblas pthread++  if flag(OpenBLAS)&& os(OSX)+    extra-lib-dirs: /usr/local/lib ++  if os(OSX)&&!flag(OpenBLAS)       frameworks: Accelerate-      -- probably dont need these extra lib dirs-      -- extra-libraries: cblas clapack+      -- extra-libraries: cblas lapack -  if os(windows)+  if os(windows) && !flag(OpenBLAS)       extra-libraries: blas lapack -  if !os(windows)&& !os(OSX)+  if !os(windows)&& !os(OSX) && !flag(OpenBLAS)     extra-libraries: blas lapack    - --  if os(OSX)-        -- where macports installs libs-        --- extra-lib-dirs: /opt/local/lib/-        --- include-dirs: /opt/local/include/-        -- where homebrew installs libs-   ---      extra-lib-dirs: /usr/local/lib/-     --   include-dirs: /usr/local/include/ --        -- by default brew doesn't link openblas-        -- extra-lib-dirs: /usr/local/opt/openblas/lib-        -- include-dirs: /usr/local/opt/openblas/include-   -- Directories containing source files.   hs-source-dirs:      src      -- Base language which the package is written in.   default-language:    Haskell2010-  -- ld-options: OpenBLAS/libopenblas.a  -Test-Suite simple-test-    type:                exitcode-stdio-1.0-    main-is:             testing/Test.hs-    build-depends:       base, vector -    default-language:    Haskell2010-    build-depends:       hblas---    extra-libraries: openblas+  ++++Test-suite unit-testsuite +   default-language: Haskell2010+   type: exitcode-stdio-1.0+   build-depends: base >=4 && < 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 +   main-is: MainUnit.hs+
readme.md view
@@ -1,4 +1,4 @@-[![Wellposed](http://www.wellposed.com/mini.png)](http://www.wellposed.com)® +[![Wellposed](http://www.wellposed.com/mini.png)](http://www.wellposed.com)™   # About hblas @@ -24,7 +24,11 @@ high performance linear algebra routines.   -## how to install+## how to install (using openblas)+By default, hblas will assume you have OpenBLAS built and installed, including+the LAPACKE interfaces for LAPACK, somewhere in your standard library path.++ * On OS X systems, things will just work. * On linux and bsd systems, the equivalent of  ```@@ -34,3 +38,11 @@  ## getting involved patches, bug reports, tests,  and other contributions welcome.++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.+++# 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.
src/Numerical/HBLAS/BLAS.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE BangPatterns , RankNTypes, GADTs #-}+{-# LANGUAGE BangPatterns , RankNTypes, GADTs, DataKinds #-}   module Numerical.HBLAS.BLAS where@@ -8,19 +8,23 @@ --    ,cgemm --    ,zgemm)  -+import Numerical.HBLAS.Constants import Numerical.HBLAS.UtilsFFI     import Numerical.HBLAS.BLAS.FFI  import Numerical.HBLAS.MatrixTypes import Control.Monad.Primitive import Data.Complex  import qualified Data.Vector.Storable.Mutable as SM+import Numerical.HBLAS.Constants+import Data.Int  ----flopsThreshold = 10000-gemmComplexity a b c = a * b * c  -- this will be wrong by some constant factor, albeit a small one+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  + -- this covers the ~6 cases for checking the dimensions for GEMM quite nicely isBadGemm tra trb  ax ay bx by cx cy = isBadGemmHelper (cds tra (ax,ay)) (cds trb (bx,by) )  (cx,cy)     where @@ -34,35 +38,72 @@ coordSwapper Transpose (a,b) = (b,a) coordSwapper ConjTranspose (a,b) = (b,a) +-- / 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 tr ax ay bdim cdim = isBadGemvHelper (cds tr (ax,ay))+    where +    cds = coordSwapper+    isBadGemvHelper (realX,realY)  =  +            minimum [realY,realX,bdim,cdim] <= 0 ||  not (realX == bdim && realY == cdim ) + encodeNiceOrder :: SOrientation x  -> CBLAS_ORDERT encodeNiceOrder SRow= encodeOrder  BLASRowMajor encodeNiceOrder SColumn= encodeOrder BLASColMajor  -encodeFFITranpose :: Transpose -> CBLAS_TRANSPOSET-encodeFFITranpose  x=  encodeTranpose $ encodeNiceTranpose x +encodeFFITranspose :: Transpose -> CBLAS_TRANSPOSET+encodeFFITranspose  x=  encodeTranspose $ encodeNiceTranspose x  -encodeNiceTranpose :: Transpose -> BLAS_Transpose-encodeNiceTranpose x = case x of +encodeNiceTranspose :: Transpose -> BLAS_Transpose+encodeNiceTranspose x = case x of          NoTranspose -> BlasNoTranspose         Transpose -> BlasTranspose         ConjTranspose -> BlasConjTranspose         ConjNoTranspose -> BlasConjNoTranspose ---data BLAS_Tranpose = BlasNoTranspose | BlasTranpose | BlasConjTranspose | BlasConjNoTranpose ---data Tranpose = NoTranpose | Tranpose | ConjTranpose | ConjNoTranpose+encodeFFIMatrixHalf :: MatUpLo -> CBLAS_UPLOT+encodeFFIMatrixHalf x = encodeUPLO $ encodeNiceUPLO x +encodeNiceUPLO :: MatUpLo -> BLASUplo+encodeNiceUPLO x = case x of+                    MatUpper  -> BUpper+                    MatLower  -> BLower -type GemmFun el orient s m = Transpose ->Transpose ->  el -> el  -> MutDenseMatrix s orient el-  ->   MutDenseMatrix s orient el  ->  MutDenseMatrix s orient el -> m ()+encodeFFITriangleSort :: MatDiag -> CBLAS_DIAGT+encodeFFITriangleSort x = encodeDiag $ encodeNiceDIAG x +encodeNiceDIAG :: MatDiag -> BlasDiag+encodeNiceDIAG x = case x of+                    MatUnit    -> BlasUnit+                    MatNonUnit -> BlasNonUnit +--data BLAS_Transpose = BlasNoTranspose | BlasTranspose | BlasConjTranspose | BlasConjNoTranspose +--data Transpose = NoTranspose | Transpose | ConjTranspose | ConjNoTranspose+++type GemmFun el orient s m = Transpose ->Transpose ->  el -> el  -> MDenseMatrix s orient el+  ->   MDenseMatrix s orient el  ->  MDenseMatrix s orient el -> m ()++type GemvFun el orient s m = Transpose -> el -> el+  -> MDenseMatrix s orient el -> MDenseVector s Direct el -> MDenseVector s Direct el -> m ()+++type TrsvFun el orient s m =+      MatUpLo -> Transpose -> MatDiag+   -> MDenseMatrix s orient el -> MDenseVector s Direct el -> m ()  ++ {- 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.  + -}  @@ -98,8 +139,8 @@                         --- 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 $  encodeFFITranpose tra -                        rawTrb <- return $   encodeFFITranpose trb+                        rawTra <- return $  encodeFFITranspose tra +                        rawTrb <- return $   encodeFFITranspose trb                                  -- example of why i want to switch to singletones                         unsafePrimToPrim $!  (if shouldCallFast cy cx ax then gemmUnsafeFFI  else gemmSafeFFI )                              blasOrder rawTra rawTrb (fromIntegral cy) (fromIntegral cx) (fromIntegral ax) @@ -114,25 +155,108 @@   ->   DenseMatrix  orient el  ->  DenseMatrix  orient el  -}  sgemm :: PrimMonad m=> -     Transpose ->Transpose ->  Float -> Float  -> MutDenseMatrix (PrimState m) orient Float-  ->   MutDenseMatrix (PrimState m) orient Float  ->  MutDenseMatrix (PrimState m) orient Float -> m ()-sgemm =  gemmAbstraction "sgemm" cblas_sgemm_unsafe cblas_sgemm_safe (\x f -> f x )                                 +     Transpose ->Transpose ->  Float -> Float  -> MDenseMatrix (PrimState m) orient Float+  ->   MDenseMatrix (PrimState m) orient Float  ->  MDenseMatrix (PrimState m) orient Float -> m ()+sgemm =  gemmAbstraction "sgemm"  cblas_sgemm_safe cblas_sgemm_unsafe (\x f -> f x )                                                            dgemm :: PrimMonad m=> -     Transpose ->Transpose ->  Double -> Double -> MutDenseMatrix (PrimState m) orient Double-  ->   MutDenseMatrix (PrimState m) orient Double   ->  MutDenseMatrix (PrimState m) orient Double -> m ()-dgemm = gemmAbstraction "dgemm" cblas_dgemm_unsafe cblas_dgemm_safe (\x f -> f x )    +     Transpose ->Transpose ->  Double -> Double -> MDenseMatrix (PrimState m) orient Double+  ->   MDenseMatrix (PrimState m) orient Double   ->  MDenseMatrix (PrimState m) orient Double -> m ()+dgemm = gemmAbstraction "dgemm"  cblas_dgemm_safe cblas_dgemm_unsafe  (\x f -> f x )        cgemm :: PrimMonad m=>  Transpose ->Transpose ->  (Complex Float) -> (Complex Float)  -> -        MutDenseMatrix (PrimState m) orient (Complex Float)  ->   -        MutDenseMatrix (PrimState m) orient (Complex Float)  ->  -        MutDenseMatrix (PrimState m) orient (Complex Float) -> m ()-cgemm = gemmAbstraction "cgemm" cblas_cgemm_unsafe cblas_cgemm_safe withRStorable_                                +        MDenseMatrix (PrimState m) orient (Complex Float)  ->   +        MDenseMatrix (PrimState m) orient (Complex Float)  ->  +        MDenseMatrix (PrimState m) orient (Complex Float) -> m ()+cgemm = gemmAbstraction "cgemm" cblas_cgemm_safe cblas_cgemm_unsafe  withRStorable_                                  zgemm :: PrimMonad m=>  Transpose ->Transpose ->  (Complex Double) -> (Complex Double )  -> -        MutDenseMatrix (PrimState m) orient (Complex Double )  ->   -        MutDenseMatrix (PrimState m) orient (Complex Double)  ->  -        MutDenseMatrix (PrimState m) orient (Complex Double) -> m ()-zgemm = gemmAbstraction "zgemm" cblas_zgemm_unsafe cblas_zgemm_safe withRStorable_  +        MDenseMatrix (PrimState m) orient (Complex Double )  ->   +        MDenseMatrix (PrimState m) orient (Complex Double)  ->  +        MDenseMatrix (PrimState m) orient (Complex Double) -> m ()+zgemm = gemmAbstraction "zgemm"  cblas_zgemm_safe cblas_zgemm_unsafe withRStorable_  +++{-# NOINLINE gemvAbstraction #-}+gemvAbstraction :: (SM.Storable el, PrimMonad m)+                => String+                -> GemvFunFFI scale el+                -> GemvFunFFI scale el+                -> (el -> (scale -> m ())-> m ())+                -> forall orient . GemvFun el orient (PrimState m) m+gemvAbstraction gemvName gemvSafeFFI gemvUnsafeFFI constHandler = gemv+  where+    shouldCallFast :: Int -> Int  -> Bool+    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." +        | otherwise = call+            where+              (newx,newy) = coordSwapper tr (ax,ay)+              call = unsafeWithPrim abuff $ \ap ->+                     unsafeWithPrim bbuff $ \bp ->+                     unsafeWithPrim cbuff $ \cp ->+                     constHandler alpha $ \alphaPtr ->+                     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 bstride) betaPtr cp (fromIntegral cstride)++sgemv :: PrimMonad m => GemvFun Float orient (PrimState m) m+sgemv = gemvAbstraction "sgemv" cblas_sgemv_safe cblas_sgemv_unsafe (flip ($))++dgemv :: PrimMonad m => GemvFun Double orient (PrimState m) m+dgemv = gemvAbstraction "dgemv" cblas_dgemv_safe cblas_dgemv_unsafe (flip ($))++cgemv :: PrimMonad m => GemvFun (Complex Float) orient (PrimState m) m+cgemv = gemvAbstraction "cgemv" cblas_cgemv_safe cblas_cgemv_unsafe withRStorable_++zgemv :: PrimMonad m => GemvFun (Complex Double) orient (PrimState m) m+zgemv = gemvAbstraction "zgemv" cblas_zgemv_safe cblas_zgemv_unsafe withRStorable_++{-# NOINLINE trsvAbstraction #-}+trsvAbstraction :: (SM.Storable el, PrimMonad m)+                => String+                -> TrsvFunFFI el -> TrsvFunFFI el+                -> forall orient . TrsvFun el orient (PrimState m) m+trsvAbstraction trsvName trsvSafeFFI trsvUnsafeFFI = trsv+  where+    shouldCallFast :: Int -> Bool+    shouldCallFast n = flopsThreshold >= (fromIntegral n)^2++    isBadTrsv :: Int -> Int -> Int -> Bool+    isBadTrsv nx ny vdim = nx < 0 || nx /= ny || nx /= vdim++    trsv uplo tra diag+      (MutableDenseMatrix ornt x y mstride mbuff)+      (MutableDenseVector _ vdim vstride vbuff)+        | isBadTrsv x y vdim =+            error $! "Bad dimension args to TRSV: x y vdim: " ++ show [x,y,vdim]+        | SM.overlaps vbuff mbuff =+            error $! "The read and write inputs for: " ++ trsvName ++ " overlap. This is a programmer error. Please fix."+        | otherwise = unsafeWithPrim mbuff $ \mp ->+                      unsafeWithPrim vbuff $ \vp ->+                        unsafePrimToPrim $! (if shouldCallFast x then trsvUnsafeFFI else trsvSafeFFI)+                          (encodeNiceOrder ornt) (encodeFFIMatrixHalf uplo) (encodeFFITranspose tra)+                          (encodeFFITriangleSort diag) (fromIntegral x) mp (fromIntegral mstride) vp+                          (fromIntegral vstride)++strsv :: PrimMonad m => TrsvFun Float orient (PrimState m) m+strsv = trsvAbstraction "strsv" cblas_strsv_safe cblas_strsv_unsafe++dtrsv :: PrimMonad m => TrsvFun Double orient (PrimState m) m+dtrsv = trsvAbstraction "dtrsv" cblas_dtrsv_safe cblas_dtrsv_unsafe++ctrsv :: PrimMonad m => TrsvFun (Complex Float) orient (PrimState m) m+ctrsv = trsvAbstraction "ctrsv" cblas_ctrsv_safe cblas_ctrsv_unsafe++ztrsv :: PrimMonad m => TrsvFun (Complex Double) orient (PrimState m) m+ztrsv = trsvAbstraction "ztrsv" cblas_ztrsv_safe cblas_ztrsv_unsafe
src/Numerical/HBLAS/BLAS/FFI.hs view
@@ -44,11 +44,11 @@  data BLAS_Transpose = BlasNoTranspose | BlasTranspose | BlasConjTranspose | BlasConjNoTranspose  -encodeTranpose :: BLAS_Transpose -> CBLAS_TRANSPOSET-encodeTranpose  BlasNoTranspose = CBLAS_TransposeT 111-encodeTranpose  BlasTranspose = CBLAS_TransposeT 112-encodeTranpose  BlasConjTranspose =  CBLAS_TransposeT 113-encodeTranpose  BlasConjNoTranspose = CBLAS_TransposeT 114+encodeTranspose :: BLAS_Transpose -> CBLAS_TRANSPOSET+encodeTranspose  BlasNoTranspose = CBLAS_TransposeT 111+encodeTranspose  BlasTranspose = CBLAS_TransposeT 112+encodeTranspose  BlasConjTranspose =  CBLAS_TransposeT 113+encodeTranspose  BlasConjNoTranspose = CBLAS_TransposeT 114  newtype CBLAS_UPLOT = CBlasUPLO CInt      deriving (Eq,Show)@@ -186,13 +186,13 @@   -foreign import ccall unsafe "cblas_scopy" cblas_sswap_unsafe :: +foreign import ccall unsafe "cblas_sswap" cblas_sswap_unsafe ::      CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO ()   -foreign import ccall unsafe "cblas_dcopy" cblas_dswap_unsafe :: +foreign import ccall unsafe "cblas_dswap" cblas_dswap_unsafe ::      CInt -> Ptr Double-> CInt -> Ptr Double -> CInt -> IO ()-foreign import ccall unsafe "cblas_ccopy" 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_zcopy" cblas_zswap_unsafe :: +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);@@ -234,6 +234,30 @@       alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, -} +type GemvFunFFI sc el =+       CBLAS_ORDERT -> CBLAS_TRANSPOSET -> CInt -> CInt+    -> sc -> Ptr el -> CInt -> Ptr el -> CInt -> sc -> Ptr el -> CInt -> IO ()++foreign import ccall unsafe "cblas_sgemv"+    cblas_sgemv_unsafe :: GemvFunFFI Float Float+foreign import ccall safe   "cblas_sgemv"+    cblas_sgemv_safe   :: GemvFunFFI Float Float++foreign import ccall unsafe "cblas_dgemv"+    cblas_dgemv_unsafe :: GemvFunFFI Double Double+foreign import ccall safe   "cblas_dgemv"+    cblas_dgemv_safe   :: GemvFunFFI Double Double++foreign import ccall unsafe "cblas_cgemv"+    cblas_cgemv_unsafe :: GemvFunFFI (Ptr (Complex Float)) (Complex Float)+foreign import ccall safe   "cblas_cgemv"+    cblas_cgemv_safe   :: GemvFunFFI (Ptr (Complex Float)) (Complex Float)++foreign import ccall unsafe "cblas_zgemv"+    cblas_zgemv_unsafe :: GemvFunFFI (Ptr (Complex Double)) (Complex Double)+foreign import ccall safe   "cblas_zgemv"+    cblas_zgemv_safe   :: GemvFunFFI (Ptr (Complex Double)) (Complex Double)+ --void cblas_sgemv(  enum CBLAS_ORDER order,    enum CBLAS_TRANSPOSE trans,    CInt m,   CInt n, --           Float alpha,   Float  *a,   CInt lda,    Float  *x,   CInt incx,    Float beta,  Float  *y,   CInt incy); --void cblas_dgemv(  enum CBLAS_ORDER order,    enum CBLAS_TRANSPOSE trans,    CInt m,   CInt n,@@ -256,6 +280,32 @@  --STRSV - solve one of the systems of equations   A*x = b, or A'*x = b, where A is a (non)unit upper(/lower) triangular matrix +--STRSV - solve one of the systems of equations   A*x = b, or A'*x = b, where A is a (non)unit upper(/lower) triangular matrix++type TrsvFunFFI el =+       CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT+    -> CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()++foreign import ccall unsafe "cblas_strsv"+  cblas_strsv_unsafe :: TrsvFunFFI Float+foreign import ccall safe   "cblas_strsv"+  cblas_strsv_safe   :: TrsvFunFFI Float++foreign import ccall unsafe "cblas_dtrsv"+  cblas_dtrsv_unsafe :: TrsvFunFFI Double+foreign import ccall safe   "cblas_dtrsv"+  cblas_dtrsv_safe   :: TrsvFunFFI Double++foreign import ccall unsafe "cblas_ctrsv"+  cblas_ctrsv_unsafe :: TrsvFunFFI (Complex Float)+foreign import ccall safe   "cblas_ctrsv"+  cblas_ctrsv_safe   :: TrsvFunFFI (Complex Float)++foreign import ccall unsafe "cblas_ztrsv"+  cblas_ztrsv_unsafe :: TrsvFunFFI (Complex Double)+foreign import ccall safe   "cblas_ztrsv"+  cblas_ztrsv_safe   :: TrsvFunFFI (Complex Double)+ --void cblas_strsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Float *A,   CInt lda, Float *X,   CInt incX); --void cblas_dtrsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Double *A,   CInt lda, Double *X,   CInt incX); --void cblas_ctrsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Float *A,   CInt lda, Float *X,   CInt incX);@@ -600,6 +650,12 @@ --  |  triangular solvers  ----------------------- ++-- +--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 ->       CInt->CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO () foreign  import ccall unsafe "cblas_strsm"
+ src/Numerical/HBLAS/Constants.lhs view
@@ -0,0 +1,7 @@+\begin{code}+module Numerical.HBLAS.Constants where +import Data.Int ++flopsThreshold :: Int64+flopsThreshold = 10000+\end{code}
+ src/Numerical/HBLAS/Lapack.lhs view
@@ -0,0 +1,13 @@+\begin{code}+module  Numerical.HBLAS.Lapack where +    ++import Numerical.HBLAS.UtilsFFI    +import Numerical.HBLAS.Lapack.FFI +import Numerical.HBLAS.MatrixTypes+import Control.Monad.Primitive+import Data.Complex +import qualified Data.Vector.Storable.Mutable as SM+++\end{code}
+ src/Numerical/HBLAS/Lapack/FFI.hs view
@@ -0,0 +1,129 @@+++module Numerical.HBLAS.Lapack.FFI where+import Foreign.Ptr+import Foreign()+import Foreign.C.Types+import Data.Complex +import Data.Int ++++{-++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 +++-}++++{-+void LAPACK_dgesvx( char* fact, char* trans, lapack_int* n, lapack_int* nrhs,+                    double* a, lapack_int* lda, double* af, lapack_int* ldaf,+                    lapack_int* ipiv, char* equed, double* r, double* c,+                    double* b, lapack_int* ldb, double* x, lapack_int* ldx,+                    double* rcond, double* ferr, double* berr, double* work,+                    lapack_int* iwork, lapack_int *info );++++-}+++{-+    fortran FFI conventions!+-}++--type Stride_C =++newtype Fact_C = Fact_C CChar+newtype Trans_C = Trans_C CChar+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 -}-> +    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 {-rcond -}-> Ptr el {- ferr-} -> Ptr el {-berr-} -> Ptr el {-work-}->+    Ptr Int32 {-iwork -}-> Ptr Int32 {-info  -} -> IO () ++++{-++the prefixes mean s=single,d=double,c=complex float,d=complex double++++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 ++-}+++{-+Xgesvx  is the s -sing+++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 +-}++--need to get around to wrapping these, but thats for another day+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)++++--lapack_int ?syev_(  char *jobz, char *uplo, lapack_int *n, ?* a, lapack_int * lda, ?* w );+-- ? is Double or Float ++newtype JobTy = JBT CChar +newtype UploTy = UPLT CChar+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 ++{-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 ++--lapack_int LAPACKE_<?>gesv( int matrix_order, lapack_int n, lapack_int nrhs, <datatype>* a, lapack_int lda, lapack_int* ipiv, <datatype>* b, lapack_int ldb );+--call sgesv( n, nrhs, a, lda, ipiv, b, ldb, info )++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 +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 "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)++{- only provide  -}++++
src/Numerical/HBLAS/MatrixTypes.hs view
@@ -7,20 +7,7 @@   ----module Numerical.HBLAS.MatrixTypes where--import qualified Data.Vector.Storable as S -import qualified Data.Vector.Storable.Mutable as SM-import Control.Monad.Primitive  ---import Data.Singletons-import Control.Monad.ST.Safe -import Data.Typeable --- import Control.Monad.Primitive--{-| PSA, the matrix data types used in the hOpenBLAS binding+{-| 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 @@ -35,10 +22,27 @@ A guiding rule of thumb for this package is that there are no generic abstractions provided, merely machinery to ensure all uses of BLAS and LAPACK operations can be used in their full generality in a human friendly type safe fashion.-It is the role of a higher leve library to provide any generic operations.+It is the role of a higher level library to provide any generic operations. --}    +One such higher level lib you can interface with easily is Numerical.+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.Mutable as SM+import Control.Monad.Primitive  +--import Data.Singletons+import Control.Monad.ST.Safe +import Data.Typeable +-- import Control.Monad.Primitive++ + {- what I really want is this, but its not possible till datakinds works on types that aren't kind *,@@ -96,11 +100,28 @@ sTranpose SRow = SColumn  -+--- | May Blas and Lapack routines allow you to implicitly tranpose your argumments data Transpose = NoTranspose | Transpose | ConjTranspose | ConjNoTranspose     deriving(Typeable,Eq,Show)++-- | For Symmetric, Hermetian or Triangular matrices, which part is modeled.+---  Applies to both Padded and Packed variants+data MatUpLo = MatUpper | MatLower +    deriving(Typeable,Eq,Show)+++-- | 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    +    deriving(Typeable,Eq,Show) {--should think long and hard before adding implicit tranposition to the internal data model+should think long and hard before adding implicit transposition to the internal data model -}  type family TransposeF (x :: Orientation) :: Orientation@@ -108,17 +129,64 @@ type instance TransposeF Row = Column type instance TransposeF Column = Row  +++++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 +-- packed triangular row major matrix wrt both their logical size and manifest sizes+-- 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 ++data DenseVector :: Variant -> * -> * where+    DenseVector :: { _VariantDenseVect ::  !(SVariant varnt)+                    ,_LogicalDimDenseVector :: {-#UNPACK#-}!Int +                    ,_StrideDenseVector :: {-#UNPACK#-} ! Int +                    ,_bufferDenseVector :: !(S.Vector elem) +                      } -> DenseVector varnt elem +#if defined(__GLASGOW_HASKELL_) && (__GLASGOW_HASKELL__ >= 707)+  deriving (Typeable)+#endif                      ++data MDenseVector :: * -> Variant -> * -> * where+    MutableDenseVector :: { _VariantMutDenseVect ::  !(SVariant varnt)+                        ,_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                      + -- | 'DenseMatrix' is for dense row or column major matrices data DenseMatrix :: Orientation -> * -> *  where      DenseMatrix ::{ _OrientationMat :: SOrientation ornt ,                     _XdimDenMat :: {-# UNPACK #-}!Int,                      _YdimDenMat :: {-# UNPACK #-}!Int ,                     _StrideDenMat :: {-# UNPACK #-} !Int , -                    _bufferDenMat :: !(S.Vector elem) }-> DenseMatrix ornt  elem +                    _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. :) +mutableVectorToList :: (PrimMonad m, S.Storable a) => S.MVector (PrimState m) a -> m [a]+mutableVectorToList mv =  do+        v <- S.unsafeFreeze mv +        return (S.toList v )+{-# NOINLINE mutableVectorToList #-}+ {- need to handle rendering a slice differently than a direct matrix  -}@@ -133,12 +201,12 @@              | otherwise = show $ mapDenseMatrix id mat                           -- | 'MDenseMatrix' -data MutDenseMatrix :: *  ->Orientation  -> * -> *  where +data MDenseMatrix :: *  ->Orientation  -> * -> *  where      MutableDenseMatrix :: { _OrientationMutMat :: SOrientation ornt ,                             _XdimDenMutMat :: {-# UNPACK #-}!Int ,                              _YdimDenMutMat ::  {-# UNPACK #-}!Int,                             _StrideDenMutMat :: {-# UNPACK #-} !Int,-                            _bufferDenMutMat :: {-# UNPACK #-} !(SM.MVector s  elem) } -> MutDenseMatrix s ornt elem+                            _bufferDenMutMat :: {-# UNPACK #-} !(SM.MVector s  elem) } -> MDenseMatrix s ornt elem  --instance (Show el,SM.Storable el, PrimMonad m )=> Show (DenseMatrix (PrimState m) Row el) where --    show mat@(DenseMatrix SRow xdim ydim stride buffer)@@ -150,7 +218,7 @@ --             |  stride == xdim = "DenseMatrix SColumn " ++ " " ++show xdim ++ " " ++ show ydim ++ " " ++ show stride ++ "(" ++ show buffer ++ ")" --             | otherwise = show $ mapDenseMatrix id mat     -type IODenseMatrix = MutDenseMatrix RealWorld +type IODenseMatrix = MDenseMatrix RealWorld  --type MutDenseMatrixIO  or elem  =   -- data PaddedSymmetricMatrix@@ -163,14 +231,14 @@ --data TriangularMatrix --data BandedMatrix {-#NOINLINE unsafeFreezeDenseMatrix #-}-unsafeFreezeDenseMatrix :: (SM.Storable elem, PrimMonad m)=> MutDenseMatrix (PrimState m) or elem -> m (DenseMatrix or elem)+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                             {-# NOINLINE unsafeThawDenseMatrix #-}-unsafeThawDenseMatrix :: (SM.Storable elem, PrimMonad m)=> DenseMatrix or elem-> m (MutDenseMatrix (PrimState m) or elem)  +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 @@ -207,15 +275,15 @@ 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) -uncheckedMutDenseMatrixIndexM :: (PrimMonad m ,S.Storable elem )=>  MutDenseMatrix (PrimState m) or elem -> (Int,Int) -> m elem -uncheckedMutDenseMatrixIndexM (MutableDenseMatrix SRow _ _ ystride arr) =  \ (x,y)->  arr `SM.unsafeRead` (x + y * ystride)-uncheckedMutDenseMatrixIndexM (MutableDenseMatrix SColumn _ _ xstride arr) = \ (x,y)->   arr `SM.unsafeRead` (y + x* xstride)+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)  swap :: (a,b)->(b,a) swap = \ (!x,!y)-> (y,x) {-# INLINE swap #-} -+-- | `map f matrix` 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 $!@@ -240,7 +308,7 @@   -+-- | 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 SRow (xdim,ydim) f = DenseMatrix SRow  xdim ydim xdim $!              S.generate (xdim * ydim) (\ix -> let !ixtup@(!_,!_) = swap $ quotRem ix xdim in @@ -249,13 +317,23 @@          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  (MutDenseMatrix (PrimState m) x a) +    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  +{-#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   --- this (uncheckedMatrixSlice) will need to have its inlining quality checked
− testing/Test.hs 
@@ -1,30 +0,0 @@-{-# LANGUAGE ScopedTypeVariables, DataKinds, CPP #-}-import Numerical.HBLAS.BLAS.FFI-import Numerical.HBLAS.BLAS-import Numerical.HBLAS.MatrixTypes -import Data.Vector.Storable.Mutable as M -import qualified Data.Vector.Storable as S --main :: IO ()-main = do-  -- Just test that the symbol resolves-  --openblas_set_num_threads_unsafe 7  -  v  :: IOVector Double <- M.replicate 10 1.0-  res <- unsafeWith v (\ptr-> cblas_ddot_unsafe 10 ptr 1   ptr 1) -#if defined(__GLASGOW_HASKELL_) && (__GLASGOW_HASKELL__ <= 704)-  --this makes 7.4 panic! -  (leftMat :: IODenseMatrix Row Float) <-  generateMutableDenseMatrix SRow (5,5) (\_ -> 1.0 ::Float )-  (rightMat :: IODenseMatrix Row Float) <-  generateMutableDenseMatrix SRow (5,5) (\_ -> 1.0 ::Float )-  (resMat :: IODenseMatrix Row Float) <-  generateMutableDenseMatrix SRow (5,5) (\_ -> 1.0 ::Float )-  sgemm NoTranspose NoTranspose 1.0 1.0 leftMat rightMat resMat-  --(MutableDenseMatrix _ _ _ _ buff) <- return resMat -  theRestMat <- unsafeFreezeDenseMatrix resMat-  putStrLn $ show theRestMat-#endif---DenseMatrix SRow  5 5 5(fromList [11.0,11.0,11.0,11.0,11.0,11.0,11.0,11.0,11.0,11.0,51.0,51.0,51.0,51.0,51.0,51.0,51.0,51.0,51.0,51.0,251.0,251.0,251.0,251.0,251.0])--- THAT IS WRONG ---Need to figure what whats going on here   --  putStrLn $ show res -  putStrLn "it works!"-
+ tests/MainUnit.hs view
@@ -0,0 +1,29 @@+module Main where++ +import   UnitBLAS.Level1+import   UnitBLAS.Level2+import   UnitBLAS.Level3++import Test.Tasty+--import Test.Tasty.SmallCheck as SC+--import Test.Tasty.QuickCheck as QC+import Test.Tasty.HUnit++import Data.List+import Data.Ord++main = defaultMain tests++tests :: TestTree+tests = testGroup "BLAS Unit Tests"    [unitTestLevel1BLAS, unitTestLevel2BLAS, unitTestLevel3BLAS] -- [unitTestShape] -- , unitTestLayout ]+++--unitTests = testGroup "Unit tests"+--  [ testCase "List comparison (different length)" $+--      [1, 2, 3] `compare` [1,2] @?= GT++--  -- the following test does not hold+--  , testCase "List comparison (same length)" $+--      [1, 2, 3] `compare` [1,2,2] @?= LT+--  ]