From 5992d92357cfd911c8f2e9f5faaa4fd8a323fd9a Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 8 May 2014 12:16:42 +0200 Subject: Data.Packed -> base (I) --- packages/base/src/C/lapack-aux.c | 1489 ++++++++++++++++++++++++++++++++++++++ packages/base/src/C/lapack-aux.h | 60 ++ 2 files changed, 1549 insertions(+) create mode 100644 packages/base/src/C/lapack-aux.c create mode 100644 packages/base/src/C/lapack-aux.h (limited to 'packages/base/src/C') diff --git a/packages/base/src/C/lapack-aux.c b/packages/base/src/C/lapack-aux.c new file mode 100644 index 0000000..e5e45ef --- /dev/null +++ b/packages/base/src/C/lapack-aux.c @@ -0,0 +1,1489 @@ +#include +#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 stepD(DVEC(x),DVEC(y)) { + DEBUGMSG("stepD") + int k; + for(k=0;k0; + } + OK +} + +//////////////////// cond ///////////////////////// + +int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) { + REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE); + DEBUGMSG("condF") + int k; + for(k=0;kyp[k]?gtp[k]:eqp[k]); + } + OK +} + +int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) { + REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE); + DEBUGMSG("condD") + int k; + for(k=0;kyp[k]?gtp[k]:eqp[k]); + } + OK +} + diff --git a/packages/base/src/C/lapack-aux.h b/packages/base/src/C/lapack-aux.h new file mode 100644 index 0000000..a3f1899 --- /dev/null +++ b/packages/base/src/C/lapack-aux.h @@ -0,0 +1,60 @@ +/* + * We have copied the definitions in f2c.h required + * to compile clapack.h, modified to support both + * 32 and 64 bit + + http://opengrok.creo.hu/dragonfly/xref/src/contrib/gcc-3.4/libf2c/readme.netlib + http://www.ibm.com/developerworks/library/l-port64.html + */ + +#ifdef _LP64 +typedef int integer; +typedef unsigned int uinteger; +typedef int logical; +typedef long longint; /* system-dependent */ +typedef unsigned long ulongint; /* system-dependent */ +#else +typedef long int integer; +typedef unsigned long int uinteger; +typedef long int logical; +typedef long long longint; /* system-dependent */ +typedef unsigned long long ulongint; /* system-dependent */ +#endif + +typedef char *address; +typedef short int shortint; +typedef float real; +typedef double doublereal; +typedef struct { real r, i; } complex; +typedef struct { doublereal r, i; } doublecomplex; +typedef short int shortlogical; +typedef char logical1; +typedef char integer1; + +typedef logical (*L_fp)(); +typedef short ftnlen; + +/********************************************************/ + +#define FVEC(A) int A##n, float*A##p +#define DVEC(A) int A##n, double*A##p +#define QVEC(A) int A##n, complex*A##p +#define CVEC(A) int A##n, doublecomplex*A##p +#define PVEC(A) int A##n, void* A##p, int A##s +#define FMAT(A) int A##r, int A##c, float* A##p +#define DMAT(A) int A##r, int A##c, double* A##p +#define QMAT(A) int A##r, int A##c, complex* A##p +#define CMAT(A) int A##r, int A##c, doublecomplex* A##p +#define PMAT(A) int A##r, int A##c, void* A##p, int A##s + +#define KFVEC(A) int A##n, const float*A##p +#define KDVEC(A) int A##n, const double*A##p +#define KQVEC(A) int A##n, const complex*A##p +#define KCVEC(A) int A##n, const doublecomplex*A##p +#define KPVEC(A) int A##n, const void* A##p, int A##s +#define KFMAT(A) int A##r, int A##c, const float* A##p +#define KDMAT(A) int A##r, int A##c, const double* A##p +#define KQMAT(A) int A##r, int A##c, const complex* A##p +#define KCMAT(A) int A##r, int A##c, const doublecomplex* A##p +#define KPMAT(A) int A##r, int A##c, const void* A##p, int A##s + -- cgit v1.2.3