From b7e57952112eae61054fc9171ddb311fbaca0841 Mon Sep 17 00:00:00 2001
From: Alberto Ruiz <aruiz@um.es>
Date: Thu, 4 Jun 2015 12:39:33 +0200
Subject: move c sources

---
 packages/base/src/Internal/C/lapack-aux.c | 1686 +++++++++++++++++++++++++++++
 packages/base/src/Internal/C/lapack-aux.h |   82 ++
 packages/base/src/Internal/C/vector-aux.c | 1134 +++++++++++++++++++
 3 files changed, 2902 insertions(+)
 create mode 100644 packages/base/src/Internal/C/lapack-aux.c
 create mode 100644 packages/base/src/Internal/C/lapack-aux.h
 create mode 100644 packages/base/src/Internal/C/vector-aux.c

(limited to 'packages/base/src/Internal/C')

diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
new file mode 100644
index 0000000..1402050
--- /dev/null
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -0,0 +1,1686 @@
+#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
+}
+
+
+int multiplyI(KOIMAT(a), KOIMAT(b), OIMAT(r)) {
+    { TRAV(r,i,j) {
+        int k;
+        AT(r,i,j) = 0;
+        for (k=0;k<ac;k++) {
+            AT(r,i,j) += AT(a,i,k) * AT(b,k,j);
+        }
+      }
+    }
+    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
+}
+
+int transI(KIMAT(x),IMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("transI");
+    int i,j;
+    for (i=0; i<tr; i++) {
+        for (j=0; j<tc; j++) {
+            tp[i*tc+j] = xp[j*xc+i];
+        }
+    }
+    OK
+}
+
+
+//////////////////// constant /////////////////////////
+
+int constantF(float * pval, FVEC(r)) {
+    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
+}
+
+
+int constantI(int * pval, IVEC(r)) {
+    DEBUGMSG("constantI")
+    int k;
+    int val = *pval;
+    for(k=0;k<rn;k++) {
+        rp[k]=val;
+    }
+    OK
+}
+
+
+//////////////////// float-double conversion /////////////////////////
+
+int float2double(FVEC(x),DVEC(y)) {
+    DEBUGMSG("float2double")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k];
+    }
+    OK
+}
+
+int float2int(KFVEC(x),IVEC(y)) {
+    DEBUGMSG("float2int")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k];
+    }
+    OK
+}
+
+
+int double2float(DVEC(x),FVEC(y)) {
+    DEBUGMSG("double2float")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k];
+    }
+    OK
+}
+
+
+int double2int(KDVEC(x),IVEC(y)) {
+    DEBUGMSG("double2int")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k];
+    }
+    OK
+}
+
+
+int int2float(KIVEC(x),FVEC(y)) {
+    DEBUGMSG("int2float")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k];
+    }
+    OK
+}
+
+
+int int2double(KIVEC(x),DVEC(y)) {
+    DEBUGMSG("int2double")
+    int k;
+    for(k=0;k<xn;k++) {
+        yp[k]=xp[k];
+    }
+    OK
+}
+
+
+//////////////////// conjugate /////////////////////////
+
+int conjugateQ(KQVEC(x),QVEC(t)) {
+    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 /////////////////////////
+
+#define STEP_IMP         \
+    int k;               \
+    for(k=0;k<xn;k++) {  \
+        yp[k]=xp[k]>0;   \
+    }                    \
+    OK
+
+int stepF(KFVEC(x),FVEC(y)) {
+    STEP_IMP
+}
+
+int stepD(KDVEC(x),DVEC(y)) {
+    STEP_IMP
+}
+
+int stepI(KIVEC(x),IVEC(y)) {
+    STEP_IMP
+}
+
+//////////////////// cond /////////////////////////
+
+#define COMPARE_IMP                               \
+    REQUIRES(xn==yn && xn==rn ,BAD_SIZE);         \
+    int k;                                        \
+    for(k=0;k<xn;k++) {                           \
+        rp[k] = xp[k]<yp[k]?-1:(xp[k]>yp[k]?1:0); \
+    }                                             \
+    OK
+
+
+int compareF(KFVEC(x),KFVEC(y),IVEC(r)) {
+    COMPARE_IMP
+}
+
+int compareD(KDVEC(x),KDVEC(y),IVEC(r)) {
+    COMPARE_IMP
+}
+
+int compareI(KIVEC(x),KIVEC(y),IVEC(r)) {
+    COMPARE_IMP
+}
+
+
+#define COND_IMP                                                            \
+    REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE);  \
+    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 condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) {
+    COND_IMP
+}
+
+int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) {
+    COND_IMP
+}
+
+int condI(KIVEC(x),KIVEC(y),KIVEC(lt),KIVEC(eq),KIVEC(gt),IVEC(r)) {
+    COND_IMP
+}
+
+
+#define CHOOSE_IMP                                                      \
+    REQUIRES(condn==ltn && ltn==eqn && ltn==gtn && ltn==rn ,BAD_SIZE);  \
+    int k;                                                              \
+    for(k=0;k<condn;k++) {                                              \
+        rp[k] = condp[k]<0?ltp[k]:(condp[k]>0?gtp[k]:eqp[k]);           \
+    }                                                                   \
+    OK
+
+int chooseF(KIVEC(cond),KFVEC(lt),KFVEC(eq),KFVEC(gt),FVEC(r)) {
+    CHOOSE_IMP
+}
+
+int chooseD(KIVEC(cond),KDVEC(lt),KDVEC(eq),KDVEC(gt),DVEC(r)) {
+    CHOOSE_IMP
+}
+
+int chooseI(KIVEC(cond),KIVEC(lt),KIVEC(eq),KIVEC(gt),IVEC(r)) {
+    CHOOSE_IMP
+}
+
+int chooseC(KIVEC(cond),KCVEC(lt),KCVEC(eq),KCVEC(gt),CVEC(r)) {
+    CHOOSE_IMP
+}
+
+int chooseQ(KIVEC(cond),KQVEC(lt),KQVEC(eq),KQVEC(gt),QVEC(r)) {
+    CHOOSE_IMP
+}
+
+//////////////////////// extract /////////////////////////////////
+
+#define EXTRACT_IMP                        \
+    int i,j,si,sj,ni,nj;                   \
+    ni = modei ? in : ip[1]-ip[0]+1;       \
+    nj = modej ? jn : jp[1]-jp[0]+1;       \
+                                           \
+    for (i=0; i<ni; i++) {                 \
+        si = modei ? ip[i] : i+ip[0];      \
+                                           \
+        for (j=0; j<nj; j++) {             \
+            sj = modej ? jp[j] : j+jp[0];  \
+                                           \
+            AT(r,i,j) = AT(m,si,sj);       \
+        }                                  \
+    }                                      \
+    OK
+
+int extractD(int modei, int modej, KIVEC(i), KIVEC(j), KODMAT(m), ODMAT(r)) {
+    EXTRACT_IMP
+}
+
+int extractF(int modei, int modej, KIVEC(i), KIVEC(j), KOFMAT(m), OFMAT(r)) {
+    EXTRACT_IMP
+}
+
+int extractC(int modei, int modej, KIVEC(i), KIVEC(j), KOCMAT(m), OCMAT(r)) {
+    EXTRACT_IMP
+}
+
+int extractQ(int modei, int modej, KIVEC(i), KIVEC(j), KOQMAT(m), OQMAT(r)) {
+    EXTRACT_IMP
+}
+
+int extractI(int modei, int modej, KIVEC(i), KIVEC(j), KOIMAT(m), OIMAT(r)) {
+    EXTRACT_IMP
+}
+
+//////////////////////// remap /////////////////////////////////
+
+#define REMAP_IMP                                               \
+    REQUIRES(ir==jr && ic==jc && ir==rr && ic==rc ,BAD_SIZE);   \
+    { TRAV(r,a,b) { AT(r,a,b) = AT(m,AT(i,a,b),AT(j,a,b)); }    \
+    }                                                           \
+    OK
+
+int remapD(KOIMAT(i), KOIMAT(j), KODMAT(m), ODMAT(r)) {
+    REMAP_IMP
+}
+
+int remapF(KOIMAT(i), KOIMAT(j), KOFMAT(m), OFMAT(r)) {
+    REMAP_IMP
+}
+
+int remapI(KOIMAT(i), KOIMAT(j), KOIMAT(m), OIMAT(r)) {
+    REMAP_IMP
+}
+
+int remapC(KOIMAT(i), KOIMAT(j), KOCMAT(m), OCMAT(r)) {
+    REMAP_IMP
+}
+
+int remapQ(KOIMAT(i), KOIMAT(j), KOQMAT(m), OQMAT(r)) {
+    REMAP_IMP
+}
+
diff --git a/packages/base/src/Internal/C/lapack-aux.h b/packages/base/src/Internal/C/lapack-aux.h
new file mode 100644
index 0000000..6ffbef1
--- /dev/null
+++ b/packages/base/src/Internal/C/lapack-aux.h
@@ -0,0 +1,82 @@
+/*
+ * 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 IVEC(A) int A##n, int*A##p
+#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 IMAT(A) int A##r, int A##c, int* A##p
+#define FMAT(A) int A##r, int A##c, float* A##p
+#define DMAT(A) int A##r, int A##c, double* A##p
+#define QMAT(A) int A##r, int A##c, complex* A##p
+#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 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 KIVEC(A) int A##n, const int*A##p
+#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 KIMAT(A) int A##r, int A##c, const int* A##p
+#define KFMAT(A) int A##r, int A##c, const float* A##p
+#define KDMAT(A) int A##r, int A##c, const double* A##p
+#define KQMAT(A) int A##r, int A##c, const complex* A##p
+#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 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 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/C/vector-aux.c b/packages/base/src/Internal/C/vector-aux.c
new file mode 100644
index 0000000..5662697
--- /dev/null
+++ b/packages/base/src/Internal/C/vector-aux.c
@@ -0,0 +1,1134 @@
+#include <complex.h>
+
+typedef double complex TCD;
+typedef float  complex TCF;
+
+#undef complex
+
+#include "lapack-aux.h"
+
+#define V(x) x##n,x##p
+
+#include <string.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+
+#define MACRO(B) do {B} while (0)
+#define ERROR(CODE) MACRO(return CODE;)
+#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
+#define OK return 0;
+
+#define MIN(A,B) ((A)<(B)?(A):(B))
+#define MAX(A,B) ((A)>(B)?(A):(B))
+
+#ifdef DBG
+#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M);
+#else
+#define DEBUGMSG(M)
+#endif
+
+#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
+
+#define BAD_SIZE 2000
+#define BAD_CODE 2001
+#define MEM      2002
+#define BAD_FILE 2003
+
+
+int sumF(KFVEC(x),FVEC(r)) {
+    DEBUGMSG("sumF");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    float res = 0;
+    for (i = 0; i < xn; i++) res += xp[i];
+    rp[0] = res;
+    OK
+}
+
+int sumR(KDVEC(x),DVEC(r)) {
+    DEBUGMSG("sumR");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    double res = 0;
+    for (i = 0; i < xn; i++) res += xp[i];
+    rp[0] = res;
+    OK
+}
+
+int sumI(KIVEC(x),IVEC(r)) {
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    int res = 0;
+    for (i = 0; i < xn; i++) res += xp[i];
+    rp[0] = res;
+    OK
+}
+
+
+int sumQ(KQVEC(x),QVEC(r)) {
+    DEBUGMSG("sumQ");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    complex res;
+    res.r = 0;
+    res.i = 0;
+    for (i = 0; i < xn; i++) {
+      res.r += xp[i].r;
+      res.i += xp[i].i;
+    }
+    rp[0] = res;
+    OK
+}
+
+int sumC(KCVEC(x),CVEC(r)) {
+    DEBUGMSG("sumC");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    doublecomplex res;
+    res.r = 0;
+    res.i = 0;
+    for (i = 0; i < xn; i++)  {
+      res.r += xp[i].r;
+      res.i += xp[i].i;
+    }
+    rp[0] = res;
+    OK
+}
+
+
+int prodF(KFVEC(x),FVEC(r)) {
+    DEBUGMSG("prodF");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    float res = 1;
+    for (i = 0; i < xn; i++) res *= xp[i];
+    rp[0] = res;
+    OK
+}
+
+int prodR(KDVEC(x),DVEC(r)) {
+    DEBUGMSG("prodR");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    double res = 1;
+    for (i = 0; i < xn; i++) res *= xp[i];
+    rp[0] = res;
+    OK
+}
+
+int prodI(KIVEC(x),IVEC(r)) {
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    int res = 1;
+    for (i = 0; i < xn; i++) res *= xp[i];
+    rp[0] = res;
+    OK
+}
+
+
+
+int prodQ(KQVEC(x),QVEC(r)) {
+    DEBUGMSG("prodQ");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    complex res;
+    float temp;
+    res.r = 1;
+    res.i = 0;
+    for (i = 0; i < xn; i++) {
+      temp  = res.r * xp[i].r - res.i * xp[i].i;
+      res.i = res.r * xp[i].i + res.i * xp[i].r;
+      res.r = temp;
+    }
+    rp[0] = res;
+    OK
+}
+
+int prodC(KCVEC(x),CVEC(r)) {
+    DEBUGMSG("prodC");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    doublecomplex res;
+    double temp;
+    res.r = 1;
+    res.i = 0;
+    for (i = 0; i < xn; i++)  {
+      temp  = res.r * xp[i].r - res.i * xp[i].i;
+      res.i = res.r * xp[i].i + res.i * xp[i].r;
+      res.r = temp;
+    }
+    rp[0] = res;
+    OK
+}
+
+
+double dnrm2_(integer*, const double*, integer*);
+double dasum_(integer*, const double*, integer*);
+
+double vector_max(KDVEC(x)) {
+    double r = xp[0];
+    int k;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]>r) {
+            r = xp[k];
+        }
+    }
+    return r;
+}
+
+double vector_min(KDVEC(x)) {
+    double r = xp[0];
+    int k;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]<r) {
+            r = xp[k];
+        }
+    }
+    return r;
+}
+
+double vector_max_index(KDVEC(x)) {
+    int k, r = 0;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]>xp[r]) {
+            r = k;
+        }
+    }
+    return r;
+}
+
+double vector_min_index(KDVEC(x)) {
+    int k, r = 0;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]<xp[r]) {
+            r = k;
+        }
+    }
+    return r;
+}
+
+int toScalarR(int code, KDVEC(x), DVEC(r)) {
+    REQUIRES(rn==1,BAD_SIZE);
+    DEBUGMSG("toScalarR");
+    double res;
+    integer one = 1;
+    integer n = xn;
+    switch(code) {
+        case 0: { res = dnrm2_(&n,xp,&one); break; }
+        case 1: { res = dasum_(&n,xp,&one);  break; }
+        case 2: { res = vector_max_index(V(x));  break; }
+        case 3: { res = vector_max(V(x));  break; }
+        case 4: { res = vector_min_index(V(x)); break; }
+        case 5: { res = vector_min(V(x)); break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+
+float snrm2_(integer*, const float*, integer*);
+float sasum_(integer*, const float*, integer*);
+
+float vector_max_f(KFVEC(x)) {
+    float r = xp[0];
+    int k;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]>r) {
+            r = xp[k];
+        }
+    }
+    return r;
+}
+
+float vector_min_f(KFVEC(x)) {
+    float r = xp[0];
+    int k;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]<r) {
+            r = xp[k];
+        }
+    }
+    return r;
+}
+
+float vector_max_index_f(KFVEC(x)) {
+    int k, r = 0;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]>xp[r]) {
+            r = k;
+        }
+    }
+    return r;
+}
+
+float vector_min_index_f(KFVEC(x)) {
+    int k, r = 0;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]<xp[r]) {
+            r = k;
+        }
+    }
+    return r;
+}
+
+
+int toScalarF(int code, KFVEC(x), FVEC(r)) {
+    REQUIRES(rn==1,BAD_SIZE);
+    DEBUGMSG("toScalarF");
+    float res;
+    integer one = 1;
+    integer n = xn;
+    switch(code) {
+        case 0: { res = snrm2_(&n,xp,&one); break; }
+        case 1: { res = sasum_(&n,xp,&one);  break; }
+        case 2: { res = vector_max_index_f(V(x));  break; }
+        case 3: { res = vector_max_f(V(x));  break; }
+        case 4: { res = vector_min_index_f(V(x)); break; }
+        case 5: { res = vector_min_f(V(x)); break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+int vector_max_i(KIVEC(x)) {
+    int r = xp[0];
+    int k;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]>r) {
+            r = xp[k];
+        }
+    }
+    return r;
+}
+
+int vector_min_i(KIVEC(x)) {
+    float r = xp[0];
+    int k;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]<r) {
+            r = xp[k];
+        }
+    }
+    return r;
+}
+
+int vector_max_index_i(KIVEC(x)) {
+    int k, r = 0;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]>xp[r]) {
+            r = k;
+        }
+    }
+    return r;
+}
+
+int vector_min_index_i(KIVEC(x)) {
+    int k, r = 0;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]<xp[r]) {
+            r = k;
+        }
+    }
+    return r;
+}
+
+
+int toScalarI(int code, KIVEC(x), IVEC(r)) {
+    REQUIRES(rn==1,BAD_SIZE);
+    int res;
+    switch(code) {
+        case 2: { res = vector_max_index_i(V(x));  break; }
+        case 3: { res = vector_max_i(V(x));  break; }
+        case 4: { res = vector_min_index_i(V(x)); break; }
+        case 5: { res = vector_min_i(V(x)); break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+
+double dznrm2_(integer*, const doublecomplex*, integer*);
+double dzasum_(integer*, const doublecomplex*, integer*);
+
+int toScalarC(int code, KCVEC(x), DVEC(r)) {
+    REQUIRES(rn==1,BAD_SIZE);
+    DEBUGMSG("toScalarC");
+    double res;
+    integer one = 1;
+    integer n = xn;
+    switch(code) {
+        case 0: { res = dznrm2_(&n,xp,&one); break; }
+        case 1: { res = dzasum_(&n,xp,&one);  break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+
+double scnrm2_(integer*, const complex*, integer*);
+double scasum_(integer*, const complex*, integer*);
+
+int toScalarQ(int code, KQVEC(x), FVEC(r)) {
+    REQUIRES(rn==1,BAD_SIZE);
+    DEBUGMSG("toScalarQ");
+    float res;
+    integer one = 1;
+    integer n = xn;
+    switch(code) {
+        case 0: { res = scnrm2_(&n,xp,&one); break; }
+        case 1: { res = scasum_(&n,xp,&one);  break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+
+inline double sign(double x) {
+    if(x>0) {
+        return +1.0;
+    } else if (x<0) {
+        return -1.0;
+    } else {
+        return 0.0;
+    }
+}
+
+inline float float_sign(float x) {
+    if(x>0) {
+        return +1.0;
+    } else if (x<0) {
+        return -1.0;
+    } else {
+        return 0.0;
+    }
+}
+
+
+#define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK }
+#define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK }
+int mapR(int code, KDVEC(x), DVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapR");
+    switch (code) {
+        OP(0,sin)
+        OP(1,cos)
+        OP(2,tan)
+        OP(3,fabs)
+        OP(4,asin)
+        OP(5,acos)
+        OP(6,atan)
+        OP(7,sinh)
+        OP(8,cosh)
+        OP(9,tanh)
+        OP(10,asinh)
+        OP(11,acosh)
+        OP(12,atanh)
+        OP(13,exp)
+        OP(14,log)
+        OP(15,sign)
+        OP(16,sqrt)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapF(int code, KFVEC(x), FVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapF");
+    switch (code) {
+        OP(0,sin)
+        OP(1,cos)
+        OP(2,tan)
+        OP(3,fabs)
+        OP(4,asin)
+        OP(5,acos)
+        OP(6,atan)
+        OP(7,sinh)
+        OP(8,cosh)
+        OP(9,tanh)
+        OP(10,asinh)
+        OP(11,acosh)
+        OP(12,atanh)
+        OP(13,exp)
+        OP(14,log)
+        OP(15,sign)
+        OP(16,sqrt)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int mapI(int code, KIVEC(x), IVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    switch (code) {
+        OP(3,abs)
+        OP(15,sign)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+
+inline double abs_complex(doublecomplex z) {
+    return sqrt(z.r*z.r + z.i*z.i);
+}
+
+inline doublecomplex complex_abs_complex(doublecomplex z) {
+    doublecomplex r;
+    r.r = abs_complex(z);
+    r.i = 0;
+    return r;
+}
+
+inline doublecomplex complex_signum_complex(doublecomplex z) {
+    doublecomplex r;
+    double mag;
+    if (z.r == 0 && z.i == 0) {
+        r.r = 0;
+        r.i = 0;
+    } else {
+        mag = abs_complex(z);
+        r.r = z.r/mag;
+        r.i = z.i/mag;
+    }
+    return r;
+}
+
+#define OPb(C,F) case C: { for(k=0;k<xn;k++) r2p[k] = F(x2p[k]); OK }
+int mapC(int code, KCVEC(x), CVEC(r)) {
+    TCD* x2p = (TCD*)xp;
+    TCD* r2p = (TCD*)rp;
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapC");
+    switch (code) {
+        OPb(0,csin)
+        OPb(1,ccos)
+        OPb(2,ctan)
+        OP(3,complex_abs_complex)
+        OPb(4,casin)
+        OPb(5,cacos)
+        OPb(6,catan)
+        OPb(7,csinh)
+        OPb(8,ccosh)
+        OPb(9,ctanh)
+        OPb(10,casinh)
+        OPb(11,cacosh)
+        OPb(12,catanh)
+        OPb(13,cexp)
+        OPb(14,clog)
+        OP(15,complex_signum_complex)
+        OPb(16,csqrt)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+
+inline complex complex_f_math_fun(doublecomplex (*cf)(doublecomplex), complex a)
+{
+  doublecomplex c;
+  doublecomplex r;
+
+  complex float_r;
+
+  c.r = a.r;
+  c.i = a.i;
+
+  r = (*cf)(c);
+
+  float_r.r = r.r;
+  float_r.i = r.i;
+
+  return float_r;
+}
+
+
+#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_f_math_fun(&F,xp[k]); OK }
+int mapQ(int code, KQVEC(x), QVEC(r)) {
+    TCF* x2p = (TCF*)xp;
+    TCF* r2p = (TCF*)rp;
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapQ");
+    switch (code) {
+        OPb(0,csinf)
+        OPb(1,ccosf)
+        OPb(2,ctanf)
+        OPC(3,complex_abs_complex)
+        OPb(4,casinf)
+        OPb(5,cacosf)
+        OPb(6,catanf)
+        OPb(7,csinhf)
+        OPb(8,ccoshf)
+        OPb(9,ctanhf)
+        OPb(10,casinhf)
+        OPb(11,cacoshf)
+        OPb(12,catanhf)
+        OPb(13,cexpf)
+        OPb(14,clogf)
+        OPC(15,complex_signum_complex)
+        OPb(16,csqrtf)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int mapValR(int code, double* pval, KDVEC(x), DVEC(r)) {
+    int k;
+    double val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValR");
+    switch (code) {
+        OPV(0,val*xp[k])
+        OPV(1,val/xp[k])
+        OPV(2,val+xp[k])
+        OPV(3,val-xp[k])
+        OPV(4,pow(val,xp[k]))
+        OPV(5,pow(xp[k],val))
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) {
+    int k;
+    float val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValF");
+    switch (code) {
+        OPV(0,val*xp[k])
+        OPV(1,val/xp[k])
+        OPV(2,val+xp[k])
+        OPV(3,val-xp[k])
+        OPV(4,pow(val,xp[k]))
+        OPV(5,pow(xp[k],val))
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mod (int a, int b) {
+    int m = a % b;
+    if (b>0) {
+        return m >=0 ? m : m+b;
+    } else {
+        return m <=0 ? m : m+b;
+    }
+}
+
+int mapValI(int code, int* pval, KIVEC(x), IVEC(r)) {
+    int k;
+    int val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValI");
+    switch (code) {
+        OPV(0,val*xp[k])
+        OPV(1,val/xp[k])
+        OPV(2,val+xp[k])
+        OPV(3,val-xp[k])
+        OPV(6,mod(val,xp[k]))
+        OPV(7,mod(xp[k],val))
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+
+inline doublecomplex complex_add(doublecomplex a, doublecomplex b) {
+    doublecomplex r;
+    r.r = a.r+b.r;
+    r.i = a.i+b.i;
+    return r;
+}
+
+#define OPVb(C,E) case C: { for(k=0;k<xn;k++) r2p[k] = E; OK }
+int mapValC(int code, doublecomplex* pval, KCVEC(x), CVEC(r)) {
+    TCD* x2p = (TCD*)xp;
+    TCD* r2p = (TCD*)rp;
+    int k;
+    TCD val = * (TCD*)pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValC");
+    switch (code) {
+        OPVb(0,val*x2p[k])
+        OPVb(1,val/x2p[k])
+        OPVb(2,val+x2p[k])
+        OPVb(3,val-x2p[k])
+        OPVb(4,cpow(val,x2p[k]))
+        OPVb(5,cpow(x2p[k],val))
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int mapValQ(int code, complex* pval, KQVEC(x), QVEC(r)) {
+    TCF* x2p = (TCF*)xp;
+    TCF* r2p = (TCF*)rp;
+    int k;
+    TCF val = *(TCF*)pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValQ");
+    switch (code) {
+        OPVb(0,val*x2p[k])
+        OPVb(1,val/x2p[k])
+        OPVb(2,val+x2p[k])
+        OPVb(3,val-x2p[k])
+        OPVb(4,cpow(val,x2p[k]))
+        OPVb(5,cpow(x2p[k],val))
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+
+#define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK }
+#define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK }
+#define OPZO(C,msg,O) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = ap[k] O bp[k]; OK }
+
+int zipR(int code, KDVEC(a), KDVEC(b), DVEC(r)) {
+REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPZO(0,"zipR Add",+)
+        OPZO(1,"zipR Sub",-)
+        OPZO(2,"zipR Mul",*)
+        OPZO(3,"zipR Div",/)
+        OPZE(4,"zipR Pow",  pow)
+        OPZE(5,"zipR ATan2",atan2)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int zipF(int code, KFVEC(a), KFVEC(b), FVEC(r)) {
+REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPZO(0,"zipR Add",+)
+        OPZO(1,"zipR Sub",-)
+        OPZO(2,"zipR Mul",*)
+        OPZO(3,"zipR Div",/)
+        OPZE(4,"zipR Pow",  pow)
+        OPZE(5,"zipR ATan2",atan2)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int zipI(int code, KIVEC(a), KIVEC(b), IVEC(r)) {
+REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPZO(0,"zipI Add",+)
+        OPZO(1,"zipI Sub",-)
+        OPZO(2,"zipI Mul",*)
+        OPZO(3,"zipI Div",/)
+        OPZO(6,"zipI Mod",%)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+
+#define OPZOb(C,msg,O) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) r2p[k] = a2p[k] O b2p[k]; OK }
+#define OPZEb(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) r2p[k] = E(a2p[k],b2p[k]); OK }
+int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
+    REQUIRES(an == bn && an == rn, BAD_SIZE);
+    TCD* a2p = (TCD*)ap;
+    TCD* b2p = (TCD*)bp;
+    TCD* r2p = (TCD*)rp;
+    int k;
+    switch(code) {
+        OPZOb(0,"zipC Add",+)
+        OPZOb(1,"zipC Sub",-)
+        OPZOb(2,"zipC Mul",*)
+        OPZOb(3,"zipC Div",/)
+        OPZEb(4,"zipC Pow",cpow)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+
+
+
+#define OPCZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = complex_f_math_op(&E,ap[k],bp[k]); OK }
+
+int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) {
+    REQUIRES(an == bn && an == rn, BAD_SIZE);
+    TCF* a2p = (TCF*)ap;
+    TCF* b2p = (TCF*)bp;
+    TCF* r2p = (TCF*)rp;
+
+    int k;
+    switch(code) {
+        OPZOb(0,"zipC Add",+)
+        OPZOb(1,"zipC Sub",-)
+        OPZOb(2,"zipC Mul",*)
+        OPZOb(3,"zipC Div",/)
+        OPZEb(4,"zipC Pow",cpowf)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+int vectorScan(char * file, int* n, double**pp){
+    FILE * fp;
+    fp = fopen (file, "r");
+    if(!fp) {
+        ERROR(BAD_FILE);
+    }
+    int nbuf = 100*100;
+    double * p = (double*)malloc(nbuf*sizeof(double));
+    int k=0;
+    double d;
+    int ok;
+    for (;;) {
+        ok = fscanf(fp,"%lf",&d);
+        if (ok<1) {
+            break;
+        }
+        if (k==nbuf) {
+            nbuf = nbuf * 2;
+            p = (double*)realloc(p,nbuf*sizeof(double));
+            // printf("R\n");
+        }
+        p[k++] = d;
+    }
+    *n = k;
+    *pp = p;
+    fclose(fp);
+    OK
+}
+
+int saveMatrix(char * file, char * format, KDMAT(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]);
+            if (c<ac-1) {
+                fprintf(fp," ");
+            } else {
+                fprintf(fp,"\n");
+            }
+        }
+    }
+    fclose(fp);
+    OK
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+#if defined (__APPLE__) || (__FreeBSD__)
+/* FreeBSD and Mac OS X do not provide random_r(), thread safety cannot be
+   guaranteed.
+   For FreeBSD and Mac OS X, nrand48() is much better than random().
+   See: http://www.evanjones.ca/random-thread-safe.html
+*/
+#pragma message "randomVector is not thread-safe in OSX and FreeBSD"
+
+inline double urandom() {
+    /* the probalility of matching will be theoretically p^3(in fact, it is not)
+       p is matching probalility of random().
+       using the test there, only 3 matches, using random(), 13783 matches
+    */
+    unsigned short state[3];
+    state[0] = random();
+    state[1] = random();
+    state[2] = random();
+
+    const long max_random = 2147483647; // 2**31 - 1
+    return (double)nrand48(state) / (double)max_random;
+}
+
+double gaussrand(int *phase, double *pV1, double *pV2, double *pS)
+{
+	double V1=*pV1, V2=*pV2, S=*pS;
+	double X;
+
+	if(*phase == 0) {
+		do {
+            double U1 = urandom();
+			double U2 = urandom();
+
+			V1 = 2 * U1 - 1;
+			V2 = 2 * U2 - 1;
+			S = V1 * V1 + V2 * V2;
+			} while(S >= 1 || S == 0);
+
+		X = V1 * sqrt(-2 * log(S) / S);
+	} else
+		X = V2 * sqrt(-2 * log(S) / S);
+
+	*phase = 1 - *phase;
+    *pV1=V1; *pV2=V2; *pS=S;
+
+	return X;
+
+}
+
+int random_vector(unsigned int seed, int code, DVEC(r)) {
+    int phase = 0;
+    double V1,V2,S;
+
+    srandom(seed);
+
+    int k;
+    switch (code) {
+      case 0: { // uniform
+        for (k=0; k<rn; k++) {
+            rp[k] = urandom();
+        }
+        OK
+      }
+      case 1: { // gaussian
+        for (k=0; k<rn; k++) {
+            rp[k] = gaussrand(&phase,&V1,&V2,&S);
+        }
+        OK
+      }
+
+      default: ERROR(BAD_CODE);
+    }
+}
+
+#else
+
+inline double urandom(struct random_data * buffer) {
+    int32_t res;
+    random_r(buffer,&res);
+    return (double)res/RAND_MAX;
+}
+
+
+// http://c-faq.com/lib/gaussian.html
+double gaussrand(struct random_data *buffer,
+                 int *phase, double *pV1, double *pV2, double *pS)
+{
+	double V1=*pV1, V2=*pV2, S=*pS;
+	double X;
+
+	if(*phase == 0) {
+		do {
+            double U1 = urandom(buffer);
+			double U2 = urandom(buffer);
+
+			V1 = 2 * U1 - 1;
+			V2 = 2 * U2 - 1;
+			S = V1 * V1 + V2 * V2;
+			} while(S >= 1 || S == 0);
+
+		X = V1 * sqrt(-2 * log(S) / S);
+	} else
+		X = V2 * sqrt(-2 * log(S) / S);
+
+	*phase = 1 - *phase;
+    *pV1=V1; *pV2=V2; *pS=S;
+
+	return X;
+
+}
+
+int random_vector(unsigned int seed, int code, DVEC(r)) {
+    struct random_data buffer;
+    char   random_state[128];
+    memset(&buffer, 0, sizeof(struct random_data));
+    memset(random_state, 0, sizeof(random_state));
+
+    initstate_r(seed,random_state,sizeof(random_state),&buffer);
+    // setstate_r(random_state,&buffer);
+    // srandom_r(seed,&buffer);
+
+    int phase = 0;
+    double V1,V2,S;
+
+    int k;
+    switch (code) {
+      case 0: { // uniform
+        for (k=0; k<rn; k++) {
+            rp[k] = urandom(&buffer);
+        }
+        OK
+      }
+      case 1: { // gaussian
+        for (k=0; k<rn; k++) {
+            rp[k] = gaussrand(&buffer,&phase,&V1,&V2,&S);
+        }
+        OK
+      }
+
+      default: ERROR(BAD_CODE);
+    }
+}
+
+#endif
+
+////////////////////////////////////////////////////////////////////////////////
+
+int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
+    int r, c;
+    for (r = 0; r < rowsn - 1; r++) {
+        rp[r] = 0;
+        for (c = rowsp[r]; c < rowsp[r+1]; c++) {
+            rp[r] += valsp[c-1] * xp[colsp[c-1]-1];
+        }
+    }
+    OK
+}
+
+int smTXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
+    int r,c;
+    for (c = 0; c < rn; c++) {
+        rp[c] = 0;
+    }
+    for (r = 0; r < rowsn - 1; r++) {
+        for (c = rowsp[r]; c < rowsp[r+1]; c++) {
+            rp[colsp[c-1]-1] += valsp[c-1] * xp[r];
+        }
+    }
+    OK
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+int
+compare_doubles (const void *a, const void *b) {
+  return *(double*)a > *(double*)b;
+}
+
+int sort_valuesD(KDVEC(v),DVEC(r)) {
+    memcpy(rp,vp,vn*sizeof(double));
+    qsort(rp,rn,sizeof(double),compare_doubles);
+    OK
+}
+
+int
+compare_floats (const void *a, const void *b) {
+  return *(float*)a > *(float*)b;
+}
+
+int sort_valuesF(KFVEC(v),FVEC(r)) {
+    memcpy(rp,vp,vn*sizeof(float));
+    qsort(rp,rn,sizeof(float),compare_floats);
+    OK
+}
+
+int
+compare_ints(const void *a, const void *b) {
+  return *(int*)a > *(int*)b;
+}
+
+int sort_valuesI(KIVEC(v),IVEC(r)) {
+    memcpy(rp,vp,vn*sizeof(int));
+    qsort(rp,rn,sizeof(int),compare_ints);
+    OK
+}
+
+////////////////////////////////////////
+
+
+#define SORTIDX_IMP(T,C)                   \
+    T* x = (T*)malloc(sizeof(T)*vn);       \
+    int k;                                 \
+    for (k=0;k<vn;k++) {                   \
+        x[k].pos = k;                      \
+        x[k].val = vp[k];                  \
+    }                                      \
+                                           \
+    qsort(x,vn,sizeof(T),C);               \
+                                           \
+    for (k=0;k<vn;k++) {                   \
+        rp[k] = x[k].pos;                  \
+    }                                      \
+    free(x);                               \
+    OK
+
+
+typedef struct SDI { int pos; double val;} DI;
+
+int compare_doubles_i (const void *a, const void *b) {
+  return ((DI*)a)->val > ((DI*)b)->val;
+}
+
+int sort_indexD(KDVEC(v),IVEC(r)) {
+    SORTIDX_IMP(DI,compare_doubles_i)
+}
+
+
+typedef struct FI { int pos; float  val;} FI;
+
+int compare_floats_i (const void *a, const void *b) {
+  return ((FI*)a)->val > ((FI*)b)->val;
+}
+
+int sort_indexF(KFVEC(v),IVEC(r)) {
+    SORTIDX_IMP(FI,compare_floats_i)
+}
+
+
+typedef struct II { int pos; int    val;} II;
+
+int compare_ints_i (const void *a, const void *b) {
+  return ((II*)a)->val > ((II*)b)->val;
+}
+
+int sort_indexI(KIVEC(v),IVEC(r)) {
+    SORTIDX_IMP(II,compare_ints_i)
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+
+int round_vector(KDVEC(v),DVEC(r)) {
+    int k;
+    for(k=0; k<vn; k++) {
+        rp[k] = round(vp[k]);
+    }
+    OK
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+int round_vector_i(KDVEC(v),IVEC(r)) {
+    int k;
+    for(k=0; k<vn; k++) {
+        rp[k] = round(vp[k]);
+    }
+    OK
+}
+
+
+int mod_vector(int m, KIVEC(v), IVEC(r)) {
+    int k;
+    for(k=0; k<vn; k++) {
+        rp[k] = vp[k] % m;
+    }
+    OK
+}
+
+int div_vector(int m, KIVEC(v), IVEC(r)) {
+    int k;
+    for(k=0; k<vn; k++) {
+        rp[k] = vp[k] / m;
+    }
+    OK
+}
+
+int range_vector(IVEC(r)) {
+    int k;
+    for(k=0; k<rn; k++) {
+        rp[k] = k;
+    }
+    OK
+}
+
-- 
cgit v1.2.3