diff options
-rw-r--r-- | lib/Data/Packed/Internal/Matrix.hs | 33 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c | 22 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h | 4 |
3 files changed, 57 insertions, 2 deletions
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index 7a17ef0..090826d 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs | |||
@@ -252,10 +252,11 @@ class (Storable a, Floating a) => Element a where | |||
252 | -> Matrix a -> Matrix a | 252 | -> Matrix a -> Matrix a |
253 | subMatrixD = subMatrix' | 253 | subMatrixD = subMatrix' |
254 | transdata :: Int -> Vector a -> Int -> Vector a | 254 | transdata :: Int -> Vector a -> Int -> Vector a |
255 | transdata = transdata' | 255 | transdata = transdataP -- transdata' |
256 | constantD :: a -> Int -> Vector a | 256 | constantD :: a -> Int -> Vector a |
257 | constantD = constant' | 257 | constantD = constantP -- constant' |
258 | conjugateD :: Vector a -> Vector a | 258 | conjugateD :: Vector a -> Vector a |
259 | conjugateD = id | ||
259 | 260 | ||
260 | instance Element Float where | 261 | instance Element Float where |
261 | transdata = transdataAux ctransF | 262 | transdata = transdataAux ctransF |
@@ -320,10 +321,27 @@ transdataAux fun c1 d c2 = | |||
320 | r2 = dim d `div` c2 | 321 | r2 = dim d `div` c2 |
321 | noneed = r1 == 1 || c1 == 1 | 322 | noneed = r1 == 1 || c1 == 1 |
322 | 323 | ||
324 | transdataP :: Storable a => Int -> Vector a -> Int -> Vector a | ||
325 | transdataP c1 d c2 = | ||
326 | if noneed | ||
327 | then d | ||
328 | else unsafePerformIO $ do | ||
329 | v <- createVector (dim d) | ||
330 | unsafeWith d $ \pd -> | ||
331 | unsafeWith v $ \pv -> | ||
332 | ctransP (fi r1) (fi c1) (castPtr pd) (fi sz) (fi r2) (fi c2) (castPtr pv) (fi sz) // check "transdataStorable" | ||
333 | return v | ||
334 | where r1 = dim d `div` c1 | ||
335 | r2 = dim d `div` c2 | ||
336 | sz = sizeOf (d @> 0) | ||
337 | noneed = r1 == 1 || c1 == 1 | ||
338 | |||
323 | foreign import ccall "transF" ctransF :: TFMFM | 339 | foreign import ccall "transF" ctransF :: TFMFM |
324 | foreign import ccall "transR" ctransR :: TMM | 340 | foreign import ccall "transR" ctransR :: TMM |
325 | foreign import ccall "transQ" ctransQ :: TQMQM | 341 | foreign import ccall "transQ" ctransQ :: TQMQM |
326 | foreign import ccall "transC" ctransC :: TCMCM | 342 | foreign import ccall "transC" ctransC :: TCMCM |
343 | foreign import ccall "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt | ||
344 | |||
327 | ---------------------------------------------------------------------- | 345 | ---------------------------------------------------------------------- |
328 | 346 | ||
329 | constant' v n = unsafePerformIO $ do | 347 | constant' v n = unsafePerformIO $ do |
@@ -359,6 +377,17 @@ constantC :: Complex Double -> Int -> Vector (Complex Double) | |||
359 | constantC = constantAux cconstantC | 377 | constantC = constantAux cconstantC |
360 | foreign import ccall "constantC" cconstantC :: Ptr (Complex Double) -> TCV | 378 | foreign import ccall "constantC" cconstantC :: Ptr (Complex Double) -> TCV |
361 | 379 | ||
380 | constantP :: Storable a => a -> Int -> Vector a | ||
381 | constantP a n = unsafePerformIO $ do | ||
382 | let sz = sizeOf a | ||
383 | v <- createVector n | ||
384 | unsafeWith v $ \p -> do | ||
385 | alloca $ \k -> do | ||
386 | poke k a | ||
387 | cconstantP (castPtr k) (fi n) (castPtr p) (fi sz) // check "constantP" | ||
388 | return v | ||
389 | foreign import ccall "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt | ||
390 | |||
362 | --------------------------------------- | 391 | --------------------------------------- |
363 | 392 | ||
364 | conjugateAux fun x = unsafePerformIO $ do | 393 | conjugateAux fun x = unsafePerformIO $ do |
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c index 09979cd..e8bbbdb 100644 --- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c +++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c | |||
@@ -1139,6 +1139,19 @@ int transC(KCMAT(x),CMAT(t)) { | |||
1139 | OK | 1139 | OK |
1140 | } | 1140 | } |
1141 | 1141 | ||
1142 | int transP(KPMAT(x), PMAT(t)) { | ||
1143 | REQUIRES(xr==tc && xc==tr,BAD_SIZE); | ||
1144 | REQUIRES(xs==ts,NOCONVER); | ||
1145 | DEBUGMSG("transP"); | ||
1146 | int i,j; | ||
1147 | for (i=0; i<tr; i++) { | ||
1148 | for (j=0; j<tc; j++) { | ||
1149 | memcpy(tp+(i*tc+j)*xs,xp +(j*xc+i)*xs,xs); | ||
1150 | } | ||
1151 | } | ||
1152 | OK | ||
1153 | } | ||
1154 | |||
1142 | //////////////////// constant ///////////////////////// | 1155 | //////////////////// constant ///////////////////////// |
1143 | 1156 | ||
1144 | int constantF(float * pval, FVEC(r)) { | 1157 | int constantF(float * pval, FVEC(r)) { |
@@ -1181,6 +1194,15 @@ int constantC(doublecomplex* pval, CVEC(r)) { | |||
1181 | OK | 1194 | OK |
1182 | } | 1195 | } |
1183 | 1196 | ||
1197 | int constantP(void* pval, PVEC(r)) { | ||
1198 | DEBUGMSG("constantP") | ||
1199 | int k; | ||
1200 | for(k=0;k<rn;k++) { | ||
1201 | memcpy(rp+k*rs,pval,rs); | ||
1202 | } | ||
1203 | OK | ||
1204 | } | ||
1205 | |||
1184 | //////////////////// float-double conversion ///////////////////////// | 1206 | //////////////////// float-double conversion ///////////////////////// |
1185 | 1207 | ||
1186 | int float2double(FVEC(x),DVEC(y)) { | 1208 | int float2double(FVEC(x),DVEC(y)) { |
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h index 3de5b19..12dcde4 100644 --- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h +++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h | |||
@@ -44,19 +44,23 @@ typedef short ftnlen; | |||
44 | #define DVEC(A) int A##n, double*A##p | 44 | #define DVEC(A) int A##n, double*A##p |
45 | #define QVEC(A) int A##n, complex*A##p | 45 | #define QVEC(A) int A##n, complex*A##p |
46 | #define CVEC(A) int A##n, doublecomplex*A##p | 46 | #define CVEC(A) int A##n, doublecomplex*A##p |
47 | #define PVEC(A) int A##n, void* A##p, int A##s | ||
47 | #define FMAT(A) int A##r, int A##c, float* A##p | 48 | #define FMAT(A) int A##r, int A##c, float* A##p |
48 | #define DMAT(A) int A##r, int A##c, double* A##p | 49 | #define DMAT(A) int A##r, int A##c, double* A##p |
49 | #define QMAT(A) int A##r, int A##c, complex* A##p | 50 | #define QMAT(A) int A##r, int A##c, complex* A##p |
50 | #define CMAT(A) int A##r, int A##c, doublecomplex* A##p | 51 | #define CMAT(A) int A##r, int A##c, doublecomplex* A##p |
52 | #define PMAT(A) int A##r, int A##c, void* A##p, int A##s | ||
51 | 53 | ||
52 | #define KFVEC(A) int A##n, const float*A##p | 54 | #define KFVEC(A) int A##n, const float*A##p |
53 | #define KDVEC(A) int A##n, const double*A##p | 55 | #define KDVEC(A) int A##n, const double*A##p |
54 | #define KQVEC(A) int A##n, const complex*A##p | 56 | #define KQVEC(A) int A##n, const complex*A##p |
55 | #define KCVEC(A) int A##n, const doublecomplex*A##p | 57 | #define KCVEC(A) int A##n, const doublecomplex*A##p |
58 | #define KPVEC(A) int A##n, const void* A##p, int A##s | ||
56 | #define KFMAT(A) int A##r, int A##c, const float* A##p | 59 | #define KFMAT(A) int A##r, int A##c, const float* A##p |
57 | #define KDMAT(A) int A##r, int A##c, const double* A##p | 60 | #define KDMAT(A) int A##r, int A##c, const double* A##p |
58 | #define KQMAT(A) int A##r, int A##c, const complex* A##p | 61 | #define KQMAT(A) int A##r, int A##c, const complex* A##p |
59 | #define KCMAT(A) int A##r, int A##c, const doublecomplex* A##p | 62 | #define KCMAT(A) int A##r, int A##c, const doublecomplex* A##p |
63 | #define KPMAT(A) int A##r, int A##c, const void* A##p, int A##s | ||
60 | 64 | ||
61 | /********************************************************/ | 65 | /********************************************************/ |
62 | 66 | ||