From bbaadc1070c13f97f1a79befdc342a73cb4b2d5b Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Mon, 27 Apr 2015 11:21:55 +0200 Subject: changelogs --- packages/base/CHANGELOG | 7 +++++++ 1 file changed, 7 insertions(+) (limited to 'packages/base/CHANGELOG') diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG index c137285..ce19c97 100644 --- a/packages/base/CHANGELOG +++ b/packages/base/CHANGELOG @@ -1,3 +1,10 @@ +0.17.0.0 +-------- + + * old compatibility modules removed + + * added "unitary" and "pairwiseD2" + 0.16.1.0 -------- -- cgit v1.2.3 From 79fd8424425087ac411f3b916582ceb4ac198142 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 22 May 2015 19:54:45 +0200 Subject: changelog --- packages/base/CHANGELOG | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'packages/base/CHANGELOG') diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG index ce19c97..ccac516 100644 --- a/packages/base/CHANGELOG +++ b/packages/base/CHANGELOG @@ -1,6 +1,10 @@ 0.17.0.0 -------- + * initial support of (C)Int elements + + * improved matrix extraction using vectors of indexes + * old compatibility modules removed * added "unitary" and "pairwiseD2" -- cgit v1.2.3 From 27334e233edd3d891d2837330289a26e54d05fd1 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 3 Jun 2015 18:45:18 +0200 Subject: thanks, changelog --- packages/base/CHANGELOG | 8 +++++--- packages/base/THANKS.md | 2 ++ packages/gsl/THANKS.md | 4 +++- 3 files changed, 10 insertions(+), 4 deletions(-) (limited to 'packages/base/CHANGELOG') diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG index ccac516..5e00f93 100644 --- a/packages/base/CHANGELOG +++ b/packages/base/CHANGELOG @@ -2,12 +2,14 @@ -------- * initial support of (C)Int elements - - * improved matrix extraction using vectors of indexes + + * improved matrix extraction using vectors of indexes (??) + + * experimental support of type safe modular arithmetic * old compatibility modules removed - * added "unitary" and "pairwiseD2" + * unitary, pairwiseD2, tr' 0.16.1.0 -------- diff --git a/packages/base/THANKS.md b/packages/base/THANKS.md index 133bd31..6600329 100644 --- a/packages/base/THANKS.md +++ b/packages/base/THANKS.md @@ -192,3 +192,5 @@ module reorganization, monadic mapVectorM, and many other improvements. - "maxc01" solved uninstallability in FreeBSD and improved urandom +- "ntfrgl" added {take,drop}Last{Rows,Columns} + diff --git a/packages/gsl/THANKS.md b/packages/gsl/THANKS.md index b066185..7c9727e 100644 --- a/packages/gsl/THANKS.md +++ b/packages/gsl/THANKS.md @@ -1,4 +1,6 @@ -Matt Peddie wrote the interpolation interface. +- Matt Peddie wrote the interpolation interface. + +- "ntfrgl" added odeSolveVWith with generalized step control function (see also the THANKS file of the hmatrix package) -- cgit v1.2.3 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/base/CHANGELOG') 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 Date: Tue, 9 Jun 2015 13:46:25 +0200 Subject: fix gaussElim, changelog --- packages/base/CHANGELOG | 6 +++--- packages/base/src/Internal/Modular.hs | 6 +++--- packages/base/src/Internal/Util.hs | 31 ++++++++++++++++++------------- 3 files changed, 24 insertions(+), 19 deletions(-) (limited to 'packages/base/CHANGELOG') diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG index b1d85fb..27c4d31 100644 --- a/packages/base/CHANGELOG +++ b/packages/base/CHANGELOG @@ -1,11 +1,11 @@ 0.17.0.0 -------- - * initial support of (C)Int elements - * improved matrix extraction using vectors of indexes (??) + + * basic support of Int32 and Int64 elements - * remap, ccompare + * remap, ccompare, sortIndex * experimental support of type safe modular arithmetic diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs index 6b34010..1289a21 100644 --- a/packages/base/src/Internal/Modular.hs +++ b/packages/base/src/Internal/Modular.hs @@ -137,7 +137,7 @@ instance forall m . KnownNat m => Container Vector (Mod m I) scalar' x = fromList [x] konst' x = i2f . konst (unMod x) build' n f = build n (fromIntegral . f) - cmap' = cmap + cmap' = mapVector atIndex' x k = fromIntegral (atIndex (f2i x) k) minIndex' = minIndex . f2i maxIndex' = maxIndex . f2i @@ -177,7 +177,7 @@ instance forall m . KnownNat m => Container Vector (Mod m Z) scalar' x = fromList [x] konst' x = i2f . konst (unMod x) build' n f = build n (fromIntegral . f) - cmap' = cmap + cmap' = mapVector atIndex' x k = fromIntegral (atIndex (f2i x) k) minIndex' = minIndex . f2i maxIndex' = maxIndex . f2i @@ -311,7 +311,7 @@ test = (ok, info) print $ am <> gaussElim am bm - bm print $ ad <> gaussElim ad bd - bd - + print g print $ g <> g print gm diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs index 1b639e1..b1fb800 100644 --- a/packages/base/src/Internal/Util.hs +++ b/packages/base/src/Internal/Util.hs @@ -512,7 +512,7 @@ block3x3 r nr c nc m = [[m ?? (er !! i, ec !! j) | j <- [0..2] ] | i <- [0..2] ] er = [ Range 0 1 (r-1), Range r 1 (r+nr-1), Drop (nr+r) ] ec = [ Range 0 1 (c-1), Range c 1 (c+nc-1), Drop (nc+c) ] -view1 :: Numeric t => Matrix t -> Maybe (t, Vector t, Vector t, Matrix t) +view1 :: Numeric t => Matrix t -> Maybe (View1 t) view1 m | rows m > 0 && cols m > 0 = Just (e, flatten m12, flatten m21 , m22) | otherwise = Nothing @@ -520,21 +520,25 @@ view1 m [[m11,m12],[m21,m22]] = block2x2 1 1 m e = m11 `atIndex` (0, 0) -unView1 :: Numeric t => (t, Vector t, Vector t, Matrix t) -> Matrix t +unView1 :: Numeric t => View1 t -> Matrix t unView1 (e,r,c,m) = fromBlocks [[scalar e, asRow r],[asColumn c, m]] +type View1 t = (t, Vector t, Vector t, Matrix t) +foldMatrix :: Numeric t => (Matrix t -> Matrix t) -> (View1 t -> View1 t) -> (Matrix t -> Matrix t) foldMatrix g f ( (f <$>) . view1 . g -> Just (e,r,c,m)) = unView1 (e, r, c, foldMatrix g f m) foldMatrix _ _ m = m -sortRowsBy h j m = m ?? (Pos (sortIndex (h (tr m ! j))), All) - -splitColsAt n = (takeColumns n &&& dropColumns n) +swapMax k m + | rows m > 0 && j>0 = (j, m ?? (Pos (idxs swapped), All)) + | otherwise = (0,m) + where + j = maxIndex $ abs (tr m ! k) + swapped = j:[1..j-1] ++ 0:[j+1..rows m-1] -down a = foldMatrix g f a +down g a = foldMatrix g f a where - g = sortRowsBy (negate.abs) 0 f (e,r,c,m) | e /= 0 = (1, r', 0, m - outer c r') | otherwise = error "singular!" @@ -547,15 +551,16 @@ down a = foldMatrix g f a -- @a <> gaussElim a b = b@ -- gaussElim - :: (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t) + :: (Eq t, Fractional t, Num (Vector t), Numeric t) => Matrix t -> Matrix t -> Matrix t -gaussElim a b = r +gaussElim a b = flipudrl r where - go x y = splitColsAt (cols a) (down $ fromBlocks [[x,y]]) - (a1,b1) = go a b - ( _, r) = go (flipud . fliprl $ a1) (flipud . fliprl $ b1) - + flipudrl = flipud . fliprl + splitColsAt n = (takeColumns n &&& dropColumns n) + go f x y = splitColsAt (cols a) (down f $ fromBlocks [[x,y]]) + (a1,b1) = go (snd . swapMax 0) a b + ( _, r) = go id (flipudrl $ a1) (flipudrl $ b1) -------------------------------------------------------------------------------- -- cgit v1.2.3 From 57487d828065ea219cdb33c9dc177b67c60b34c7 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Sun, 14 Jun 2015 19:49:10 +0200 Subject: minor changes --- packages/base/CHANGELOG | 5 +- packages/base/src/Internal/Matrix.hs | 10 +--- packages/base/src/Internal/ST.hs | 59 ++++++++++++++++++---- packages/base/src/Internal/Static.hs | 3 ++ packages/base/src/Internal/Util.hs | 44 ++++++++-------- packages/base/src/Numeric/LinearAlgebra.hs | 10 ++-- packages/base/src/Numeric/LinearAlgebra/Data.hs | 3 +- packages/base/src/Numeric/LinearAlgebra/Devel.hs | 7 ++- packages/base/src/Numeric/LinearAlgebra/HMatrix.hs | 4 +- packages/base/src/Numeric/LinearAlgebra/Static.hs | 6 +-- packages/tests/src/Numeric/LinearAlgebra/Tests.hs | 1 + 11 files changed, 96 insertions(+), 56 deletions(-) (limited to 'packages/base/CHANGELOG') diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG index 27c4d31..93b2594 100644 --- a/packages/base/CHANGELOG +++ b/packages/base/CHANGELOG @@ -7,7 +7,10 @@ * remap, ccompare, sortIndex - * experimental support of type safe modular arithmetic + * experimental support of type safe modular arithmetic, including linear + systems and lu factorization + + * elementary row operations in ST monad * old compatibility modules removed diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs index e0f5ed2..e4b1226 100644 --- a/packages/base/src/Internal/Matrix.hs +++ b/packages/base/src/Internal/Matrix.hs @@ -262,15 +262,7 @@ compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2 ------------------------------------------------------------------ -{- | Supported matrix elements. - - This class provides optimized internal - operations for selected element types. - It provides unoptimised defaults for any 'Storable' type, - so you can create instances simply as: - - >instance Element Foo --} +-- | Supported matrix elements. class (Storable a) => Element a where transdata :: Int -> Vector a -> Int -> Vector a constantD :: a -> Int -> Vector a diff --git a/packages/base/src/Internal/ST.hs b/packages/base/src/Internal/ST.hs index a84ca25..434fe63 100644 --- a/packages/base/src/Internal/ST.hs +++ b/packages/base/src/Internal/ST.hs @@ -10,7 +10,7 @@ -- Stability : provisional -- -- In-place manipulation inside the ST monad. --- See examples/inplace.hs in the distribution. +-- See @examples/inplace.hs@ in the repository. -- ----------------------------------------------------------------------------- @@ -21,8 +21,8 @@ module Internal.ST ( -- * Mutable Matrices STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix, readMatrix, writeMatrix, modifyMatrix, liftSTMatrix, - axpy, scal, swap, extractMatrix, setMatrix, rowOpST, - mutable, +-- axpy, scal, swap, rowOp, + mutable, extractMatrix, setMatrix, rowOper, RowOper(..), RowRange(..), ColRange(..), -- * Unsafe functions newUndefinedVector, unsafeReadVector, unsafeWriteVector, @@ -178,16 +178,55 @@ newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c) -------------------------------------------------------------------------------- -rowOpST :: Element t => Int -> t -> Int -> Int -> Int -> Int -> STMatrix s t -> ST s () -rowOpST c x i1 i2 j1 j2 (STMatrix m) = unsafeIOToST (rowOp c x i1 i2 j1 j2 m) +data ColRange = AllCols + | ColRange Int Int + | Col Int + | FromCol Int -axpy (STMatrix m) a i j = rowOpST 0 a i j 0 (cols m -1) (STMatrix m) -scal (STMatrix m) a i = rowOpST 1 a i i 0 (cols m -1) (STMatrix m) -swap (STMatrix m) i j = rowOpST 2 0 i j 0 (cols m -1) (STMatrix m) +getColRange c AllCols = (0,c-1) +getColRange c (ColRange a b) = (a `mod` c, b `mod` c) +getColRange c (Col a) = (a `mod` c, a `mod` c) +getColRange c (FromCol a) = (a `mod` c, c-1) -extractMatrix (STMatrix m) i1 i2 j1 j2 = unsafeIOToST (extractR m 0 (idxs[i1,i2]) 0 (idxs[j1,j2])) +data RowRange = AllRows + | RowRange Int Int + | Row Int + | FromRow Int + +getRowRange r AllRows = (0,r-1) +getRowRange r (RowRange a b) = (a `mod` r, b `mod` r) +getRowRange r (Row a) = (a `mod` r, a `mod` r) +getRowRange r (FromRow a) = (a `mod` r, r-1) + +data RowOper t = AXPY t Int Int ColRange + | SCAL t RowRange ColRange + | SWAP Int Int ColRange + +rowOper :: (Num t, Element t) => RowOper t -> STMatrix s t -> ST s () + +rowOper (AXPY x i1 i2 r) (STMatrix m) = unsafeIOToST $ rowOp 0 x i1' i2' j1 j2 m + where + (j1,j2) = getColRange (cols m) r + i1' = i1 `mod` (rows m) + i2' = i2 `mod` (rows m) + +rowOper (SCAL x rr rc) (STMatrix m) = unsafeIOToST $ rowOp 1 x i1 i2 j1 j2 m + where + (i1,i2) = getRowRange (rows m) rr + (j1,j2) = getColRange (cols m) rc + +rowOper (SWAP i1 i2 r) (STMatrix m) = unsafeIOToST $ rowOp 2 0 i1' i2' j1 j2 m + where + (j1,j2) = getColRange (cols m) r + i1' = i1 `mod` (rows m) + i2' = i2 `mod` (rows m) + + +extractMatrix (STMatrix m) rr rc = unsafeIOToST (extractR m 0 (idxs[i1,i2]) 0 (idxs[j1,j2])) + where + (i1,i2) = getRowRange (rows m) rr + (j1,j2) = getColRange (cols m) rc --------------------------------------------------------------------------------- mutable :: Storable t => (forall s . (Int, Int) -> STMatrix s t -> ST s u) -> Matrix t -> (Matrix t,u) mutable f a = runST $ do diff --git a/packages/base/src/Internal/Static.hs b/packages/base/src/Internal/Static.hs index 01c2205..0068313 100644 --- a/packages/base/src/Internal/Static.hs +++ b/packages/base/src/Internal/Static.hs @@ -34,6 +34,9 @@ import Text.Printf -------------------------------------------------------------------------------- +type ℝ = Double +type ℂ = Complex Double + newtype Dim (n :: Nat) t = Dim t deriving Show diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs index 2650ac8..09ba21c 100644 --- a/packages/base/src/Internal/Util.hs +++ b/packages/base/src/Internal/Util.hs @@ -65,7 +65,7 @@ import Internal.Element import Internal.Container import Internal.Vectorized import Internal.IO -import Internal.Algorithms hiding (i,Normed,swap) +import Internal.Algorithms hiding (i,Normed,swap,linearSolve') import Numeric.Matrix() import Numeric.Vector() import Internal.Random @@ -155,7 +155,7 @@ infixl 3 & (&) :: Vector Double -> Vector Double -> Vector Double a & b = vjoin [a,b] -{- | horizontal concatenation of real matrices +{- | horizontal concatenation >>> ident 3 ||| konst 7 (3,4) (3><7) @@ -165,7 +165,7 @@ a & b = vjoin [a,b] -} infixl 3 ||| -(|||) :: Matrix Double -> Matrix Double -> Matrix Double +(|||) :: Element t => Matrix t -> Matrix t -> Matrix t a ||| b = fromBlocks [[a,b]] -- | a synonym for ('|||') (unicode 0x00a6, broken bar) @@ -174,9 +174,9 @@ infixl 3 ¦ (¦) = (|||) --- | vertical concatenation of real matrices +-- | vertical concatenation -- -(===) :: Matrix Double -> Matrix Double -> Matrix Double +(===) :: Element t => Matrix t -> Matrix t -> Matrix t infixl 2 === a === b = fromBlocks [[a],[b]] @@ -588,7 +588,7 @@ gaussElim_2 a b = flipudrl r where flipudrl = flipud . fliprl splitColsAt n = (takeColumns n &&& dropColumns n) - go f x y = splitColsAt (cols a) (down f $ fromBlocks [[x,y]]) + go f x y = splitColsAt (cols a) (down f $ x ||| y) (a1,b1) = go (snd . swapMax 0) a b ( _, r) = go id (flipudrl $ a1) (flipudrl $ b1) @@ -600,7 +600,7 @@ gaussElim_1 gaussElim_1 x y = dropColumns (rows x) (flipud $ fromRows s2) where - rs = toRows $ fromBlocks [[x , y]] + rs = toRows $ x ||| y s1 = fromRows $ pivotDown (rows x) 0 rs -- interesting s2 = pivotUp (rows x-1) (toRows $ flipud s1) @@ -637,12 +637,15 @@ pivotUp n xs -------------------------------------------------------------------------------- -gaussElim a b = dropColumns (rows a) $ fst $ mutable gaussST (fromBlocks [[a,b]]) +gaussElim a b = dropColumns (rows a) $ fst $ mutable gaussST (a ||| b) gaussST (r,_) x = do let n = r-1 + axpy m a i j = rowOper (AXPY a i j AllCols) m + swap m i j = rowOper (SWAP i j AllCols) m + scal m a i = rowOper (SCAL a (Row i) AllCols) m forM_ [0..n] $ \i -> do - c <- maxIndex . abs . flatten <$> extractMatrix x i n i i + c <- maxIndex . abs . flatten <$> extractMatrix x (FromRow i) (Col i) swap x i (i+c) a <- readMatrix x i i when (a == 0) $ error "singular!" @@ -656,22 +659,23 @@ gaussST (r,_) x = do axpy x (-b) i j -luST ok (r,c) x = do - let n = r-1 - axpy' m a i j = rowOpST 0 a i j (i+1) (c-1) m - p <- thawMatrix . asColumn . range $ r - forM_ [0..n] $ \i -> do - k <- maxIndex . abs . flatten <$> extractMatrix x i n i i - writeMatrix p i 0 (fi (k+i)) + +luST ok (r,_) x = do + let axpy m a i j = rowOper (AXPY a i j (FromCol (i+1))) m + swap m i j = rowOper (SWAP i j AllCols) m + p <- newUndefinedVector r + forM_ [0..r-1] $ \i -> do + k <- maxIndex . abs . flatten <$> extractMatrix x (FromRow i) (Col i) + writeVector p i (k+i) swap x i (i+k) a <- readMatrix x i i when (ok a) $ do - forM_ [i+1..n] $ \j -> do + forM_ [i+1..r-1] $ \j -> do b <- (/a) <$> readMatrix x j i - axpy' x (-b) i j + axpy x (-b) i j writeMatrix x j i b - v <- unsafeFreezeMatrix p - return (map ti $ toList $ flatten v) + v <- unsafeFreezeVector p + return (toList v) -------------------------------------------------------------------------------- diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs index 0f8efa4..fe524cc 100644 --- a/packages/base/src/Numeric/LinearAlgebra.hs +++ b/packages/base/src/Numeric/LinearAlgebra.hs @@ -80,6 +80,7 @@ module Numeric.LinearAlgebra ( cholSolve, cgSolve, cgSolve', + linearSolve', -- * Inverse and pseudoinverse inv, pinv, pinvTol, @@ -136,8 +137,9 @@ module Numeric.LinearAlgebra ( Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample, -- * Misc - meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, gaussElim, luST, magnit, - ℝ,ℂ,iC, + meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, magnit, + haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, + iC, -- * Auxiliary classes Element, Container, Product, Numeric, LSDiv, Complexable, RealElement, @@ -156,7 +158,7 @@ import Numeric.Vector() import Internal.Matrix import Internal.Container hiding ((<>)) import Internal.Numeric hiding (mul) -import Internal.Algorithms hiding (linearSolve,Normed,orth,luPacked') +import Internal.Algorithms hiding (linearSolve,Normed,orth,luPacked',linearSolve') import qualified Internal.Algorithms as A import Internal.Util import Internal.Random @@ -240,3 +242,5 @@ orth m = orthSVD (Left (1*eps)) m (leftSV m) luPacked' x = mutable (luST (magnit 0)) x +linearSolve' x y = gaussElim x y + diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs index 1c9bb68..fffc2bd 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Data.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs @@ -53,8 +53,7 @@ module Numeric.LinearAlgebra.Data( -- * Matrix extraction Extractor(..), (??), - takeRows, takeLastRows, dropRows, dropLastRows, - takeColumns, takeLastColumns, dropColumns, dropLastColumns, + takeRows, dropRows, takeColumns, dropColumns, subMatrix, (?), (¿), fliprl, flipud, remap, -- * Block matrix diff --git a/packages/base/src/Numeric/LinearAlgebra/Devel.hs b/packages/base/src/Numeric/LinearAlgebra/Devel.hs index f572656..36c5f03 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Devel.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Devel.hs @@ -20,8 +20,7 @@ module Numeric.LinearAlgebra.Devel( module Internal.Foreign, -- * FFI tools - -- | Illustrative usage examples can be found - -- in the @examples\/devel@ folder included in the package. + -- | See @examples/devel@ in the repository. createVector, createMatrix, vec, mat, omat, @@ -36,7 +35,7 @@ module Numeric.LinearAlgebra.Devel( -- * ST -- | In-place manipulation inside the ST monad. - -- See examples\/inplace.hs in the distribution. + -- See @examples/inplace.hs@ in the repository. -- ** Mutable Vectors STVector, newVector, thawVector, freezeVector, runSTVector, @@ -44,7 +43,7 @@ module Numeric.LinearAlgebra.Devel( -- ** Mutable Matrices STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix, readMatrix, writeMatrix, modifyMatrix, liftSTMatrix, - axpy,scal,swap, extractMatrix, setMatrix, mutable, rowOpST, + mutable, extractMatrix, setMatrix, rowOper, RowOper(..), RowRange(..), ColRange(..), -- ** Unsafe functions newUndefinedVector, unsafeReadVector, unsafeWriteVector, diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs index 327f284..11c2487 100644 --- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs +++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs @@ -13,10 +13,10 @@ compatibility with previous version, to be removed module Numeric.LinearAlgebra.HMatrix ( module Numeric.LinearAlgebra, - (¦),(——) + (¦),(——),ℝ,ℂ, ) where import Numeric.LinearAlgebra import Internal.Util - + diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs index dee5b2c..a657bd0 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Static.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs @@ -28,7 +28,7 @@ This module is under active development and the interface is subject to changes. module Numeric.LinearAlgebra.Static( -- * Vector - ℝ, R, + ℝ, R, vec2, vec3, vec4, (&), (#), split, headTail, vector, linspace, range, dim, @@ -71,10 +71,6 @@ import Data.Proxy(Proxy) import Internal.Static import Control.Arrow((***)) - - - - ud1 :: R n -> Vector ℝ ud1 (R (Dim v)) = v diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs index ffa45e7..148bbb9 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs @@ -772,6 +772,7 @@ luBenchN f n x msg = do time (msg ++ " "++ show n) (f m) luBench = do + putStrLn "" luBenchN luPacked 1000 (5::R) "luPacked Double " luBenchN luPacked' 1000 (5::R) "luPacked' Double " luBenchN luPacked' 1000 (5::Mod 9973 I) "luPacked' I mod 9973" -- cgit v1.2.3 From 7376d022b12a27db5a396f89806a709555c1c522 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Mon, 15 Jun 2015 12:52:46 +0200 Subject: documentation, more general cond, remove some unicode, minor changes --- packages/base/CHANGELOG | 2 +- packages/base/src/Internal/CG.hs | 8 +- packages/base/src/Internal/Container.hs | 85 ++++++++++++++++++---- packages/base/src/Internal/Element.hs | 38 ++++++++-- packages/base/src/Internal/Matrix.hs | 1 - packages/base/src/Internal/Modular.hs | 9 ++- packages/base/src/Internal/Numeric.hs | 12 +-- packages/base/src/Internal/Util.hs | 38 ++++++---- packages/base/src/Numeric/LinearAlgebra.hs | 48 +++++++++++- packages/base/src/Numeric/LinearAlgebra/Data.hs | 19 +++-- packages/base/src/Numeric/LinearAlgebra/HMatrix.hs | 5 +- packages/base/src/Numeric/LinearAlgebra/Static.hs | 16 ++-- 12 files changed, 217 insertions(+), 64 deletions(-) (limited to 'packages/base/CHANGELOG') diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG index 93b2594..c3e118c 100644 --- a/packages/base/CHANGELOG +++ b/packages/base/CHANGELOG @@ -5,7 +5,7 @@ * basic support of Int32 and Int64 elements - * remap, ccompare, sortIndex + * remap, more general cond * experimental support of type safe modular arithmetic, including linear systems and lu factorization diff --git a/packages/base/src/Internal/CG.hs b/packages/base/src/Internal/CG.hs index fd14212..758d130 100644 --- a/packages/base/src/Internal/CG.hs +++ b/packages/base/src/Internal/CG.hs @@ -45,13 +45,13 @@ cg sym at a (CGState p r r2 x _) = CGState p' r' r'2 x' rdx ap1 = a p ap | sym = ap1 | otherwise = at ap1 - pap | sym = p <·> ap1 + pap | sym = p <.> ap1 | otherwise = norm2 ap1 ** 2 alpha = r2 / pap dx = scale alpha p x' = x + dx r' = r - scale alpha ap - r'2 = r' <·> r' + r'2 = r' <.> r' beta = r'2 / r2 p' = r' + scale beta p @@ -75,9 +75,9 @@ solveG mat ma meth rawb x0' ϵb ϵx b = mat rawb x0 = if x0' == 0 then konst 0 (dim b) else x0' r0 = b - a x0 - r20 = r0 <·> r0 + r20 = r0 <.> r0 p0 = r0 - nb2 = b <·> b + nb2 = b <.> b ok CGState {..} = cgr2 --- | An infix synonym for 'dot' -(<.>) :: Numeric t => Vector t -> Vector t -> t -(<.>) = dot +infixr 8 <.> +{- | An infix synonym for 'dot' +>>> vector [1,2,3,4] <.> vector [-2,0,1,1] +5.0 -infixr 8 <·>, #> +>>> let 𝑖 = 0:+1 :: C +>>> fromList [1+𝑖,1] <.> fromList [1,1+𝑖] +2.0 :+ 0.0 -{- | infix synonym for 'dot' +-} ->>> vector [1,2,3,4] <·> vector [-2,0,1,1] -5.0 +(<.>) :: Numeric t => Vector t -> Vector t -> t +(<.>) = dot ->>> let 𝑖 = 0:+1 :: ℂ ->>> fromList [1+𝑖,1] <·> fromList [1,1+𝑖] -2.0 :+ 0.0 -(the dot symbol "·" is obtained by Alt-Gr .) --} -(<·>) :: Numeric t => Vector t -> Vector t -> t -(<·>) = dot {- | infix synonym for 'app' @@ -91,6 +86,7 @@ infixr 8 <·>, #> fromList [140.0,320.0] -} +infixr 8 #> (#>) :: Numeric t => Matrix t -> Vector t -> Vector t (#>) = mXv @@ -264,6 +260,24 @@ instance Numeric Z sortVector :: (Ord t, Element t) => Vector t -> Vector t sortVector = sortV +{- | + +>>> m <- randn 4 10 +>>> disp 2 m +4x10 +-0.31 0.41 0.43 -0.19 -0.17 -0.23 -0.17 -1.04 -0.07 -1.24 + 0.26 0.19 0.14 0.83 -1.54 -0.09 0.37 -0.63 0.71 -0.50 +-0.11 -0.10 -1.29 -1.40 -1.04 -0.89 -0.68 0.35 -1.46 1.86 + 1.04 -0.29 0.19 -0.75 -2.20 -0.01 1.06 0.11 -2.09 -1.58 + +>>> disp 2 $ m ?? (All, Pos $ sortIndex (m!1)) +4x10 +-0.17 -1.04 -1.24 -0.23 0.43 0.41 -0.31 -0.17 -0.07 -0.19 +-1.54 -0.63 -0.50 -0.09 0.14 0.19 0.26 0.37 0.71 0.83 +-1.04 0.35 1.86 -0.89 -1.29 -0.10 -0.11 -0.68 -1.46 -1.40 +-2.20 0.11 -1.58 -0.01 0.19 -0.29 1.04 1.06 -2.09 -0.75 + +-} sortIndex :: (Ord t, Element t) => Vector t -> Vector I sortIndex = sortI @@ -273,11 +287,50 @@ ccompare = ccompare' cselect :: (Container c t) => c I -> c t -> c t -> c t -> c t cselect = cselect' +{- | Extract elements from positions given in matrices of rows and columns. + +>>> r +(3><3) + [ 1, 1, 1 + , 1, 2, 2 + , 1, 2, 3 ] +>>> c +(3><3) + [ 0, 1, 5 + , 2, 2, 1 + , 4, 4, 1 ] +>>> m +(4><6) + [ 0, 1, 2, 3, 4, 5 + , 6, 7, 8, 9, 10, 11 + , 12, 13, 14, 15, 16, 17 + , 18, 19, 20, 21, 22, 23 ] + +>>> remap r c m +(3><3) + [ 6, 7, 11 + , 8, 14, 13 + , 10, 16, 19 ] + +The indexes are autoconformable. + +>>> c' +(3><1) + [ 1 + , 2 + , 4 ] +>>> remap r c' m +(3><3) + [ 7, 7, 7 + , 8, 14, 14 + , 10, 16, 22 ] + +-} remap :: Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t remap i j m | minElement i >= 0 && maxElement i < fromIntegral (rows m) && minElement j >= 0 && maxElement j < fromIntegral (cols m) = remapM i' j' m - | otherwise = error $ "out of range index in rmap" + | otherwise = error $ "out of range index in remap" where [i',j'] = conformMs [i,j] diff --git a/packages/base/src/Internal/Element.hs b/packages/base/src/Internal/Element.hs index 4007491..51d5686 100644 --- a/packages/base/src/Internal/Element.hs +++ b/packages/base/src/Internal/Element.hs @@ -80,7 +80,7 @@ breakAt c l = (a++[c],tail b) where (a,b) = break (==c) l -------------------------------------------------------------------------------- - +-- | Specification of indexes for the operator '??'. data Extractor = All | Range Int Int Int @@ -102,7 +102,32 @@ ppext (Drop n) = printf "Drop %d" n ppext (TakeLast n) = printf "TakeLast %d" n ppext (DropLast n) = printf "DropLast %d" n +{- | General matrix slicing. + +>>> m +(4><5) + [ 0, 1, 2, 3, 4 + , 5, 6, 7, 8, 9 + , 10, 11, 12, 13, 14 + , 15, 16, 17, 18, 19 ] + +>>> m ?? (Take 3, DropLast 2) +(3><3) + [ 0, 1, 2 + , 5, 6, 7 + , 10, 11, 12 ] + +>>> m ?? (Pos (idxs[2,1]), All) +(2><5) + [ 10, 11, 12, 13, 14 + , 5, 6, 7, 8, 9 ] + +>>> m ?? (PosCyc (idxs[-7,80]), Range 4 (-2) 0) +(2><3) + [ 9, 7, 5 + , 4, 2, 0 ] +-} infixl 9 ?? (??) :: Element t => Matrix t -> (Extractor,Extractor) -> Matrix t @@ -328,27 +353,30 @@ r >< c = f where ---------------------------------------------------------------- --- | Creates a matrix with the first n rows of another matrix takeRows :: Element t => Int -> Matrix t -> Matrix t takeRows n mt = subMatrix (0,0) (n, cols mt) mt + -- | Creates a matrix with the last n rows of another matrix takeLastRows :: Element t => Int -> Matrix t -> Matrix t takeLastRows n mt = subMatrix (rows mt - n, 0) (n, cols mt) mt --- | Creates a copy of a matrix without the first n rows + dropRows :: Element t => Int -> Matrix t -> Matrix t dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt + -- | Creates a copy of a matrix without the last n rows dropLastRows :: Element t => Int -> Matrix t -> Matrix t dropLastRows n mt = subMatrix (0,0) (rows mt - n, cols mt) mt --- |Creates a matrix with the first n columns of another matrix + takeColumns :: Element t => Int -> Matrix t -> Matrix t takeColumns n mt = subMatrix (0,0) (rows mt, n) mt + -- |Creates a matrix with the last n columns of another matrix takeLastColumns :: Element t => Int -> Matrix t -> Matrix t takeLastColumns n mt = subMatrix (0, cols mt - n) (rows mt, n) mt --- | Creates a copy of a matrix without the first n columns + dropColumns :: Element t => Int -> Matrix t -> Matrix t dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt + -- | Creates a copy of a matrix without the last n columns dropLastColumns :: Element t => Int -> Matrix t -> Matrix t dropLastColumns n mt = subMatrix (0,0) (rows mt, cols mt - n) mt diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs index e4b1226..75e92a5 100644 --- a/packages/base/src/Internal/Matrix.hs +++ b/packages/base/src/Internal/Matrix.hs @@ -378,7 +378,6 @@ foreign import ccall unsafe "transL" ctransL :: TMM Z ---------------------------------------------------------------------- --- | Extracts a submatrix from a matrix. subMatrix :: Element a => (Int,Int) -- ^ (r0,c0) starting position -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs index 3b27310..0d906bb 100644 --- a/packages/base/src/Internal/Modular.hs +++ b/packages/base/src/Internal/Modular.hs @@ -9,8 +9,8 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE GADTs #-} -{-# LANGUAGE TypeFamilies #-} - +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} {- | Module : Internal.Modular @@ -23,7 +23,7 @@ Proof of concept of statically checked modular arithmetic. -} module Internal.Modular( - Mod + Mod, type (./.) ) where import Internal.Vector @@ -49,6 +49,9 @@ import Data.Complex newtype Mod (n :: Nat) t = Mod {unMod:: t} deriving (Storable) +infixr 5 ./. +type (./.) x n = Mod n x + instance (Integral t, Enum t, KnownNat m) => Enum (Mod m t) where toEnum = l0 (\m x -> fromIntegral $ x `mod` (fromIntegral m)) diff --git a/packages/base/src/Internal/Numeric.hs b/packages/base/src/Internal/Numeric.hs index ca17c23..efcde2c 100644 --- a/packages/base/src/Internal/Numeric.hs +++ b/packages/base/src/Internal/Numeric.hs @@ -481,14 +481,16 @@ step = step' -- , 0.0, 100.0, 7.0, 8.0 -- , 0.0, 0.0, 100.0, 12.0 ] -- +-- >>> let chop x = cond (abs x) 1E-6 0 0 x +-- cond - :: (Ord e, Container c e) + :: (Ord e, Container c e, Container c x) => c e -- ^ a -> c e -- ^ b - -> c e -- ^ l - -> c e -- ^ e - -> c e -- ^ g - -> c e -- ^ result + -> c x -- ^ l + -> c x -- ^ e + -> c x -- ^ g + -> c x -- ^ result cond a b l e g = cselect' (ccompare' a b) l e g diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs index 09ba21c..f08f710 100644 --- a/packages/base/src/Internal/Util.hs +++ b/packages/base/src/Internal/Util.hs @@ -85,7 +85,7 @@ type ℤ = Int type ℂ = Complex Double -- | imaginary unit -iC :: ℂ +iC :: C iC = 0:+1 {- | Create a real vector. @@ -94,7 +94,7 @@ iC = 0:+1 fromList [1.0,2.0,3.0,4.0,5.0] -} -vector :: [ℝ] -> Vector ℝ +vector :: [R] -> Vector R vector = fromList {- | Create a real matrix. @@ -108,8 +108,8 @@ vector = fromList -} matrix :: Int -- ^ number of columns - -> [ℝ] -- ^ elements in row order - -> Matrix ℝ + -> [R] -- ^ elements in row order + -> Matrix R matrix c = reshape c . fromList @@ -260,34 +260,34 @@ norm = pnorm PNorm2 class Normed a where - norm_0 :: a -> ℝ - norm_1 :: a -> ℝ - norm_2 :: a -> ℝ - norm_Inf :: a -> ℝ + norm_0 :: a -> R + norm_1 :: a -> R + norm_2 :: a -> R + norm_Inf :: a -> R -instance Normed (Vector ℝ) +instance Normed (Vector R) where norm_0 v = sumElements (step (abs v - scalar (eps*normInf v))) norm_1 = pnorm PNorm1 norm_2 = pnorm PNorm2 norm_Inf = pnorm Infinity -instance Normed (Vector ℂ) +instance Normed (Vector C) where norm_0 v = sumElements (step (fst (fromComplex (abs v)) - scalar (eps*normInf v))) norm_1 = pnorm PNorm1 norm_2 = pnorm PNorm2 norm_Inf = pnorm Infinity -instance Normed (Matrix ℝ) +instance Normed (Matrix R) where norm_0 = norm_0 . flatten norm_1 = pnorm PNorm1 norm_2 = pnorm PNorm2 norm_Inf = pnorm Infinity -instance Normed (Matrix ℂ) +instance Normed (Matrix C) where norm_0 = norm_0 . flatten norm_1 = pnorm PNorm1 @@ -323,12 +323,22 @@ instance Normed (Vector (Complex Float)) norm_Inf = norm_Inf . double -norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> ℝ +norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R norm_Frob = norm_2 . flatten -norm_nuclear :: Field t => Matrix t -> ℝ +norm_nuclear :: Field t => Matrix t -> R norm_nuclear = sumElements . singularValues +{- | Check if the absolute value or complex magnitude is greater than a given threshold + +>>> magnit 1E-6 (1E-12 :: R) +False +>>> magnit 1E-6 (3+iC :: C) +True +>>> magnit 0 (3 :: I ./. 5) +True + +-} magnit :: (Element t, Normed (Vector t)) => R -> t -> Bool magnit e x = norm_1 (fromList [x]) > e diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs index fe524cc..2e6b8ca 100644 --- a/packages/base/src/Numeric/LinearAlgebra.hs +++ b/packages/base/src/Numeric/LinearAlgebra.hs @@ -47,7 +47,7 @@ module Numeric.LinearAlgebra ( -- * Products -- ** dot - dot, (<·>), + dot, (<.>), -- ** matrix-vector app, (#>), (<#), (!#>), -- ** matrix-matrix @@ -240,7 +240,53 @@ nullspace m = nullspaceSVD (Left (1*eps)) m (rightSV m) -- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'. orth m = orthSVD (Left (1*eps)) m (leftSV m) +{- | Experimental implementation of LU factorization + working on any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'. + +>>> let m = ident 5 + (5><5) [0..] :: Matrix (Z ./. 17) +(5><5) + [ 1, 1, 2, 3, 4 + , 5, 7, 7, 8, 9 + , 10, 11, 13, 13, 14 + , 15, 16, 0, 2, 2 + , 3, 4, 5, 6, 8 ] + +>>> let (l,u,p,s) = luFact $ luPacked' m +>>> l +(5><5) + [ 1, 0, 0, 0, 0 + , 6, 1, 0, 0, 0 + , 12, 7, 1, 0, 0 + , 7, 10, 7, 1, 0 + , 8, 2, 6, 11, 1 ] +>>> u +(5><5) + [ 15, 16, 0, 2, 2 + , 0, 13, 7, 13, 14 + , 0, 0, 15, 0, 11 + , 0, 0, 0, 15, 15 + , 0, 0, 0, 0, 1 ] + +-} luPacked' x = mutable (luST (magnit 0)) x +{- | Experimental implementation of gaussian elimination + working on any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'. + +>>> let a = (2><2) [1,2,3,5] :: Matrix (Z ./. 13) +(2><2) + [ 1, 2 + , 3, 5 ] +>>> b +(2><3) + [ 5, 1, 3 + , 8, 6, 3 ] + +>>> let x = linearSolve' a b +(2><3) + [ 4, 7, 4 + , 7, 10, 6 ] + +-} linearSolve' x y = gaussElim x y diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs index fffc2bd..2a45771 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Data.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE TypeOperators #-} + -------------------------------------------------------------------------------- {- | Module : Numeric.LinearAlgebra.Data @@ -17,11 +19,11 @@ module Numeric.LinearAlgebra.Data( -- | 1D arrays are storable vectors from the vector package. There is no distinction -- between row and column vectors. - fromList, toList, vector, (|>), + fromList, toList, (|>), vector, range, idxs, -- * Matrix - matrix, (><), tr, tr', + (><), matrix, tr, tr', -- * Dimensions @@ -43,7 +45,7 @@ module Numeric.LinearAlgebra.Data( Indexable(..), -- * Construction - scalar, Konst(..), Build(..), assoc, accum, linspace, range, idxs, -- ones, zeros, + scalar, Konst(..), Build(..), assoc, accum, linspace, -- ones, zeros, -- * Diagonal ident, diag, diagl, diagRect, takeDiag, @@ -53,16 +55,19 @@ module Numeric.LinearAlgebra.Data( -- * Matrix extraction Extractor(..), (??), + takeRows, dropRows, takeColumns, dropColumns, - subMatrix, (?), (¿), fliprl, flipud, remap, + subMatrix, (?), (¿), fliprl, flipud, + + remap, -- * Block matrix fromBlocks, (|||), (===), diagBlock, repmat, toBlocks, toBlocksEvery, -- * Mapping functions conj, cmap, cmod, - - step, cond, ccompare, cselect, + + step, cond, -- * Find elements find, maxIndex, minIndex, maxElement, minElement, @@ -87,7 +92,7 @@ module Numeric.LinearAlgebra.Data( separable, fromArray2D, module Data.Complex, - R,C,I,Z,Mod, + R,C,I,Z,Mod, type(./.), Vector, Matrix, GMatrix, nRows, nCols ) where diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs index 11c2487..8adaaaf 100644 --- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs +++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs @@ -13,10 +13,13 @@ compatibility with previous version, to be removed module Numeric.LinearAlgebra.HMatrix ( module Numeric.LinearAlgebra, - (¦),(——),ℝ,ℂ, + (¦),(——),ℝ,ℂ,(<·>) ) where import Numeric.LinearAlgebra import Internal.Util +infixr 8 <·> +(<·>) :: Numeric t => Vector t -> Vector t -> t +(<·>) = dot diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs index a657bd0..d0a790d 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Static.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs @@ -44,7 +44,7 @@ module Numeric.LinearAlgebra.Static( -- * Complex C, M, Her, her, 𝑖, -- * Products - (<>),(#>),(<·>), + (<>),(#>),(<.>), -- * Linear Systems linSolve, (<\>), -- * Factorizations @@ -55,13 +55,13 @@ module Numeric.LinearAlgebra.Static( Disp(..), Domain(..), withVector, withMatrix, toRows, toColumns, - Sized(..), Diag(..), Sym, sym, mTm, unSym + Sized(..), Diag(..), Sym, sym, mTm, unSym, (<·>) ) where import GHC.TypeLits import Numeric.LinearAlgebra hiding ( - (<>),(#>),(<·>),Konst(..),diag, disp,(===),(|||), + (<>),(#>),(<.>),Konst(..),diag, disp,(===),(|||), row,col,vector,matrix,linspace,toRows,toColumns, (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH', eigenvalues,eigenvaluesSH,eigenvaluesSH',build, @@ -207,6 +207,10 @@ infixr 8 <·> (<·>) :: R n -> R n -> ℝ (<·>) = dotR +infixr 8 <.> +(<.>) :: R n -> R n -> ℝ +(<.>) = dotR + -------------------------------------------------------------------------------- class Diag m d | m -> d @@ -496,7 +500,7 @@ appC m v = mkC (extract m LA.#> extract v) dotC :: KnownNat n => C n -> C n -> ℂ dotC (unwrap -> u) (unwrap -> v) | singleV u || singleV v = sumElements (conj u * v) - | otherwise = u LA.<·> v + | otherwise = u LA.<.> v crossC :: C 3 -> C 3 -> C 3 @@ -584,12 +588,12 @@ test = (ok,info) where q = tm :: L 10 3 - thingD = vjoin [ud1 u, 1] LA.<·> tr m LA.#> m LA.#> ud1 v + thingD = vjoin [ud1 u, 1] LA.<.> tr m LA.#> m LA.#> ud1 v where m = LA.matrix 3 [1..30] precS = (1::Double) + (2::Double) * ((1 :: R 3) * (u & 6)) <·> konst 2 #> v - precD = 1 + 2 * vjoin[ud1 u, 6] LA.<·> LA.konst 2 (LA.size (ud1 u) +1, LA.size (ud1 v)) LA.#> ud1 v + precD = 1 + 2 * vjoin[ud1 u, 6] LA.<.> LA.konst 2 (LA.size (ud1 u) +1, LA.size (ud1 v)) LA.#> ud1 v splittest -- cgit v1.2.3 From b4873dbd201e0e887fb9cb5b5fe55774fa6fbe78 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Sun, 12 Jul 2015 14:10:51 +0200 Subject: documentation --- packages/base/CHANGELOG | 12 ++++++++---- packages/base/hmatrix.cabal | 6 +++++- packages/base/src/Internal/Algorithms.hs | 6 +++--- packages/base/src/Internal/Container.hs | 2 +- packages/base/src/Internal/Element.hs | 4 ++-- packages/base/src/Internal/Vector.hs | 2 +- packages/base/src/Numeric/LinearAlgebra.hs | 13 +++++++++---- packages/base/src/Numeric/LinearAlgebra/Data.hs | 21 +++++++++++++++------ packages/base/src/Numeric/LinearAlgebra/HMatrix.hs | 6 +++++- packages/base/src/Numeric/LinearAlgebra/Static.hs | 4 ++-- 10 files changed, 51 insertions(+), 25 deletions(-) (limited to 'packages/base/CHANGELOG') diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG index c3e118c..0336a28 100644 --- a/packages/base/CHANGELOG +++ b/packages/base/CHANGELOG @@ -1,21 +1,25 @@ 0.17.0.0 -------- - * improved matrix extraction using vectors of indexes (??) + * improved matrix extraction (??) and rectangular matrix slices without data copy * basic support of Int32 and Int64 elements - * remap, more general cond + * remap, more general cond, sortIndex * experimental support of type safe modular arithmetic, including linear systems and lu factorization - * elementary row operations in ST monad + * elementary row operations and inplace matrix slice products in the ST monad - * old compatibility modules removed + * Improved development tools. + + * old compatibility modules removed, simpler organization of internal modules * unitary, pairwiseD2, tr' + * ldlPacked, ldlSolve + 0.16.1.0 -------- diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal index f725341..7b25349 100644 --- a/packages/base/hmatrix.cabal +++ b/packages/base/hmatrix.cabal @@ -7,7 +7,11 @@ Maintainer: Alberto Ruiz Stability: provisional Homepage: https://github.com/albertoruiz/hmatrix Synopsis: Numeric Linear Algebra -Description: Linear algebra based on BLAS and LAPACK. +Description: Linear systems, matrix decompositions, and other numerical computations based on BLAS and LAPACK. + . + The standard interface is provided by the module "Numeric.LinearAlgebra". + . + A safer interface with statically checked dimensions is provided by "Numeric.LinearAlgebra.Static". . Code examples: diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs index c8b2d3e..99c90aa 100644 --- a/packages/base/src/Internal/Algorithms.hs +++ b/packages/base/src/Internal/Algorithms.hs @@ -470,14 +470,14 @@ rq m = {-# SCC "rq" #-} (r,q) where -- | Hessenberg factorization. -- --- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary +-- If @(p,h) = hess m@ then @m == p \<> h \<> tr p@, where p is unitary -- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal). hess :: Field t => Matrix t -> (Matrix t, Matrix t) hess = hess' -- | Schur factorization. -- --- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary +-- If @(u,s) = schur m@ then @m == u \<> s \<> tr u@, where u is unitary -- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is -- upper triangular in 2x2 blocks. -- @@ -497,7 +497,7 @@ cholSH = {-# SCC "cholSH" #-} cholSH' -- | Cholesky factorization of a positive definite hermitian or symmetric matrix. -- --- If @c = chol m@ then @c@ is upper triangular and @m == ctrans c \<> c@. +-- If @c = chol m@ then @c@ is upper triangular and @m == tr c \<> c@. chol :: Field t => Matrix t -> Matrix t chol m | exactHermitian m = cholSH m | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" diff --git a/packages/base/src/Internal/Container.hs b/packages/base/src/Internal/Container.hs index 8926fac..307c6a8 100644 --- a/packages/base/src/Internal/Container.hs +++ b/packages/base/src/Internal/Container.hs @@ -72,7 +72,7 @@ infixr 8 <.> -{- | infix synonym for 'app' +{- | dense matrix-vector product >>> let m = (2><3) [1..] >>> m diff --git a/packages/base/src/Internal/Element.hs b/packages/base/src/Internal/Element.hs index 6d86f3d..a459678 100644 --- a/packages/base/src/Internal/Element.hs +++ b/packages/base/src/Internal/Element.hs @@ -325,9 +325,9 @@ takeDiag m = fromList [flatten m @> (k*cols m+k) | k <- [0 .. min (rows m) (cols ------------------------------------------------------------ -{- | create a general matrix +{- | Create a matrix from a list of elements ->>> (2><3) [2, 4, 7+2*𝑖, -3, 11, 0] +>>> (2><3) [2, 4, 7+2*iC, -3, 11, 0] (2><3) [ 2.0 :+ 0.0, 4.0 :+ 0.0, 7.0 :+ 2.0 , (-3.0) :+ (-0.0), 11.0 :+ 0.0, 0.0 :+ 0.0 ] diff --git a/packages/base/src/Internal/Vector.hs b/packages/base/src/Internal/Vector.hs index 29b6797..de4e670 100644 --- a/packages/base/src/Internal/Vector.hs +++ b/packages/base/src/Internal/Vector.hs @@ -127,7 +127,7 @@ n |> l l' = take n l --- | Create a vector of indexes, useful for matrix extraction using '??' +-- | Create a vector of indexes, useful for matrix extraction using '(??)' idxs :: [Int] -> Vector I idxs js = fromList (map fromIntegral js) :: Vector I diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs index dd4cc67..9cab635 100644 --- a/packages/base/src/Numeric/LinearAlgebra.hs +++ b/packages/base/src/Numeric/LinearAlgebra.hs @@ -8,11 +8,16 @@ License : BSD3 Maintainer : Alberto Ruiz Stability : provisional + -} ----------------------------------------------------------------------------- module Numeric.LinearAlgebra ( - -- * Basic types and data processing + -- * Basic Types and data manipulation + -- | This package works with 2D ('Matrix') and 1D ('Vector') + -- arrays of real ('R') or complex ('C') double precision numbers. + -- Single precision and machine integers are also supported for + -- basic arithmetic and data manipulation. module Numeric.LinearAlgebra.Data, -- * Numeric classes @@ -51,9 +56,9 @@ module Numeric.LinearAlgebra ( -- ** dot dot, (<.>), -- ** matrix-vector - app, (#>), (<#), (!#>), + (#>), (<#), (!#>), -- ** matrix-matrix - mul, (<>), + (<>), -- | The matrix product is also implemented in the "Data.Monoid" instance, where -- single-element matrices (created from numeric literals or using 'scalar') -- are used for scaling. @@ -172,7 +177,7 @@ import Internal.Sparse((!#>)) import Internal.CG import Internal.Conversion -{- | infix synonym of 'mul' +{- | dense matrix product >>> let a = (3><5) [1..] >>> a diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs index d2843c2..a389aac 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Data.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs @@ -8,21 +8,29 @@ License : BSD3 Maintainer : Alberto Ruiz Stability : provisional -Basic data processing. +This module provides functions for creation and manipulation of vectors and matrices, IO, and other utilities. -} -------------------------------------------------------------------------------- module Numeric.LinearAlgebra.Data( + -- * Elements + R,C,I,Z,type(./.), + -- * Vector - -- | 1D arrays are storable vectors from the vector package. There is no distinction - -- between row and column vectors. + {- | 1D arrays are storable vectors directly reexported from the vector package. + -} fromList, toList, (|>), vector, range, idxs, -- * Matrix + {- | The main data type of hmatrix is a 2D dense array defined on top of + a storable vector. The internal representation is suitable for direct + interface with standard numeric libraries. + -} + (><), matrix, tr, tr', -- * Dimensions @@ -56,8 +64,9 @@ module Numeric.LinearAlgebra.Data( -- * Matrix extraction Extractor(..), (??), - takeRows, dropRows, takeColumns, dropColumns, - subMatrix, (?), (¿), fliprl, flipud, + (?), (¿), fliprl, flipud, + + subMatrix, takeRows, dropRows, takeColumns, dropColumns, remap, @@ -92,7 +101,7 @@ module Numeric.LinearAlgebra.Data( separable, fromArray2D, module Data.Complex, - R,C,I,Z,Mod, type(./.), + Mod, Vector, Matrix, GMatrix, nRows, nCols ) where diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs index 8adaaaf..bac1c0c 100644 --- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs +++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs @@ -13,7 +13,7 @@ compatibility with previous version, to be removed module Numeric.LinearAlgebra.HMatrix ( module Numeric.LinearAlgebra, - (¦),(——),ℝ,ℂ,(<·>) + (¦),(——),ℝ,ℂ,(<·>),app,mul ) where import Numeric.LinearAlgebra @@ -23,3 +23,7 @@ infixr 8 <·> (<·>) :: Numeric t => Vector t -> Vector t -> t (<·>) = dot +app m v = m #> v + +mul a b = a <> b + diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs index d0a790d..0dab0e6 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Static.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs @@ -22,7 +22,7 @@ Stability : experimental Experimental interface with statically checked dimensions. -This module is under active development and the interface is subject to changes. +See code examples at http://dis.um.es/~alberto/hmatrix/static.html. -} @@ -65,7 +65,7 @@ import Numeric.LinearAlgebra hiding ( row,col,vector,matrix,linspace,toRows,toColumns, (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH', eigenvalues,eigenvaluesSH,eigenvaluesSH',build, - qr,size,app,mul,dot,chol,range,R,C) + qr,size,dot,chol,range,R,C) import qualified Numeric.LinearAlgebra as LA import Data.Proxy(Proxy) import Internal.Static -- cgit v1.2.3 From 39f2bbe937ccbf786af0a326e7aa065890ee331e Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 17 Jul 2015 21:03:13 +0200 Subject: documentation --- packages/base/CHANGELOG | 26 +++++++++++++++++--------- packages/base/hmatrix.cabal | 4 ++-- packages/base/src/Numeric/LinearAlgebra.hs | 16 +++++++++------- 3 files changed, 28 insertions(+), 18 deletions(-) (limited to 'packages/base/CHANGELOG') diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG index 0336a28..5b0adc1 100644 --- a/packages/base/CHANGELOG +++ b/packages/base/CHANGELOG @@ -1,24 +1,32 @@ 0.17.0.0 -------- - * improved matrix extraction (??) and rectangular matrix slices without data copy + * eigSH, chol, and other functions that work with Hermitian or symmetric matrices + now take a special "Sym" argument that can be created by means of "sym" + or "mTm". The unchecked versions of those functions have been removed and we + use "trustSym" to create the Sym type when the matrix is known to be Hermitian/symmetric. + + * Improved matrix extraction (??) and rectangular matrix slices without data copy + + * Basic support of Int32 and Int64 elements - * basic support of Int32 and Int64 elements - * remap, more general cond, sortIndex - * experimental support of type safe modular arithmetic, including linear - systems and lu factorization - - * elementary row operations and inplace matrix slice products in the ST monad + * Experimental support of type safe modular arithmetic, including linear + system solver and LU factorization + + * Elementary row operations and inplace matrix slice products in the ST monad * Improved development tools. - * old compatibility modules removed, simpler organization of internal modules + * Old compatibility modules removed, simpler organization of internal modules * unitary, pairwiseD2, tr' - * ldlPacked, ldlSolve + * ldlPacked, ldlSolve for indefinite symmetric systems (apparently not faster + than the general solver based on the LU) + + * LU, LDL, and QR types for these compact decompositions. 0.16.1.0 -------- diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal index 7b25349..31bea3e 100644 --- a/packages/base/hmatrix.cabal +++ b/packages/base/hmatrix.cabal @@ -9,9 +9,9 @@ Homepage: https://github.com/albertoruiz/hmatrix Synopsis: Numeric Linear Algebra Description: Linear systems, matrix decompositions, and other numerical computations based on BLAS and LAPACK. . - The standard interface is provided by the module "Numeric.LinearAlgebra". + Standard interface: "Numeric.LinearAlgebra". . - A safer interface with statically checked dimensions is provided by "Numeric.LinearAlgebra.Static". + Safer interface with statically checked dimensions: "Numeric.LinearAlgebra.Static". . Code examples: diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs index a1c0158..ea54932 100644 --- a/packages/base/src/Numeric/LinearAlgebra.hs +++ b/packages/base/src/Numeric/LinearAlgebra.hs @@ -35,8 +35,9 @@ module Numeric.LinearAlgebra ( -- * Autoconformable dimensions -- | - -- In arithmetic operations single-element vectors and matrices - -- (created from numeric literals or using 'scalar') automatically + -- In most operations, single-element vectors and matrices + -- (created from numeric literals or using 'scalar'), and matrices + -- with just one row or column, automatically -- expand to match the dimensions of the other operand: -- -- >>> 5 + 2*ident 3 :: Matrix Double @@ -45,11 +46,12 @@ module Numeric.LinearAlgebra ( -- , 5.0, 7.0, 5.0 -- , 5.0, 5.0, 7.0 ] -- - -- >>> matrix 3 [1..9] + matrix 1 [10,20,30] - -- (3><3) - -- [ 11.0, 12.0, 13.0 - -- , 24.0, 25.0, 26.0 - -- , 37.0, 38.0, 39.0 ] + -- >>> (4><3) [1..] + row [10,20,30] + -- (4><3) + -- [ 11.0, 22.0, 33.0 + -- , 14.0, 25.0, 36.0 + -- , 17.0, 28.0, 39.0 + -- , 20.0, 31.0, 42.0 ] -- -- * Products -- cgit v1.2.3 From 792864d3ec95d198a751591256c302aed11d8392 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Tue, 21 Jul 2015 11:11:49 +0200 Subject: change names: Herm, unSym, mTm, and rnf instances --- packages/base/CHANGELOG | 4 +- packages/base/src/Internal/Algorithms.hs | 71 +++++++++++----------- packages/base/src/Internal/Util.hs | 2 +- packages/base/src/Numeric/LinearAlgebra.hs | 4 +- packages/base/src/Numeric/LinearAlgebra/Static.hs | 2 +- .../src/Numeric/LinearAlgebra/Tests/Instances.hs | 13 ++-- .../src/Numeric/LinearAlgebra/Tests/Properties.hs | 2 +- 7 files changed, 50 insertions(+), 48 deletions(-) (limited to 'packages/base/CHANGELOG') diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG index 5b0adc1..581d2ac 100644 --- a/packages/base/CHANGELOG +++ b/packages/base/CHANGELOG @@ -2,9 +2,9 @@ -------- * eigSH, chol, and other functions that work with Hermitian or symmetric matrices - now take a special "Sym" argument that can be created by means of "sym" + now take a special "Herm" argument that can be created by means of "sym" or "mTm". The unchecked versions of those functions have been removed and we - use "trustSym" to create the Sym type when the matrix is known to be Hermitian/symmetric. + use "trustSym" to create the Herm type when the matrix is known to be Hermitian/symmetric. * Improved matrix extraction (??) and rectangular matrix slices without data copy diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs index ee3ddff..c4f1a60 100644 --- a/packages/base/src/Internal/Algorithms.hs +++ b/packages/base/src/Internal/Algorithms.hs @@ -4,12 +4,6 @@ {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE FlexibleContexts, FlexibleInstances #-} -{-# LANGUAGE CPP #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE TypeFamilies #-} - ----------------------------------------------------------------------------- {- | Module : Internal.Algorithms @@ -376,8 +370,8 @@ ldlPackedSH x = {-# SCC "ldlPacked" #-} LDL m p (m,p) = ldlPacked' x -- | Obtains the LDL decomposition of a matrix in a compact data structure suitable for 'ldlSolve'. -ldlPacked :: Field t => Her t -> LDL t -ldlPacked (Her m) = ldlPackedSH m +ldlPacked :: Field t => Herm t -> LDL t +ldlPacked (Herm m) = ldlPackedSH m -- | Solution of a linear system (for several right hand sides) from a precomputed LDL factorization obtained by 'ldlPacked'. -- @@ -461,18 +455,23 @@ fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575] 3.000 5.000 6.000 -} -eigSH :: Field t => Her t -> (Vector Double, Matrix t) -eigSH (Her m) = eigSH' m +eigSH :: Field t => Herm t -> (Vector Double, Matrix t) +eigSH (Herm m) = eigSH' m -- | Eigenvalues (in descending order) of a complex hermitian or real symmetric matrix. -eigenvaluesSH :: Field t => Her t -> Vector Double -eigenvaluesSH (Her m) = eigenvaluesSH' m +eigenvaluesSH :: Field t => Herm t -> Vector Double +eigenvaluesSH (Herm m) = eigenvaluesSH' m -------------------------------------------------------------- -- | QR decomposition of a matrix in compact form. (The orthogonal matrix is not explicitly formed.) data QR t = QR (Matrix t) (Vector t) +instance (NFData t, Numeric t) => NFData (QR t) + where + rnf (QR m _) = rnf m + + -- | QR factorization. -- -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular. @@ -533,12 +532,12 @@ cholSH = cholSH' -- | Cholesky factorization of a positive definite hermitian or symmetric matrix. -- -- If @c = chol m@ then @c@ is upper triangular and @m == tr c \<> c@. -chol :: Field t => Her t -> Matrix t -chol (Her m) = {-# SCC "chol" #-} cholSH' m +chol :: Field t => Herm t -> Matrix t +chol (Herm m) = {-# SCC "chol" #-} cholSH' m -- | Similar to 'chol', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'. -mbChol :: Field t => Her t -> Maybe (Matrix t) -mbChol (Her m) = {-# SCC "mbChol" #-} mbCholSH' m +mbChol :: Field t => Herm t -> Maybe (Matrix t) +mbChol (Herm m) = {-# SCC "mbChol" #-} mbCholSH' m @@ -977,10 +976,10 @@ relativeError norm a b = r -- | Generalized symmetric positive definite eigensystem Av = lBv, -- for A and B symmetric, B positive definite. geigSH :: Field t - => Her t -- ^ A - -> Her t -- ^ B + => Herm t -- ^ A + -> Herm t -- ^ B -> (Vector Double, Matrix t) -geigSH (Her a) (Her b) = geigSH' a b +geigSH (Herm a) (Herm b) = geigSH' a b geigSH' :: Field t => Matrix t -- ^ A @@ -999,29 +998,33 @@ geigSH' a b = (l,v') -- | A matrix that, by construction, it is known to be complex Hermitian or real symmetric. -- --- It can be created using 'sym', 'xTx', or 'trustSym', and the matrix can be extracted using 'her'. -data Her t = Her (Matrix t) deriving Show +-- It can be created using 'sym', 'mTm', or 'trustSym', and the matrix can be extracted using 'unSym'. +newtype Herm t = Herm (Matrix t) deriving Show + +instance (NFData t, Numeric t) => NFData (Herm t) + where + rnf (Herm m) = rnf m --- | Extract the general matrix from a 'Her' structure, forgetting its symmetric or Hermitian property. -her :: Her t -> Matrix t -her (Her x) = x +-- | Extract the general matrix from a 'Herm' structure, forgetting its symmetric or Hermitian property. +unSym :: Herm t -> Matrix t +unSym (Herm x) = x -- | Compute the complex Hermitian or real symmetric part of a square matrix (@(x + tr x)/2@). -sym :: Field t => Matrix t -> Her t -sym x = Her (scale 0.5 (tr x `add` x)) +sym :: Field t => Matrix t -> Herm t +sym x = Herm (scale 0.5 (tr x `add` x)) -- | Compute the contraction @tr x <> x@ of a general matrix. -xTx :: Numeric t => Matrix t -> Her t -xTx x = Her (tr x `mXm` x) +mTm :: Numeric t => Matrix t -> Herm t +mTm x = Herm (tr x `mXm` x) -instance Field t => Linear t Her where - scale x (Her m) = Her (scale x m) +instance Field t => Linear t Herm where + scale x (Herm m) = Herm (scale x m) -instance Field t => Additive (Her t) where - add (Her a) (Her b) = Her (a `add` b) +instance Field t => Additive (Herm t) where + add (Herm a) (Herm b) = Herm (a `add` b) -- | At your own risk, declare that a matrix is complex Hermitian or real symmetric -- for usage in 'chol', 'eigSH', etc. Only a triangular part of the matrix will be used. -trustSym :: Matrix t -> Her t -trustSym x = (Her x) +trustSym :: Matrix t -> Herm t +trustSym x = (Herm x) diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs index 36b7855..cf42961 100644 --- a/packages/base/src/Internal/Util.hs +++ b/packages/base/src/Internal/Util.hs @@ -462,7 +462,7 @@ null1 :: Matrix R -> Vector R null1 = last . toColumns . snd . rightSV -- | solution of overconstrained homogeneous symmetric linear system -null1sym :: Her R -> Vector R +null1sym :: Herm R -> Vector R null1sym = last . toColumns . snd . eigSH -------------------------------------------------------------------------------- diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs index ea54932..6a9c33a 100644 --- a/packages/base/src/Numeric/LinearAlgebra.hs +++ b/packages/base/src/Numeric/LinearAlgebra.hs @@ -154,9 +154,9 @@ module Numeric.LinearAlgebra ( -- * Misc meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, magnit, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, - iC, sym, xTx, trustSym, her, + iC, sym, mTm, trustSym, unSym, -- * Auxiliary classes - Element, Container, Product, Numeric, LSDiv, Her, + Element, Container, Product, Numeric, LSDiv, Herm, Complexable, RealElement, RealOf, ComplexOf, SingleOf, DoubleOf, IndexOf, diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs index ded69fa..843c727 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Static.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs @@ -65,7 +65,7 @@ import Numeric.LinearAlgebra hiding ( row,col,vector,matrix,linspace,toRows,toColumns, (<\>),fromList,takeDiag,svd,eig,eigSH, eigenvalues,eigenvaluesSH,build, - qr,size,dot,chol,range,R,C,Her,her,sym) + qr,size,dot,chol,range,R,C,sym,mTm,unSym) import qualified Numeric.LinearAlgebra as LA import Data.Proxy(Proxy) import Internal.Static diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs index 4704989..3d5441d 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs @@ -81,8 +81,7 @@ instance (Field a, Arbitrary a) => Arbitrary (Rot a) where -- a complex hermitian or real symmetric matrix ---newtype (Her a) = Her (Matrix a) deriving Show -instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where +instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Herm a) where arbitrary = do Sq m <- arbitrary let m' = m/2 @@ -127,7 +126,7 @@ instance (Numeric a, ArbitraryField a, Num (Vector a)) arbitrary = do m <- arbitrary let (_,v) = eigSH m - n = rows (her m) + n = rows (unSym m) l <- replicateM n (choose (0,100)) let s = diag (fromList l) p = v <> real s <> tr v @@ -161,8 +160,8 @@ fM m = m :: FM zM m = m :: ZM -rHer m = her m :: RM -cHer m = her m :: CM +rHer m = unSym m :: RM +cHer m = unSym m :: CM rRot (Rot m) = m :: RM cRot (Rot m) = m :: CM @@ -176,8 +175,8 @@ cWC (WC m) = m :: CM rSqWC (SqWC m) = m :: RM cSqWC (SqWC m) = m :: CM -rSymWC (SqWC m) = sym m :: Her R -cSymWC (SqWC m) = sym m :: Her C +rSymWC (SqWC m) = sym m :: Herm R +cSymWC (SqWC m) = sym m :: Herm C rPosDef (PosDef m) = m :: RM cPosDef (PosDef m) = m :: CM diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs index 2ac3588..720b7bd 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs @@ -263,7 +263,7 @@ multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a) linearSolveProp f m = f m m |~| ident (rows m) -linearSolvePropH f m = f m (her m) |~| ident (rows (her m)) +linearSolvePropH f m = f m (unSym m) |~| ident (rows (unSym m)) linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b) where q = min (rows a) (cols a) -- cgit v1.2.3