hblas (empty) → 0.1.0.0
raw patch · 9 files changed
+1411/−0 lines, 9 filesdep +basedep +hblasdep +primitivebuild-type:Customsetup-changed
Dependencies added: base, hblas, primitive, storable-complex, vector
Files
- LICENSE +30/−0
- Setup.hs +5/−0
- hblas.cabal +131/−0
- readme.md +36/−0
- src/Numerical/HBLAS/BLAS.hs +138/−0
- src/Numerical/HBLAS/BLAS/FFI.hs +666/−0
- src/Numerical/HBLAS/MatrixTypes.hs +295/−0
- src/Numerical/HBLAS/UtilsFFI.hs +80/−0
- testing/Test.hs +30/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2013, Carter Tazio Schonwald++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Carter Tazio Schonwald nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,5 @@+import Distribution.Simple++++main = defaultMain
+ hblas.cabal view
@@ -0,0 +1,131 @@+-- 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) +-- 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.1.0.0++-- A short (one-line) description of the package.+synopsis: BLAS and Lapack bindings for OpenBLAS++-- A longer description of the package.+description: User friendly, simple bindings to BLAS and Lapack, favoring OpenBLAS as the substrate.++-- The license under which the package is released.+license: BSD3++-- The file containing the license text.+license-file: LICENSE++-- The package author(s).+author: Carter Tazio Schonwald++-- An email address to which users can send suggestions, bug reports, and +-- patches.+maintainer: carter at wellposed dot com++-- A copyright notice.+-- copyright: +cabal-version: >=1.10+category: Math++build-type: Custom+++extra-source-files: + readme.md +++source-repository head+ type: git+ location: git://github.com/wellposed/hOpenBLAS.git+ +-- Extra files to be distributed with the package, such as examples or a +-- README.+++-- Constraint on the version of Cabal needed to build this package.+-- flag lib-Werror+-- default: False+-- manual: True+-- flag lib-Wall +-- default: False+-- manual: True+++library++ -- Modules exported by the library.+ exposed-modules: Numerical.HBLAS.BLAS.FFI,+ Numerical.HBLAS.MatrixTypes,+ Numerical.HBLAS.UtilsFFI+ Numerical.HBLAS.BLAS+ -- Modules included in this library but not exported.+ -- other-modules: + + -- LANGUAGE extensions used by modules in this package.+ -- other-extensions: + --- if flag(lib-Werror)+ --- ghc-options: -Werror+ --- if flag(lib-Wall)+ --- ghc-options: -Wall + -- -Wall no wall for now+ -- 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)++ frameworks: Accelerate+ -- probably dont need these extra lib dirs+ -- extra-libraries: cblas clapack++ if os(windows)+ extra-libraries: blas lapack++ if !os(windows)&& !os(OSX)+ 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+++++ + +
+ readme.md view
@@ -0,0 +1,36 @@+[](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). +++[](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. ++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. +++## how to install+* On OS X systems, things will just work.+* On linux and bsd systems, the equivalent of +```+sudo apt-get install libblas liblapack+```+is all you should have to do before hand++## getting involved+patches, bug reports, tests, and other contributions welcome.
+ src/Numerical/HBLAS/BLAS.hs view
@@ -0,0 +1,138 @@+{-# LANGUAGE BangPatterns , RankNTypes, GADTs #-}+++module Numerical.HBLAS.BLAS where+--(+-- dgemm+-- ,sgemm+-- ,cgemm+-- ,zgemm) +++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++--+flopsThreshold = 10000+gemmComplexity a b c = a * b * c -- this will be wrong by some constant factor, albeit a small one+++-- 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 + 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 Transpose (a,b) = (b,a)+coordSwapper ConjTranspose (a,b) = (b,a)+++encodeNiceOrder :: SOrientation x -> CBLAS_ORDERT+encodeNiceOrder SRow= encodeOrder BLASRowMajor+encodeNiceOrder SColumn= encodeOrder BLASColMajor+++encodeFFITranpose :: Transpose -> CBLAS_TRANSPOSET+encodeFFITranpose x= encodeTranpose $ encodeNiceTranpose x ++encodeNiceTranpose :: Transpose -> BLAS_Transpose+encodeNiceTranpose 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+++type GemmFun el orient s m = Transpose ->Transpose -> el -> el -> MutDenseMatrix s orient el+ -> MutDenseMatrix s orient el -> MutDenseMatrix s orient 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++-}+++---- | 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-> + --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 + 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) + | 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 = + {- 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 -> + do (ax,ay) <- 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 $ encodeFFITranpose tra + rawTrb <- return $ encodeFFITranpose 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) + alphaPtr ap (fromIntegral astride) bp (fromIntegral bstride) betaPtr cp (fromIntegral cstride)+++{-pureGemm :: PrimMonad m=>+(Transpose ->Transpose -> el -> el -> MutDenseMatrix (PrimState m) orient el+ -> MutDenseMatrix (PrimState m) orient el -> + MutDenseMatrix (PrimState m) orient el -> m ())->+ Transpose ->Transpose -> el -> el -> DenseMatrix orient el+ -> 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 ) + ++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 ) + ++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_ ++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_
+ src/Numerical/HBLAS/BLAS/FFI.hs view
@@ -0,0 +1,666 @@++{-# LANGUAGE ForeignFunctionInterface, GeneralizedNewtypeDeriving #-}+++module Numerical.HBLAS.BLAS.FFI where ++import Foreign.Ptr+import Foreign()+import Foreign.C.Types+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, +++-}+{-+++typedef enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102} CBLAS_ORDER;+typedef enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113, CblasConjNoTrans=114} CBLAS_TRANSPOSE;+typedef enum CBLAS_UPLO {CblasUpper=121, CblasLower=122} CBLAS_UPLO;+typedef enum CBLAS_DIAG {CblasNonUnit=131, CblasUnit=132} CBLAS_DIAG;+typedef enum CBLAS_SIDE {CblasLeft=141, CblasRight=142} CBLAS_SIDE;++-}++newtype CBLAS_INDEX = CBIndex CSize + deriving (Eq,Show)+newtype CBLAS_ORDERT = CBOInt CInt+ deriving (Eq,Show)+data BLASOrder = BLASRowMajor | BLASColMajor + deriving (Eq,Show)++encodeOrder :: BLASOrder -> CBLAS_ORDERT+encodeOrder BLASRowMajor = CBOInt 101+encodeOrder BLASColMajor = CBOInt 102 ++newtype CBLAS_TRANSPOSET = CBLAS_TransposeT{ unCBLAS_TransposeT :: CInt } deriving (Eq, Show)++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++newtype CBLAS_UPLOT = CBlasUPLO CInt + deriving (Eq,Show)+data BLASUplo = BUpper | BLower+ deriving (Eq,Show)++encodeUPLO :: BLASUplo -> CBLAS_UPLOT +encodeUPLO BUpper = CBlasUPLO 121 +encodeUPLO BLower = CBlasUPLO 122++newtype CBLAS_DIAGT = CBLAS_DiagT CUChar + deriving (Show,Eq)+++data BlasDiag = BlasNonUnit | BlasUnit + deriving (Eq,Show )+ +encodeDiag :: BlasDiag -> CBLAS_DIAGT +encodeDiag BlasNonUnit = CBLAS_DiagT 131+encodeDiag BlasUnit = CBLAS_DiagT 132++newtype CBLAS_SIDET = CBLAS_SideT { unCBLAS_SideT :: CUChar } + deriving (Eq, Show)+data BlasSide = BlasLeft | BlasRight + deriving (Eq,Show)+encodeSide :: BlasSide -> CBLAS_SIDET +encodeSide BlasLeft = CBLAS_SideT 141+encodeSide BlasRight = CBLAS_SideT 142+++--typedef enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102} CBLAS_ORDER;+--typedef enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113, CblasConjNoTrans=114} CBLAS_TRANSPOSE;+--typedef enum CBLAS_UPLO {CblasUpper=121, CblasLower=122} CBLAS_UPLO;+--typedef enum CBLAS_DIAG {CblasNonUnit=131, CblasUnit=132} CBLAS_DIAG;+--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_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+--Float cblas_sdsdot( CInt n, Float alpha, Float *x, CInt incx, Float *y, CInt incy);+--Double cblas_dsdot ( CInt n, Float *x, CInt incx, Float *y, CInt incy);+--Float cblas_sdot( CInt n, Float *x, CInt incx, Float *y, CInt incy);+--Double cblas_ddot( CInt n, Double *x, CInt incx, Double *y, CInt incy);++++{-+not doing these right now, because requires handling return value as a complex number,+we can only handle pointers to complex numbers right now+-}+--openblas_complex_Float cblas_cdotu( CInt n, Float *x, CInt incx, Float *y, CInt incy);+--openblas_complex_Float cblas_cdotc( CInt n, Float *x, CInt incx, Float *y, CInt incy);+--openblas_complex_Double cblas_zdotu( CInt n, Double *x, CInt incx, Double *y, CInt incy);+--openblas_complex_Double cblas_zdotc( CInt n, Double *x, CInt incx, Double *y, CInt incy);+++++--- 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:: + CInt -> Ptr Float -> CInt -> IO Float+foreign import ccall unsafe "cblas_dasum" cblas_dasum_unsafe :: + CInt -> Ptr Double -> CInt -> IO Double+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 :: + 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);+--Float cblas_scasum( CInt n, Float *x, CInt incx);+--Double cblas_dzasum( CInt n, Double *x, CInt incx);+++foreign import ccall unsafe "cblas_snrm2" cblas_snrm2_unsafe :: + CInt -> Ptr Float -> CInt -> IO Float+foreign import ccall unsafe "cblas_dnrm2" cblas_dnrm2_unsafe :: + CInt -> Ptr Double -> CInt -> IO Double+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 :: + CInt -> Ptr (Complex Double) -> CInt -> IO Double+++--Float cblas_snrm2 ( CInt N, Float *X, CInt incX);+--Double cblas_dnrm2 ( CInt N, Double *X, CInt incX);+--Float cblas_scnrm2( CInt N, Float *X, CInt incX);+--Double cblas_dznrm2( CInt N, Double *X, CInt incX);++++foreign import ccall unsafe "cblas_isamax" cblas_isamax_unsafe :: + CInt -> Ptr Float -> CInt -> IO CInt+foreign import ccall unsafe "cblas_idamax" cblas_idamax_unsafe :: + CInt -> Ptr Float -> CInt -> IO CInt+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 +--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);+--CBLAS_INDEX cblas_izamax( CInt n, Double *x, CInt incx);+++++--void cblas_saxpy( CInt n, Float alpha, Float *x, CInt incx, Float *y, CInt incy);+--void cblas_daxpy( CInt n, Double alpha, Double *x, CInt incx, Double *y, CInt incy);+--void cblas_caxpy( CInt n, Float *alpha, Float *x, CInt incx, Float *y, CInt incy);+--void cblas_zaxpy( CInt n, Double *alpha, Double *x, CInt incx, Double *y, CInt incy);++++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 :: + 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 () ++--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);+--void cblas_ccopy( CInt n, Float *x, CInt incx, Float *y, CInt incy);+--void cblas_zcopy( CInt n, Double *x, CInt incx, Double *y, CInt incy);++++foreign import ccall unsafe "cblas_scopy" cblas_sswap_unsafe :: + CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO () +foreign import ccall unsafe "cblas_dcopy" cblas_dswap_unsafe :: + CInt -> Ptr Double-> CInt -> Ptr Double -> CInt -> IO ()+foreign import ccall unsafe "cblas_ccopy" cblas_cswap_unsafe :: + CInt -> Ptr (Complex Float) -> CInt -> Ptr (Complex Float) -> CInt -> IO ()+foreign import ccall unsafe "cblas_zcopy" 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);+--void cblas_cswap( CInt n, Float *x, CInt incx, Float *y, CInt incy);+--void cblas_zswap( CInt n, Double *x, CInt incx, Double *y, CInt incy);++++--void cblas_srot( CInt N, Float *X, CInt incX, Float *Y, CInt incY, Float c, Float s);+--void cblas_drot( CInt N, Double *X, CInt incX, Double *Y, CInt incY, Double c, Double s);++--void cblas_srotg(Float *a, Float *b, Float *c, Float *s);+--void cblas_drotg(Double *a, Double *b, Double *c, Double *s);++--void cblas_srotm( CInt N, Float *X, CInt incX, Float *Y, CInt incY, Float *P);+--void cblas_drotm( CInt N, Double *X, CInt incX, Double *Y, CInt incY, Double *P);++--void cblas_srotmg(Float *d1, Float *d2, Float *b1, Float b2, Float *P);+--void cblas_drotmg(Double *d1, Double *d2, Double *b1, Double b2, Double *P);++--void cblas_sscal( CInt N, Float alpha, Float *X, CInt incX);+--void cblas_dscal( CInt N, Double alpha, Double *X, CInt incX);+--void cblas_cscal( CInt N, Float *alpha, Float *X, CInt incX);+--void cblas_zscal( CInt N, Double *alpha, Double *X, CInt incX);+--void cblas_csscal( CInt N, Float alpha, Float *X, CInt incX);+--void cblas_zdscal( CInt N, Double alpha, Double *X, CInt incX);+++++---------------------+----- BLAS LEVEL 2+---------------------++{-+matrix vector product for general matrices+ perform one of the matrix-vector operations y :=+ alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,+-}++--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,+-- Double alpha, Double *a, CInt lda, Double *x, CInt incx, Double beta, Double *y, CInt incy);+--void cblas_cgemv( 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_zgemv( enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE trans, CInt m, CInt n,+-- Double *alpha, Double *a, CInt lda, Double *x, CInt incx, Double *beta, Double *y, CInt incy);+++-- perform the rank 1 operation A := alpha*x*y' + A,++--void cblas_sger ( enum CBLAS_ORDER order, CInt M, CInt N, Float alpha, Float *X, CInt incX, Float *Y, CInt incY, Float *A, CInt lda);+--void cblas_dger ( enum CBLAS_ORDER order, CInt M, CInt N, Double alpha, Double *X, CInt incX, Double *Y, CInt incY, Double *A, CInt lda);+--void cblas_cgeru( enum CBLAS_ORDER order, CInt M, CInt N, Float *alpha, Float *X, CInt incX, Float *Y, CInt incY, Float *A, CInt lda);+--void cblas_cgerc( enum CBLAS_ORDER order, CInt M, CInt N, Float *alpha, Float *X, CInt incX, Float *Y, CInt incY, Float *A, CInt lda);+--void cblas_zgeru( enum CBLAS_ORDER order, CInt M, CInt N, Double *alpha, Double *X, CInt incX, Double *Y, CInt incY, Double *A, CInt lda);+--void cblas_zgerc( enum CBLAS_ORDER order, CInt M, CInt N, Double *alpha, Double *X, CInt incX, Double *Y, CInt incY, Double *A, CInt lda);+++--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++--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);+--void cblas_ztrsv( 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_strmv( 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_dtrmv( 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_ctrmv( 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_ztrmv( 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_ssyr( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Float alpha, Float *X, CInt incX, Float *A, CInt lda);+--void cblas_dsyr( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Double alpha, Double *X, CInt incX, Double *A, CInt lda);+--void cblas_cher( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Float alpha, Float *X, CInt incX, Float *A, CInt lda);+--void cblas_zher( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Double alpha, Double *X, CInt incX, Double *A, CInt lda);++++--void cblas_ssyr2( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Float alpha, Float *X,+-- CInt incX, Float *Y, CInt incY, Float *A, CInt lda);+--void cblas_dsyr2( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Double alpha, Double *X,+-- CInt incX, Double *Y, CInt incY, Double *A, CInt lda);+--void cblas_cher2( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Float *alpha, Float *X, CInt incX,+-- Float *Y, CInt incY, Float *A, CInt lda);+--void cblas_zher2( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Double *alpha, Double *X, CInt incX,+-- Double *Y, CInt incY, Double *A, CInt lda);+++++--void cblas_sgbmv( enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE TransA, CInt M, CInt N,+-- CInt KL, CInt KU, Float alpha, Float *A, CInt lda, Float *X, CInt incX, Float beta, Float *Y, CInt incY);+--void cblas_dgbmv( enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE TransA, CInt M, CInt N,+-- CInt KL, CInt KU, Double alpha, Double *A, CInt lda, Double *X, CInt incX, Double beta, Double *Y, CInt incY);+--void cblas_cgbmv( enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE TransA, CInt M, CInt N,+-- CInt KL, CInt KU, Float *alpha, Float *A, CInt lda, Float *X, CInt incX, Float *beta, Float *Y, CInt incY);+--void cblas_zgbmv( enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE TransA, CInt M, CInt N,+-- CInt KL, CInt KU, Double *alpha, Double *A, CInt lda, Double *X, CInt incX, Double *beta, Double *Y, CInt incY);+++--void cblas_ssbmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, CInt K, Float alpha, Float *A,+-- CInt lda, Float *X, CInt incX, Float beta, Float *Y, CInt incY);+--void cblas_dsbmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, CInt K, Double alpha, Double *A,+-- CInt lda, Double *X, CInt incX, Double beta, Double *Y, CInt incY);++++--void cblas_stbmv( 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);+--void cblas_dtbmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag,+-- CInt N, CInt K, Double *A, CInt lda, Double *X, CInt incX);+--void cblas_ctbmv( 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);+--void cblas_ztbmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag,+-- 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 +----------------+--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);+--void cblas_dtbsv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag,+-- CInt N, CInt K, Double *A, CInt lda, Double *X, CInt incX);+--void cblas_ctbsv( 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);+--void cblas_ztbsv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag,+-- CInt N, CInt K, Double *A, CInt lda, Double *X, CInt incX);++-------------------------------------------------------------------------+-- | matrix vector product Av, writes result into v, where A is a packed triangular nxn matrix+-------------------------------------------------------------------------+--void cblas_stpmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag,+-- CInt N, Float *Ap, Float *X, CInt incX);+--void cblas_dtpmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag,+-- CInt N, Double *Ap, Double *X, CInt incX);+--void cblas_ctpmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag,+-- CInt N, Float *Ap, Float *X, CInt incX);+--void cblas_ztpmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag,+-- 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. +--------------------------------------------------+--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);+--void cblas_dtpsv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag,+-- CInt N, Double *Ap, Double *X, CInt incX);+--void cblas_ctpsv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag,+-- CInt N, Float *Ap, Float *X, CInt incX);+--void cblas_ztpsv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA, enum CBLAS_DIAG Diag,+-- CInt N, Double *Ap, Double *X, CInt incX);+++----------------------------------+---- | (unpacked) symmetric matrix vector product x:=Av, writes result x into v+---------------------------------++type SymvFunFFI el = CBLAS_ORDERT -> CBLAS_UPLOT -> CInt -> el -> Ptr el -> CInt ->+ 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);+--void cblas_dsymv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Double alpha, Double *A,+-- CInt lda, Double *X, CInt incX, Double beta, Double *Y, CInt incY);++--------------------------------+---- | hermitian matrix vector product x:=Av, writes result x into v+--------------------------------+--void cblas_chemv( 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);+--void cblas_zhemv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Double *alpha, Double *A,+-- CInt lda, Double *X, CInt incX, Double *beta, Double *Y, CInt incY);++++---------------+--- | packed symmetric matrix * vector product y:= alpha * Av + beta * y+---------------++--void cblas_sspmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Float alpha, Float *Ap,+-- Float *X, CInt incX, Float beta, Float *Y, CInt incY);+--void cblas_dspmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Double alpha, Double *Ap,+-- Double *X, CInt incX, Double beta, Double *Y, CInt incY);+++++--void cblas_sspr( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Float alpha, Float *X, CInt incX, Float *Ap);+--void cblas_dspr( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Double alpha, Double *X, CInt incX, Double *Ap);++--void cblas_chpr( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Float alpha, Float *X, CInt incX, Float *A);+--void cblas_zhpr( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Double alpha, Double *X, CInt incX, Double *A);++++--void cblas_sspr2( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Float alpha, Float *X, CInt incX, Float *Y, CInt incY, Float *A);+--void cblas_dspr2( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Double alpha, Double *X, CInt incX, Double *Y, CInt incY, Double *A);+--void cblas_chpr2( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Float *alpha, Float *X, CInt incX, Float *Y, CInt incY, Float *Ap);+--void cblas_zhpr2( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, Double *alpha, Double *X, CInt incX, Double *Y, CInt incY, Double *Ap);+++--void cblas_chbmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, CInt K,+-- Float *alpha, Float *A, CInt lda, Float *X, CInt incX, Float *beta, Float *Y, CInt incY);+--void cblas_zhbmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N, CInt K,+-- Double *alpha, Double *A, CInt lda, Double *X, CInt incX, Double *beta, Double *Y, CInt incY);+++--void cblas_chpmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N,+-- Float *alpha, Float *Ap, Float *X, CInt incX, Float *beta, Float *Y, CInt incY);+--void cblas_zhpmv( enum CBLAS_ORDER order, enum CBLAS_UPLO Uplo, CInt N,+-- Double *alpha, Double *Ap, Double *X, CInt incX, Double *beta, Double *Y, CInt incY);+++{-+++-}++--------------------------------------------------------------------------------+--------------------------------------------------------------------------------+------------------------------ | BLAS LEVEL 3 ROUTINES +--------------------------------------------------------------------------------+----------------------- | Level 3 ops are faster than Levels 1 or 2 +--------------------------------------------------------------------------------+++--void cblas_sgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const blasint M, const blasint N, const blasint K,+-- const float alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float beta, float *C, const blasint ldc);++--void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const blasint M, 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);+--void cblas_cgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const blasint M, const blasint N, const blasint K,+-- const float *alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float *beta, float *C, const blasint ldc);+--void cblas_zgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const blasint M, 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);++-- | 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-> + scale -> {- C -} Ptr el -> CInt -> IO ()++{- C := alpha*op( A )*op( B ) + beta*C , -}++-- matrix mult!+foreign import ccall unsafe "cblas_sgemm" + cblas_sgemm_unsafe :: GemmFunFFI Float Float++foreign import ccall unsafe "cblas_dgemm" + cblas_dgemm_unsafe :: GemmFunFFI Double Double++foreign import ccall unsafe "cblas_cgemm" + cblas_cgemm_unsafe :: GemmFunFFI (Ptr(Complex Float)) (Complex Float)++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" + cblas_sgemm_safe :: GemmFunFFI Float Float++foreign import ccall "cblas_dgemm" + cblas_dgemm_safe :: GemmFunFFI Double Double++foreign import ccall "cblas_cgemm" + cblas_cgemm_safe :: GemmFunFFI (Ptr(Complex Float)) (Complex Float)++foreign import ccall "cblas_zgemm" + cblas_zgemm_safe :: GemmFunFFI (Ptr (Complex Double)) (Complex Double)++-----------------------------------------+----- | Matrix mult for Symmetric Matrices+-----------------------------------------+++--void cblas_ssymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side, const enum CBLAS_UPLO Uplo, const blasint M, const blasint N,+-- const float alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float beta, float *C, const blasint ldc);+--void cblas_dsymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side, const enum CBLAS_UPLO Uplo, const blasint M, const blasint N,+-- const double alpha, const double *A, const blasint lda, const double *B, const blasint ldb, const double beta, double *C, const blasint ldc);+--void cblas_csymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side, const enum CBLAS_UPLO Uplo, const blasint M, const blasint N,+-- const float *alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float *beta, float *C, const blasint ldc);+--void cblas_zsymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side, const enum CBLAS_UPLO Uplo, const blasint M, const blasint N,+-- const double *alpha, const double *A, const blasint lda, const double *B, const blasint ldb, const double *beta, double *C, const blasint ldc);++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_dsymm" + cblas_dsymm_unsafe :: SymmFunFFI Double Double ++foreign import ccall unsafe "cblas_csymm" + cblas_csymm_unsafe :: SymmFunFFI (Ptr (Complex Float )) (Complex Float)++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 ++foreign import ccall "cblas_dsymm" + cblas_dsymm_safe :: SymmFunFFI Double Double ++foreign import ccall "cblas_csymm" + cblas_csymm_safe :: SymmFunFFI (Ptr (Complex Float )) (Complex Float)++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 +------------------------------------+++--void cblas_ssyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,+-- const blasint N, const blasint K, const float alpha, const float *A, const blasint lda, const float beta, float *C, const blasint ldc);+--void cblas_dsyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,+-- const blasint N, const blasint K, const double alpha, const double *A, const blasint lda, const double beta, double *C, const blasint ldc);+--void cblas_csyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,+-- const blasint N, const blasint K, const float *alpha, const float *A, const blasint lda, const float *beta, float *C, const blasint ldc);+--void cblas_zsyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,+-- const blasint N, const blasint K, const double *alpha, const double *A, const blasint lda, const double *beta, double *C, const blasint ldc);++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" + cblas_dsyrk_unsafe :: SyrkFunFFI Double Double+foreign import ccall unsafe "cblas_csyrk" + cblas_csyrk_unsafe :: SyrkFunFFI (Ptr(Complex Float)) (Complex Float)+foreign import ccall unsafe "cblas_zsyrk" + cblas_zsyrk_unsafe :: SyrkFunFFI (Ptr(Complex Double)) (Complex Double)++----------------------+----- | Symmetric Rank 2k matrix update, C= alpha* A*B' + alpha* B*A' + beta * C+----- or C= alpha* A'*B + alpha* B'*A + beta * C+-------------------+++--void cblas_ssyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,+-- const blasint N, const blasint K, const float alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float beta, float *C, const blasint ldc);+--void cblas_dsyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,+-- 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);+--void cblas_csyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,+-- const blasint N, const blasint K, const float *alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float *beta, float *C, const blasint ldc);+--void cblas_zsyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,+ --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 -> + scale ->Ptr el -> CInt -> IO ()++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 ++++-------------------------------+-------- | matrix matrix product for triangular matrices+------------------------------++++++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) ++--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);+--void cblas_dtrmm( enum CBLAS_ORDER Order, enum CBLAS_SIDE Side, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA,+-- enum CBLAS_DIAG Diag, CInt M, CInt N, Double alpha, Double *A, CInt lda, Double *B, CInt ldb);+--void cblas_ctrmm( 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);+--void cblas_ztrmm( enum CBLAS_ORDER Order, enum CBLAS_SIDE Side, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA,+-- enum CBLAS_DIAG Diag, CInt M, CInt N, Double *alpha, Double *A, CInt lda, Double *B, CInt ldb);++------------------------+-- | triangular solvers +-----------------------++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) ++--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);+--void cblas_dtrsm( enum CBLAS_ORDER Order, enum CBLAS_SIDE Side, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA,+-- enum CBLAS_DIAG Diag, CInt M, CInt N, Double alpha, Double *A, CInt lda, Double *B, CInt ldb);+--void cblas_ctrsm( 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);+--void cblas_ztrsm( enum CBLAS_ORDER Order, enum CBLAS_SIDE Side, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE TransA,+-- enum CBLAS_DIAG Diag, CInt M, CInt N, Double *alpha, Double *A, CInt lda, Double *B, CInt ldb);++-------------------------+-- | hermitian matrix mult+------------------------++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) ++--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);+--void cblas_zhemm( enum CBLAS_ORDER Order, enum CBLAS_SIDE Side, enum CBLAS_UPLO Uplo, CInt M, CInt N,+-- Double *alpha, Double *A, CInt lda, Double *B, CInt ldb, Double *beta, Double *C, CInt ldc);++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) +--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,+-- Double alpha, Double *A, CInt lda, Double beta, Double *C, CInt ldc);++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" + cblas_zher2k_unsafe :: Her2kFunFFI Double (Complex Double)+++--void cblas_cher2k( enum CBLAS_ORDER Order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE Trans, CInt N, CInt K,+-- Float *alpha, Float *A, CInt lda, Float *B, CInt ldb, Float beta, Float *C, CInt ldc);+--void cblas_zher2k( enum CBLAS_ORDER Order, enum CBLAS_UPLO Uplo, enum CBLAS_TRANSPOSE Trans, CInt N, CInt K,+-- Double *alpha, Double *A, CInt lda, Double *B, CInt ldb, Double beta, Double *C, CInt ldc);++----void cblas_xerbla(CInt p, char *rout, char *form, ...);
+ src/Numerical/HBLAS/MatrixTypes.hs view
@@ -0,0 +1,295 @@++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE GADTs,ScopedTypeVariables, PolyKinds,FlexibleInstances,DeriveDataTypeable #-}+++++++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+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 +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.++The primary and hence default format is Dense Row and Column Major Matrices,+but support will be added for other formats that BLAS and LAPACK provide operations for.++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.++-} ++{-+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 ++-}++--data Eff s where+-- 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 Eff = Pure | Mut ++--data EVector :: Eff -> * -> * -> * where +-- PureVector :: S.Vector el -> EVector Pure () el +-- MutVector :: SM.MVector s el -> EVector Mut s el +++data Orientation = Row | Column + deriving (Eq,Show,Typeable)++data SOrientation :: Orientation -> * where+ SRow :: SOrientation Row + SColumn :: SOrientation Column +#if defined(__GLASGOW_HASKELL_) && (__GLASGOW_HASKELL__ >= 707)+ deriving (Typeable)+#endif ++instance Show (SOrientation Row) where+ show !a = "SRow"+instance Show (SOrientation Column) where+ show !a = "SColumn" +instance Eq (SOrientation Row) where+ (==) !a !b = True +instance Eq (SOrientation Column) where+ (==) !a !b = True +++sTranpose :: (x~ TransposeF y, y~TransposeF x ) =>SOrientation x -> SOrientation y +sTranpose SColumn = SRow+sTranpose SRow = SColumn++++data Transpose = NoTranspose | Transpose | ConjTranspose | ConjNoTranspose+ deriving(Typeable,Eq,Show)+{-+should think long and hard before adding implicit tranposition to the internal data model+-}++type family TransposeF (x :: Orientation) :: Orientation++type instance TransposeF Row = Column+type instance TransposeF Column = Row ++-- | '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 +#if defined(__GLASGOW_HASKELL_) && (__GLASGOW_HASKELL__ >= 707)+ deriving (Typeable)+#endif++{-+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 ++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 ++-- | 'MDenseMatrix' +data MutDenseMatrix :: * ->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++--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 ++--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 ++type IODenseMatrix = MutDenseMatrix RealWorld +--type MutDenseMatrixIO or elem = ++-- data PaddedSymmetricMatrix+-- data PaddedHermetianMatrix+--data PaddedTriangularMatrix +--- these three may just be wrappers for general dense matrices ++--data SymmetricMatrix+--data HermitianMatrix -- may just be a wrapper for symmetric?+--data TriangularMatrix+--data BandedMatrix+{-#NOINLINE unsafeFreezeDenseMatrix #-}+unsafeFreezeDenseMatrix :: (SM.Storable elem, PrimMonad m)=> MutDenseMatrix (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 (DenseMatrix ornt a b c v) = do + mv <- S.unsafeThaw v+ return $! MutableDenseMatrix ornt a b c mv +++--freezeDenseMatrix + ++getDenseMatrixRow :: DenseMatrix or elem -> Int+getDenseMatrixRow (DenseMatrix _ _ ydim _ _)= ydim+++getDenseMatrixColumn :: DenseMatrix or elem -> Int+getDenseMatrixColumn (DenseMatrix _ xdim _ _ _)= xdim+++getDenseMatrixLeadingDimStride :: DenseMatrix or elem -> Int +getDenseMatrixLeadingDimStride (DenseMatrix _ _ _ stride _ ) = stride+++getDenseMatrixArray :: DenseMatrix or elem -> S.Vector elem +getDenseMatrixArray (DenseMatrix _ _ _ _ arr) = arr++getDenseMatrixOrientation :: DenseMatrix or elem -> SOrientation or +getDenseMatrixOrientation m = _OrientationMat m ++++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 (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)++swap :: (a,b)->(b,a)+swap = \ (!x,!y)-> (y,x)+{-# INLINE swap #-}+++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 _ _) = + 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 _ _) = + 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 + --- dont need the swap for column major+++++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 ) +generateDenseMatrix SColumn (xdim,ydim) f = DenseMatrix SColumn xdim ydim ydim $!+ S.generate (xdim * ydim ) (\ix -> let ixtup@(!_,!_) = ( quotRem ix ydim ) in + f ixtup ) ++{-# NOINLINE generateMutableDenseMatrix #-}+generateMutableDenseMatrix :: (S.Storable a,PrimMonad m)=> + SOrientation x -> (Int,Int)->((Int,Int)-> a) -> m (MutDenseMatrix (PrimState m) x a) +generateMutableDenseMatrix sor dims fun = do+ x <- unsafeThawDenseMatrix $! generateDenseMatrix sor dims fun + return x +++--- 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 (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) +uncheckedDenseMatrixSlice (DenseMatrix SColumn _ ydim xstride arr) (xstart,ystart) (xend,yend) = res+ 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. +-- 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 (DenseMatrix orient x y stride arr)= (DenseMatrix (sTranpose orient) y x stride arr)++++{-+need an init with index/ map with index, etc utils≤++-}+++
+ src/Numerical/HBLAS/UtilsFFI.hs view
@@ -0,0 +1,80 @@+{-# LANGUAGE Trustworthy #-}+{- VERY TRUST WORTHY :) -}+module Numerical.HBLAS.UtilsFFI where++++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 Foreign.Ptr++{-+the IO version of these various utils is in Base.+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 #-} +++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 ++ 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 #-} ++touchForeignPtrPrim ::PrimMonad m => ForeignPtr a -> m ()+touchForeignPtrPrim fp = unsafePrimToPrim $! touchForeignPtr fp+{-# NOINLINE touchForeignPtrPrim #-}+++unsafeWithPrim ::( Storable a, PrimMonad m )=> MVector (PrimState m) a -> (Ptr a -> m b) -> m b+{-# INLINE unsafeWithPrim #-}+unsafeWithPrim (MVector _ fp) fun = withForeignPtrPrim fp fun+++unsafeWithPrimLen ::( Storable a, PrimMonad m )=> MVector (PrimState m) a -> ((Ptr a, Int )-> m b) -> m b+{-# INLINE unsafeWithPrimLen #-}+unsafeWithPrimLen (MVector n fp ) fun = withForeignPtrPrim fp (\x -> fun (x,n))+++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 + res <- withForeignPtrPrim fp (\x -> fun x)+ 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 + res <- withForeignPtrPrim fp (\x -> fun (x,n))+ touchForeignPtrPrim fp + return res +
+ testing/Test.hs view
@@ -0,0 +1,30 @@+{-# 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!"+