wellposed / hblas

haskell bindings for blas and lapack
www.wellposed.com
BSD 3-Clause "New" or "Revised" License
49 stars 19 forks source link

add rank 1 real matrix update / outer product #14

Open cartazio opened 10 years ago

cartazio commented 10 years ago

sger and dger

http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-BD2E87B3-5FA7-4E0C-88E2-1982AB0773A2.htm

https://github.com/wellposed/hblas/blob/master/src/Numerical/HBLAS/BLAS/FFI.hs#L247-L250

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

the haskelly high level version would look like (assuming you punt on dealing with modelling the vectors properly because i've not added strided vector support to hblas yet, and you assume inx and incy are set to 1)

sger ::( PrimMonad s) =>Float -> Float 
   -> Mvector (PrimState s) Float -> MVector (PrimState s) Float 
  ->  MutDenseMatrix (PrimState m) orient Float -> m () 

note this would actually be done slightly differently once I add strided vector support

a mini explanation of some of the API stuff I do currently is on this other ticket https://github.com/wellposed/hblas/issues/13

don't hesitate to ask for help

cartazio commented 10 years ago

could adapt the very dumb lacking test program https://github.com/wellposed/hblas/blob/master/testing/Test.hs to test this

adrian-nitu-92 commented 10 years ago

Small mistake : the function signature should be :

sger ::( PrimMonad s) =>
Float -> Mvector (PrimState s) Float -> MVector (PrimState s) Float 
  ->  MutDenseMatrix (PrimState m) orient Float -> m () 
cartazio commented 10 years ago

good catch! As specced in http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-BD2E87B3-5FA7-4E0C-88E2-1982AB0773A2.htm it should be A := alpha*x*y'+ A,

cartazio commented 9 years ago

note that this may not be implemented correctly in terms of certain corner cases, reopening for now and this ticket + associated audit needs to be done before i do the next release