4 files changed, 98 insertions, 104 deletions

diff --git a/packages/base/src/C/lapack-aux.c b/packages/base/src/C/lapack-aux.c
index e76d31e..c2cdc62 100644
--- a/packages/base/src/C/lapack-aux.c
+++ b/packages/base/src/C/lapack-aux.c
@@ -1290,19 +1290,18 @@ int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) {
 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;
-    }
+    int ai,ak,bk,bj;
+
+    n = ta ? ar : ac;
+    
+    if (ta==0) { ai =  1; ak = ar; } else { ai = ar; ak =  1; }
+    if (tb==0) { bk =  1; bj = br; } else { bk = br; bj =  1; }
+    
     for (i=0;i<rr;i++) {
         for (j=0;j<rc;j++) {
-            rp[i*rc+j] = 0;
+            rp[i+rr*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;
+                rp[i+rr*j] += ap[ai*i+ak*k] * bp[bk*k+bj*j];
             }
         }
     }
@@ -1622,65 +1621,44 @@ int chooseD(KIVEC(cond),KDVEC(lt),KDVEC(eq),KDVEC(gt),DVEC(r)) {
 
 //////////////////////// 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) {                                                               \
-        if (mode==0) {                           \
-            for (k=0; k<jp[1]-jp[0]+1; k++) {                                                   \
-                s = k + jp[0];                                                         \
-                printf("%d\n",s); \
-                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<mc; i++) {                                             \
-                    rp[rc*k+i] = mp[mc*s+i];                                       \
-                }                                                                  \
-            }  \
-        }                                                               \
-    } else {                                                                   \
-        if (mode==0) {                           \
-            for (k=0; k<jp[1]-jp[0]+1; k++) {                                                   \
-                s = k + jp[0];                                                         \
-                printf("%d\n",s); \
-                for (i=0; i<mr; i++) {                                             \
-                    rp[rc*k+i] = mp[mc*i+s];                                       \
-                }                                                                  \
-            }                           \
-        } else {        \
-            for (k=0;k<jn;k++) {                                                   \
-                s = jp[k];                                                         \
-                for (i=0; i<mr; i++) {                                             \
-                    rp[rc*k+i] = mp[mc*i+s];                                       \
-                }                                                                  \
-            }                                                                      \
-        }  \
-    }                                                                          \
+#define EXTRACT_IMP                        \
+    int i,j,si,sj,ni,nj,ai,aj;             \
+    if (tm==0) {                           \
+        ai=mc; aj=1;                       \
+    } else {                               \
+        ai=1, aj=mc;                       \
+    }                                      \
+    ni = modei ? in : ip[1]-ip[0]+1;       \
+    nj = modej ? jn : jp[1]-jp[0]+1;       \
+                                           \
+    for (i=0; i<ni; i++) {                 \
+        si = modei ? ip[i] : i+ip[0];      \
+                                           \
+        for (j=0; j<nj; j++) {             \
+            sj = modej ? jp[j] : j+jp[0];  \
+                                           \
+            rp[rc*i+j] = mp[ai*si+aj*sj];  \
+        }                                  \
+    }                                      \
     OK
 
-
-int extractRD(int mode, int tm, KIVEC(j), KDMAT(m), DMAT(r)) {
+int extractD(int modei, int modej, int tm, KIVEC(i), KIVEC(j), KDMAT(m), DMAT(r)) {
     EXTRACT_IMP
 }
 
-int extractRF(int mode, int tm, KIVEC(j), KFMAT(m), FMAT(r)) {
+int extractF(int modei, int modej, int tm, KIVEC(i), KIVEC(j), KFMAT(m), FMAT(r)) {
     EXTRACT_IMP
 }
 
-int extractRC(int mode, int tm, KIVEC(j), KCMAT(m), CMAT(r)) {
+int extractC(int modei, int modej, int tm, KIVEC(i), KIVEC(j), KCMAT(m), CMAT(r)) {
     EXTRACT_IMP
 }
 
-int extractRQ(int mode, int tm, KIVEC(j), KQMAT(m), QMAT(r)) {
+int extractQ(int modei, int modej, int tm, KIVEC(i), KIVEC(j), KQMAT(m), QMAT(r)) {
     EXTRACT_IMP
 }
 
-int extractRI(int mode, int tm, KIVEC(j), KIMAT(m), IMAT(r)) {
+int extractI(int modei, int modej, int tm, KIVEC(i), KIVEC(j), KIMAT(m), IMAT(r)) {
     EXTRACT_IMP
 }
 
diff --git a/packages/base/src/Data/Packed/Internal/Matrix.hs b/packages/base/src/Data/Packed/Internal/Matrix.hs
index 1aee7d3..76d2204 100644
--- a/packages/base/src/Data/Packed/Internal/Matrix.hs
+++ b/packages/base/src/Data/Packed/Internal/Matrix.hs
@@ -265,32 +265,32 @@ class (Storable a) => Element a where
     transdata = transdataP -- transdata'
     constantD  :: a -> Int -> Vector a
     constantD = constantP -- constant'
-    extractR :: Matrix a -> CInt -> Idxs -> Matrix a
+    extractR :: Matrix a -> CInt -> Idxs -> CInt -> Idxs -> Matrix a
 
 instance Element Float where
     transdata  = transdataAux ctransF
     constantD  = constantAux cconstantF
-    extractR   = extractAux c_extractRF
+    extractR   = extractAux c_extractF
 
 instance Element Double where
     transdata  = transdataAux ctransR
     constantD  = constantAux cconstantR
-    extractR   = extractAux c_extractRD
+    extractR   = extractAux c_extractD
 
 instance Element (Complex Float) where
     transdata  = transdataAux ctransQ
     constantD  = constantAux cconstantQ
-    extractR   = extractAux c_extractRQ
+    extractR   = extractAux c_extractQ
 
 instance Element (Complex Double) where
     transdata  = transdataAux ctransC
     constantD  = constantAux cconstantC
-    extractR   = extractAux c_extractRC
+    extractR   = extractAux c_extractC
     
 instance Element (CInt) where
     transdata  = transdataAux ctransI
     constantD  = constantAux cconstantI
-    extractR   = extractAux c_extractRI
+    extractR   = extractAux c_extractI
 
 
 -------------------------------------------------------------------
@@ -443,26 +443,26 @@ 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 mode v = unsafePerformIO $ do
-    let nr | mode == 0 = fromIntegral $ max 0 (v@>1 - v@>0 + 1)
-           | otherwise = dim v
-    r <- createMatrix RowMajor nr (cols m)
-    app3 (f mode (isT m)) vec v mat (tt m) mat r "extractAux"
+
+extractAux f m moder vr modec vc = unsafePerformIO $ do
+    let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr
+        nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc
+    r <- createMatrix RowMajor nr nc
+    app4 (f moder modec (isT m)) vec vr vec vc mat (tt m) mat r "extractAux"
     return r
 
-foreign import ccall unsafe "extractRD" c_extractRD
-    :: CInt -> CInt -> CIdxs (CM Double (CM Double (IO CInt)))
+foreign import ccall unsafe "extractD" c_extractD
+    :: CInt -> CInt -> CInt -> CIdxs (CIdxs (CM Double (CM Double (IO CInt))))
 
-foreign import ccall unsafe "extractRF" c_extractRF
-    :: CInt -> CInt -> CIdxs (CM Float (CM Float (IO CInt)))
+foreign import ccall unsafe "extractF" c_extractF
+    :: CInt -> CInt -> CInt -> CIdxs (CIdxs (CM Float (CM Float (IO CInt))))
 
-foreign import ccall unsafe "extractRC" c_extractRC
-    :: CInt -> CInt -> CIdxs (CM (Complex Double) (CM (Complex Double) (IO CInt)))
+foreign import ccall unsafe "extractC" c_extractC
+    :: CInt -> CInt -> CInt -> CIdxs (CIdxs (CM (Complex Double) (CM (Complex Double) (IO CInt))))
 
-foreign import ccall unsafe "extractRQ" c_extractRQ
-    :: CInt -> CInt -> CIdxs (CM (Complex Float) (CM (Complex Float) (IO CInt)))
+foreign import ccall unsafe "extractQ" c_extractQ
+    :: CInt -> CInt -> CInt -> CIdxs (CIdxs (CM (Complex Float) (CM (Complex Float) (IO CInt))))
 
-foreign import ccall unsafe "extractRI" c_extractRI
-    :: CInt -> CInt -> CIdxs (CM CInt (CM CInt (IO CInt)))
+foreign import ccall unsafe "extractI" c_extractI
+    :: CInt -> CInt -> CInt -> CIdxs (CIdxs (CM CInt (CM CInt (IO CInt))))
 
diff --git a/packages/base/src/Data/Packed/Internal/Numeric.hs b/packages/base/src/Data/Packed/Internal/Numeric.hs
index 00ec70c..353877a 100644
--- a/packages/base/src/Data/Packed/Internal/Numeric.hs
+++ b/packages/base/src/Data/Packed/Internal/Numeric.hs
@@ -39,7 +39,7 @@ module Data.Packed.Internal.Numeric (
     roundVector,
     RealOf, ComplexOf, SingleOf, DoubleOf,
     IndexOf,
-    CInt, Extractor(..), (??),(¿¿),
+    CInt, Extractor(..), (??),
     module Data.Complex
 ) where
 
@@ -69,27 +69,50 @@ type instance ArgOf Matrix a = a -> a -> a
 
 --------------------------------------------------------------------------
 
-data Extractor = All | Range Int Int | At [Int] | AtCyc [Int] | Take Int | Drop Int
+data Extractor
+    = All
+    | Range Int Int
+    | At [Int]
+    | AtCyc [Int]
+    | Take Int
+    | Drop Int
+   deriving Show
 
 idxs js = fromList (map fromIntegral js) :: Idxs
 
-infixl 9 ??, ¿¿
-(??),(¿¿) :: Element t => Matrix t -> Extractor -> Matrix t
-
-m ?? All = m
-m ?? Take 0 = (0><cols m) []
-m ?? Take n | abs n >= rows m = m
-m ?? Drop 0 = m
-m ?? Drop n | abs n >= rows m = (0><cols m) []
-m ?? Range a b | a > b = m ?? Take 0
-m ?? Range a b | a < 0 || b >= cols m = error $
-    printf "can't extract rows %d to %d from matrix %dx%d" a b (rows m) (cols m)
-m ?? At ps | minimum ps < 0 || maximum ps >= rows m = error $
-    printf "can't extract rows %s from matrix %dx%d" (show ps) (rows m) (cols m)
-
-m ?? er = extractR m mode js
+infixl 9 ??
+(??)  :: Element t => Matrix t -> (Extractor,Extractor) -> Matrix t
+
+
+extractError m e = error $ printf "can't extract %s from matrix %dx%d" (show e) (rows m) (cols m)
+
+m ?? e@(Range a b,_) | a < 0 || b >= rows m = extractError m e
+m ?? e@(_,Range a b) | a < 0 || b >= cols m = extractError m e
+m ?? e@(At ps,_) | minimum ps < 0 || maximum ps >= rows m = extractError m e
+m ?? e@(_,At ps) | minimum ps < 0 || maximum ps >= cols m = extractError m e
+
+m ?? (All,All) = m
+
+m ?? (Range a b,e) | a > b = m ?? (Take 0,e)
+m ?? (e,Range a b) | a > b = m ?? (e,Take 0)
+
+m ?? (Take 0,e) = (0><cols m) [] ?? (All,e)
+m ?? (e,Take 0) = (rows m><0) [] ?? (e,All)
+
+m ?? (Take n,e) | abs n > rows m = m ?? (All,e)
+m ?? (e,Take n) | abs n > cols m = m ?? (e,All)
+
+m ?? (Drop 0,e) = m ?? (All,e)
+m ?? (e,Drop 0) = m ?? (e,All)
+
+m ?? (Drop n,e) | abs n > rows m = m ?? (Take 0,e)
+m ?? (e,Drop n) | abs n > cols m = m ?? (e,Take 0)
+
+
+m ?? (er,ec) = extractR m moder rs modec cs
   where
-    (mode,js) = mkExt (rows m) er
+    (moder,rs) = mkExt (rows m) er
+    (modec,cs) = mkExt (cols m) ec
     ran a b = (0, idxs [a,b])
     pos ks  = (1, idxs ks)
     mkExt _ (At ks)       = pos ks
@@ -104,13 +127,6 @@ m ?? er = extractR m mode js
         | otherwise       = mkExt n (Take (n+k))
 
 
-m ¿¿ Range a b | a < 0 || b > cols m -1 = error $
-    printf "can't extract columns %d to %d from matrix %dx%d" a b (rows m) (cols m)
-
-m ¿¿ At ps | minimum ps < 0 || maximum ps >= cols m = error $
-    printf "can't extract columns %s from matrix %dx%d" (show ps) (rows m) (cols m)
-m ¿¿ ec = trans (trans m ?? ec)
-
 -------------------------------------------------------------------
 
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs
index e4485fc..c94350f 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Data.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs
@@ -37,7 +37,7 @@ module Numeric.LinearAlgebra.Data(
     fromList, toList, subVector, takesV, vjoin,
     flatten, reshape, asRow, asColumn, row, col,
     fromRows, toRows, fromColumns, toColumns, fromLists, toLists, fromArray2D,
-    Extractor(..), (??), (¿¿),
+    Extractor(..), (??),
     takeRows, dropRows, takeColumns, dropColumns, subMatrix, (?), (¿), fliprl, flipud,
   
     -- * Block matrix
@@ -81,6 +81,6 @@ import Numeric.LinearAlgebra.Util hiding ((&),(#))
 import Data.Complex
 import Numeric.Sparse
 import Data.Packed.Internal.Vector(Idxs)
-import Data.Packed.Internal.Numeric(CInt,Extractor(..),(??),(¿¿))
+import Data.Packed.Internal.Numeric(CInt,Extractor(..),(??))