summaryrefslogtreecommitdiff
path: root/packages/base
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base')
-rw-r--r--packages/base/src/C/lapack-aux.c177
-rw-r--r--packages/base/src/C/lapack-aux.h2
-rw-r--r--packages/base/src/C/vector-aux.c27
-rw-r--r--packages/base/src/Data/Packed/Internal/Matrix.hs45
-rw-r--r--packages/base/src/Data/Packed/Internal/Signatures.hs7
-rw-r--r--packages/base/src/Data/Packed/Internal/Vector.hs6
6 files changed, 261 insertions, 3 deletions
diff --git a/packages/base/src/C/lapack-aux.c b/packages/base/src/C/lapack-aux.c
index e5e45ef..d56d466 100644
--- a/packages/base/src/C/lapack-aux.c
+++ b/packages/base/src/C/lapack-aux.c
@@ -1287,6 +1287,29 @@ int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) {
1287 OK 1287 OK
1288} 1288}
1289 1289
1290int multiplyI(int ta, int tb, KIMAT(a), KIMAT(b), IMAT(r)) {
1291 int i,j,k;
1292 int n;
1293 int u, v;
1294 if (ta==0) {
1295 n = ac;
1296 } else {
1297 n = ar;
1298 }
1299 for (i=0;i<rr;i++) {
1300 for (j=0;j<rc;j++) {
1301 rp[i*rc+j] = 0;
1302 for (k=0; k<n; k++) {
1303 u = ta==0 ? ap[i*ac+k] : ap[k*ac+i];
1304 v = tb==0 ? bp[k*bc+j] : bp[j*bc+k];
1305 rp[i*rc+j] += u*v;
1306 }
1307 }
1308 }
1309 OK
1310}
1311
1312
1290//////////////////// transpose ///////////////////////// 1313//////////////////// transpose /////////////////////////
1291 1314
1292int transF(KFMAT(x),FMAT(t)) { 1315int transF(KFMAT(x),FMAT(t)) {
@@ -1350,6 +1373,19 @@ int transP(KPMAT(x), PMAT(t)) {
1350 OK 1373 OK
1351} 1374}
1352 1375
1376int transI(KIMAT(x),IMAT(t)) {
1377 REQUIRES(xr==tc && xc==tr,BAD_SIZE);
1378 DEBUGMSG("transI");
1379 int i,j;
1380 for (i=0; i<tr; i++) {
1381 for (j=0; j<tc; j++) {
1382 tp[i*tc+j] = xp[j*xc+i];
1383 }
1384 }
1385 OK
1386}
1387
1388
1353//////////////////// constant ///////////////////////// 1389//////////////////// constant /////////////////////////
1354 1390
1355int constantF(float * pval, FVEC(r)) { 1391int constantF(float * pval, FVEC(r)) {
@@ -1401,6 +1437,18 @@ int constantP(void* pval, PVEC(r)) {
1401 OK 1437 OK
1402} 1438}
1403 1439
1440
1441int constantI(int * pval, IVEC(r)) {
1442 DEBUGMSG("constantI")
1443 int k;
1444 int val = *pval;
1445 for(k=0;k<rn;k++) {
1446 rp[k]=val;
1447 }
1448 OK
1449}
1450
1451
1404//////////////////// float-double conversion ///////////////////////// 1452//////////////////// float-double conversion /////////////////////////
1405 1453
1406int float2double(FVEC(x),DVEC(y)) { 1454int float2double(FVEC(x),DVEC(y)) {
@@ -1412,6 +1460,16 @@ int float2double(FVEC(x),DVEC(y)) {
1412 OK 1460 OK
1413} 1461}
1414 1462
1463int float2int(KFVEC(x),IVEC(y)) {
1464 DEBUGMSG("float2int")
1465 int k;
1466 for(k=0;k<xn;k++) {
1467 yp[k]=xp[k];
1468 }
1469 OK
1470}
1471
1472
1415int double2float(DVEC(x),FVEC(y)) { 1473int double2float(DVEC(x),FVEC(y)) {
1416 DEBUGMSG("double2float") 1474 DEBUGMSG("double2float")
1417 int k; 1475 int k;
@@ -1421,6 +1479,37 @@ int double2float(DVEC(x),FVEC(y)) {
1421 OK 1479 OK
1422} 1480}
1423 1481
1482
1483int double2int(KDVEC(x),IVEC(y)) {
1484 DEBUGMSG("double2int")
1485 int k;
1486 for(k=0;k<xn;k++) {
1487 yp[k]=xp[k];
1488 }
1489 OK
1490}
1491
1492
1493int int2float(KIVEC(x),FVEC(y)) {
1494 DEBUGMSG("int2float")
1495 int k;
1496 for(k=0;k<xn;k++) {
1497 yp[k]=xp[k];
1498 }
1499 OK
1500}
1501
1502
1503int int2double(KIVEC(x),DVEC(y)) {
1504 DEBUGMSG("int2double")
1505 int k;
1506 for(k=0;k<xn;k++) {
1507 yp[k]=xp[k];
1508 }
1509 OK
1510}
1511
1512
1424//////////////////// conjugate ///////////////////////// 1513//////////////////// conjugate /////////////////////////
1425 1514
1426int conjugateQ(KQVEC(x),QVEC(t)) { 1515int conjugateQ(KQVEC(x),QVEC(t)) {
@@ -1465,8 +1554,30 @@ int stepD(DVEC(x),DVEC(y)) {
1465 OK 1554 OK
1466} 1555}
1467 1556
1557
1468//////////////////// cond ///////////////////////// 1558//////////////////// cond /////////////////////////
1469 1559
1560int compareF(KFVEC(x),KFVEC(y),IVEC(r)) {
1561 REQUIRES(xn==yn && xn==rn ,BAD_SIZE);
1562 DEBUGMSG("compareF")
1563 int k;
1564 for(k=0;k<xn;k++) {
1565 rp[k] = xp[k]<yp[k]?-1:(xp[k]>yp[k]?1:0);
1566 }
1567 OK
1568}
1569
1570int compareD(KDVEC(x),KDVEC(y),IVEC(r)) {
1571 REQUIRES(xn==yn && xn==rn ,BAD_SIZE);
1572 DEBUGMSG("compareD")
1573 int k;
1574 for(k=0;k<xn;k++) {
1575 rp[k] = xp[k]<yp[k]?-1:(xp[k]>yp[k]?1:0);
1576 }
1577 OK
1578}
1579
1580
1470int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) { 1581int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) {
1471 REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE); 1582 REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE);
1472 DEBUGMSG("condF") 1583 DEBUGMSG("condF")
@@ -1487,3 +1598,69 @@ int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) {
1487 OK 1598 OK
1488} 1599}
1489 1600
1601
1602int chooseF(KIVEC(cond),KFVEC(lt),KFVEC(eq),KFVEC(gt),FVEC(r)) {
1603 REQUIRES(condn==ltn && ltn==eqn && ltn==gtn && ltn==rn ,BAD_SIZE);
1604 DEBUGMSG("chooseF")
1605 int k;
1606 for(k=0;k<condn;k++) {
1607 rp[k] = condp[k]<0?ltp[k]:(condp[k]>0?gtp[k]:eqp[k]);
1608 }
1609 OK
1610}
1611
1612
1613int chooseD(KIVEC(cond),KDVEC(lt),KDVEC(eq),KDVEC(gt),DVEC(r)) {
1614 REQUIRES(condn==ltn && ltn==eqn && ltn==gtn && ltn==rn ,BAD_SIZE);
1615 DEBUGMSG("chooseD")
1616 int k;
1617 for(k=0;k<condn;k++) {
1618 rp[k] = condp[k]<0?ltp[k]:(condp[k]>0?gtp[k]:eqp[k]);
1619 }
1620 OK
1621}
1622
1623//////////////////////// extract /////////////////////////////////
1624
1625#define EXTRACT_IMP \
1626 REQUIRES((tm == 0 && jn==rr && mc==rc) || (jn==rr && mr==rc) ,BAD_SIZE); \
1627 DEBUGMSG("extractRD") \
1628 int k,i,s; \
1629 if (tm==0) { \
1630 for (k=0;k<jn;k++) { \
1631 s = jp[k]; \
1632 for (i=0; i<mc; i++) { \
1633 rp[rc*k+i] = mp[mc*s+i]; \
1634 } \
1635 } \
1636 } else { \
1637 for (k=0;k<jn;k++) { \
1638 s = jp[k]; \
1639 for (i=0; i<mr; i++) { \
1640 rp[rc*k+i] = mp[mc*i+s]; \
1641 } \
1642 } \
1643 } \
1644 OK
1645
1646
1647int extractRD(int tm, KIVEC(j), KDMAT(m), DMAT(r)) {
1648 EXTRACT_IMP
1649}
1650
1651int extractRF(int tm, KIVEC(j), KFMAT(m), FMAT(r)) {
1652 EXTRACT_IMP
1653}
1654
1655int extractRC(int tm, KIVEC(j), KCMAT(m), CMAT(r)) {
1656 EXTRACT_IMP
1657}
1658
1659int extractRQ(int tm, KIVEC(j), KQMAT(m), QMAT(r)) {
1660 EXTRACT_IMP
1661}
1662
1663int extractRI(int tm, KIVEC(j), KIMAT(m), IMAT(r)) {
1664 EXTRACT_IMP
1665}
1666
diff --git a/packages/base/src/C/lapack-aux.h b/packages/base/src/C/lapack-aux.h
index c95a2a3..b49c7c9 100644
--- a/packages/base/src/C/lapack-aux.h
+++ b/packages/base/src/C/lapack-aux.h
@@ -42,6 +42,7 @@ typedef short ftnlen;
42#define QVEC(A) int A##n, complex*A##p 42#define QVEC(A) int A##n, complex*A##p
43#define CVEC(A) int A##n, doublecomplex*A##p 43#define CVEC(A) int A##n, doublecomplex*A##p
44#define PVEC(A) int A##n, void* A##p, int A##s 44#define PVEC(A) int A##n, void* A##p, int A##s
45#define IMAT(A) int A##r, int A##c, int* A##p
45#define FMAT(A) int A##r, int A##c, float* A##p 46#define FMAT(A) int A##r, int A##c, float* A##p
46#define DMAT(A) int A##r, int A##c, double* A##p 47#define DMAT(A) int A##r, int A##c, double* A##p
47#define QMAT(A) int A##r, int A##c, complex* A##p 48#define QMAT(A) int A##r, int A##c, complex* A##p
@@ -54,6 +55,7 @@ typedef short ftnlen;
54#define KQVEC(A) int A##n, const complex*A##p 55#define KQVEC(A) int A##n, const complex*A##p
55#define KCVEC(A) int A##n, const doublecomplex*A##p 56#define KCVEC(A) int A##n, const doublecomplex*A##p
56#define KPVEC(A) int A##n, const void* A##p, int A##s 57#define KPVEC(A) int A##n, const void* A##p, int A##s
58#define KIMAT(A) int A##r, int A##c, const int* A##p
57#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
58#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
59#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
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index abeba76..58afc49 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -894,3 +894,30 @@ int round_vector(KDVEC(v),DVEC(r)) {
894 OK 894 OK
895} 895}
896 896
897////////////////////////////////////////////////////////////////////////////////
898
899int round_vector_i(KDVEC(v),IVEC(r)) {
900 int k;
901 for(k=0; k<vn; k++) {
902 rp[k] = round(vp[k]);
903 }
904 OK
905}
906
907
908int mod_vector(int m, KIVEC(v), IVEC(r)) {
909 int k;
910 for(k=0; k<vn; k++) {
911 rp[k] = vp[k] % m;
912 }
913 OK
914}
915
916int div_vector(int m, KIVEC(v), IVEC(r)) {
917 int k;
918 for(k=0; k<vn; k++) {
919 rp[k] = vp[k] / m;
920 }
921 OK
922}
923
diff --git a/packages/base/src/Data/Packed/Internal/Matrix.hs b/packages/base/src/Data/Packed/Internal/Matrix.hs
index 150b978..be5fb03 100644
--- a/packages/base/src/Data/Packed/Internal/Matrix.hs
+++ b/packages/base/src/Data/Packed/Internal/Matrix.hs
@@ -265,23 +265,33 @@ class (Storable a) => Element a where
265 transdata = transdataP -- transdata' 265 transdata = transdataP -- transdata'
266 constantD :: a -> Int -> Vector a 266 constantD :: a -> Int -> Vector a
267 constantD = constantP -- constant' 267 constantD = constantP -- constant'
268 268 extractR :: Matrix a -> Idxs -> Matrix a
269 269
270instance Element Float where 270instance Element Float where
271 transdata = transdataAux ctransF 271 transdata = transdataAux ctransF
272 constantD = constantAux cconstantF 272 constantD = constantAux cconstantF
273 extractR = extractAux c_extractRF
273 274
274instance Element Double where 275instance Element Double where
275 transdata = transdataAux ctransR 276 transdata = transdataAux ctransR
276 constantD = constantAux cconstantR 277 constantD = constantAux cconstantR
278 extractR = extractAux c_extractRD
277 279
278instance Element (Complex Float) where 280instance Element (Complex Float) where
279 transdata = transdataAux ctransQ 281 transdata = transdataAux ctransQ
280 constantD = constantAux cconstantQ 282 constantD = constantAux cconstantQ
283 extractR = extractAux c_extractRQ
281 284
282instance Element (Complex Double) where 285instance Element (Complex Double) where
283 transdata = transdataAux ctransC 286 transdata = transdataAux ctransC
284 constantD = constantAux cconstantC 287 constantD = constantAux cconstantC
288 extractR = extractAux c_extractRC
289
290instance Element (CInt) where
291 transdata = transdataAux ctransI
292 constantD = constantAux cconstantI
293 extractR = extractAux c_extractRI
294
285 295
286------------------------------------------------------------------- 296-------------------------------------------------------------------
287 297
@@ -289,6 +299,7 @@ transdataAux fun c1 d c2 =
289 if noneed 299 if noneed
290 then d 300 then d
291 else unsafePerformIO $ do 301 else unsafePerformIO $ do
302 -- putStrLn "T"
292 v <- createVector (dim d) 303 v <- createVector (dim d)
293 unsafeWith d $ \pd -> 304 unsafeWith d $ \pd ->
294 unsafeWith v $ \pv -> 305 unsafeWith v $ \pv ->
@@ -317,6 +328,7 @@ foreign import ccall unsafe "transF" ctransF :: TFMFM
317foreign import ccall unsafe "transR" ctransR :: TMM 328foreign import ccall unsafe "transR" ctransR :: TMM
318foreign import ccall unsafe "transQ" ctransQ :: TQMQM 329foreign import ccall unsafe "transQ" ctransQ :: TQMQM
319foreign import ccall unsafe "transC" ctransC :: TCMCM 330foreign import ccall unsafe "transC" ctransC :: TCMCM
331foreign import ccall unsafe "transI" ctransI :: CM CInt (CM CInt (IO CInt))
320foreign import ccall unsafe "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt 332foreign import ccall unsafe "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt
321 333
322---------------------------------------------------------------------- 334----------------------------------------------------------------------
@@ -336,6 +348,8 @@ foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV
336 348
337foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV 349foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV
338 350
351foreign import ccall unsafe "constantI" cconstantI :: Ptr CInt -> CV CInt (IO CInt)
352
339constantP :: Storable a => a -> Int -> Vector a 353constantP :: Storable a => a -> Int -> Vector a
340constantP a n = unsafePerformIO $ do 354constantP a n = unsafePerformIO $ do
341 let sz = sizeOf a 355 let sz = sizeOf a
@@ -421,3 +435,32 @@ instance (Storable t, NFData t) => NFData (Matrix t)
421 d = dim v 435 d = dim v
422 v = xdat m 436 v = xdat m
423 437
438---------------------------------------------------------------
439
440isT Matrix{order = ColumnMajor} = 1
441isT Matrix{order = RowMajor} = 0
442
443tt x@Matrix{order = ColumnMajor} = trans x
444tt x@Matrix{order = RowMajor} = x
445
446--extractAux :: Matrix Double -> Idxs -> Matrix Double
447extractAux f m v = unsafePerformIO $ do
448 r <- createMatrix RowMajor (dim v) (cols m)
449 app3 (f (isT m)) vec v mat (tt m) mat r "extractAux"
450 return r
451
452foreign import ccall unsafe "extractRD" c_extractRD
453 :: CInt -> CIdxs (CM Double (CM Double (IO CInt)))
454
455foreign import ccall unsafe "extractRF" c_extractRF
456 :: CInt -> CIdxs (CM Float (CM Float (IO CInt)))
457
458foreign import ccall unsafe "extractRC" c_extractRC
459 :: CInt -> CIdxs (CM (Complex Double) (CM (Complex Double) (IO CInt)))
460
461foreign import ccall unsafe "extractRQ" c_extractRQ
462 :: CInt -> CIdxs (CM (Complex Float) (CM (Complex Float) (IO CInt)))
463
464foreign import ccall unsafe "extractRI" c_extractRI
465 :: CInt -> CIdxs (CM CInt (CM CInt (IO CInt)))
466
diff --git a/packages/base/src/Data/Packed/Internal/Signatures.hs b/packages/base/src/Data/Packed/Internal/Signatures.hs
index acc3070..37dac16 100644
--- a/packages/base/src/Data/Packed/Internal/Signatures.hs
+++ b/packages/base/src/Data/Packed/Internal/Signatures.hs
@@ -1,6 +1,6 @@
1-- | 1-- |
2-- Module : Data.Packed.Internal.Signatures 2-- Module : Data.Packed.Internal.Signatures
3-- Copyright : (c) Alberto Ruiz 2009 3-- Copyright : (c) Alberto Ruiz 2009-15
4-- License : BSD3 4-- License : BSD3
5-- Maintainer : Alberto Ruiz 5-- Maintainer : Alberto Ruiz
6-- Stability : provisional 6-- Stability : provisional
@@ -68,3 +68,8 @@ type TCVM = CInt -> PC -> TM --
68type TMCVM = CInt -> CInt -> PD -> TCVM -- 68type TMCVM = CInt -> CInt -> PD -> TCVM --
69type TMMCVM = CInt -> CInt -> PD -> TMCVM -- 69type TMMCVM = CInt -> CInt -> PD -> TMCVM --
70 70
71type CM b r = CInt -> CInt -> Ptr b -> r
72type CV b r = CInt -> Ptr b -> r
73
74type CIdxs r = CV CInt r
75
diff --git a/packages/base/src/Data/Packed/Internal/Vector.hs b/packages/base/src/Data/Packed/Internal/Vector.hs
index b49f379..2a6ed2c 100644
--- a/packages/base/src/Data/Packed/Internal/Vector.hs
+++ b/packages/base/src/Data/Packed/Internal/Vector.hs
@@ -24,7 +24,8 @@ module Data.Packed.Internal.Vector (
24 cloneVector, 24 cloneVector,
25 unsafeToForeignPtr, 25 unsafeToForeignPtr,
26 unsafeFromForeignPtr, 26 unsafeFromForeignPtr,
27 unsafeWith 27 unsafeWith,
28 Idxs
28) where 29) where
29 30
30import Data.Packed.Internal.Common 31import Data.Packed.Internal.Common
@@ -56,6 +57,8 @@ import Data.Vector.Storable(Vector,
56 unsafeWith) 57 unsafeWith)
57 58
58 59
60type Idxs = Vector CInt
61
59-- | Number of elements 62-- | Number of elements
60dim :: (Storable t) => Vector t -> Int 63dim :: (Storable t) => Vector t -> Int
61dim = Vector.length 64dim = Vector.length
@@ -243,6 +246,7 @@ double2FloatV v = unsafePerformIO $ do
243foreign import ccall unsafe "float2double" c_float2double:: TFV 246foreign import ccall unsafe "float2double" c_float2double:: TFV
244foreign import ccall unsafe "double2float" c_double2float:: TVF 247foreign import ccall unsafe "double2float" c_double2float:: TVF
245 248
249
246--------------------------------------------------------------- 250---------------------------------------------------------------
247 251
248stepF :: Vector Float -> Vector Float 252stepF :: Vector Float -> Vector Float