use slice interface for lapack funcs (wip)

6 files changed, 163 insertions, 143 deletions

diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index cdbaab9..baa0570 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -38,6 +38,9 @@ typedef float  complex TCF;
 // #define OK return 0;
 // #endif
 
+
+// printf("%dx%d %d:%d\n",ar,ac,aXr,aXc);
+
 #define TRACEMAT(M) {int q; printf(" %d x %d: ",M##r,M##c); \
                      for(q=0;q<M##r*M##c;q++) printf("%.1f ",M##p[q]); printf("\n");}
 
@@ -56,7 +59,7 @@ inline int mod (int a, int b);
 
 inline int64_t mod_l (int64_t a, int64_t b);
 
-//---------------------------------------
+////////////////////////////////////////////////////////////////////////////////
 void asm_finit() {
 #ifdef i386
 
@@ -78,8 +81,6 @@ void asm_finit() {
 #endif
 }
 
-//---------------------------------------
-
 #if NANDEBUG
 
 #define CHECKNANR(M,msg)                     \
@@ -109,16 +110,16 @@ for(k=0; k<(M##r * M##c); k++) {             \
 #define CHECKNANR(M,msg)
 #endif
 
-//---------------------------------------
 
-//////////////////// real svd ////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+//////////////////// real svd ///////////////////////////////////////////////////
 
 /* Subroutine */ int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
         doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
         ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
         integer *info);
 
-int svd_l_R(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) {
+int svd_l_R(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
     integer m = ar;
     integer n = ac;
     integer q = MIN(m,n);
@@ -181,7 +182,7 @@ int svd_l_R(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) {
         doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
         integer *iwork, integer *info);
 
-int svd_l_Rdd(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) {
+int svd_l_Rdd(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
     integer m = ar;
     integer n = ac;
     integer q = MIN(m,n);
@@ -231,7 +232,7 @@ int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
     integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
     integer *lwork, doublereal *rwork, integer *info);
 
-int svd_l_C(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
+int svd_l_C(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
     integer m = ar;
     integer n = ac;
     integer q = MIN(m,n);
@@ -297,7 +298,7 @@ int zgesdd_ (char *jobz, integer *m, integer *n,
     integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
     integer *lwork, doublereal *rwork, integer* iwork, integer *info);
 
-int svd_l_Cdd(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
+int svd_l_Cdd(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
     //printf("entro\n");
     integer m = ar;
     integer n = ac;
@@ -358,7 +359,7 @@ int svd_l_Cdd(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
         integer *ldvl, doublecomplex *vr, integer *ldvr, doublecomplex *work,
         integer *lwork, doublereal *rwork, integer *info);
 
-int eig_l_C(KCMAT(a), CMAT(u), CVEC(s),CMAT(v)) {
+int eig_l_C(KOCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
     integer n = ar;
     REQUIRES(ac==n && sn==n, BAD_SIZE);
     REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
@@ -413,7 +414,7 @@ int eig_l_C(KCMAT(a), CMAT(u), CVEC(s),CMAT(v)) {
         integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work,
         integer *lwork, integer *info);
 
-int eig_l_R(KDMAT(a),DMAT(u), CVEC(s),DMAT(v)) {
+int eig_l_R(KODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) {
     integer n = ar;
     REQUIRES(ac==n && sn==n, BAD_SIZE);
     REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
@@ -461,7 +462,7 @@ int eig_l_R(KDMAT(a),DMAT(u), CVEC(s),DMAT(v)) {
          integer *lda, doublereal *w, doublereal *work, integer *lwork,
         integer *info);
 
-int eig_l_S(int wantV,KDMAT(a),DVEC(s),DMAT(v)) {
+int eig_l_S(int wantV,KODMAT(a),DVEC(s),ODMAT(v)) {
     integer n = ar;
     REQUIRES(ac==n && sn==n, BAD_SIZE);
     REQUIRES(vr==n && vc==n, BAD_SIZE);
@@ -499,7 +500,7 @@ int eig_l_S(int wantV,KDMAT(a),DVEC(s),DMAT(v)) {
         *a, integer *lda, doublereal *w, doublecomplex *work, integer *lwork,
         doublereal *rwork, integer *info);
 
-int eig_l_H(int wantV,KCMAT(a),DVEC(s),CMAT(v)) {
+int eig_l_H(int wantV,KOCMAT(a),DVEC(s),OCMAT(v)) {
     integer n = ar;
     REQUIRES(ac==n && sn==n, BAD_SIZE);
     REQUIRES(vr==n && vc==n, BAD_SIZE);
@@ -541,7 +542,7 @@ int eig_l_H(int wantV,KCMAT(a),DVEC(s),CMAT(v)) {
 /* Subroutine */ int dgesv_(integer *n, integer *nrhs, doublereal *a, integer
         *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info);
 
-int linearSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
+int linearSolveR_l(KODMAT(a),KODMAT(b),ODMAT(x)) {
     integer n = ar;
     integer nhrs = bc;
     REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
@@ -571,7 +572,7 @@ int linearSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
         integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer *
         info);
 
-int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
+int linearSolveC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) {
     integer n = ar;
     integer nhrs = bc;
     REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
@@ -601,7 +602,7 @@ int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
         doublereal *a, integer *lda, doublereal *b, integer *ldb, integer *
         info);
 
-int cholSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
+int cholSolveR_l(KODMAT(a),KODMAT(b),ODMAT(x)) {
     integer n = ar;
     integer nhrs = bc;
     REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
@@ -623,7 +624,7 @@ int cholSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
         doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
         integer *info);
 
-int cholSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
+int cholSolveC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) {
     integer n = ar;
     integer nhrs = bc;
     REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
@@ -645,7 +646,7 @@ int cholSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
         nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb,
         doublereal *work, integer *lwork, integer *info);
 
-int linearSolveLSR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
+int linearSolveLSR_l(KODMAT(a),KODMAT(b),ODMAT(x)) {
     integer m = ar;
     integer n = ac;
     integer nrhs = bc;
@@ -693,7 +694,7 @@ int linearSolveLSR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
         nrhs, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
         doublecomplex *work, integer *lwork, integer *info);
 
-int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
+int linearSolveLSC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) {
     integer m = ar;
     integer n = ac;
     integer nrhs = bc;
@@ -742,7 +743,7 @@ int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
         s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork,
          integer *info);
 
-int linearSolveSVDR_l(double rcond,KDMAT(a),KDMAT(b),DMAT(x)) {
+int linearSolveSVDR_l(double rcond,KODMAT(a),KODMAT(b),ODMAT(x)) {
     integer m = ar;
     integer n = ac;
     integer nrhs = bc;
@@ -801,7 +802,7 @@ int zgelss_(integer *m, integer *n, integer *nhrs,
     doublecomplex *work, integer* lwork, doublereal* rwork,
     integer *info);
 
-int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) {
+int linearSolveSVDC_l(double rcond, KOCMAT(a),KOCMAT(b),OCMAT(x)) {
     integer m = ar;
     integer n = ac;
     integer nrhs = bc;
@@ -859,7 +860,7 @@ int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) {
 /* Subroutine */ int zpotrf_(char *uplo, integer *n, doublecomplex *a,
         integer *lda, integer *info);
 
-int chol_l_H(KCMAT(a),CMAT(l)) {
+int chol_l_H(KOCMAT(a),OCMAT(l)) {
     integer n = ar;
     REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE);
     DEBUGMSG("chol_l_H");
@@ -871,9 +872,9 @@ int chol_l_H(KCMAT(a),CMAT(l)) {
     CHECK(res,res);
     doublecomplex zero = {0.,0.};
     int r,c;
-    for (r=0; r<lr-1; r++) {
-        for(c=r+1; c<lc; c++) {
-            lp[r*lc+c] = zero;
+    for (r=0; r<lr; r++) {
+        for(c=0; c<r; c++) {
+            AT(l,r,c) = zero;
         }
     }
     OK
@@ -883,7 +884,7 @@ int chol_l_H(KCMAT(a),CMAT(l)) {
 /* Subroutine */ int dpotrf_(char *uplo, integer *n, doublereal *a, integer *
         lda, integer *info);
 
-int chol_l_S(KDMAT(a),DMAT(l)) {
+int chol_l_S(KODMAT(a),ODMAT(l)) {
     integer n = ar;
     REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE);
     DEBUGMSG("chol_l_S");
@@ -894,9 +895,9 @@ int chol_l_S(KDMAT(a),DMAT(l)) {
     CHECK(res>0,NODEFPOS);
     CHECK(res,res);
     int r,c;
-    for (r=0; r<lr-1; r++) {
-        for(c=r+1; c<lc; c++) {
-            lp[r*lc+c] = 0.;
+    for (r=0; r<lr; r++) {
+        for(c=0; c<r; c++) {
+            AT(l,r,c) = 0.;
         }
     }
     OK
@@ -907,7 +908,7 @@ int chol_l_S(KDMAT(a),DMAT(l)) {
 /* Subroutine */ int dgeqr2_(integer *m, integer *n, doublereal *a, integer *
         lda, doublereal *tau, doublereal *work, integer *info);
 
-int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r)) {
+int qr_l_R(KODMAT(a), DVEC(tau), ODMAT(r)) {
     integer m = ar;
     integer n = ac;
     integer mn = MIN(m,n);
@@ -926,7 +927,7 @@ int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r)) {
 /* Subroutine */ int zgeqr2_(integer *m, integer *n, doublecomplex *a,
         integer *lda, doublecomplex *tau, doublecomplex *work, integer *info);
 
-int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
+int qr_l_C(KOCMAT(a), CVEC(tau), OCMAT(r)) {
     integer m = ar;
     integer n = ac;
     integer mn = MIN(m,n);
@@ -946,7 +947,7 @@ int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
         a, integer *lda, doublereal *tau, doublereal *work, integer *lwork,
         integer *info);
 
-int c_dorgqr(KDMAT(a), KDVEC(tau), DMAT(r)) {
+int c_dorgqr(KODMAT(a), KDVEC(tau), ODMAT(r)) {
     integer m = ar;
     integer n = MIN(ac,ar);
     integer k = taun;
@@ -966,7 +967,7 @@ int c_dorgqr(KDMAT(a), KDVEC(tau), DMAT(r)) {
         doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *
         work, integer *lwork, integer *info);
 
-int c_zungqr(KCMAT(a), KCVEC(tau), CMAT(r)) {
+int c_zungqr(KOCMAT(a), KCVEC(tau), OCMAT(r)) {
     integer m = ar;
     integer n = MIN(ac,ar);
     integer k = taun;
@@ -989,7 +990,7 @@ int c_zungqr(KCMAT(a), KCVEC(tau), CMAT(r)) {
         doublereal *a, integer *lda, doublereal *tau, doublereal *work,
         integer *lwork, integer *info);
 
-int hess_l_R(KDMAT(a), DVEC(tau), DMAT(r)) {
+int hess_l_R(KODMAT(a), DVEC(tau), ODMAT(r)) {
     integer m = ar;
     integer n = ac;
     integer mn = MIN(m,n);
@@ -1012,7 +1013,7 @@ int hess_l_R(KDMAT(a), DVEC(tau), DMAT(r)) {
         doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *
         work, integer *lwork, integer *info);
 
-int hess_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
+int hess_l_C(KOCMAT(a), CVEC(tau), OCMAT(r)) {
     integer m = ar;
     integer n = ac;
     integer mn = MIN(m,n);
@@ -1037,7 +1038,7 @@ int hess_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
         doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work,
         integer *lwork, logical *bwork, integer *info);
 
-int schur_l_R(KDMAT(a), DMAT(u), DMAT(s)) {
+int schur_l_R(KODMAT(a), ODMAT(u), ODMAT(s)) {
     integer m = ar;
     integer n = ac;
     REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE);
@@ -1077,7 +1078,7 @@ int schur_l_R(KDMAT(a), DMAT(u), DMAT(s)) {
         doublecomplex *vs, integer *ldvs, doublecomplex *work, integer *lwork,
          doublereal *rwork, logical *bwork, integer *info);
 
-int schur_l_C(KCMAT(a), CMAT(u), CMAT(s)) {
+int schur_l_C(KOCMAT(a), OCMAT(u), OCMAT(s)) {
     integer m = ar;
     integer n = ac;
     REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE);
@@ -1109,7 +1110,7 @@ int schur_l_C(KCMAT(a), CMAT(u), CMAT(s)) {
 /* Subroutine */ int dgetrf_(integer *m, integer *n, doublereal *a, integer *
         lda, integer *ipiv, integer *info);
 
-int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r)) {
+int lu_l_R(KODMAT(a), DVEC(ipiv), ODMAT(r)) {
     integer m = ar;
     integer n = ac;
     integer mn = MIN(m,n);
@@ -1135,7 +1136,7 @@ int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r)) {
 /* Subroutine */ int zgetrf_(integer *m, integer *n, doublecomplex *a,
         integer *lda, integer *ipiv, integer *info);
 
-int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r)) {
+int lu_l_C(KOCMAT(a), DVEC(ipiv), OCMAT(r)) {
     integer m = ar;
     integer n = ac;
     integer mn = MIN(m,n);
@@ -1164,7 +1165,7 @@ int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r)) {
         doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *
         ldb, integer *info);
 
-int luS_l_R(KDMAT(a), KDVEC(ipiv), KDMAT(b), DMAT(x)) {
+int luS_l_R(KODMAT(a), KDVEC(ipiv), KODMAT(b), ODMAT(x)) {
   integer m = ar;
   integer n = ac;
   integer mrhs = br;
@@ -1189,7 +1190,7 @@ int luS_l_R(KDMAT(a), KDVEC(ipiv), KDMAT(b), DMAT(x)) {
         doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *b,
         integer *ldb, integer *info);
 
-int luS_l_C(KCMAT(a), KDVEC(ipiv), KCMAT(b), CMAT(x)) {
+int luS_l_C(KOCMAT(a), KDVEC(ipiv), KOCMAT(b), OCMAT(x)) {
     integer m = ar;
     integer n = ac;
     integer mrhs = br;
@@ -1215,7 +1216,7 @@ void dgemm_(char *, char *, integer *, integer *, integer *,
            double *, const double *, integer *, const double *,
            integer *, double *, double *, integer *);
 
-int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)) {
+int multiplyR(int ta, int tb, KODMAT(a),KODMAT(b),ODMAT(r)) {
     //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
     DEBUGMSG("dgemm_");
     CHECKNANR(a,"NaN multR Input\n")
@@ -1237,7 +1238,7 @@ void zgemm_(char *, char *, integer *, integer *, integer *,
            doublecomplex *, const doublecomplex *, integer *, const doublecomplex *,
            integer *, doublecomplex *, doublecomplex *, integer *);
 
-int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) {
+int multiplyC(int ta, int tb, KOCMAT(a),KOCMAT(b),OCMAT(r)) {
     //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
     DEBUGMSG("zgemm_");
     CHECKNANC(a,"NaN multC Input\n")
@@ -1262,7 +1263,7 @@ void sgemm_(char *, char *, integer *, integer *, integer *,
             float *, const float *, integer *, const float *,
            integer *, float *, float *, integer *);
 
-int multiplyF(int ta, int tb, KFMAT(a),KFMAT(b),FMAT(r)) {
+int multiplyF(int ta, int tb, KOFMAT(a),KOFMAT(b),OFMAT(r)) {
     //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
     DEBUGMSG("sgemm_");
     integer m = ta?ac:ar;
@@ -1281,7 +1282,7 @@ void cgemm_(char *, char *, integer *, integer *, integer *,
            complex *, const complex *, integer *, const complex *,
            integer *, complex *, complex *, integer *);
 
-int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) {
+int multiplyQ(int ta, int tb, KOQMAT(a),KOQMAT(b),OQMAT(r)) {
     //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
     DEBUGMSG("cgemm_");
     integer m = ta?ac:ar;
@@ -1564,13 +1565,13 @@ int remapQ(KOIMAT(i), KOIMAT(j), KOQMAT(m), OQMAT(r)) {
 
 ////////////////////////////////////////////////////////////////////////////////
 
-int saveMatrix(char * file, char * format, KDMAT(a)){
+int saveMatrix(char * file, char * format, KODMAT(a)){
     FILE * fp;
     fp = fopen (file, "w");
     int r, c;
     for (r=0;r<ar; r++) {
         for (c=0; c<ac; c++) {
-            fprintf(fp,format,ap[r*ac+c]);
+            fprintf(fp,format,AT(a,r,c));
             if (c<ac-1) {
                 fprintf(fp," ");
             } else {
diff --git a/packages/base/src/Internal/C/lapack-aux.h b/packages/base/src/Internal/C/lapack-aux.h
index bf8c5e9..b38ca7a 100644
--- a/packages/base/src/Internal/C/lapack-aux.h
+++ b/packages/base/src/Internal/C/lapack-aux.h
@@ -52,16 +52,6 @@ typedef short ftnlen;
 #define CMAT(A) int A##r, int A##c, doublecomplex* A##p
 #define PMAT(A) int A##r, int A##c, void* A##p, int A##s
 
-#define OIMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, int* A##p
-#define OLMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, int64_t* A##p
-#define OFMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, float* A##p
-#define ODMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, double* A##p
-#define OQMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, complex* A##p
-#define OCMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, doublecomplex* A##p
-
-#define VECG(T,A) int A##n, T* A##p
-#define MATG(T,A) int A##r, int A##c, int A##Xr, int A##Xc, T* A##p
-
 #define KIVEC(A) int A##n, const int*A##p
 #define KLVEC(A) int A##n, const int64_t*A##p
 #define KFVEC(A) int A##n, const float*A##p
@@ -78,12 +68,22 @@ typedef short ftnlen;
 #define KCMAT(A) int A##r, int A##c, const doublecomplex* A##p
 #define KPMAT(A) int A##r, int A##c, const void* A##p, int A##s
 
-#define KOIMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, const int* A##p
-#define KOLMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, const int64_t* A##p
-#define KOFMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, const float* A##p
-#define KODMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, const double* A##p
-#define KOQMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, const complex* A##p
-#define KOCMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, const doublecomplex* A##p
+#define VECG(T,A) int A##n, T* A##p
+#define MATG(T,A) int A##r, int A##c, int A##Xr, int A##Xc, T* A##p
+
+#define OIMAT(A) MATG(int,A)
+#define OLMAT(A) MATG(int64_t,A)
+#define OFMAT(A) MATG(float,A)
+#define ODMAT(A) MATG(double,A)
+#define OQMAT(A) MATG(complex,A)
+#define OCMAT(A) MATG(doublecomplex,A)
+
+#define KOIMAT(A) MATG(const int,A)
+#define KOLMAT(A) MATG(const int64_t,A)
+#define KOFMAT(A) MATG(const float,A)
+#define KODMAT(A) MATG(const double,A)
+#define KOQMAT(A) MATG(const complex,A)
+#define KOCMAT(A) MATG(const doublecomplex,A)
 
 #define AT(m,i,j) (m##p[(i)*m##Xr + (j)*m##Xc])
 #define TRAV(m,i,j) int i,j; for (i=0;i<m##r;i++) for (j=0;j<m##c;j++)
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index 5319e95..2c7148b 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -29,16 +29,12 @@ import System.IO.Unsafe(unsafePerformIO)
 -----------------------------------------------------------------------------------
 
 infixl 1 #
-a # b = applyRaw a b
+a # b = apply a b
 {-# INLINE (#) #-}
 
-infixl 1 #!
-a #! b = apply a b
-{-# INLINE (#!) #-}
-
 -----------------------------------------------------------------------------------
 
-type TMMM t = t ..> t ..> t ..> Ok
+type TMMM t = t ::> t ::> t ::> Ok
 
 type F = Float
 type Q = Complex Float
@@ -47,8 +43,8 @@ foreign import ccall unsafe "multiplyR" dgemmc :: CInt -> CInt -> TMMM R
 foreign import ccall unsafe "multiplyC" zgemmc :: CInt -> CInt -> TMMM C
 foreign import ccall unsafe "multiplyF" sgemmc :: CInt -> CInt -> TMMM F
 foreign import ccall unsafe "multiplyQ" cgemmc :: CInt -> CInt -> TMMM Q
-foreign import ccall unsafe "multiplyI" c_multiplyI :: I -> I ::> I ::> I ::> Ok
-foreign import ccall unsafe "multiplyL" c_multiplyL :: Z -> Z ::> Z ::> Z ::> Ok
+foreign import ccall unsafe "multiplyI" c_multiplyI :: I -> TMMM I
+foreign import ccall unsafe "multiplyL" c_multiplyL :: Z -> TMMM Z
 
 isT (rowOrder -> False) = 0
 isT _                   = 1
@@ -84,7 +80,7 @@ multiplyI m a b = unsafePerformIO $ do
     when (cols a /= rows b) $ error $
         "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b
     s <- createMatrix ColumnMajor (rows a) (cols b)
-    c_multiplyI m #! a #! b #! s #|"c_multiplyI"
+    c_multiplyI m # a # b # s #|"c_multiplyI"
     return s
 
 multiplyL :: Z -> Matrix Z -> Matrix Z -> Matrix Z
@@ -92,12 +88,12 @@ multiplyL m a b = unsafePerformIO $ do
     when (cols a /= rows b) $ error $
         "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b
     s <- createMatrix ColumnMajor (rows a) (cols b)
-    c_multiplyL m #! a #! b #! s #|"c_multiplyL"
+    c_multiplyL m # a # b # s #|"c_multiplyL"
     return s
 
 -----------------------------------------------------------------------------
 
-type TSVD t = t ..> t ..> R :> t ..> Ok
+type TSVD t = t ::> t ::> R :> t ::> Ok
 
 foreign import ccall unsafe "svd_l_R" dgesvd :: TSVD R
 foreign import ccall unsafe "svd_l_C" zgesvd :: TSVD C
@@ -126,8 +122,9 @@ svdAux f st x = unsafePerformIO $ do
     v <- createMatrix ColumnMajor c c
     f # x # u # s # v #| st
     return (u,s,v)
-  where r = rows x
-        c = cols x
+  where
+    r = rows x
+    c = cols x
 
 
 -- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'.
@@ -152,9 +149,10 @@ thinSVDAux f st x = unsafePerformIO $ do
     v <- createMatrix ColumnMajor q c
     f # x # u # s # v #| st
     return (u,s,v)
-  where r = rows x
-        c = cols x
-        q = min r c
+  where
+    r = rows x
+    c = cols x
+    q = min r c
 
 
 -- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'.
@@ -177,10 +175,11 @@ svAux f st x = unsafePerformIO $ do
     s <- createVector q
     g # x # s #| st
     return s
-  where r = rows x
-        c = cols x
-        q = min r c
-        g ra ca pa nb pb = f ra ca pa 0 0 nullPtr nb pb 0 0 nullPtr
+  where
+    r = rows x
+    c = cols x
+    q = min r c
+    g ra ca xra xca pa nb pb = f ra ca xra xca pa 0 0 0 0 nullPtr nb pb 0 0 0 0 nullPtr
 
 
 -- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'.
@@ -196,10 +195,11 @@ rightSVAux f st x = unsafePerformIO $ do
     v <- createMatrix ColumnMajor c c
     g # x # s # v #| st
     return (s,v)
-  where r = rows x
-        c = cols x
-        q = min r c
-        g ra ca pa = f ra ca pa 0 0 nullPtr
+  where
+    r = rows x
+    c = cols x
+    q = min r c
+    g ra ca xra xca pa = f ra ca xra xca pa 0 0 0 0 nullPtr
 
 
 -- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/  with jobu == \'A\' and jobvt == \'N\'.
@@ -215,25 +215,27 @@ leftSVAux f st x = unsafePerformIO $ do
     s <- createVector q
     g # x # u # s #| st
     return (u,s)
-  where r = rows x
-        c = cols x
-        q = min r c
-        g ra ca pa ru cu pu nb pb = f ra ca pa ru cu pu nb pb 0 0 nullPtr
+  where
+    r = rows x
+    c = cols x
+    q = min r c
+    g ra ca xra xca pa ru cu xru xcu pu nb pb = f ra ca xra xca pa ru cu xru xcu pu nb pb 0 0 0 0 nullPtr
 
 -----------------------------------------------------------------------------
 
-foreign import ccall unsafe "eig_l_R" dgeev :: R ..> R ..> C :> R ..> Ok
-foreign import ccall unsafe "eig_l_C" zgeev :: C ..> C ..> C :> C ..> Ok
-foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> R ..> R :> R ..> Ok
-foreign import ccall unsafe "eig_l_H" zheev :: CInt -> C ..> R :> C ..> Ok
+foreign import ccall unsafe "eig_l_R" dgeev :: R ::> R ::> C :> R ::> Ok
+foreign import ccall unsafe "eig_l_C" zgeev :: C ::> C ::> C :> C ::> Ok
+foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> R ::> R :> R ::> Ok
+foreign import ccall unsafe "eig_l_H" zheev :: CInt -> C ::> R :> C ::> Ok
 
 eigAux f st m = unsafePerformIO $ do
-        l <- createVector r
-        v <- createMatrix ColumnMajor r r
-        g # m # l # v #| st
-        return (l,v)
-  where r = rows m
-        g ra ca pa = f ra ca pa 0 0 nullPtr
+    l <- createVector r
+    v <- createMatrix ColumnMajor r r
+    g # m # l # v #| st
+    return (l,v)
+  where
+    r = rows m
+    g ra ca xra xca pa = f ra ca xra xca pa 0 0 0 0 nullPtr
 
 
 -- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/.
@@ -242,11 +244,12 @@ eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Dou
 eigC = eigAux zgeev "eigC" . fmat
 
 eigOnlyAux f st m = unsafePerformIO $ do
-        l <- createVector r
-        g # m # l #| st
-        return l
-  where r = rows m
-        g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr
+    l <- createVector r
+    g # m # l #| st
+    return l
+  where
+    r = rows m
+    g ra ca xra xca pa nl pl = f ra ca xra xca pa 0 0 0 0 nullPtr nl pl 0 0 0 0 nullPtr
 
 -- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'.
 -- The eigenvalues are not sorted.
@@ -264,12 +267,13 @@ eigR m = (s', v'')
 
 eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
 eigRaux m = unsafePerformIO $ do
-        l <- createVector r
-        v <- createMatrix ColumnMajor r r
-        g # m # l # v #| "eigR"
-        return (l,v)
-  where r = rows m
-        g ra ca pa = dgeev ra ca pa 0 0 nullPtr
+    l <- createVector r
+    v <- createMatrix ColumnMajor r r
+    g # m # l # v #| "eigR"
+    return (l,v)
+  where
+    r = rows m
+    g ra ca xra xca pa = dgeev ra ca xra xca pa 0 0 0 0 nullPtr
 
 fixeig1 s = toComplex' (subVector 0 r (asReal s), subVector r r (asReal s))
     where r = dim s
@@ -291,11 +295,12 @@ eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat
 -----------------------------------------------------------------------------
 
 eigSHAux f st m = unsafePerformIO $ do
-        l <- createVector r
-        v <- createMatrix ColumnMajor r r
-        f # m # l # v #| st
-        return (l,v)
-  where r = rows m
+    l <- createVector r
+    v <- createMatrix ColumnMajor r r
+    f # m # l # v #| st
+    return (l,v)
+  where
+    r = rows m
 
 -- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/.
 -- The eigenvectors are the columns of v.
@@ -314,8 +319,9 @@ eigS' = eigSHAux (dsyev 1) "eigS'" . fmat
 -- The eigenvalues are sorted in descending order (use 'eigH'' for ascending order).
 eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
 eigH m = (s', fliprl v)
-    where (s,v) = eigH' (fmat m)
-          s' = fromList . reverse . toList $  s
+  where
+    (s,v) = eigH' (fmat m)
+    s' = fromList . reverse . toList $  s
 
 -- | 'eigH' in ascending order
 eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
@@ -346,10 +352,11 @@ linearSolveSQAux g f st a b
         f # a # b # s #| st
         return s
     | otherwise = error $ st ++ " of nonsquare matrix"
-  where n1 = rows a
-        n2 = cols a
-        r  = rows b
-        c  = cols b
+  where
+    n1 = rows a
+    n2 = cols a
+    r  = rows b
+    c  = cols b
 
 -- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'.
 linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
@@ -375,6 +382,7 @@ cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Comp
 cholSolveC a b = linearSolveSQAux id zpotrs "cholSolveC" (fmat a) (fmat b)
 
 -----------------------------------------------------------------------------------
+
 foreign import ccall unsafe "linearSolveLSR_l" dgels :: TMMM R
 foreign import ccall unsafe "linearSolveLSC_l" zgels :: TMMM C
 foreign import ccall unsafe "linearSolveSVDR_l" dgelss :: Double -> TMMM R
@@ -384,9 +392,10 @@ linearSolveAux f st a b = unsafePerformIO $ do
     r <- createMatrix ColumnMajor (max m n) nrhs
     f # a # b # r #| st
     return r
-  where m = rows a
-        n = cols a
-        nrhs = cols b
+  where
+    m = rows a
+    n = cols a
+    nrhs = cols b
 
 -- | Least squared error solution of an overconstrained real linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /dgels/. For rank-deficient systems use 'linearSolveSVDR'.
 linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
@@ -418,7 +427,7 @@ linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b)
 
 -----------------------------------------------------------------------------------
 
-type TMM t = t ..> t ..> Ok
+type TMM t = t ::> t ::> Ok
 
 foreign import ccall unsafe "chol_l_H" zpotrf :: TMM C
 foreign import ccall unsafe "chol_l_S" dpotrf :: TMM R
@@ -427,7 +436,8 @@ cholAux f st a = do
     r <- createMatrix ColumnMajor n n
     f # a # r #| st
     return r
-  where n = rows a
+  where
+    n = rows a
 
 -- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.
 cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
@@ -447,7 +457,7 @@ mbCholS =  unsafePerformIO . mbCatch . cholAux dpotrf "cholS" . fmat
 
 -----------------------------------------------------------------------------------
 
-type TMVM t = t ..> t :> t ..> Ok
+type TMVM t = t ::> t :> t ::> Ok
 
 foreign import ccall unsafe "qr_l_R" dgeqr2 :: TMVM R
 foreign import ccall unsafe "qr_l_C" zgeqr2 :: TMVM C
@@ -504,13 +514,14 @@ hessAux f st a = unsafePerformIO $ do
     tau <- createVector (mn-1)
     f # a # tau # r #| st
     return (r,tau)
-  where m = rows a
-        n = cols a
-        mn = min m n
+  where
+    m = rows a
+    n = cols a
+    mn = min m n
 
 -----------------------------------------------------------------------------------
-foreign import ccall unsafe "schur_l_R" dgees :: TMMM R
-foreign import ccall unsafe "schur_l_C" zgees :: TMMM C
+foreign import ccall unsafe "schur_l_R" dgees :: R ::> R ::> R ::> Ok
+foreign import ccall unsafe "schur_l_C" zgees :: C ::> C ::> C ::> Ok
 
 -- | Schur factorization of a square real matrix, using LAPACK's /dgees/.
 schurR :: Matrix Double -> (Matrix Double, Matrix Double)
@@ -525,11 +536,12 @@ schurAux f st a = unsafePerformIO $ do
     s <- createMatrix ColumnMajor n n
     f # a # u # s #| st
     return (u,s)
-  where n = rows a
+  where
+    n = rows a
 
 -----------------------------------------------------------------------------------
 foreign import ccall unsafe "lu_l_R" dgetrf :: TMVM R
-foreign import ccall unsafe "lu_l_C" zgetrf :: C ..> R :> C ..> Ok
+foreign import ccall unsafe "lu_l_C" zgetrf :: C ::> R :> C ::> Ok
 
 -- | LU factorization of a general real matrix, using LAPACK's /dgetrf/.
 luR :: Matrix Double -> (Matrix Double, [Int])
@@ -544,12 +556,13 @@ luAux f st a = unsafePerformIO $ do
     piv <- createVector (min n m)
     f # a # piv # lu #| st
     return (lu, map (pred.round) (toList piv))
-  where n = rows a
-        m = cols a
+  where
+    n = rows a
+    m = cols a
 
 -----------------------------------------------------------------------------------
 
-type Tlus t = t ..> Double :> t ..> t ..> Ok
+type Tlus t = t ::> Double :> t ::> t ::> Ok
 
 foreign import ccall unsafe "luS_l_R" dgetrs :: Tlus R
 foreign import ccall unsafe "luS_l_C" zgetrs :: Tlus C
diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs
index c4f95d8..54d9cb8 100644
--- a/packages/base/src/Internal/Modular.hs
+++ b/packages/base/src/Internal/Modular.hs
@@ -398,6 +398,7 @@ test = (ok, info)
         print lgm
         print $ lgm <> lgm
         
+        putStrLn "checkGen"
         print (checkGen (gen 5 :: Matrix R))
         print (checkGen (gen 5 :: Matrix Float))
         print (checkGen (cgen 5 :: Matrix C))
@@ -408,6 +409,7 @@ test = (ok, info)
         print $ mutable (luST (const True)) (gen 5 :: Matrix R)
         print $ mutable (luST (const True)) (gen 5 :: Matrix (Mod 11 Z))
 
+        putStrLn "checkLU"
         print $ checkLU (magnit 0) (gen 5 :: Matrix R)
         print $ checkLU (magnit 0) (gen 5 :: Matrix Float)
         print $ checkLU (magnit 0) (cgen 5 :: Matrix C)
@@ -415,6 +417,7 @@ test = (ok, info)
         print $ checkLU (magnit 0) (gen 5 :: Matrix (Mod 7 I))
         print $ checkLU (magnit 0) (gen 5 :: Matrix (Mod 7 Z))
 
+        putStrLn "checkSolve"
         print $ checkSolve (gen 5 :: Matrix R)
         print $ checkSolve (gen 5 :: Matrix Float)
         print $ checkSolve (cgen 5 :: Matrix C)
@@ -422,6 +425,7 @@ test = (ok, info)
         print $ checkSolve (gen 5 :: Matrix (Mod 7 I))
         print $ checkSolve (gen 5 :: Matrix (Mod 7 Z))
         
+        putStrLn "luSolve'"
         print $ luSolve' (luPacked' tmm) (ident (rows tmm))
         print $ invershur tmm
 
diff --git a/packages/base/src/Internal/ST.hs b/packages/base/src/Internal/ST.hs
index 73cdf0c..23fda99 100644
--- a/packages/base/src/Internal/ST.hs
+++ b/packages/base/src/Internal/ST.hs
@@ -150,6 +150,7 @@ unsafeFreezeMatrix (STMatrix x) = unsafeIOToST . return $ x
 freezeMatrix :: (Storable t) => STMatrix s t -> ST s (Matrix t)
 freezeMatrix m = liftSTMatrix id m
 
+-- FIXME
 cloneMatrix m = cloneVector (xdat m) >>= return . (\d' -> m{xdat = d'})
 
 {-# INLINE safeIndexM #-}
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index 98eb4ef..258c3a3 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -845,10 +845,11 @@ viewBlock' r c m
     | otherwise        = Block m11 m12 m21 m22
   where
     (rt,ct) = size m
-    m11 = sliceMatrix (0,0) (r,c)       m
-    m12 = sliceMatrix (0,c) (r,ct-c)    m
-    m21 = sliceMatrix (r,0) (rt-r,c)    m
-    m22 = sliceMatrix (r,c) (rt-r,ct-c) m
+    m11 = subm (0,0) (r,c)       m
+    m12 = subm (0,c) (r,ct-c)    m
+    m21 = subm (r,0) (rt-r,c)    m
+    m22 = subm (r,c) (rt-r,ct-c) m
+    subm = sliceMatrix
 
 viewBlock m = viewBlock' n n m
   where