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)) {
     OK
 }
 
+int multiplyI(int ta, int tb, KIMAT(a), KIMAT(b), IMAT(r)) {
+    int i,j,k;
+    int n;
+    int u, v;
+    if (ta==0) {
+        n = ac;
+    } else {
+        n = ar;
+    }
+    for (i=0;i<rr;i++) {
+        for (j=0;j<rc;j++) {
+            rp[i*rc+j] = 0;
+            for (k=0; k<n; k++) {
+                u = ta==0 ? ap[i*ac+k] : ap[k*ac+i];
+                v = tb==0 ? bp[k*bc+j] : bp[j*bc+k];
+                rp[i*rc+j] += u*v;
+            }
+        }
+    }
+    OK
+}
+
+
 //////////////////// transpose /////////////////////////
 
 int transF(KFMAT(x),FMAT(t)) {
@@ -1350,6 +1373,19 @@ int transP(KPMAT(x), PMAT(t)) {
     OK
 }
 
+int transI(KIMAT(x),IMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("transI");
+    int i,j;
+    for (i=0; i<tr; i++) {
+        for (j=0; j<tc; j++) {
+            tp[i*tc+j] = xp[j*xc+i];
+        }
+    }
+    OK
+}
+
+
 //////////////////// constant /////////////////////////
 
 int constantF(float * pval, FVEC(r)) {
@@ -1401,6 +1437,18 @@ int constantP(void* pval, PVEC(r)) {
     OK
 }
 
+
+int constantI(int * pval, IVEC(r)) {
+    DEBUGMSG("constantI")
+    int k;
+    int val = *pval;
+    for(k=0;k<rn;k++) {
+        rp[k]=val;
+    }
+    OK
+}
+
+
 //////////////////// float-double conversion /////////////////////////
 
 int float2double(FVEC(x),DVEC(y)) {
@@ -1412,6 +1460,16 @@ int float2double(FVEC(x),DVEC(y)) {
     OK
 }
 
+int float2int(KFVEC(x),IVEC(y)) {
+    DEBUGMSG("float2int")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k];
+    }
+    OK
+}
+
+
 int double2float(DVEC(x),FVEC(y)) {
     DEBUGMSG("double2float")
     int k;
@@ -1421,6 +1479,37 @@ int double2float(DVEC(x),FVEC(y)) {
     OK
 }
 
+
+int double2int(KDVEC(x),IVEC(y)) {
+    DEBUGMSG("double2int")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k];
+    }
+    OK
+}
+
+
+int int2float(KIVEC(x),FVEC(y)) {
+    DEBUGMSG("int2float")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k];
+    }
+    OK
+}
+
+
+int int2double(KIVEC(x),DVEC(y)) {
+    DEBUGMSG("int2double")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k];
+    }
+    OK
+}
+
+
 //////////////////// conjugate /////////////////////////
 
 int conjugateQ(KQVEC(x),QVEC(t)) {
@@ -1465,8 +1554,30 @@ int stepD(DVEC(x),DVEC(y)) {
     OK
 }
 
+
 //////////////////// cond /////////////////////////
 
+int compareF(KFVEC(x),KFVEC(y),IVEC(r)) {
+    REQUIRES(xn==yn && xn==rn ,BAD_SIZE);
+    DEBUGMSG("compareF")
+    int k;
+    for(k=0;k<xn;k++) {
+        rp[k] = xp[k]<yp[k]?-1:(xp[k]>yp[k]?1:0);
+    }
+    OK
+}
+
+int compareD(KDVEC(x),KDVEC(y),IVEC(r)) {
+    REQUIRES(xn==yn && xn==rn ,BAD_SIZE);
+    DEBUGMSG("compareD")
+    int k;
+    for(k=0;k<xn;k++) {
+        rp[k] = xp[k]<yp[k]?-1:(xp[k]>yp[k]?1:0);
+    }
+    OK
+}
+
+
 int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) {
     REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE);
     DEBUGMSG("condF")
@@ -1487,3 +1598,69 @@ int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) {
     OK
 }
 
+
+int chooseF(KIVEC(cond),KFVEC(lt),KFVEC(eq),KFVEC(gt),FVEC(r)) {
+    REQUIRES(condn==ltn && ltn==eqn && ltn==gtn && ltn==rn ,BAD_SIZE);
+    DEBUGMSG("chooseF")
+    int k;
+    for(k=0;k<condn;k++) {
+        rp[k] = condp[k]<0?ltp[k]:(condp[k]>0?gtp[k]:eqp[k]);
+    }
+    OK
+}
+
+
+int chooseD(KIVEC(cond),KDVEC(lt),KDVEC(eq),KDVEC(gt),DVEC(r)) {
+    REQUIRES(condn==ltn && ltn==eqn && ltn==gtn && ltn==rn ,BAD_SIZE);
+    DEBUGMSG("chooseD")
+    int k;
+    for(k=0;k<condn;k++) {
+        rp[k] = condp[k]<0?ltp[k]:(condp[k]>0?gtp[k]:eqp[k]);
+    }
+    OK
+}
+
+//////////////////////// extract /////////////////////////////////
+
+#define EXTRACT_IMP \
+    REQUIRES((tm == 0 && jn==rr && mc==rc) || (jn==rr && mr==rc) ,BAD_SIZE);   \
+    DEBUGMSG("extractRD")                                                      \
+    int k,i,s;                                                                 \
+    if (tm==0) {                                                               \
+        for (k=0;k<jn;k++) {                                                   \
+            s = jp[k];                                                         \
+            for (i=0; i<mc; i++) {                                             \
+                rp[rc*k+i] = mp[mc*s+i];                                       \
+            }                                                                  \
+        }                                                                      \
+    } else {                                                                   \
+        for (k=0;k<jn;k++) {                                                   \
+            s = jp[k];                                                         \
+            for (i=0; i<mr; i++) {                                             \
+                rp[rc*k+i] = mp[mc*i+s];                                       \
+            }                                                                  \
+        }                                                                      \
+    }                                                                          \
+    OK
+
+
+int extractRD(int tm, KIVEC(j), KDMAT(m), DMAT(r)) {
+    EXTRACT_IMP
+}
+
+int extractRF(int tm, KIVEC(j), KFMAT(m), FMAT(r)) {
+    EXTRACT_IMP
+}
+
+int extractRC(int tm, KIVEC(j), KCMAT(m), CMAT(r)) {
+    EXTRACT_IMP
+}
+
+int extractRQ(int tm, KIVEC(j), KQMAT(m), QMAT(r)) {
+    EXTRACT_IMP
+}
+
+int extractRI(int tm, KIVEC(j), KIMAT(m), IMAT(r)) {
+    EXTRACT_IMP
+}
+
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;
 #define QVEC(A) int A##n, complex*A##p
 #define CVEC(A) int A##n, doublecomplex*A##p
 #define PVEC(A) int A##n, void* A##p, int A##s
+#define IMAT(A) int A##r, int A##c, int* A##p
 #define FMAT(A) int A##r, int A##c, float* A##p
 #define DMAT(A) int A##r, int A##c, double* A##p
 #define QMAT(A) int A##r, int A##c, complex* A##p
@@ -54,6 +55,7 @@ typedef short ftnlen;
 #define KQVEC(A) int A##n, const complex*A##p
 #define KCVEC(A) int A##n, const doublecomplex*A##p
 #define KPVEC(A) int A##n, const void* A##p, int A##s
+#define KIMAT(A) int A##r, int A##c, const int* A##p
 #define KFMAT(A) int A##r, int A##c, const float* A##p
 #define KDMAT(A) int A##r, int A##c, const double* A##p
 #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)) {
     OK
 }
 
+////////////////////////////////////////////////////////////////////////////////
+
+int round_vector_i(KDVEC(v),IVEC(r)) {
+    int k;
+    for(k=0; k<vn; k++) {
+        rp[k] = round(vp[k]);
+    }
+    OK
+}
+
+
+int mod_vector(int m, KIVEC(v), IVEC(r)) {
+    int k;
+    for(k=0; k<vn; k++) {
+        rp[k] = vp[k] % m;
+    }
+    OK
+}
+
+int div_vector(int m, KIVEC(v), IVEC(r)) {
+    int k;
+    for(k=0; k<vn; k++) {
+        rp[k] = vp[k] / m;
+    }
+    OK
+}
+
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
     transdata = transdataP -- transdata'
     constantD  :: a -> Int -> Vector a
     constantD = constantP -- constant'
-
+    extractR :: Matrix a -> Idxs -> Matrix a
 
 instance Element Float where
     transdata  = transdataAux ctransF
     constantD  = constantAux cconstantF
+    extractR   = extractAux c_extractRF
 
 instance Element Double where
     transdata  = transdataAux ctransR
     constantD  = constantAux cconstantR
+    extractR   = extractAux c_extractRD
 
 instance Element (Complex Float) where
     transdata  = transdataAux ctransQ
     constantD  = constantAux cconstantQ
+    extractR   = extractAux c_extractRQ
 
 instance Element (Complex Double) where
     transdata  = transdataAux ctransC
     constantD  = constantAux cconstantC
+    extractR   = extractAux c_extractRC
+    
+instance Element (CInt) where
+    transdata  = transdataAux ctransI
+    constantD  = constantAux cconstantI
+    extractR   = extractAux c_extractRI
+
 
 -------------------------------------------------------------------
 
@@ -289,6 +299,7 @@ transdataAux fun c1 d c2 =
     if noneed
         then d
         else unsafePerformIO $ do
+            -- putStrLn "T"
             v <- createVector (dim d)
             unsafeWith d $ \pd ->
                 unsafeWith v $ \pv ->
@@ -317,6 +328,7 @@ foreign import ccall unsafe "transF" ctransF :: TFMFM
 foreign import ccall unsafe "transR" ctransR :: TMM
 foreign import ccall unsafe "transQ" ctransQ :: TQMQM
 foreign import ccall unsafe "transC" ctransC :: TCMCM
+foreign import ccall unsafe "transI" ctransI :: CM CInt (CM CInt (IO CInt))
 foreign import ccall unsafe "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt
 
 ----------------------------------------------------------------------
@@ -336,6 +348,8 @@ foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV
 
 foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV
 
+foreign import ccall unsafe "constantI" cconstantI :: Ptr CInt -> CV CInt (IO CInt)
+
 constantP :: Storable a => a -> Int -> Vector a
 constantP a n = unsafePerformIO $ do
     let sz = sizeOf a
@@ -421,3 +435,32 @@ instance (Storable t, NFData t) => NFData (Matrix t)
         d = dim v
         v = xdat m
 
+---------------------------------------------------------------
+
+isT Matrix{order = ColumnMajor} = 1
+isT Matrix{order = RowMajor} = 0
+
+tt x@Matrix{order = ColumnMajor} = trans x
+tt x@Matrix{order = RowMajor} = x
+
+--extractAux :: Matrix Double -> Idxs -> Matrix Double
+extractAux f m v = unsafePerformIO $ do
+    r <- createMatrix RowMajor (dim v) (cols m)
+    app3 (f (isT m)) vec v mat (tt m) mat r "extractAux"
+    return r
+
+foreign import ccall unsafe "extractRD" c_extractRD
+    :: CInt -> CIdxs (CM Double (CM Double (IO CInt)))
+
+foreign import ccall unsafe "extractRF" c_extractRF
+    :: CInt -> CIdxs (CM Float (CM Float (IO CInt)))
+
+foreign import ccall unsafe "extractRC" c_extractRC
+    :: CInt -> CIdxs (CM (Complex Double) (CM (Complex Double) (IO CInt)))
+
+foreign import ccall unsafe "extractRQ" c_extractRQ
+    :: CInt -> CIdxs (CM (Complex Float) (CM (Complex Float) (IO CInt)))
+
+foreign import ccall unsafe "extractRI" c_extractRI
+    :: CInt -> CIdxs (CM CInt (CM CInt (IO CInt)))
+
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 @@
 -- |
 -- Module      :  Data.Packed.Internal.Signatures
--- Copyright   :  (c) Alberto Ruiz 2009
+-- Copyright   :  (c) Alberto Ruiz 2009-15
 -- License     :  BSD3
 -- Maintainer  :  Alberto Ruiz
 -- Stability   :  provisional
@@ -68,3 +68,8 @@ type TCVM = CInt -> PC -> TM                    --
 type TMCVM = CInt -> CInt -> PD -> TCVM         --
 type TMMCVM = CInt -> CInt -> PD -> TMCVM       --
 
+type CM b r = CInt -> CInt -> Ptr b -> r
+type CV b r = CInt -> Ptr b -> r
+
+type CIdxs r = CV CInt r
+
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 (
     cloneVector,
     unsafeToForeignPtr,
     unsafeFromForeignPtr,
-    unsafeWith
+    unsafeWith,
+    Idxs
 ) where
 
 import Data.Packed.Internal.Common
@@ -56,6 +57,8 @@ import Data.Vector.Storable(Vector,
                             unsafeWith)
 
 
+type Idxs = Vector CInt
+
 -- | Number of elements
 dim :: (Storable t) => Vector t -> Int
 dim = Vector.length
@@ -243,6 +246,7 @@ double2FloatV v = unsafePerformIO $ do
 foreign import ccall unsafe "float2double" c_float2double:: TFV
 foreign import ccall unsafe "double2float" c_double2float:: TVF
 
+
 ---------------------------------------------------------------
 
 stepF :: Vector Float -> Vector Float