diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,5 @@ +import Distribution.Simple + + + +main = defaultMain diff --git a/hblas.cabal b/hblas.cabal new file mode 100644 --- /dev/null +++ b/hblas.cabal @@ -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 + + + + + + + diff --git a/readme.md b/readme.md new file mode 100644 --- /dev/null +++ b/readme.md @@ -0,0 +1,36 @@ +[![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). + + +[![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. + +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. diff --git a/src/Numerical/HBLAS/BLAS.hs b/src/Numerical/HBLAS/BLAS.hs new file mode 100644 --- /dev/null +++ b/src/Numerical/HBLAS/BLAS.hs @@ -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_ diff --git a/src/Numerical/HBLAS/BLAS/FFI.hs b/src/Numerical/HBLAS/BLAS/FFI.hs new file mode 100644 --- /dev/null +++ b/src/Numerical/HBLAS/BLAS/FFI.hs @@ -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, ...); diff --git a/src/Numerical/HBLAS/MatrixTypes.hs b/src/Numerical/HBLAS/MatrixTypes.hs new file mode 100644 --- /dev/null +++ b/src/Numerical/HBLAS/MatrixTypes.hs @@ -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≤ + +-} + + + diff --git a/src/Numerical/HBLAS/UtilsFFI.hs b/src/Numerical/HBLAS/UtilsFFI.hs new file mode 100644 --- /dev/null +++ b/src/Numerical/HBLAS/UtilsFFI.hs @@ -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 + diff --git a/testing/Test.hs b/testing/Test.hs new file mode 100644 --- /dev/null +++ b/testing/Test.hs @@ -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!" +