Data.Packed -> base (I)

14 files changed, 1606 insertions, 176 deletions

diff --git a/packages/base/hmatrix-base.cabal b/packages/base/hmatrix-base.cabal
index 148c8f3..96f90e2 100644
--- a/packages/base/hmatrix-base.cabal
+++ b/packages/base/hmatrix-base.cabal
@@ -15,22 +15,42 @@ cabal-version:      >=1.8
 
 build-type:         Simple
 
-extra-source-files:
+extra-source-files: src/C/lapack-aux.h
 
 library
 
-    Build-Depends:      base >= 4 && < 5
+    Build-Depends:      base,
+                        binary,
+                        array,
+                        deepseq,
+                        storable-complex,
+                        vector >= 0.8
 
     hs-source-dirs:     src
 
-    exposed-modules:    
-
-    other-modules:      
+    exposed-modules:    Data.Packed,
+                        Data.Packed.Vector,
+                        Data.Packed.Matrix,
+                        Data.Packed.Foreign,
+                        Data.Packed.ST,
+                        Data.Packed.Development
+                        
+    other-modules:      Data.Packed.Internal,
+                        Data.Packed.Internal.Common,
+                        Data.Packed.Internal.Signatures,
+                        Data.Packed.Internal.Vector,
+                        Data.Packed.Internal.Matrix
+
+    C-sources:          src/C/lapack-aux.c
+                        
 
     extensions:         ForeignFunctionInterface,
                         CPP
 
     ghc-options:        -Wall
+                        -fno-warn-missing-signatures
+                        -fno-warn-orphans
+--                        -fno-warn-unused-binds
 
     cc-options:         -O4 -msse2 -Wall
 
diff --git a/packages/base/src/C/lapack-aux.c b/packages/base/src/C/lapack-aux.c
new file mode 100644
index 0000000..e5e45ef
--- /dev/null
+++ b/packages/base/src/C/lapack-aux.c
@@ -0,0 +1,1489 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <time.h>
+#include "lapack-aux.h"
+
+#define MACRO(B) do {B} while (0)
+#define ERROR(CODE) MACRO(return CODE;)
+#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
+
+#define MIN(A,B) ((A)<(B)?(A):(B))
+#define MAX(A,B) ((A)>(B)?(A):(B))
+
+// #define DBGL
+
+#ifdef DBGL
+#define DEBUGMSG(M) printf("\nLAPACK "M"\n");
+#else
+#define DEBUGMSG(M)
+#endif
+
+#define OK return 0;
+
+// #ifdef DBGL
+// #define DEBUGMSG(M) printf("LAPACK Wrapper "M"\n: "); size_t t0 = time(NULL);
+// #define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;);
+// #else
+// #define DEBUGMSG(M)
+// #define OK return 0;
+// #endif
+
+#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");}
+
+#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
+
+#define BAD_SIZE 2000
+#define BAD_CODE 2001
+#define MEM      2002
+#define BAD_FILE 2003
+#define SINGULAR 2004
+#define NOCONVER 2005
+#define NODEFPOS 2006
+#define NOSPRTD  2007
+
+//---------------------------------------
+void asm_finit() {
+#ifdef i386
+
+//  asm("finit");
+
+    static unsigned char buf[108];
+    asm("FSAVE %0":"=m" (buf));
+
+    #if FPUDEBUG
+    if(buf[8]!=255 || buf[9]!=255) {  // print warning in red
+        printf("%c[;31mWarning: FPU TAG = %x %x\%c[0m\n",0x1B,buf[8],buf[9],0x1B);
+    }
+    #endif
+
+    #if NANDEBUG
+    asm("FRSTOR %0":"=m" (buf));
+    #endif
+
+#endif
+}
+
+//---------------------------------------
+
+#if NANDEBUG
+
+#define CHECKNANR(M,msg)                     \
+{ int k;                                     \
+for(k=0; k<(M##r * M##c); k++) {             \
+    if(M##p[k] != M##p[k]) {                 \
+        printf(msg);                         \
+        TRACEMAT(M)                          \
+        /*exit(1);*/                         \
+    }                                        \
+}                                            \
+}
+
+#define CHECKNANC(M,msg)                     \
+{ int k;                                     \
+for(k=0; k<(M##r * M##c); k++) {             \
+    if(  M##p[k].r != M##p[k].r              \
+      || M##p[k].i != M##p[k].i) {           \
+        printf(msg);                         \
+        /*exit(1);*/                         \
+    }                                        \
+}                                            \
+}
+
+#else
+#define CHECKNANC(M,msg)
+#define CHECKNANR(M,msg)
+#endif
+
+//---------------------------------------
+
+//////////////////// 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)) {
+    integer m = ar;
+    integer n = ac;
+    integer q = MIN(m,n);
+    REQUIRES(sn==q,BAD_SIZE);
+    REQUIRES(up==NULL || (ur==m && (uc==m || uc==q)),BAD_SIZE);
+    char* jobu  = "A";
+    if (up==NULL) {
+        jobu = "N";
+    } else {
+        if (uc==q) {
+            jobu = "S";
+        }
+    }
+    REQUIRES(vp==NULL || (vc==n && (vr==n || vr==q)),BAD_SIZE);
+    char* jobvt  = "A";
+    integer ldvt = n;
+    if (vp==NULL) {
+        jobvt = "N";
+    } else {
+        if (vr==q) {
+            jobvt = "S";
+            ldvt = q;
+        }
+    }
+    DEBUGMSG("svd_l_R");
+    double *B = (double*)malloc(m*n*sizeof(double));
+    CHECK(!B,MEM);
+    memcpy(B,ap,m*n*sizeof(double));
+    integer lwork = -1;
+    integer res;
+    // ask for optimal lwork
+    double ans;
+    dgesvd_ (jobu,jobvt,
+             &m,&n,B,&m,
+             sp,
+             up,&m,
+             vp,&ldvt,
+             &ans, &lwork,
+             &res);
+    lwork = ceil(ans);
+    double * work = (double*)malloc(lwork*sizeof(double));
+    CHECK(!work,MEM);
+    dgesvd_ (jobu,jobvt,
+             &m,&n,B,&m,
+             sp,
+             up,&m,
+             vp,&ldvt,
+             work, &lwork,
+             &res);
+    CHECK(res,res);
+    free(work);
+    free(B);
+    OK
+}
+
+// (alternative version)
+
+/* Subroutine */ int dgesdd_(char *jobz, integer *m, integer *n, doublereal *
+        a, integer *lda, doublereal *s, doublereal *u, integer *ldu,
+        doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
+        integer *iwork, integer *info);
+
+int svd_l_Rdd(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) {
+    integer m = ar;
+    integer n = ac;
+    integer q = MIN(m,n);
+    REQUIRES(sn==q,BAD_SIZE);
+    REQUIRES((up == NULL && vp == NULL)
+             || (ur==m && vc==n
+                &&   ((uc == q && vr == q)
+                   || (uc == m && vc==n))),BAD_SIZE);
+    char* jobz  = "A";
+    integer ldvt = n;
+    if (up==NULL) {
+        jobz = "N";
+    } else {
+        if (uc==q && vr == q) {
+            jobz = "S";
+            ldvt = q;
+        }
+    }
+    DEBUGMSG("svd_l_Rdd");
+    double *B = (double*)malloc(m*n*sizeof(double));
+    CHECK(!B,MEM);
+    memcpy(B,ap,m*n*sizeof(double));
+    integer* iwk = (integer*) malloc(8*q*sizeof(integer));
+    CHECK(!iwk,MEM);
+    integer lwk = -1;
+    integer res;
+    // ask for optimal lwk
+    double ans;
+    dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res);
+    lwk = ans;
+    double * workv = (double*)malloc(lwk*sizeof(double));
+    CHECK(!workv,MEM);
+    dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res);
+    CHECK(res,res);
+    free(iwk);
+    free(workv);
+    free(B);
+    OK
+}
+
+//////////////////// complex svd ////////////////////////////////////
+
+// not in clapack.h
+
+int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
+    doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
+    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)) {
+    integer m = ar;
+    integer n = ac;
+    integer q = MIN(m,n);
+    REQUIRES(sn==q,BAD_SIZE);
+    REQUIRES(up==NULL || (ur==m && (uc==m || uc==q)),BAD_SIZE);
+    char* jobu  = "A";
+    if (up==NULL) {
+        jobu = "N";
+    } else {
+        if (uc==q) {
+            jobu = "S";
+        }
+    }
+    REQUIRES(vp==NULL || (vc==n && (vr==n || vr==q)),BAD_SIZE);
+    char* jobvt  = "A";
+    integer ldvt = n;
+    if (vp==NULL) {
+        jobvt = "N";
+    } else {
+        if (vr==q) {
+            jobvt = "S";
+            ldvt = q;
+        }
+    }DEBUGMSG("svd_l_C");
+    doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
+    CHECK(!B,MEM);
+    memcpy(B,ap,m*n*sizeof(doublecomplex));
+
+    double *rwork = (double*) malloc(5*q*sizeof(double));
+    CHECK(!rwork,MEM);
+    integer lwork = -1;
+    integer res;
+    // ask for optimal lwork
+    doublecomplex ans;
+    zgesvd_ (jobu,jobvt,
+             &m,&n,B,&m,
+             sp,
+             up,&m,
+             vp,&ldvt,
+             &ans, &lwork,
+             rwork,
+             &res);
+    lwork = ceil(ans.r);
+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    CHECK(!work,MEM);
+    zgesvd_ (jobu,jobvt,
+             &m,&n,B,&m,
+             sp,
+             up,&m,
+             vp,&ldvt,
+             work, &lwork,
+             rwork,
+             &res);
+    CHECK(res,res);
+    free(work);
+    free(rwork);
+    free(B);
+    OK
+}
+
+int zgesdd_ (char *jobz, integer *m, integer *n,
+    doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
+    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)) {
+    //printf("entro\n");
+    integer m = ar;
+    integer n = ac;
+    integer q = MIN(m,n);
+    REQUIRES(sn==q,BAD_SIZE);
+    REQUIRES((up == NULL && vp == NULL)
+             || (ur==m && vc==n
+                &&   ((uc == q && vr == q)
+                   || (uc == m && vc==n))),BAD_SIZE);
+    char* jobz  = "A";
+    integer ldvt = n;
+    if (up==NULL) {
+        jobz = "N";
+    } else {
+        if (uc==q && vr == q) {
+            jobz = "S";
+            ldvt = q;
+        }
+    }
+    DEBUGMSG("svd_l_Cdd");
+    doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
+    CHECK(!B,MEM);
+    memcpy(B,ap,m*n*sizeof(doublecomplex));
+    integer* iwk = (integer*) malloc(8*q*sizeof(integer));
+    CHECK(!iwk,MEM);
+    int lrwk;
+    if (0 && *jobz == 'N') {
+        lrwk = 5*q; // does not work, crash at free below
+    } else {
+        lrwk = 5*q*q + 7*q;
+    }
+    double *rwk = (double*)malloc(lrwk*sizeof(double));;
+    CHECK(!rwk,MEM);
+    //printf("%s %ld %d\n",jobz,q,lrwk);
+    integer lwk = -1;
+    integer res;
+    // ask for optimal lwk
+    doublecomplex ans;
+    zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res);
+    lwk = ans.r;
+    //printf("lwk = %ld\n",lwk);
+    doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex));
+    CHECK(!workv,MEM);
+    zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res);
+    //printf("res = %ld\n",res);
+    CHECK(res,res);
+    free(workv); // printf("freed workv\n");
+    free(rwk);   // printf("freed rwk\n");
+    free(iwk);   // printf("freed iwk\n");
+    free(B);     // printf("freed B, salgo\n");
+    OK
+}
+
+//////////////////// general complex eigensystem ////////////
+
+/* Subroutine */ int zgeev_(char *jobvl, char *jobvr, integer *n,
+        doublecomplex *a, integer *lda, doublecomplex *w, doublecomplex *vl,
+        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)) {
+    integer n = ar;
+    REQUIRES(ac==n && sn==n, BAD_SIZE);
+    REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
+    char jobvl = up==NULL?'N':'V';
+    REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE);
+    char jobvr = vp==NULL?'N':'V';
+    DEBUGMSG("eig_l_C");
+    doublecomplex *B = (doublecomplex*)malloc(n*n*sizeof(doublecomplex));
+    CHECK(!B,MEM);
+    memcpy(B,ap,n*n*sizeof(doublecomplex));
+    double *rwork = (double*) malloc(2*n*sizeof(double));
+    CHECK(!rwork,MEM);
+    integer lwork = -1;
+    integer res;
+    // ask for optimal lwork
+    doublecomplex ans;
+    //printf("ask zgeev\n");
+    zgeev_  (&jobvl,&jobvr,
+             &n,B,&n,
+             sp,
+             up,&n,
+             vp,&n,
+             &ans, &lwork,
+             rwork,
+             &res);
+    lwork = ceil(ans.r);
+    //printf("ans = %d\n",lwork);
+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    CHECK(!work,MEM);
+    //printf("zgeev\n");
+    zgeev_  (&jobvl,&jobvr,
+             &n,B,&n,
+             sp,
+             up,&n,
+             vp,&n,
+             work, &lwork,
+             rwork,
+             &res);
+    CHECK(res,res);
+    free(work);
+    free(rwork);
+    free(B);
+    OK
+}
+
+
+
+//////////////////// general real eigensystem ////////////
+
+/* Subroutine */ int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *
+        a, integer *lda, doublereal *wr, doublereal *wi, doublereal *vl,
+        integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work,
+        integer *lwork, integer *info);
+
+int eig_l_R(KDMAT(a),DMAT(u), CVEC(s),DMAT(v)) {
+    integer n = ar;
+    REQUIRES(ac==n && sn==n, BAD_SIZE);
+    REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
+    char jobvl = up==NULL?'N':'V';
+    REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE);
+    char jobvr = vp==NULL?'N':'V';
+    DEBUGMSG("eig_l_R");
+    double *B = (double*)malloc(n*n*sizeof(double));
+    CHECK(!B,MEM);
+    memcpy(B,ap,n*n*sizeof(double));
+    integer lwork = -1;
+    integer res;
+    // ask for optimal lwork
+    double ans;
+    //printf("ask dgeev\n");
+    dgeev_  (&jobvl,&jobvr,
+             &n,B,&n,
+             (double*)sp, (double*)sp+n,
+             up,&n,
+             vp,&n,
+             &ans, &lwork,
+             &res);
+    lwork = ceil(ans);
+    //printf("ans = %d\n",lwork);
+    double * work = (double*)malloc(lwork*sizeof(double));
+    CHECK(!work,MEM);
+    //printf("dgeev\n");
+    dgeev_  (&jobvl,&jobvr,
+             &n,B,&n,
+             (double*)sp, (double*)sp+n,
+             up,&n,
+             vp,&n,
+             work, &lwork,
+             &res);
+    CHECK(res,res);
+    free(work);
+    free(B);
+    OK
+}
+
+
+//////////////////// symmetric real eigensystem ////////////
+
+/* Subroutine */ int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a,
+         integer *lda, doublereal *w, doublereal *work, integer *lwork,
+        integer *info);
+
+int eig_l_S(int wantV,KDMAT(a),DVEC(s),DMAT(v)) {
+    integer n = ar;
+    REQUIRES(ac==n && sn==n, BAD_SIZE);
+    REQUIRES(vr==n && vc==n, BAD_SIZE);
+    char jobz = wantV?'V':'N';
+    DEBUGMSG("eig_l_S");
+    memcpy(vp,ap,n*n*sizeof(double));
+    integer lwork = -1;
+    char uplo = 'U';
+    integer res;
+    // ask for optimal lwork
+    double ans;
+    //printf("ask dsyev\n");
+    dsyev_  (&jobz,&uplo,
+             &n,vp,&n,
+             sp,
+             &ans, &lwork,
+             &res);
+    lwork = ceil(ans);
+    //printf("ans = %d\n",lwork);
+    double * work = (double*)malloc(lwork*sizeof(double));
+    CHECK(!work,MEM);
+    dsyev_  (&jobz,&uplo,
+             &n,vp,&n,
+             sp,
+             work, &lwork,
+             &res);
+    CHECK(res,res);
+    free(work);
+    OK
+}
+
+//////////////////// hermitian complex eigensystem ////////////
+
+/* Subroutine */ int zheev_(char *jobz, char *uplo, integer *n, doublecomplex
+        *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)) {
+    integer n = ar;
+    REQUIRES(ac==n && sn==n, BAD_SIZE);
+    REQUIRES(vr==n && vc==n, BAD_SIZE);
+    char jobz = wantV?'V':'N';
+    DEBUGMSG("eig_l_H");
+    memcpy(vp,ap,2*n*n*sizeof(double));
+    double *rwork = (double*) malloc((3*n-2)*sizeof(double));
+    CHECK(!rwork,MEM);
+    integer lwork = -1;
+    char uplo = 'U';
+    integer res;
+    // ask for optimal lwork
+    doublecomplex ans;
+    //printf("ask zheev\n");
+    zheev_  (&jobz,&uplo,
+             &n,vp,&n,
+             sp,
+             &ans, &lwork,
+             rwork,
+             &res);
+    lwork = ceil(ans.r);
+    //printf("ans = %d\n",lwork);
+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    CHECK(!work,MEM);
+    zheev_  (&jobz,&uplo,
+             &n,vp,&n,
+             sp,
+             work, &lwork,
+             rwork,
+             &res);
+    CHECK(res,res);
+    free(work);
+    free(rwork);
+    OK
+}
+
+//////////////////// general real linear system ////////////
+
+/* 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)) {
+    integer n = ar;
+    integer nhrs = bc;
+    REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
+    DEBUGMSG("linearSolveR_l");
+    double*AC = (double*)malloc(n*n*sizeof(double));
+    memcpy(AC,ap,n*n*sizeof(double));
+    memcpy(xp,bp,n*nhrs*sizeof(double));
+    integer * ipiv = (integer*)malloc(n*sizeof(integer));
+    integer res;
+    dgesv_  (&n,&nhrs,
+             AC, &n,
+             ipiv,
+             xp, &n,
+             &res);
+    if(res>0) {
+        return SINGULAR;
+    }
+    CHECK(res,res);
+    free(ipiv);
+    free(AC);
+    OK
+}
+
+//////////////////// general complex linear system ////////////
+
+/* Subroutine */ int zgesv_(integer *n, integer *nrhs, doublecomplex *a,
+        integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer *
+        info);
+
+int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
+    integer n = ar;
+    integer nhrs = bc;
+    REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
+    DEBUGMSG("linearSolveC_l");
+    doublecomplex*AC = (doublecomplex*)malloc(n*n*sizeof(doublecomplex));
+    memcpy(AC,ap,n*n*sizeof(doublecomplex));
+    memcpy(xp,bp,n*nhrs*sizeof(doublecomplex));
+    integer * ipiv = (integer*)malloc(n*sizeof(integer));
+    integer res;
+    zgesv_  (&n,&nhrs,
+             AC, &n,
+             ipiv,
+             xp, &n,
+             &res);
+    if(res>0) {
+        return SINGULAR;
+    }
+    CHECK(res,res);
+    free(ipiv);
+    free(AC);
+    OK
+}
+
+//////// symmetric positive definite real linear system using Cholesky ////////////
+
+/* Subroutine */ int dpotrs_(char *uplo, integer *n, integer *nrhs,
+        doublereal *a, integer *lda, doublereal *b, integer *ldb, integer *
+        info);
+
+int cholSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
+    integer n = ar;
+    integer nhrs = bc;
+    REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
+    DEBUGMSG("cholSolveR_l");
+    memcpy(xp,bp,n*nhrs*sizeof(double));
+    integer res;
+    dpotrs_ ("U",
+             &n,&nhrs,
+             (double*)ap, &n,
+             xp, &n,
+             &res);
+    CHECK(res,res);
+    OK
+}
+
+//////// Hermitian positive definite real linear system using Cholesky ////////////
+
+/* Subroutine */ int zpotrs_(char *uplo, integer *n, integer *nrhs,
+        doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
+        integer *info);
+
+int cholSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
+    integer n = ar;
+    integer nhrs = bc;
+    REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
+    DEBUGMSG("cholSolveC_l");
+    memcpy(xp,bp,n*nhrs*sizeof(doublecomplex));
+    integer res;
+    zpotrs_  ("U",
+             &n,&nhrs,
+             (doublecomplex*)ap, &n,
+             xp, &n,
+             &res);
+    CHECK(res,res);
+    OK
+}
+
+//////////////////// least squares real linear system ////////////
+
+/* Subroutine */ int dgels_(char *trans, integer *m, integer *n, integer *
+        nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb,
+        doublereal *work, integer *lwork, integer *info);
+
+int linearSolveLSR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
+    integer m = ar;
+    integer n = ac;
+    integer nrhs = bc;
+    integer ldb = xr;
+    REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);
+    DEBUGMSG("linearSolveLSR_l");
+    double*AC = (double*)malloc(m*n*sizeof(double));
+    memcpy(AC,ap,m*n*sizeof(double));
+    if (m>=n) {
+        memcpy(xp,bp,m*nrhs*sizeof(double));
+    } else {
+        int k;
+        for(k = 0; k<nrhs; k++) {
+            memcpy(xp+ldb*k,bp+m*k,m*sizeof(double));
+        }
+    }
+    integer res;
+    integer lwork = -1;
+    double ans;
+    dgels_  ("N",&m,&n,&nrhs,
+             AC,&m,
+             xp,&ldb,
+             &ans,&lwork,
+             &res);
+    lwork = ceil(ans);
+    //printf("ans = %d\n",lwork);
+    double * work = (double*)malloc(lwork*sizeof(double));
+    dgels_  ("N",&m,&n,&nrhs,
+             AC,&m,
+             xp,&ldb,
+             work,&lwork,
+             &res);
+    if(res>0) {
+        return SINGULAR;
+    }
+    CHECK(res,res);
+    free(work);
+    free(AC);
+    OK
+}
+
+//////////////////// least squares complex linear system ////////////
+
+/* Subroutine */ int zgels_(char *trans, integer *m, integer *n, integer *
+        nrhs, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
+        doublecomplex *work, integer *lwork, integer *info);
+
+int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
+    integer m = ar;
+    integer n = ac;
+    integer nrhs = bc;
+    integer ldb = xr;
+    REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);
+    DEBUGMSG("linearSolveLSC_l");
+    doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
+    memcpy(AC,ap,m*n*sizeof(doublecomplex));
+    if (m>=n) {
+        memcpy(xp,bp,m*nrhs*sizeof(doublecomplex));
+    } else {
+        int k;
+        for(k = 0; k<nrhs; k++) {
+            memcpy(xp+ldb*k,bp+m*k,m*sizeof(doublecomplex));
+        }
+    }
+    integer res;
+    integer lwork = -1;
+    doublecomplex ans;
+    zgels_  ("N",&m,&n,&nrhs,
+             AC,&m,
+             xp,&ldb,
+             &ans,&lwork,
+             &res);
+    lwork = ceil(ans.r);
+    //printf("ans = %d\n",lwork);
+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    zgels_  ("N",&m,&n,&nrhs,
+             AC,&m,
+             xp,&ldb,
+             work,&lwork,
+             &res);
+    if(res>0) {
+        return SINGULAR;
+    }
+    CHECK(res,res);
+    free(work);
+    free(AC);
+    OK
+}
+
+//////////////////// least squares real linear system using SVD ////////////
+
+/* Subroutine */ int dgelss_(integer *m, integer *n, integer *nrhs,
+        doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *
+        s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork,
+         integer *info);
+
+int linearSolveSVDR_l(double rcond,KDMAT(a),KDMAT(b),DMAT(x)) {
+    integer m = ar;
+    integer n = ac;
+    integer nrhs = bc;
+    integer ldb = xr;
+    REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);
+    DEBUGMSG("linearSolveSVDR_l");
+    double*AC = (double*)malloc(m*n*sizeof(double));
+    double*S = (double*)malloc(MIN(m,n)*sizeof(double));
+    memcpy(AC,ap,m*n*sizeof(double));
+    if (m>=n) {
+        memcpy(xp,bp,m*nrhs*sizeof(double));
+    } else {
+        int k;
+        for(k = 0; k<nrhs; k++) {
+            memcpy(xp+ldb*k,bp+m*k,m*sizeof(double));
+        }
+    }
+    integer res;
+    integer lwork = -1;
+    integer rank;
+    double ans;
+    dgelss_  (&m,&n,&nrhs,
+             AC,&m,
+             xp,&ldb,
+             S,
+             &rcond,&rank,
+             &ans,&lwork,
+             &res);
+    lwork = ceil(ans);
+    //printf("ans = %d\n",lwork);
+    double * work = (double*)malloc(lwork*sizeof(double));
+    dgelss_  (&m,&n,&nrhs,
+             AC,&m,
+             xp,&ldb,
+             S,
+             &rcond,&rank,
+             work,&lwork,
+             &res);
+    if(res>0) {
+        return NOCONVER;
+    }
+    CHECK(res,res);
+    free(work);
+    free(S);
+    free(AC);
+    OK
+}
+
+//////////////////// least squares complex linear system using SVD ////////////
+
+// not in clapack.h
+
+int zgelss_(integer *m, integer *n, integer *nhrs,
+    doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, doublereal *s,
+    doublereal *rcond, integer* rank,
+    doublecomplex *work, integer* lwork, doublereal* rwork,
+    integer *info);
+
+int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) {
+    integer m = ar;
+    integer n = ac;
+    integer nrhs = bc;
+    integer ldb = xr;
+    REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);
+    DEBUGMSG("linearSolveSVDC_l");
+    doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
+    double*S = (double*)malloc(MIN(m,n)*sizeof(double));
+    double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double));
+    memcpy(AC,ap,m*n*sizeof(doublecomplex));
+    if (m>=n) {
+        memcpy(xp,bp,m*nrhs*sizeof(doublecomplex));
+    } else {
+        int k;
+        for(k = 0; k<nrhs; k++) {
+            memcpy(xp+ldb*k,bp+m*k,m*sizeof(doublecomplex));
+        }
+    }
+    integer res;
+    integer lwork = -1;
+    integer rank;
+    doublecomplex ans;
+    zgelss_  (&m,&n,&nrhs,
+             AC,&m,
+             xp,&ldb,
+             S,
+             &rcond,&rank,
+             &ans,&lwork,
+             RWORK,
+             &res);
+    lwork = ceil(ans.r);
+    //printf("ans = %d\n",lwork);
+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    zgelss_  (&m,&n,&nrhs,
+             AC,&m,
+             xp,&ldb,
+             S,
+             &rcond,&rank,
+             work,&lwork,
+             RWORK,
+             &res);
+    if(res>0) {
+        return NOCONVER;
+    }
+    CHECK(res,res);
+    free(work);
+    free(RWORK);
+    free(S);
+    free(AC);
+    OK
+}
+
+//////////////////// Cholesky factorization /////////////////////////
+
+/* Subroutine */ int zpotrf_(char *uplo, integer *n, doublecomplex *a,
+        integer *lda, integer *info);
+
+int chol_l_H(KCMAT(a),CMAT(l)) {
+    integer n = ar;
+    REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE);
+    DEBUGMSG("chol_l_H");
+    memcpy(lp,ap,n*n*sizeof(doublecomplex));
+    char uplo = 'U';
+    integer res;
+    zpotrf_ (&uplo,&n,lp,&n,&res);
+    CHECK(res>0,NODEFPOS);
+    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;
+        }
+    }
+    OK
+}
+
+
+/* Subroutine */ int dpotrf_(char *uplo, integer *n, doublereal *a, integer *
+        lda, integer *info);
+
+int chol_l_S(KDMAT(a),DMAT(l)) {
+    integer n = ar;
+    REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE);
+    DEBUGMSG("chol_l_S");
+    memcpy(lp,ap,n*n*sizeof(double));
+    char uplo = 'U';
+    integer res;
+    dpotrf_ (&uplo,&n,lp,&n,&res);
+    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.;
+        }
+    }
+    OK
+}
+
+//////////////////// QR factorization /////////////////////////
+
+/* 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)) {
+    integer m = ar;
+    integer n = ac;
+    integer mn = MIN(m,n);
+    REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE);
+    DEBUGMSG("qr_l_R");
+    double *WORK = (double*)malloc(n*sizeof(double));
+    CHECK(!WORK,MEM);
+    memcpy(rp,ap,m*n*sizeof(double));
+    integer res;
+    dgeqr2_ (&m,&n,rp,&m,taup,WORK,&res);
+    CHECK(res,res);
+    free(WORK);
+    OK
+}
+
+/* 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)) {
+    integer m = ar;
+    integer n = ac;
+    integer mn = MIN(m,n);
+    REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE);
+    DEBUGMSG("qr_l_C");
+    doublecomplex *WORK = (doublecomplex*)malloc(n*sizeof(doublecomplex));
+    CHECK(!WORK,MEM);
+    memcpy(rp,ap,m*n*sizeof(doublecomplex));
+    integer res;
+    zgeqr2_ (&m,&n,rp,&m,taup,WORK,&res);
+    CHECK(res,res);
+    free(WORK);
+    OK
+}
+
+/* Subroutine */ int dorgqr_(integer *m, integer *n, integer *k, doublereal *
+        a, integer *lda, doublereal *tau, doublereal *work, integer *lwork,
+        integer *info);
+
+int c_dorgqr(KDMAT(a), KDVEC(tau), DMAT(r)) {
+    integer m = ar;
+    integer n = MIN(ac,ar);
+    integer k = taun;
+    DEBUGMSG("c_dorgqr");
+    integer lwork = 8*n; // FIXME
+    double *WORK = (double*)malloc(lwork*sizeof(double));
+    CHECK(!WORK,MEM);
+    memcpy(rp,ap,m*k*sizeof(double));
+    integer res;
+    dorgqr_ (&m,&n,&k,rp,&m,(double*)taup,WORK,&lwork,&res);
+    CHECK(res,res);
+    free(WORK);
+    OK
+}
+
+/* Subroutine */ int zungqr_(integer *m, integer *n, integer *k,
+        doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *
+        work, integer *lwork, integer *info);
+
+int c_zungqr(KCMAT(a), KCVEC(tau), CMAT(r)) {
+    integer m = ar;
+    integer n = MIN(ac,ar);
+    integer k = taun;
+    DEBUGMSG("z_ungqr");
+    integer lwork = 8*n; // FIXME
+    doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    CHECK(!WORK,MEM);
+    memcpy(rp,ap,m*k*sizeof(doublecomplex));
+    integer res;
+    zungqr_ (&m,&n,&k,rp,&m,(doublecomplex*)taup,WORK,&lwork,&res);
+    CHECK(res,res);
+    free(WORK);
+    OK
+}
+
+
+//////////////////// Hessenberg factorization /////////////////////////
+
+/* Subroutine */ int dgehrd_(integer *n, integer *ilo, integer *ihi,
+        doublereal *a, integer *lda, doublereal *tau, doublereal *work,
+        integer *lwork, integer *info);
+
+int hess_l_R(KDMAT(a), DVEC(tau), DMAT(r)) {
+    integer m = ar;
+    integer n = ac;
+    integer mn = MIN(m,n);
+    REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE);
+    DEBUGMSG("hess_l_R");
+    integer lwork = 5*n; // fixme
+    double *WORK = (double*)malloc(lwork*sizeof(double));
+    CHECK(!WORK,MEM);
+    memcpy(rp,ap,m*n*sizeof(double));
+    integer res;
+    integer one = 1;
+    dgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res);
+    CHECK(res,res);
+    free(WORK);
+    OK
+}
+
+
+/* Subroutine */ int zgehrd_(integer *n, integer *ilo, integer *ihi,
+        doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *
+        work, integer *lwork, integer *info);
+
+int hess_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
+    integer m = ar;
+    integer n = ac;
+    integer mn = MIN(m,n);
+    REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE);
+    DEBUGMSG("hess_l_C");
+    integer lwork = 5*n; // fixme
+    doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    CHECK(!WORK,MEM);
+    memcpy(rp,ap,m*n*sizeof(doublecomplex));
+    integer res;
+    integer one = 1;
+    zgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res);
+    CHECK(res,res);
+    free(WORK);
+    OK
+}
+
+//////////////////// Schur factorization /////////////////////////
+
+/* Subroutine */ int dgees_(char *jobvs, char *sort, L_fp select, integer *n,
+        doublereal *a, integer *lda, integer *sdim, doublereal *wr,
+        doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work,
+        integer *lwork, logical *bwork, integer *info);
+
+int schur_l_R(KDMAT(a), DMAT(u), DMAT(s)) {
+    integer m = ar;
+    integer n = ac;
+    REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE);
+    DEBUGMSG("schur_l_R");
+    //int k;
+    //printf("---------------------------\n");
+    //printf("%p: ",ap); for(k=0;k<n*n;k++) printf("%f ",ap[k]); printf("\n");
+    //printf("%p: ",up); for(k=0;k<n*n;k++) printf("%f ",up[k]); printf("\n");
+    //printf("%p: ",sp); for(k=0;k<n*n;k++) printf("%f ",sp[k]); printf("\n");
+    memcpy(sp,ap,n*n*sizeof(double));
+    integer lwork = 6*n; // fixme
+    double *WORK = (double*)malloc(lwork*sizeof(double));
+    double *WR = (double*)malloc(n*sizeof(double));
+    double *WI = (double*)malloc(n*sizeof(double));
+    // WR and WI not really required in this call
+    logical *BWORK = (logical*)malloc(n*sizeof(logical));
+    integer res;
+    integer sdim;
+    dgees_ ("V","N",NULL,&n,sp,&n,&sdim,WR,WI,up,&n,WORK,&lwork,BWORK,&res);
+    //printf("%p: ",ap); for(k=0;k<n*n;k++) printf("%f ",ap[k]); printf("\n");
+    //printf("%p: ",up); for(k=0;k<n*n;k++) printf("%f ",up[k]); printf("\n");
+    //printf("%p: ",sp); for(k=0;k<n*n;k++) printf("%f ",sp[k]); printf("\n");
+    if(res>0) {
+        return NOCONVER;
+    }
+    CHECK(res,res);
+    free(WR);
+    free(WI);
+    free(BWORK);
+    free(WORK);
+    OK
+}
+
+
+/* Subroutine */ int zgees_(char *jobvs, char *sort, L_fp select, integer *n,
+        doublecomplex *a, integer *lda, integer *sdim, doublecomplex *w,
+        doublecomplex *vs, integer *ldvs, doublecomplex *work, integer *lwork,
+         doublereal *rwork, logical *bwork, integer *info);
+
+int schur_l_C(KCMAT(a), CMAT(u), CMAT(s)) {
+    integer m = ar;
+    integer n = ac;
+    REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE);
+    DEBUGMSG("schur_l_C");
+    memcpy(sp,ap,n*n*sizeof(doublecomplex));
+    integer lwork = 6*n; // fixme
+    doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    doublecomplex *W = (doublecomplex*)malloc(n*sizeof(doublecomplex));
+    // W not really required in this call
+    logical *BWORK = (logical*)malloc(n*sizeof(logical));
+    double *RWORK = (double*)malloc(n*sizeof(double));
+    integer res;
+    integer sdim;
+    zgees_ ("V","N",NULL,&n,sp,&n,&sdim,W,
+                            up,&n,
+                            WORK,&lwork,RWORK,BWORK,&res);
+    if(res>0) {
+        return NOCONVER;
+    }
+    CHECK(res,res);
+    free(W);
+    free(BWORK);
+    free(WORK);
+    OK
+}
+
+//////////////////// LU factorization /////////////////////////
+
+/* 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)) {
+    integer m = ar;
+    integer n = ac;
+    integer mn = MIN(m,n);
+    REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);
+    DEBUGMSG("lu_l_R");
+    integer* auxipiv = (integer*)malloc(mn*sizeof(integer));
+    memcpy(rp,ap,m*n*sizeof(double));
+    integer res;
+    dgetrf_ (&m,&n,rp,&m,auxipiv,&res);
+    if(res>0) {
+        res = 0; // fixme
+    }
+    CHECK(res,res);
+    int k;
+    for (k=0; k<mn; k++) {
+        ipivp[k] = auxipiv[k];
+    }
+    free(auxipiv);
+    OK
+}
+
+
+/* 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)) {
+    integer m = ar;
+    integer n = ac;
+    integer mn = MIN(m,n);
+    REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);
+    DEBUGMSG("lu_l_C");
+    integer* auxipiv = (integer*)malloc(mn*sizeof(integer));
+    memcpy(rp,ap,m*n*sizeof(doublecomplex));
+    integer res;
+    zgetrf_ (&m,&n,rp,&m,auxipiv,&res);
+    if(res>0) {
+        res = 0; // fixme
+    }
+    CHECK(res,res);
+    int k;
+    for (k=0; k<mn; k++) {
+        ipivp[k] = auxipiv[k];
+    }
+    free(auxipiv);
+    OK
+}
+
+
+//////////////////// LU substitution /////////////////////////
+
+/* Subroutine */ int dgetrs_(char *trans, integer *n, integer *nrhs,
+        doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *
+        ldb, integer *info);
+
+int luS_l_R(KDMAT(a), KDVEC(ipiv), KDMAT(b), DMAT(x)) {
+  integer m = ar;
+  integer n = ac;
+  integer mrhs = br;
+  integer nrhs = bc;
+
+  REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
+  integer* auxipiv = (integer*)malloc(n*sizeof(integer));
+  int k;
+  for (k=0; k<n; k++) {
+    auxipiv[k] = (integer)ipivp[k];
+  }
+  integer res;
+  memcpy(xp,bp,mrhs*nrhs*sizeof(double));
+  dgetrs_ ("N",&n,&nrhs,(/*no const (!?)*/ double*)ap,&m,auxipiv,xp,&mrhs,&res);
+  CHECK(res,res);
+  free(auxipiv);
+  OK
+}
+
+
+/* Subroutine */ int zgetrs_(char *trans, integer *n, integer *nrhs,
+        doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *b,
+        integer *ldb, integer *info);
+
+int luS_l_C(KCMAT(a), KDVEC(ipiv), KCMAT(b), CMAT(x)) {
+    integer m = ar;
+    integer n = ac;
+    integer mrhs = br;
+    integer nrhs = bc;
+
+    REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
+    integer* auxipiv = (integer*)malloc(n*sizeof(integer));
+    int k;
+    for (k=0; k<n; k++) {
+        auxipiv[k] = (integer)ipivp[k];
+    }
+    integer res;
+    memcpy(xp,bp,mrhs*nrhs*sizeof(doublecomplex));
+    zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&m,auxipiv,xp,&mrhs,&res);
+    CHECK(res,res);
+    free(auxipiv);
+    OK
+}
+
+//////////////////// Matrix Product /////////////////////////
+
+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)) {
+    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    DEBUGMSG("dgemm_");
+    CHECKNANR(a,"NaN multR Input\n")
+    CHECKNANR(b,"NaN multR Input\n")
+    integer m = ta?ac:ar;
+    integer n = tb?br:bc;
+    integer k = ta?ar:ac;
+    integer lda = ar;
+    integer ldb = br;
+    integer ldc = rr;
+    double alpha = 1;
+    double beta = 0;
+    dgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc);
+    CHECKNANR(r,"NaN multR Output\n")
+    OK
+}
+
+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)) {
+    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    DEBUGMSG("zgemm_");
+    CHECKNANC(a,"NaN multC Input\n")
+    CHECKNANC(b,"NaN multC Input\n")
+    integer m = ta?ac:ar;
+    integer n = tb?br:bc;
+    integer k = ta?ar:ac;
+    integer lda = ar;
+    integer ldb = br;
+    integer ldc = rr;
+    doublecomplex alpha = {1,0};
+    doublecomplex beta = {0,0};
+    zgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,
+           ap,&lda,
+           bp,&ldb,&beta,
+           rp,&ldc);
+    CHECKNANC(r,"NaN multC Output\n")
+    OK
+}
+
+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)) {
+    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    DEBUGMSG("sgemm_");
+    integer m = ta?ac:ar;
+    integer n = tb?br:bc;
+    integer k = ta?ar:ac;
+    integer lda = ar;
+    integer ldb = br;
+    integer ldc = rr;
+    float alpha = 1;
+    float beta = 0;
+    sgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc);
+    OK
+}
+
+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)) {
+    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    DEBUGMSG("cgemm_");
+    integer m = ta?ac:ar;
+    integer n = tb?br:bc;
+    integer k = ta?ar:ac;
+    integer lda = ar;
+    integer ldb = br;
+    integer ldc = rr;
+    complex alpha = {1,0};
+    complex beta = {0,0};
+    cgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,
+           ap,&lda,
+           bp,&ldb,&beta,
+           rp,&ldc);
+    OK
+}
+
+//////////////////// transpose /////////////////////////
+
+int transF(KFMAT(x),FMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("transF");
+    int i,j;
+    for (i=0; i<tr; i++) {
+        for (j=0; j<tc; j++) {
+        tp[i*tc+j] = xp[j*xc+i];
+        }
+    }
+    OK
+}
+
+int transR(KDMAT(x),DMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("transR");
+    int i,j;
+    for (i=0; i<tr; i++) {
+        for (j=0; j<tc; j++) {
+        tp[i*tc+j] = xp[j*xc+i];
+        }
+    }
+    OK
+}
+
+int transQ(KQMAT(x),QMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("transQ");
+    int i,j;
+    for (i=0; i<tr; i++) {
+        for (j=0; j<tc; j++) {
+        tp[i*tc+j] = xp[j*xc+i];
+        }
+    }
+    OK
+}
+
+int transC(KCMAT(x),CMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("transC");
+    int i,j;
+    for (i=0; i<tr; i++) {
+        for (j=0; j<tc; j++) {
+        tp[i*tc+j] = xp[j*xc+i];
+        }
+    }
+    OK
+}
+
+int transP(KPMAT(x), PMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    REQUIRES(xs==ts,NOCONVER);
+    DEBUGMSG("transP");
+    int i,j;
+    for (i=0; i<tr; i++) {
+        for (j=0; j<tc; j++) {
+          memcpy(tp+(i*tc+j)*xs,xp +(j*xc+i)*xs,xs);
+        }
+    }
+    OK
+}
+
+//////////////////// constant /////////////////////////
+
+int constantF(float * pval, FVEC(r)) {
+    DEBUGMSG("constantF")
+    int k;
+    double val = *pval;
+    for(k=0;k<rn;k++) {
+        rp[k]=val;
+    }
+    OK
+}
+
+int constantR(double * pval, DVEC(r)) {
+    DEBUGMSG("constantR")
+    int k;
+    double val = *pval;
+    for(k=0;k<rn;k++) {
+        rp[k]=val;
+    }
+    OK
+}
+
+int constantQ(complex* pval, QVEC(r)) {
+    DEBUGMSG("constantQ")
+    int k;
+    complex val = *pval;
+    for(k=0;k<rn;k++) {
+        rp[k]=val;
+    }
+    OK
+}
+
+int constantC(doublecomplex* pval, CVEC(r)) {
+    DEBUGMSG("constantC")
+    int k;
+    doublecomplex val = *pval;
+    for(k=0;k<rn;k++) {
+        rp[k]=val;
+    }
+    OK
+}
+
+int constantP(void* pval, PVEC(r)) {
+    DEBUGMSG("constantP")
+    int k;
+    for(k=0;k<rn;k++) {
+      memcpy(rp+k*rs,pval,rs);
+    }
+    OK
+}
+
+//////////////////// float-double conversion /////////////////////////
+
+int float2double(FVEC(x),DVEC(y)) {
+    DEBUGMSG("float2double")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k];
+    }
+    OK
+}
+
+int double2float(DVEC(x),FVEC(y)) {
+    DEBUGMSG("double2float")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k];
+    }
+    OK
+}
+
+//////////////////// conjugate /////////////////////////
+
+int conjugateQ(KQVEC(x),QVEC(t)) {
+    REQUIRES(xn==tn,BAD_SIZE);
+    DEBUGMSG("conjugateQ");
+    int k;
+    for(k=0;k<xn;k++) {
+        tp[k].r =  xp[k].r;
+        tp[k].i = -xp[k].i;
+    }
+    OK
+}
+
+int conjugateC(KCVEC(x),CVEC(t)) {
+    REQUIRES(xn==tn,BAD_SIZE);
+    DEBUGMSG("conjugateC");
+    int k;
+    for(k=0;k<xn;k++) {
+        tp[k].r =  xp[k].r;
+        tp[k].i = -xp[k].i;
+    }
+    OK
+}
+
+//////////////////// step /////////////////////////
+
+int stepF(FVEC(x),FVEC(y)) {
+    DEBUGMSG("stepF")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k]>0;
+    }
+    OK
+}
+
+int stepD(DVEC(x),DVEC(y)) {
+    DEBUGMSG("stepD")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k]>0;
+    }
+    OK
+}
+
+//////////////////// cond /////////////////////////
+
+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")
+    int k;
+    for(k=0;k<xn;k++) {
+        rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]);
+    }
+    OK
+}
+
+int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) {
+    REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE);
+    DEBUGMSG("condD")
+    int k;
+    for(k=0;k<xn;k++) {
+        rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]);
+    }
+    OK
+}
+
diff --git a/packages/base/src/C/lapack-aux.h b/packages/base/src/C/lapack-aux.h
new file mode 100644
index 0000000..a3f1899
--- /dev/null
+++ b/packages/base/src/C/lapack-aux.h
@@ -0,0 +1,60 @@
+/*
+ * We have copied the definitions in f2c.h required
+ * to compile clapack.h, modified to support both
+ * 32 and 64 bit
+
+      http://opengrok.creo.hu/dragonfly/xref/src/contrib/gcc-3.4/libf2c/readme.netlib
+      http://www.ibm.com/developerworks/library/l-port64.html
+ */
+
+#ifdef _LP64
+typedef int integer;
+typedef unsigned int uinteger;
+typedef int logical;
+typedef long longint;           /* system-dependent */
+typedef unsigned long ulongint; /* system-dependent */
+#else
+typedef long int integer;
+typedef unsigned long int uinteger;
+typedef long int logical;
+typedef long long longint;              /* system-dependent */
+typedef unsigned long long ulongint;    /* system-dependent */
+#endif
+
+typedef char *address;
+typedef short int shortint;
+typedef float real;
+typedef double doublereal;
+typedef struct { real r, i; } complex;
+typedef struct { doublereal r, i; } doublecomplex;
+typedef short int shortlogical;
+typedef char logical1;
+typedef char integer1;
+
+typedef logical (*L_fp)();
+typedef short ftnlen;
+
+/********************************************************/
+
+#define FVEC(A) int A##n, float*A##p
+#define DVEC(A) int A##n, double*A##p
+#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 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
+#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 KFVEC(A) int A##n, const float*A##p
+#define KDVEC(A) int A##n, const double*A##p
+#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 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
+#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
+
diff --git a/packages/hmatrix/src/Data/Packed.hs b/packages/base/src/Data/Packed.hs
index 957aab8..c66718a 100644
--- a/packages/hmatrix/src/Data/Packed.hs
+++ b/packages/base/src/Data/Packed.hs
@@ -1,29 +1,25 @@
 -----------------------------------------------------------------------------
 {- |
 Module      :  Data.Packed
-Copyright   :  (c) Alberto Ruiz 2006-2010
-License     :  GPL-style
-
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Copyright   :  (c) Alberto Ruiz 2006-2014
+License     :  BSD3
+Maintainer  :  Alberto Ruiz
 Stability   :  provisional
-Portability :  uses ffi
 
 Types for dense 'Vector' and 'Matrix' of 'Storable' elements.
 
 -}
 -----------------------------------------------------------------------------
-{-# OPTIONS_HADDOCK hide #-}
 
 module Data.Packed (
+    -- * Vector
+    --
+    -- | Vectors are @Data.Vector.Storable.Vector@ from the \"vector\" package.
     module Data.Packed.Vector,
+    -- * Matrix
     module Data.Packed.Matrix,
---    module Numeric.Conversion,
---    module Data.Packed.Random,
---    module Data.Complex
 ) where
 
 import Data.Packed.Vector
 import Data.Packed.Matrix
---import Data.Packed.Random
---import Data.Complex
---import Numeric.Conversion
+
diff --git a/packages/hmatrix/src/Data/Packed/Development.hs b/packages/base/src/Data/Packed/Development.hs
index 471e560..777b6c5 100644
--- a/packages/hmatrix/src/Data/Packed/Development.hs
+++ b/packages/base/src/Data/Packed/Development.hs
@@ -3,9 +3,8 @@
 -- |
 -- Module      :  Data.Packed.Development
 -- Copyright   :  (c) Alberto Ruiz 2009
--- License     :  GPL
---
--- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- License     :  BSD3
+-- Maintainer  :  Alberto Ruiz
 -- Stability   :  provisional
 -- Portability :  portable
 --
@@ -14,7 +13,6 @@
 -- in the @examples\/devel@ folder included in the package.
 --
 -----------------------------------------------------------------------------
-{-# OPTIONS_HADDOCK hide #-}
 
 module Data.Packed.Development (
     createVector, createMatrix,
@@ -30,3 +28,4 @@ module Data.Packed.Development (
 ) where
 
 import Data.Packed.Internal
+
diff --git a/packages/hmatrix/src/Data/Packed/Foreign.hs b/packages/base/src/Data/Packed/Foreign.hs
index 1ec3694..efa51ca 100644
--- a/packages/hmatrix/src/Data/Packed/Foreign.hs
+++ b/packages/base/src/Data/Packed/Foreign.hs
@@ -6,7 +6,6 @@
 -- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3) 
 -- @
 --
-{-# OPTIONS_HADDOCK hide #-}
 module Data.Packed.Foreign 
     ( app
     , appVector, appVectorLen
diff --git a/packages/hmatrix/src/Data/Packed/Internal.hs b/packages/base/src/Data/Packed/Internal.hs
index 537e51e..537e51e 100644
--- a/packages/hmatrix/src/Data/Packed/Internal.hs
+++ b/packages/base/src/Data/Packed/Internal.hs
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Common.hs b/packages/base/src/Data/Packed/Internal/Common.hs
index edef3c2..615bbdf 100644
--- a/packages/hmatrix/src/Data/Packed/Internal/Common.hs
+++ b/packages/base/src/Data/Packed/Internal/Common.hs
@@ -1,18 +1,15 @@
 {-# LANGUAGE CPP #-}
------------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Internal.Common
 -- Copyright   :  (c) Alberto Ruiz 2007
--- License     :  GPL-style
---
--- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- License     :  BSD3
+-- Maintainer  :  Alberto Ruiz
 -- Stability   :  provisional
--- Portability :  portable (uses FFI)
+--
 --
 -- Development utilities.
 --
------------------------------------------------------------------------------
--- #hide
+
 
 module Data.Packed.Internal.Common(
   Adapt,
@@ -21,12 +18,11 @@ module Data.Packed.Internal.Common(
   (//), check, mbCatch,
   splitEvery, common, compatdim,
   fi,
-  table
+  table,
+  finit
 ) where
 
-import Foreign
 import Control.Monad(when)
-import Foreign.C.String(peekCString)
 import Foreign.C.Types
 import Foreign.Storable.Complex()
 import Data.List(transpose,intersperse)
@@ -153,19 +149,12 @@ check msg f = do
     finit
 #endif
     err <- f
-    when (err/=0) $ if err > 1024
-                      then (error (msg++": "++errorCode err)) -- our errors
-                      else do                                 -- GSL errors
-                        ps <- gsl_strerror err
-                        s <- peekCString ps
-                        error (msg++": "++s)
+    when (err/=0) $ error (msg++": "++errorCode err)
     return ()
 
--- | description of GSL error codes
-foreign import ccall unsafe "gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar)
-
 -- | Error capture and conversion to Maybe
 mbCatch :: IO x -> IO (Maybe x)
 mbCatch act = E.catch (Just `fmap` act) f
     where f :: SomeException -> IO (Maybe x)
           f _ = return Nothing
+
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs b/packages/base/src/Data/Packed/Internal/Matrix.hs
index 9719fc0..9b831cc 100644
--- a/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs
+++ b/packages/base/src/Data/Packed/Internal/Matrix.hs
@@ -2,20 +2,16 @@
 {-# LANGUAGE FlexibleContexts         #-}
 {-# LANGUAGE FlexibleInstances        #-}
 {-# LANGUAGE BangPatterns             #-}
------------------------------------------------------------------------------
+
 -- |
 -- Module      :  Data.Packed.Internal.Matrix
 -- Copyright   :  (c) Alberto Ruiz 2007
--- License     :  GPL-style
---
--- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- License     :  BSD3
+-- Maintainer  :  Alberto Ruiz
 -- Stability   :  provisional
--- Portability :  portable (uses FFI)
 --
 -- Internal matrix representation
 --
------------------------------------------------------------------------------
--- #hide
 
 module Data.Packed.Internal.Matrix(
     Matrix(..), rows, cols, cdat, fdat,
@@ -30,7 +26,6 @@ module Data.Packed.Internal.Matrix(
     subMatrix,
     liftMatrix, liftMatrix2,
     (@@>), atM',
-    saveMatrix,
     singleton,
     emptyM,
     size, shSize, conformVs, conformMs, conformVTo, conformMTo
@@ -46,7 +41,6 @@ import Foreign.Ptr(Ptr, castPtr)
 import Foreign.Storable(Storable, peekElemOff, pokeElemOff, poke, sizeOf)
 import Data.Complex(Complex)
 import Foreign.C.Types
-import Foreign.C.String(newCString)
 import System.IO.Unsafe(unsafePerformIO)
 import Control.DeepSeq
 
@@ -290,34 +284,6 @@ instance Element (Complex Double) where
 
 -------------------------------------------------------------------
 
-transdata' :: Storable a => Int -> Vector a -> Int -> Vector a
-transdata' c1 v c2 =
-    if noneed
-        then v
-        else unsafePerformIO $ do
-                w <- createVector (r2*c2)
-                unsafeWith v $ \p ->
-                    unsafeWith w $ \q -> do
-                        let go (-1) _ = return ()
-                            go !i (-1) = go (i-1) (c1-1)
-                            go !i !j = do x <- peekElemOff p (i*c1+j)
-                                          pokeElemOff      q (j*c2+i) x
-                                          go i (j-1)
-                        go (r1-1) (c1-1)
-                return w
-  where r1 = dim v `div` c1
-        r2 = dim v `div` c2
-        noneed = dim v == 0 || r1 == 1 || c1 == 1
-
--- {-# SPECIALIZE transdata' :: Int -> Vector Double -> Int ->  Vector Double #-}
--- {-# SPECIALIZE transdata' :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double) #-}
-
--- I don't know how to specialize...
--- The above pragmas only seem to work on top level defs
--- Fortunately everything seems to work using the above class
-
--- C versions, still a little faster:
-
 transdataAux fun c1 d c2 =
     if noneed
         then d
@@ -354,16 +320,6 @@ foreign import ccall unsafe "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -
 
 ----------------------------------------------------------------------
 
-constant' v n = unsafePerformIO $ do
-    w <- createVector n
-    unsafeWith w $ \p -> do
-        let go (-1) = return ()
-            go !k = pokeElemOff p k v >> go (k-1)
-        go (n-1)
-    return w
-
--- C versions
-
 constantAux fun x n = unsafePerformIO $ do
     v <- createVector n
     px <- newArray [x]
@@ -371,20 +327,12 @@ constantAux fun x n = unsafePerformIO $ do
     free px
     return v
 
-constantF :: Float -> Int -> Vector Float
-constantF = constantAux cconstantF
 foreign import ccall unsafe "constantF" cconstantF :: Ptr Float -> TF
 
-constantR :: Double -> Int -> Vector Double
-constantR = constantAux cconstantR
 foreign import ccall unsafe "constantR" cconstantR :: Ptr Double -> TV
 
-constantQ :: Complex Float -> Int -> Vector (Complex Float)
-constantQ = constantAux cconstantQ
 foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV
 
-constantC :: Complex Double -> Int -> Vector (Complex Double)
-constantC = constantAux cconstantC
 foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV
 
 constantP :: Storable a => a -> Int -> Vector a
@@ -429,23 +377,6 @@ subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m)
 
 --------------------------------------------------------------------------
 
--- | Saves a matrix as 2D ASCII table.
-saveMatrix :: FilePath
-           -> String     -- ^ format (%f, %g, %e)
-           -> Matrix Double
-           -> IO ()
-saveMatrix filename fmt m = do
-    charname <- newCString filename
-    charfmt <- newCString fmt
-    let o = if orderOf m == RowMajor then 1 else 0
-    app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
-    free charname
-    free charfmt
-
-foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
-
-----------------------------------------------------------------------
-
 maxZ xs = if minimum xs == 0 then 0 else maximum xs
 
 conformMs ms = map (conformMTo (r,c)) ms
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs b/packages/base/src/Data/Packed/Internal/Signatures.hs
index 2835720..acc3070 100644
--- a/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs
+++ b/packages/base/src/Data/Packed/Internal/Signatures.hs
@@ -1,16 +1,13 @@
------------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Internal.Signatures
 -- Copyright   :  (c) Alberto Ruiz 2009
--- License     :  GPL-style
---
--- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- License     :  BSD3
+-- Maintainer  :  Alberto Ruiz
 -- Stability   :  provisional
--- Portability :  portable (uses FFI)
 --
 -- Signatures of the C functions.
 --
------------------------------------------------------------------------------
+
 
 module Data.Packed.Internal.Signatures where
 
@@ -70,3 +67,4 @@ type TVCV = CInt -> PD -> TCV                   --
 type TCVM = CInt -> PC -> TM                    --
 type TMCVM = CInt -> CInt -> PD -> TCVM         --
 type TMMCVM = CInt -> CInt -> PD -> TMCVM       --
+
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Vector.hs b/packages/base/src/Data/Packed/Internal/Vector.hs
index 6d03438..d0bc143 100644
--- a/packages/hmatrix/src/Data/Packed/Internal/Vector.hs
+++ b/packages/base/src/Data/Packed/Internal/Vector.hs
@@ -1,17 +1,14 @@
 {-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-}
------------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Internal.Vector
 -- Copyright   :  (c) Alberto Ruiz 2007
--- License     :  GPL-style
---
--- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- License     :  BSD3
+-- Maintainer  :  Alberto Ruiz
 -- Stability   :  provisional
--- Portability :  portable (uses FFI)
 --
 -- Vector implementation
 --
------------------------------------------------------------------------------
+--------------------------------------------------------------------------------
 
 module Data.Packed.Internal.Vector (
     Vector, dim,
@@ -24,7 +21,6 @@ module Data.Packed.Internal.Vector (
     asComplex, asReal, float2DoubleV, double2FloatV,
     stepF, stepD, condF, condD,
     conjugateQ, conjugateC,
-    fwriteVector, freadVector, fprintfVector, fscanfVector,
     cloneVector,
     unsafeToForeignPtr,
     unsafeFromForeignPtr,
@@ -33,12 +29,10 @@ module Data.Packed.Internal.Vector (
 
 import Data.Packed.Internal.Common
 import Data.Packed.Internal.Signatures
-import Foreign.Marshal.Alloc(free)
 import Foreign.Marshal.Array(peekArray, copyArray, advancePtr)
 import Foreign.ForeignPtr(ForeignPtr, castForeignPtr)
 import Foreign.Ptr(Ptr)
 import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf)
-import Foreign.C.String
 import Foreign.C.Types
 import Data.Complex
 import Control.Monad(when)
@@ -474,48 +468,4 @@ mapVectorWithIndex f v = unsafePerformIO $ do
     return w
 {-# INLINE mapVectorWithIndex #-}
 
--------------------------------------------------------------------
-
-
--- | Loads a vector from an ASCII file (the number of elements must be known in advance).
-fscanfVector :: FilePath -> Int -> IO (Vector Double)
-fscanfVector filename n = do
-    charname <- newCString filename
-    res <- createVector n
-    app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
-    free charname
-    return res
-
-foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV
-
--- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file.
-fprintfVector :: FilePath -> String -> Vector Double -> IO ()
-fprintfVector filename fmt v = do
-    charname <- newCString filename
-    charfmt <- newCString fmt
-    app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
-    free charname
-    free charfmt
-
-foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV
-
--- | Loads a vector from a binary file (the number of elements must be known in advance).
-freadVector :: FilePath -> Int -> IO (Vector Double)
-freadVector filename n = do
-    charname <- newCString filename
-    res <- createVector n
-    app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
-    free charname
-    return res
-
-foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
-
--- | Saves the elements of a vector to a binary file.
-fwriteVector :: FilePath -> Vector Double -> IO ()
-fwriteVector filename v = do
-    charname <- newCString filename
-    app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
-    free charname
-
-foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
 
diff --git a/packages/hmatrix/src/Data/Packed/Matrix.hs b/packages/base/src/Data/Packed/Matrix.hs
index d94d167..d94d167 100644
--- a/packages/hmatrix/src/Data/Packed/Matrix.hs
+++ b/packages/base/src/Data/Packed/Matrix.hs
diff --git a/packages/hmatrix/src/Data/Packed/ST.hs b/packages/base/src/Data/Packed/ST.hs
index 1cef296..dae457c 100644
--- a/packages/hmatrix/src/Data/Packed/ST.hs
+++ b/packages/base/src/Data/Packed/ST.hs
@@ -6,9 +6,8 @@
 -- |
 -- Module      :  Data.Packed.ST
 -- Copyright   :  (c) Alberto Ruiz 2008
--- License     :  GPL-style
---
--- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- License     :  BSD3
+-- Maintainer  :  Alberto Ruiz
 -- Stability   :  provisional
 -- Portability :  portable
 --
@@ -16,7 +15,6 @@
 -- See examples/inplace.hs in the distribution.
 --
 -----------------------------------------------------------------------------
-{-# OPTIONS_HADDOCK hide #-}
 
 module Data.Packed.ST (
     -- * Mutable Vectors
@@ -177,3 +175,4 @@ newUndefinedMatrix ord r c = unsafeIOToST $ fmap STMatrix $ createMatrix ord r c
 {-# NOINLINE newMatrix #-}
 newMatrix :: Storable t => t -> Int -> Int -> ST s (STMatrix s t)
 newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c)
+
diff --git a/packages/hmatrix/src/Data/Packed/Vector.hs b/packages/base/src/Data/Packed/Vector.hs
index b5a4318..b5a4318 100644
--- a/packages/hmatrix/src/Data/Packed/Vector.hs
+++ b/packages/base/src/Data/Packed/Vector.hs