hblas 0.3.0.1 → 0.3.1.0

raw patch · 10 files changed

+490/−432 lines, 10 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• hblas.cabal +53/−49
• readme.md +73/−22
• src/Numerical/HBLAS/BLAS.hs +31/−31
• src/Numerical/HBLAS/BLAS/FFI.hs +120/−120
• src/Numerical/HBLAS/BLAS/Internal.hs +42/−42
• src/Numerical/HBLAS/Constants.lhs +2/−2
• src/Numerical/HBLAS/Lapack.lhs +5/−5
• src/Numerical/HBLAS/Lapack/FFI.hs +22/−22
• src/Numerical/HBLAS/MatrixTypes.hs +114/−111
• src/Numerical/HBLAS/UtilsFFI.hs +28/−28

@@ -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 -    <https://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm Intel MKL BLAS and LAPACK documention>, +    Lapack, please read the+    <https://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm Intel MKL BLAS and LAPACK documention>,     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 -----  -    -  

@@ -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.

@@ -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

@@ -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)  

@@ -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)  

@@ -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

@@ -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  

@@ -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, <datatype>* a, lapack_int lda, lapack_int* ipiv, <datatype>* 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)

@@ -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)  

@@ -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