From b7e57952112eae61054fc9171ddb311fbaca0841 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 4 Jun 2015 12:39:33 +0200 Subject: move c sources --- packages/base/CHANGELOG | 2 + packages/base/hmatrix.cabal | 6 +- packages/base/src/C/lapack-aux.c | 1686 ----------------------------- packages/base/src/C/lapack-aux.h | 82 -- packages/base/src/C/vector-aux.c | 1134 ------------------- 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 +++++++++++++++++++ 8 files changed, 2907 insertions(+), 2905 deletions(-) delete mode 100644 packages/base/src/C/lapack-aux.c delete mode 100644 packages/base/src/C/lapack-aux.h delete mode 100644 packages/base/src/C/vector-aux.c 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') diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG index 5e00f93..b1d85fb 100644 --- a/packages/base/CHANGELOG +++ b/packages/base/CHANGELOG @@ -4,6 +4,8 @@ * initial support of (C)Int elements * improved matrix extraction using vectors of indexes (??) + + * remap, ccompare * experimental support of type safe modular arithmetic diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal index 6614633..6cbe8c7 100644 --- a/packages/base/hmatrix.cabal +++ b/packages/base/hmatrix.cabal @@ -20,7 +20,7 @@ build-type: Simple extra-source-files: THANKS.md CHANGELOG -extra-source-files: src/C/lapack-aux.h +extra-source-files: src/Internal/C/lapack-aux.h flag openblas description: Link with OpenBLAS (https://github.com/xianyi/OpenBLAS) optimized libraries. @@ -77,8 +77,8 @@ library Numeric.LinearAlgebra.Util Numeric.LinearAlgebra.Util.Modular - C-sources: src/C/lapack-aux.c - src/C/vector-aux.c + C-sources: src/Internal/C/lapack-aux.c + src/Internal/C/vector-aux.c extensions: ForeignFunctionInterface, diff --git a/packages/base/src/C/lapack-aux.c b/packages/base/src/C/lapack-aux.c deleted file mode 100644 index 1402050..0000000 --- a/packages/base/src/C/lapack-aux.c +++ /dev/null @@ -1,1686 +0,0 @@ -#include -#include -#include -#include -#include -#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=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; k0) { - 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; k0) { - 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; k0) { - 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; k0) { - 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=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=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;k0) { - 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=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; k0; \ - } \ - 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;kyp[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;kyp[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;k0?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 - -typedef double complex TCD; -typedef float complex TCF; - -#undef complex - -#include "lapack-aux.h" - -#define V(x) x##n,x##p - -#include -#include -#include -#include -#include - -#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; kr) { - r = xp[k]; - } - } - return r; -} - -double vector_min(KDVEC(x)) { - double r = xp[0]; - int k; - for (k = 1; kxp[r]) { - r = k; - } - } - return r; -} - -double vector_min_index(KDVEC(x)) { - int k, r = 0; - for (k = 1; kr) { - r = xp[k]; - } - } - return r; -} - -float vector_min_f(KFVEC(x)) { - float r = xp[0]; - int k; - for (k = 1; kxp[r]) { - r = k; - } - } - return r; -} - -float vector_min_index_f(KFVEC(x)) { - int k, r = 0; - for (k = 1; kr) { - r = xp[k]; - } - } - return r; -} - -int vector_min_i(KIVEC(x)) { - float r = xp[0]; - int k; - for (k = 1; kxp[r]) { - r = k; - } - } - return r; -} - -int vector_min_index_i(KIVEC(x)) { - int k, r = 0; - for (k = 1; k0) { - 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;k0) { - 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= 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= 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 *(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;kval > ((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 +#include +#include +#include +#include +#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=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; k0) { + 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; k0) { + 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; k0) { + 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; k0) { + 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=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=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;k0) { + 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=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; k0; \ + } \ + 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;kyp[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;kyp[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;k0?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 + +typedef double complex TCD; +typedef float complex TCF; + +#undef complex + +#include "lapack-aux.h" + +#define V(x) x##n,x##p + +#include +#include +#include +#include +#include + +#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; kr) { + r = xp[k]; + } + } + return r; +} + +double vector_min(KDVEC(x)) { + double r = xp[0]; + int k; + for (k = 1; kxp[r]) { + r = k; + } + } + return r; +} + +double vector_min_index(KDVEC(x)) { + int k, r = 0; + for (k = 1; kr) { + r = xp[k]; + } + } + return r; +} + +float vector_min_f(KFVEC(x)) { + float r = xp[0]; + int k; + for (k = 1; kxp[r]) { + r = k; + } + } + return r; +} + +float vector_min_index_f(KFVEC(x)) { + int k, r = 0; + for (k = 1; kr) { + r = xp[k]; + } + } + return r; +} + +int vector_min_i(KIVEC(x)) { + float r = xp[0]; + int k; + for (k = 1; kxp[r]) { + r = k; + } + } + return r; +} + +int vector_min_index_i(KIVEC(x)) { + int k, r = 0; + for (k = 1; k0) { + 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;k0) { + 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= 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= 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 *(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;kval > ((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