diff --git a/hmatrix.cabal b/hmatrix.cabal --- a/hmatrix.cabal +++ b/hmatrix.cabal @@ -1,5 +1,5 @@ Name: hmatrix -Version: 0.20.0.0 +Version: 0.20.1 License: BSD3 License-file: LICENSE Author: Alberto Ruiz @@ -18,7 +18,7 @@ Category: Math tested-with: GHC==8.2, GHC==8.4 -cabal-version: >=1.8 +cabal-version: >=1.18 build-type: Simple @@ -29,22 +29,30 @@ flag openblas description: Link with OpenBLAS (https://github.com/xianyi/OpenBLAS) optimized libraries. default: False - manual: True + manual: True flag disable-default-paths description: When enabled, don't add default hardcoded include/link dirs by default. Needed for hermetic builds like in nix. default: False - manual: True + manual: True +flag no-random_r + description: When enabled, don't depend on the random_r() C function. + default: False + manual: True + library - Build-Depends: base >= 4.9 && < 5, + default-language: Haskell2010 + + Build-Depends: base >= 4.8 && < 5, binary, array, deepseq, random, split, bytestring, + primitive, storable-complex, semigroups, vector >= 0.11 @@ -84,7 +92,7 @@ src/Internal/C/vector-aux.c - extensions: ForeignFunctionInterface + other-extensions: ForeignFunctionInterface ghc-options: -Wall -fno-warn-missing-signatures @@ -99,6 +107,9 @@ cc-options: -msse2 + if flag(no-random_r) + cc-options: -DNO_RANDOM_R + if os(OSX) if flag(openblas) if !flag(disable-default-paths) @@ -128,10 +139,11 @@ extra-lib-dirs: /usr/local/lib include-dirs: /usr/local/include extra-libraries: gfortran + extra-lib-dirs: /usr/local/lib/gcc9 /usr/local/lib/gcc8 /usr/local/lib/gcc7 if os(windows) if flag(openblas) - extra-libraries: libopenblas, libgcc_s_seh-1, libgfortran, libquadmath-0 + extra-libraries: openblas else extra-libraries: blas lapack @@ -150,4 +162,3 @@ source-repository head type: git location: https://github.com/albertoruiz/hmatrix - diff --git a/src/Internal/Algorithms.hs b/src/Internal/Algorithms.hs --- a/src/Internal/Algorithms.hs +++ b/src/Internal/Algorithms.hs @@ -1,7 +1,6 @@ {-# LANGUAGE FlexibleContexts, FlexibleInstances #-} {-# LANGUAGE CPP #-} {-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE TypeFamilies #-} {-# OPTIONS_GHC -fno-warn-missing-signatures #-} @@ -27,6 +26,10 @@ UpLo(..) ) where +#if MIN_VERSION_base(4,11,0) +import Prelude hiding ((<>)) +#endif + import Internal.Vector import Internal.Matrix import Internal.Element @@ -1159,4 +1162,3 @@ -- for usage in 'chol', 'eigSH', etc. Only a triangular part of the matrix will be used. trustSym :: Matrix t -> Herm t trustSym x = (Herm x) - diff --git a/src/Internal/C/lapack-aux.c b/src/Internal/C/lapack-aux.c --- a/src/Internal/C/lapack-aux.c +++ b/src/Internal/C/lapack-aux.c @@ -45,6 +45,9 @@ for(q=0;q 0) { ret = CODE; } else if(RES < 0) { ret = RES; }) +#define UNWIND(RES,CODE,LABEL) MACRO(if(RES) { ret = CODE; goto LABEL; }) #define BAD_SIZE 2000 #define BAD_CODE 2001 @@ -116,6 +119,7 @@ integer *info); int svd_l_R(ODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) { + integer ret = 0; integer m = ar; integer n = ac; integer q = MIN(m,n); @@ -152,9 +156,12 @@ vp,&ldvt, &ans, &lwork, &res); + CHECK(res,res); + lwork = ceil(ans); double * work = (double*)malloc(lwork*sizeof(double)); CHECK(!work,MEM); + dgesvd_ (jobu,jobvt, &m,&n,ap,&m, sp, @@ -162,9 +169,10 @@ vp,&ldvt, work, &lwork, &res); - CHECK(res,res); + + MARK(res, res); free(work); - OK + return ret; } // (alternative version) @@ -175,9 +183,10 @@ integer *iwork, integer *info); int svd_l_Rdd(ODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) { - integer m = ar; - integer n = ac; - integer q = MIN(m,n); + integer ret = 0; + integer m = ar; + integer n = ac; + integer q = MIN(m,n); REQUIRES(sn==q,BAD_SIZE); REQUIRES((up == NULL && vp == NULL) || (ur==m && vc==n @@ -195,20 +204,27 @@ } DEBUGMSG("svd_l_Rdd"); integer* iwk = (integer*) malloc(8*q*sizeof(integer)); - CHECK(!iwk,MEM); + UNWIND(!iwk,MEM,cleanup0); integer lwk = -1; integer res; // ask for optimal lwk double ans; dgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res); + UNWIND(res,res,cleanup1); + lwk = ans; double * workv = (double*)malloc(lwk*sizeof(double)); - CHECK(!workv,MEM); + UNWIND(!workv,MEM,cleanup1); + dgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res); - CHECK(res,res); - free(iwk); + UNWIND(res,res,cleanup2); + +cleanup2: free(workv); - OK +cleanup1: + free(iwk); +cleanup0: + return ret; } //////////////////// complex svd //////////////////////////////////// @@ -219,11 +235,14 @@ integer *lwork, doublereal *rwork, integer *info); int svd_l_C(OCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) { - integer m = ar; - integer n = ac; - integer q = MIN(m,n); + integer ret = 0; + integer m = ar; + integer n = ac; + integer q = MIN(m,n); REQUIRES(sn==q,BAD_SIZE); REQUIRES(up==NULL || (ur==m && (uc==m || uc==q)),BAD_SIZE); + REQUIRES(vp==NULL || (vc==n && (vr==n || vr==q)),BAD_SIZE); + char* jobu = "A"; if (up==NULL) { jobu = "N"; @@ -232,7 +251,6 @@ jobu = "S"; } } - REQUIRES(vp==NULL || (vc==n && (vr==n || vr==q)),BAD_SIZE); char* jobvt = "A"; integer ldvt = n; if (vp==NULL) { @@ -245,7 +263,8 @@ }DEBUGMSG("svd_l_C"); double *rwork = (double*) malloc(5*q*sizeof(double)); - CHECK(!rwork,MEM); + UNWIND(!rwork,MEM,cleanup0); + integer lwork = -1; integer res; // ask for optimal lwork @@ -258,9 +277,12 @@ &ans, &lwork, rwork, &res); + UNWIND(res,res,cleanup1); + lwork = ceil(ans.r); doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); - CHECK(!work,MEM); + UNWIND(!work,MEM,cleanup1); + zgesvd_ (jobu,jobvt, &m,&n,ap,&m, sp, @@ -269,10 +291,14 @@ work, &lwork, rwork, &res); - CHECK(res,res); + UNWIND(res,res,cleanup2); + +cleanup2: free(work); +cleanup1: free(rwork); - OK +cleanup0: + return ret; } int zgesdd_ (char *jobz, integer *m, integer *n, @@ -281,49 +307,68 @@ integer *lwork, doublereal *rwork, integer* iwork, integer *info); int svd_l_Cdd(OCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) { - integer m = ar; - integer n = ac; - integer q = MIN(m,n); - REQUIRES(sn==q,BAD_SIZE); + integer ret = 0; + integer m = ar; + integer n = ac; + integer mx = MAX(m,n); + integer mn = MIN(m,n); + REQUIRES(sn==mn,BAD_SIZE); REQUIRES((up == NULL && vp == NULL) || (ur==m && vc==n - && ((uc == q && vr == q) + && ((uc == mn && vr == mn) || (uc == m && vc==n))),BAD_SIZE); char* jobz = "A"; integer ldvt = n; if (up==NULL) { jobz = "N"; } else { - if (uc==q && vr == q) { + if (uc==mn && vr == mn) { jobz = "S"; - ldvt = q; + ldvt = mn; } } DEBUGMSG("svd_l_Cdd"); - integer* iwk = (integer*) malloc(8*q*sizeof(integer)); - CHECK(!iwk,MEM); + integer* iwk = (integer*) malloc(8*mn*sizeof(integer)); + UNWIND(!iwk,MEM,cleanup0); + + // Docs: http://www.netlib.org/lapack/explore-html/d8/d54/zgesdd_8f_source.html + // RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) + // Let mx = max(M,N) and mn = min(M,N). + // If JOBZ = 'N', LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn); + // else if mx >> mn, LRWORK >= 5*mn*mn + 5*mn; + // else LRWORK >= max( 5*mn*mn + 5*mn, + // 2*mx*mn + 2*mn*mn + mn ). int lrwk; - if (0 && *jobz == 'N') { - lrwk = 5*q; // does not work, crash at free below + if (*jobz == 'N') { + lrwk = 7*mn; } else { - lrwk = 5*q*q + 7*q; + lrwk = MAX(5*mn*mn + 7*mn, 2*mx*mn + 2*mn*mn + mn); } - double *rwk = (double*)malloc(lrwk*sizeof(double));; - CHECK(!rwk,MEM); + double *rwk = (double*)malloc(MAX(1, lrwk)*sizeof(double));; + UNWIND(!rwk,MEM,cleanup1); + integer lwk = -1; integer res; // ask for optimal lwk doublecomplex ans; zgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res); + UNWIND(res,res,cleanup2); + lwk = ans.r; doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex)); - CHECK(!workv,MEM); + UNWIND(!workv,MEM,cleanup2); + zgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res); - CHECK(res,res); + UNWIND(res,res,cleanup3); + +cleanup3: free(workv); +cleanup2: free(rwk); +cleanup1: free(iwk); - OK +cleanup0: + return ret; } //////////////////// general complex eigensystem //////////// @@ -334,15 +379,18 @@ integer *lwork, doublereal *rwork, integer *info); int eig_l_C(OCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) { - integer n = ar; + integer ret = 0; + integer n = ar; REQUIRES(ac==n && sn==n, BAD_SIZE); REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE); char jobvl = up==NULL?'N':'V'; REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE); char jobvr = vp==NULL?'N':'V'; DEBUGMSG("eig_l_C"); + double *rwork = (double*) malloc(2*n*sizeof(double)); - CHECK(!rwork,MEM); + UNWIND(!rwork,MEM,cleanup0); + integer lwork = -1; integer res; // ask for optimal lwork @@ -355,9 +403,13 @@ &ans, &lwork, rwork, &res); + + UNWIND(res,res,cleanup1); + lwork = ceil(ans.r); doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); - CHECK(!work,MEM); + UNWIND(!work,MEM,cleanup1); + zgeev_ (&jobvl,&jobvr, &n,ap,&n, sp, @@ -366,10 +418,15 @@ work, &lwork, rwork, &res); - CHECK(res,res); + + UNWIND(res,res,cleanup2); + +cleanup2: free(work); +cleanup1: free(rwork); - OK +cleanup0: + return ret; } @@ -382,7 +439,8 @@ integer *lwork, integer *info); int eig_l_R(ODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) { - integer n = ar; + integer ret = 0; + integer n = ar; REQUIRES(ac==n && sn==n, BAD_SIZE); REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE); char jobvl = up==NULL?'N':'V'; @@ -400,6 +458,8 @@ vp,&n, &ans, &lwork, &res); + CHECK(res,res); + lwork = ceil(ans); double * work = (double*)malloc(lwork*sizeof(double)); CHECK(!work,MEM); @@ -410,9 +470,10 @@ vp,&n, work, &lwork, &res); - CHECK(res,res); + MARK(res,res); + free(work); - OK + return ret; } //////////////////// generalized real eigensystem //////////// @@ -425,7 +486,8 @@ integer *lwork, integer *info); int eig_l_G(ODMAT(a), ODMAT(b), CVEC(alpha), DVEC(beta), ODMAT(vl), ODMAT(vr)) { - integer n = ar; + integer ret = 0; + integer n = ar; REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE); REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE); char jobvl = vlp==NULL?'N':'V'; @@ -443,9 +505,12 @@ vlp, &n, vrp, &n, &ans, &lwork, &res); + CHECK(res,res); + lwork = ceil(ans); double * work = (double*)malloc(lwork*sizeof(double)); CHECK(!work,MEM); + dggev_ (&jobvl,&jobvr, &n, ap,&n,bp,&n, @@ -453,9 +518,10 @@ vlp, &n, vrp, &n, work, &lwork, &res); - CHECK(res,res); + MARK(res,res); + free(work); - OK + return ret; } //////////////////// generalized complex eigensystem //////////// @@ -468,7 +534,8 @@ doublereal *rwork, integer *info); int eig_l_GC(OCMAT(a), OCMAT(b), CVEC(alpha), CVEC(beta), OCMAT(vl), OCMAT(vr)) { - integer n = ar; + integer ret = 0; + integer n = ar; REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE); REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE); char jobvl = vlp==NULL?'N':'V'; @@ -476,7 +543,8 @@ char jobvr = vrp==NULL?'N':'V'; DEBUGMSG("eig_l_GC"); double *rwork = (double*) malloc(8*n*sizeof(double)); - CHECK(!rwork,MEM); + UNWIND(!rwork,MEM,cleanup0); + integer lwork = -1; integer res; // ask for optimal lwork @@ -488,9 +556,12 @@ vlp, &n, vrp, &n, &ans, &lwork, rwork, &res); + UNWIND(res,res,cleanup1); + lwork = ceil(ans.r); doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); - CHECK(!work,MEM); + UNWIND(!work,MEM,cleanup1); + zggev_ (&jobvl,&jobvr, &n, ap,&n,bp,&n, @@ -498,9 +569,14 @@ vlp, &n, vrp, &n, work, &lwork, rwork, &res); - CHECK(res,res); + UNWIND(res,res,cleanup2); + +cleanup2: free(work); - OK +cleanup1: + free(rwork); +cleanup0: + return ret; } //////////////////// symmetric real eigensystem //////////// @@ -510,7 +586,8 @@ integer *info); int eig_l_S(int wantV,DVEC(s),ODMAT(v)) { - integer n = sn; + integer ret = 0; + integer n = sn; REQUIRES(vr==n && vc==n, BAD_SIZE); char jobz = wantV?'V':'N'; DEBUGMSG("eig_l_S"); @@ -524,17 +601,21 @@ sp, &ans, &lwork, &res); + CHECK(res,res); + lwork = ceil(ans); double * work = (double*)malloc(lwork*sizeof(double)); CHECK(!work,MEM); + dsyev_ (&jobz,&uplo, &n,vp,&n, sp, work, &lwork, &res); - CHECK(res,res); + MARK(res,res); + free(work); - OK + return ret; } //////////////////// hermitian complex eigensystem //////////// @@ -544,12 +625,15 @@ doublereal *rwork, integer *info); int eig_l_H(int wantV,DVEC(s),OCMAT(v)) { - integer n = sn; + integer ret = 0; + integer n = sn; + REQUIRES(vr==n && vc==n, BAD_SIZE); char jobz = wantV?'V':'N'; DEBUGMSG("eig_l_H"); double *rwork = (double*) malloc((3*n-2)*sizeof(double)); - CHECK(!rwork,MEM); + UNWIND(!rwork,MEM,cleanup0); + integer lwork = -1; char uplo = 'U'; integer res; @@ -561,19 +645,26 @@ &ans, &lwork, rwork, &res); + UNWIND(res,res,cleanup1); + lwork = ceil(ans.r); doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); - CHECK(!work,MEM); + UNWIND(!work,MEM,cleanup1); + zheev_ (&jobz,&uplo, &n,vp,&n, sp, work, &lwork, rwork, &res); - CHECK(res,res); + UNWIND(res,res,cleanup2); + +cleanup2: free(work); +cleanup1: free(rwork); - OK +cleanup0: + return ret; } //////////////////// general real linear system //////////// @@ -582,23 +673,25 @@ *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info); int linearSolveR_l(ODMAT(a),ODMAT(b)) { - integer n = ar; + integer ret = 0; + integer n = ar; integer nhrs = bc; + REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); DEBUGMSG("linearSolveR_l"); integer * ipiv = (integer*)malloc(n*sizeof(integer)); + CHECK(!ipiv,MEM); + integer res; dgesv_ (&n,&nhrs, ap, &n, ipiv, bp, &n, &res); - if(res>0) { - return SINGULAR; - } - CHECK(res,res); + CONVERGED(res,SINGULAR); + free(ipiv); - OK + return ret; } //////////////////// general complex linear system //////////// @@ -608,23 +701,25 @@ info); int linearSolveC_l(OCMAT(a),OCMAT(b)) { - integer n = ar; + integer ret = 0; + integer n = ar; integer nhrs = bc; + REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); DEBUGMSG("linearSolveC_l"); integer * ipiv = (integer*)malloc(n*sizeof(integer)); + CHECK(!ipiv,MEM); + integer res; zgesv_ (&n,&nhrs, ap, &n, ipiv, bp, &n, &res); - if(res>0) { - return SINGULAR; - } - CHECK(res,res); + CONVERGED(res,SINGULAR); + free(ipiv); - OK + return ret; } //////// symmetric positive definite real linear system using Cholesky //////////// @@ -768,27 +863,37 @@ integer *info); int triDiagSolveR_l(DVEC(dl), DVEC(d), DVEC(du), ODMAT(b)) { - integer n = dn; + integer ret = 0; + integer n = dn; integer nhrs = bc; REQUIRES(n >= 1 && dln == dn - 1 && dun == dn - 1 && br == n, BAD_SIZE); DEBUGMSG("triDiagSolveR_l"); integer res; integer* ipiv = (integer*)malloc(n*sizeof(integer)); - double* du2 = (double*)malloc((n - 2)*sizeof(double)); + UNWIND(!ipiv,MEM,cleanup0); + + double* du2 = (double*)malloc((n - 2)*sizeof(double)); + UNWIND(!du2,MEM,cleanup1); + dgttrf_ (&n, dlp, dp, dup, du2, ipiv, &res); - CHECK(res,res); + UNWIND(res,res,cleanup2); + dgttrs_ ("N", &n,&nhrs, dlp, dp, dup, du2, ipiv, bp, &n, &res); - CHECK(res,res); - free(ipiv); + UNWIND(res,res,cleanup2); + +cleanup2: free(du2); - OK +cleanup1: + free(ipiv); +cleanup0: + return ret; } //////// tridiagonal complex linear system //////////// @@ -804,27 +909,37 @@ integer *info); int triDiagSolveC_l(CVEC(dl), CVEC(d), CVEC(du), OCMAT(b)) { - integer n = dn; + integer ret = 0; + integer n = dn; integer nhrs = bc; REQUIRES(n >= 1 && dln == dn - 1 && dun == dn - 1 && br == n, BAD_SIZE); DEBUGMSG("triDiagSolveC_l"); integer res; integer* ipiv = (integer*)malloc(n*sizeof(integer)); + UNWIND(!ipiv,MEM,cleanup0); + doublecomplex* du2 = (doublecomplex*)malloc((n - 2)*sizeof(doublecomplex)); + UNWIND(!du2,MEM,cleanup1); + zgttrf_ (&n, dlp, dp, dup, du2, ipiv, &res); - CHECK(res,res); + UNWIND(res,res,cleanup2); + zgttrs_ ("N", &n,&nhrs, dlp, dp, dup, du2, ipiv, bp, &n, &res); - CHECK(res,res); - free(ipiv); + UNWIND(res,res,cleanup2); + +cleanup2: free(du2); - OK +cleanup1: + free(ipiv); +cleanup0: + return ret; } //////////////////// least squares real linear system //////////// @@ -834,10 +949,11 @@ doublereal *work, integer *lwork, integer *info); int linearSolveLSR_l(ODMAT(a),ODMAT(b)) { - integer m = ar; - integer n = ac; + integer ret = 0; + integer m = ar; + integer n = ac; integer nrhs = bc; - integer ldb = bXc; + integer ldb = bXc; REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE); DEBUGMSG("linearSolveLSR_l"); integer res; @@ -848,19 +964,21 @@ bp,&ldb, &ans,&lwork, &res); + CHECK(res,res); + lwork = ceil(ans); double * work = (double*)malloc(lwork*sizeof(double)); + CHECK(!work,MEM); + dgels_ ("N",&m,&n,&nrhs, ap,&m, bp,&ldb, work,&lwork, &res); - if(res>0) { - return SINGULAR; - } - CHECK(res,res); + CONVERGED(res,SINGULAR); + free(work); - OK + return ret; } //////////////////// least squares complex linear system //////////// @@ -870,10 +988,11 @@ doublecomplex *work, integer *lwork, integer *info); int linearSolveLSC_l(OCMAT(a),OCMAT(b)) { - integer m = ar; - integer n = ac; + integer ret = 0; + integer m = ar; + integer n = ac; integer nrhs = bc; - integer ldb = bXc; + integer ldb = bXc; REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE); DEBUGMSG("linearSolveLSC_l"); integer res; @@ -884,19 +1003,21 @@ bp,&ldb, &ans,&lwork, &res); + CHECK(res,res); + lwork = ceil(ans.r); doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); + CHECK(!work,MEM); + zgels_ ("N",&m,&n,&nrhs, ap,&m, bp,&ldb, work,&lwork, &res); - if(res>0) { - return SINGULAR; - } - CHECK(res,res); + CONVERGED(res,SINGULAR); + free(work); - OK + return ret; } //////////////////// least squares real linear system using SVD //////////// @@ -907,13 +1028,17 @@ integer *info); int linearSolveSVDR_l(double rcond,ODMAT(a),ODMAT(b)) { - integer m = ar; - integer n = ac; + integer ret = 0; + integer m = ar; + integer n = ac; integer nrhs = bc; - integer ldb = bXc; + integer ldb = bXc; REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE); DEBUGMSG("linearSolveSVDR_l"); - double*S = (double*)malloc(MIN(m,n)*sizeof(double)); + + double * S = (double*)malloc(MIN(m,n)*sizeof(double)); + UNWIND(!S,MEM,cleanup0); + integer res; integer lwork = -1; integer rank; @@ -925,8 +1050,12 @@ &rcond,&rank, &ans,&lwork, &res); + UNWIND(res,res,cleanup1); + lwork = ceil(ans); double * work = (double*)malloc(lwork*sizeof(double)); + UNWIND(!work,MEM,cleanup1); + dgelss_ (&m,&n,&nrhs, ap,&m, bp,&ldb, @@ -934,13 +1063,15 @@ &rcond,&rank, work,&lwork, &res); - if(res>0) { - return NOCONVER; - } - CHECK(res,res); + + CONVERGED(res,NOCONVER); + free(work); +cleanup1: free(S); - OK +cleanup0: + return ret; + } //////////////////// least squares complex linear system using SVD //////////// @@ -952,14 +1083,20 @@ integer *info); int linearSolveSVDC_l(double rcond, OCMAT(a),OCMAT(b)) { - integer m = ar; - integer n = ac; + integer ret = 0; + integer m = ar; + integer n = ac; integer nrhs = bc; - integer ldb = bXc; + integer ldb = bXc; REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE); DEBUGMSG("linearSolveSVDC_l"); + double*S = (double*)malloc(MIN(m,n)*sizeof(double)); + UNWIND(!S,MEM,cleanup0); + double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double)); + UNWIND(!S,MEM,cleanup1); + integer res; integer lwork = -1; integer rank; @@ -972,8 +1109,12 @@ &ans,&lwork, RWORK, &res); + UNWIND(res,res,cleanup2); + lwork = ceil(ans.r); doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); + UNWIND(!work,MEM,cleanup2); + zgelss_ (&m,&n,&nrhs, ap,&m, bp,&ldb, @@ -982,14 +1123,16 @@ work,&lwork, RWORK, &res); - if(res>0) { - return NOCONVER; - } - CHECK(res,res); + CONVERGED(res,NOCONVER); + free(work); +cleanup2: free(RWORK); +cleanup1: free(S); - OK +cleanup0: + return ret; + } //////////////////// Cholesky factorization ///////////////////////// @@ -1042,36 +1185,43 @@ lda, doublereal *tau, doublereal *work, integer *info); int qr_l_R(DVEC(tau), ODMAT(r)) { - integer m = rr; - integer n = rc; - integer mn = MIN(m,n); + integer ret = 0; + integer m = rr; + integer n = rc; + integer mn = MIN(m,n); REQUIRES(m>=1 && n >=1 && taun == mn, BAD_SIZE); DEBUGMSG("qr_l_R"); double *WORK = (double*)malloc(n*sizeof(double)); CHECK(!WORK,MEM); + integer res; dgeqr2_ (&m,&n,rp,&m,taup,WORK,&res); - CHECK(res,res); + MARK(res,res); + free(WORK); - OK + return ret; } int zgeqr2_(integer *m, integer *n, doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *work, integer *info); int qr_l_C(CVEC(tau), OCMAT(r)) { - integer m = rr; - integer n = rc; - integer mn = MIN(m,n); + integer ret = 0; + integer m = rr; + integer n = rc; + integer mn = MIN(m,n); REQUIRES(m>=1 && n >=1 && taun == mn, BAD_SIZE); DEBUGMSG("qr_l_C"); + doublecomplex *WORK = (doublecomplex*)malloc(n*sizeof(doublecomplex)); CHECK(!WORK,MEM); + integer res; zgeqr2_ (&m,&n,rp,&m,taup,WORK,&res); - CHECK(res,res); + MARK(res,res); + free(WORK); - OK + return ret; } int dorgqr_(integer *m, integer *n, integer *k, doublereal * @@ -1079,18 +1229,21 @@ integer *info); int c_dorgqr(KDVEC(tau), ODMAT(r)) { - integer m = rr; - integer n = MIN(rc,rr); - integer k = taun; + integer ret = 0; + integer m = rr; + integer n = MIN(rc,rr); + integer k = taun; DEBUGMSG("c_dorgqr"); integer lwork = 8*n; // FIXME double *WORK = (double*)malloc(lwork*sizeof(double)); CHECK(!WORK,MEM); + integer res; dorgqr_ (&m,&n,&k,rp,&m,(double*)taup,WORK,&lwork,&res); - CHECK(res,res); + MARK(res,res); + free(WORK); - OK + return ret; } int zungqr_(integer *m, integer *n, integer *k, @@ -1098,18 +1251,21 @@ work, integer *lwork, integer *info); int c_zungqr(KCVEC(tau), OCMAT(r)) { - integer m = rr; - integer n = MIN(rc,rr); - integer k = taun; + integer ret = 0; + integer m = rr; + integer n = MIN(rc,rr); + integer k = taun; DEBUGMSG("z_ungqr"); integer lwork = 8*n; // FIXME doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); CHECK(!WORK,MEM); + integer res; zungqr_ (&m,&n,&k,rp,&m,(doublecomplex*)taup,WORK,&lwork,&res); - CHECK(res,res); + MARK(res,res); + free(WORK); - OK + return ret; } @@ -1120,20 +1276,23 @@ integer *lwork, integer *info); int hess_l_R(DVEC(tau), ODMAT(r)) { - integer m = rr; - integer n = rc; - integer mn = MIN(m,n); + integer ret = 0; + integer m = rr; + integer n = rc; + integer mn = MIN(m,n); REQUIRES(m>=1 && n == m && taun == mn-1, BAD_SIZE); DEBUGMSG("hess_l_R"); integer lwork = 5*n; // FIXME double *WORK = (double*)malloc(lwork*sizeof(double)); CHECK(!WORK,MEM); + integer res; integer one = 1; dgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res); - CHECK(res,res); + MARK(res,res); + free(WORK); - OK + return ret; } @@ -1142,20 +1301,23 @@ work, integer *lwork, integer *info); int hess_l_C(CVEC(tau), OCMAT(r)) { - integer m = rr; - integer n = rc; - integer mn = MIN(m,n); + integer ret = 0; + integer m = rr; + integer n = rc; + integer mn = MIN(m,n); REQUIRES(m>=1 && n == m && taun == mn-1, BAD_SIZE); DEBUGMSG("hess_l_C"); integer lwork = 5*n; // FIXME doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); CHECK(!WORK,MEM); + integer res; integer one = 1; zgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res); - CHECK(res,res); + MARK(res,res); + free(WORK); - OK + return ret; } //////////////////// Schur factorization ///////////////////////// @@ -1166,28 +1328,35 @@ integer *lwork, logical *bwork, integer *info); int schur_l_R(ODMAT(u), ODMAT(s)) { - integer m = sr; - integer n = sc; + integer ret = 0; + integer m = sr; + integer n = sc; REQUIRES(m>=1 && n==m && ur==n && uc==n, BAD_SIZE); DEBUGMSG("schur_l_R"); integer lwork = 6*n; // FIXME double *WORK = (double*)malloc(lwork*sizeof(double)); - double *WR = (double*)malloc(n*sizeof(double)); - double *WI = (double*)malloc(n*sizeof(double)); + UNWIND(!WORK,MEM,cleanup0); + double *WR = (double*)malloc(n*sizeof(double)); + UNWIND(!WORK,MEM,cleanup1); + double *WI = (double*)malloc(n*sizeof(double)); + UNWIND(!WORK,MEM,cleanup2); // WR and WI not really required in this call logical *BWORK = (logical*)malloc(n*sizeof(logical)); + UNWIND(!BWORK,MEM,cleanup3); integer res; integer sdim; dgees_ ("V","N",NULL,&n,sp,&n,&sdim,WR,WI,up,&n,WORK,&lwork,BWORK,&res); - if(res>0) { - return NOCONVER; - } - CHECK(res,res); - free(WR); - free(WI); + CONVERGED(res,NOCONVER); + free(BWORK); +cleanup3: + free(WI); +cleanup2: + free(WR); +cleanup1: free(WORK); - OK +cleanup0: + return ret; } @@ -1197,29 +1366,40 @@ doublereal *rwork, logical *bwork, integer *info); int schur_l_C(OCMAT(u), OCMAT(s)) { - integer m = sr; - integer n = sc; + integer ret = 0; + integer m = sr; + integer n = sc; REQUIRES(m>=1 && n==m && ur==n && uc==n, BAD_SIZE); DEBUGMSG("schur_l_C"); integer lwork = 6*n; // FIXME doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); - doublecomplex *W = (doublecomplex*)malloc(n*sizeof(doublecomplex)); + UNWIND(!WORK,MEM,cleanup0); + + doublecomplex *W = (doublecomplex*)malloc(n*sizeof(doublecomplex)); + UNWIND(!W,MEM,cleanup1); + // W not really required in this call logical *BWORK = (logical*)malloc(n*sizeof(logical)); - double *RWORK = (double*)malloc(n*sizeof(double)); + UNWIND(!BWORK,MEM,cleanup2); + + double *RWORK = (double*)malloc(n*sizeof(double)); + UNWIND(!RWORK,MEM,cleanup3); 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); + CONVERGED(res,NOCONVER); + + free(RWORK); +cleanup3: free(BWORK); +cleanup2: + free(W); +cleanup1: free(WORK); - OK +cleanup0: + return ret; } //////////////////// LU factorization ///////////////////////// @@ -1228,24 +1408,30 @@ lda, integer *ipiv, integer *info); int lu_l_R(DVEC(ipiv), ODMAT(r)) { - integer m = rr; - integer n = rc; - integer mn = MIN(m,n); + integer ret = 0; + integer m = rr; + integer n = rc; + 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)); + UNWIND(!auxipiv,MEM,cleanup0); + 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)); + UNWIND(!auxipiv,MEM,cleanup0); + integer res; zgetrf_ (&m,&n,rp,&m,auxipiv,&res); if(res>0) { res = 0; // FIXME } - CHECK(res,res); - int k; - for (k=0; k=1 && rc==n && ipivn == n, BAD_SIZE); DEBUGMSG("ldl_R"); + integer* auxipiv = (integer*)malloc(n*sizeof(integer)); + UNWIND(!auxipiv,MEM,cleanup0); + integer res; integer lda = rXc; integer lwork = -1; @@ -1343,15 +1547,22 @@ dsytrf_ ("L",&n,rp,&lda,auxipiv,&ans,&lwork,&res); lwork = ceil(ans); doublereal* work = (doublereal*)malloc(lwork*sizeof(doublereal)); + UNWIND(!work,MEM,cleanup1); + dsytrf_ ("L",&n,rp,&lda,auxipiv,work,&lwork,&res); - CHECK(res,res); + UNWIND(res,res,cleanup2); + int k; for (k=0; k=1 && rc==n && ipivn == n, BAD_SIZE); DEBUGMSG("ldl_R"); integer* auxipiv = (integer*)malloc(n*sizeof(integer)); + UNWIND(!auxipiv,MEM,cleanup0); + integer res; integer lda = rXc; integer lwork = -1; @@ -1370,15 +1585,21 @@ zhetrf_ ("L",&n,rp,&lda,auxipiv,&ans,&lwork,&res); lwork = ceil(ans.r); doublecomplex* work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); + UNWIND(!work,MEM,cleanup1); + zhetrf_ ("L",&n,rp,&lda,auxipiv,work,&lwork,&res); - CHECK(res,res); + UNWIND(res,res,cleanup2); int k; for (k=0; k Build (Int,Int) (e -> e -> e) Matrix e +instance (Num e, Container Vector e) => Build (Int,Int) (e -> e -> e) Matrix e where build = build' diff --git a/src/Internal/Conversion.hs b/src/Internal/Conversion.hs --- a/src/Internal/Conversion.hs +++ b/src/Internal/Conversion.hs @@ -3,7 +3,6 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-} -{-# LANGUAGE UndecidableInstances #-} ----------------------------------------------------------------------------- -- | diff --git a/src/Internal/Devel.hs b/src/Internal/Devel.hs --- a/src/Internal/Devel.hs +++ b/src/Internal/Devel.hs @@ -16,7 +16,8 @@ import Foreign.C.Types ( CInt ) --import Foreign.Storable.Complex () import Foreign.Ptr(Ptr) -import Control.Exception as E ( SomeException, catch ) +import Control.Exception (SomeException, SomeAsyncException (..)) +import qualified Control.Exception as Exception import Internal.Vector(Vector,avec) import Foreign.Storable(Storable) @@ -59,9 +60,21 @@ -- | Error capture and conversion to Maybe mbCatch :: IO x -> IO (Maybe x) -mbCatch act = E.catch (Just `fmap` act) f - where f :: SomeException -> IO (Maybe x) - f _ = return Nothing +mbCatch act = + hush <$> + Exception.tryJust + (\e -> if isSyncException e then Just e else Nothing) + act + + where + hush :: Either a b -> Maybe b + hush = either (const Nothing) Just + + isSyncException :: SomeException -> Bool + isSyncException e = + case Exception.fromException e of + Just (SomeAsyncException _) -> False + Nothing -> True -------------------------------------------------------------------------------- diff --git a/src/Internal/Element.hs b/src/Internal/Element.hs --- a/src/Internal/Element.hs +++ b/src/Internal/Element.hs @@ -1,7 +1,6 @@ {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# OPTIONS_GHC -fno-warn-orphans #-} @@ -40,7 +39,7 @@ import Data.Binary -instance (Binary (Vector a), Element a) => Binary (Matrix a) where +instance (Binary a, Element a) => Binary (Matrix a) where put m = do put (cols m) put (flatten m) diff --git a/src/Internal/Modular.hs b/src/Internal/Modular.hs --- a/src/Internal/Modular.hs +++ b/src/Internal/Modular.hs @@ -1,20 +1,14 @@ {-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} -{-# LANGUAGE KindSignatures #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE Rank2Types #-} {-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE GADTs #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE TypeOperators #-} {-# OPTIONS_GHC -fno-warn-missing-signatures #-} -{-# OPTIONS_GHC -fno-warn-missing-methods #-} {- | Module : Internal.Modular @@ -83,11 +77,11 @@ where compare a b = compare (unMod a) (unMod b) -instance (Integral t, Real t, KnownNat m, Integral (Mod m t)) => Real (Mod m t) +instance (Integral t, Real t, KnownNat m) => Real (Mod m t) where toRational x = toInteger x % 1 -instance (Integral t, KnownNat m, Num (Mod m t)) => Integral (Mod m t) +instance (Integral t, KnownNat m) => Integral (Mod m t) where toInteger = toInteger . unMod quotRem a b = (Mod q, Mod r) @@ -95,7 +89,7 @@ (q,r) = quotRem (unMod a) (unMod b) -- | this instance is only valid for prime m -instance (Show (Mod m t), Num (Mod m t), Eq t, KnownNat m) => Fractional (Mod m t) +instance (Integral t, Show t, Eq t, KnownNat m) => Fractional (Mod m t) where recip x | x*r == 1 = r @@ -151,6 +145,7 @@ gemm u a b c = gemmg (c_gemmMI m') (f2i u) (f2iM a) (f2iM b) (f2iM c) where m' = fromIntegral . natVal $ (undefined :: Proxy m) + reorderV strides dims = i2f . reorderAux c_reorderI strides dims . f2i instance KnownNat m => Element (Mod m Z) where @@ -168,6 +163,7 @@ gemm u a b c = gemmg (c_gemmML m') (f2i u) (f2iM a) (f2iM b) (f2iM c) where m' = fromIntegral . natVal $ (undefined :: Proxy m) + reorderV strides dims = i2f . reorderAux c_reorderL strides dims . f2i instance KnownNat m => CTrans (Mod m I) diff --git a/src/Internal/Sparse.hs b/src/Internal/Sparse.hs --- a/src/Internal/Sparse.hs +++ b/src/Internal/Sparse.hs @@ -1,11 +1,11 @@ -{-# LANGUAGE RecordWildCards #-} -{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE BangPatterns #-} {-# LANGUAGE FlexibleInstances #-} - -{-# OPTIONS_GHC -fno-warn-missing-signatures #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE RecordWildCards #-} module Internal.Sparse( - GMatrix(..), CSR(..), mkCSR, fromCSR, + GMatrix(..), CSR(..), mkCSR, fromCSR, impureCSR, mkSparse, mkDiagR, mkDense, AssocMatrix, toDense, @@ -16,10 +16,12 @@ import Internal.Matrix import Internal.Numeric import qualified Data.Vector.Storable as V -import Data.Function(on) +import qualified Data.Vector.Storable.Mutable as M import Control.Arrow((***)) -import Control.Monad(when) -import Data.List(groupBy, sort) +import Control.Monad(when, foldM) +import Control.Monad.ST (runST) +import Control.Monad.Primitive (PrimMonad) +import Data.List(sort) import Foreign.C.Types(CInt(..)) import Internal.Devel @@ -27,10 +29,7 @@ import Foreign(Ptr) import Text.Printf(printf) -infixl 0 ~!~ -c ~!~ msg = when c (error msg) - -type AssocMatrix = [((Int,Int),Double)] +type AssocMatrix = [(IndexOf Matrix, Double)] data CSR = CSR { csrVals :: Vector Double @@ -49,23 +48,79 @@ } deriving Show +-- | Produce a CSR sparse matrix from a association matrix. mkCSR :: AssocMatrix -> CSR -mkCSR sm' = CSR{..} +mkCSR ms = + runST $ impureCSR runFold $ sort ms + where + runFold next initialise xtract as0 = do + i0 <- initialise + acc <- foldM next i0 as0 + xtract acc + +-- | Produce a CSR sparse matrix by applying a generic folding function. +-- +-- This allows one to build a CSR from an effectful streaming source +-- when combined with libraries like pipes, io-streams, or streaming. +-- +-- For example +-- +-- > impureCSR Pipes.Prelude.foldM :: PrimMonad m => Producer AssocEntry m () -> m CSR +-- > impureCSR Streaming.Prelude.foldM :: PrimMonad m => Stream (Of AssocEntry) m r -> m (Of CSR r) +-- +impureCSR + :: PrimMonad m + => (forall x . (x -> (IndexOf Matrix, Double) -> m x) -> m x -> (x -> m CSR) -> r) + -> r +impureCSR f = f next begin done where - sm = sort sm' - rws = map ((fromList *** fromList) - . unzip - . map ((succ.fi.snd) *** id) - ) - . groupBy ((==) `on` (fst.fst)) - $ sm - rszs = map (fi . dim . fst) rws - csrRows = fromList (scanl (+) 1 rszs) - csrVals = vjoin (map snd rws) - csrCols = vjoin (map fst rws) - csrNRows = dim csrRows - 1 - csrNCols = fromIntegral (V.maximum csrCols) + sfi = succ . fi + begin = do + mv <- M.unsafeNew 64 + mr <- M.unsafeNew 64 + mc <- M.unsafeNew 64 + return (mv, mr, mc, 0, 0, 0, -1) + next (!mv, !mr, !mc, !idxVC, !idxR, !maxC, !curRow) ((r,c),d) = do + when (r < curRow) $ + error (printf "impureCSR: row %i specified after %i" r curRow) + + let lenVC = M.length mv + lenR = M.length mr + maxC' = max maxC c + + (mv', mc') <- + if idxVC >= lenVC then do + mv' <- M.unsafeGrow mv lenVC + mc' <- M.unsafeGrow mc lenVC + return (mv', mc') + else + return (mv, mc) + + mr' <- + if idxR >= lenR - 1 then + M.unsafeGrow mr lenR + else + return mr + + M.unsafeWrite mc' idxVC (sfi c) + M.unsafeWrite mv' idxVC d + + idxR' <- + foldM + (\idxR' _ -> idxR' + 1 <$ M.unsafeWrite mr' idxR' (sfi idxVC)) + idxR [1 .. (r-curRow)] + + return (mv', mr', mc', idxVC + 1, idxR', maxC', r) + + done (!mv, !mr, !mc, !idxVC, !idxR, !maxC, !curR) = do + M.unsafeWrite mr idxR (sfi idxVC) + vv <- V.unsafeFreeze (M.unsafeTake idxVC mv) + vc <- V.unsafeFreeze (M.unsafeTake idxVC mc) + vr <- V.unsafeFreeze (M.unsafeTake (idxR + 1) mr) + return $ CSR vv vc vr (succ curR) (succ maxC) + + {- | General matrix with specialized internal representations for dense, sparse, diagonal, banded, and constant elements. @@ -129,6 +184,7 @@ nCols = csrNCols +mkDiagR :: Int -> Int -> Vector Double -> GMatrix mkDiagR r c v | dim v <= min r c = Diag{..} | otherwise = error $ printf "mkDiagR: incorrect sizes (%d,%d) [%d]" r c (dim v) @@ -144,13 +200,17 @@ gmXv :: GMatrix -> Vector Double -> Vector Double gmXv SparseR { gmCSR = CSR{..}, .. } v = unsafePerformIO $ do - dim v /= nCols ~!~ printf "gmXv (CSR): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v) + when (dim v /= nCols) $ + error (printf "gmXv (CSR): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)) + r <- createVector nRows (csrVals # csrCols # csrRows # v #! r) c_smXv #|"CSRXv" return r gmXv SparseC { gmCSC = CSC{..}, .. } v = unsafePerformIO $ do - dim v /= nCols ~!~ printf "gmXv (CSC): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v) + when (dim v /= nCols) $ + error (printf "gmXv (CSC): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)) + r <- createVector nRows (cscVals # cscRows # cscCols # v #! r) c_smTXv #|"CSCXv" return r diff --git a/src/Internal/Static.hs b/src/Internal/Static.hs --- a/src/Internal/Static.hs +++ b/src/Internal/Static.hs @@ -16,7 +16,6 @@ {-# LANGUAGE DeriveGeneric #-} {-# OPTIONS_GHC -fno-warn-missing-signatures #-} -{-# OPTIONS_GHC -fno-warn-simplifiable-class-constraints #-} {- | Module : Internal.Static @@ -131,9 +130,9 @@ where du = fromIntegral . natVal $ (undefined :: Proxy n) dv = fromIntegral . natVal $ (undefined :: Proxy m) - u' | du > 1 && LA.size u == 1 = LA.konst (u D.@> 0) du + u' | du /= 1 && LA.size u == 1 = LA.konst (u D.@> 0) du | otherwise = u - v' | dv > 1 && LA.size v == 1 = LA.konst (v D.@> 0) dv + v' | dv /= 1 && LA.size v == 1 = LA.konst (v D.@> 0) dv | otherwise = v @@ -367,7 +366,7 @@ negate = lift1F negate fromInteger x = Dim (fromInteger x) -instance (Num (Vector t), Num (Matrix t), Fractional t, Numeric t) => Fractional (Dim n (Vector t)) +instance (Num (Vector t), Fractional t, Numeric t) => Fractional (Dim n (Vector t)) where fromRational x = Dim (fromRational x) (/) = lift2F (/) @@ -392,7 +391,7 @@ pi = Dim pi -instance (Num (Matrix t), Numeric t) => Num (Dim m (Dim n (Matrix t))) +instance (Num (Vector t), Numeric t) => Num (Dim m (Dim n (Matrix t))) where (+) = (lift2F . lift2F) (+) (*) = (lift2F . lift2F) (*) @@ -402,12 +401,12 @@ negate = (lift1F . lift1F) negate fromInteger x = Dim (Dim (fromInteger x)) -instance (Num (Vector t), Num (Matrix t), Fractional t, Numeric t) => Fractional (Dim m (Dim n (Matrix t))) +instance (Num (Vector t), Fractional t, Numeric t) => Fractional (Dim m (Dim n (Matrix t))) where fromRational x = Dim (Dim (fromRational x)) (/) = (lift2F.lift2F) (/) -instance (Num (Vector t), Floating (Matrix t), Fractional t, Numeric t) => Floating (Dim m (Dim n (Matrix t))) where +instance (Floating (Vector t), Floating t, Numeric t) => Floating (Dim m (Dim n (Matrix t))) where sin = (lift1F . lift1F) sin cos = (lift1F . lift1F) cos tan = (lift1F . lift1F) tan diff --git a/src/Internal/Util.hs b/src/Internal/Util.hs --- a/src/Internal/Util.hs +++ b/src/Internal/Util.hs @@ -417,7 +417,7 @@ instance Element t => Indexable (Matrix t) (Vector t) where - m!j = subVector (j*c) c (flatten m) + m ! j = subVector (j*c) c (flatten m) where c = cols m @@ -912,4 +912,3 @@ , remap r (tr c) p == ep , tr p ?? (PosCyc (idxs[-5,13]), Pos (idxs[3,7,1])) == (2><3) [35,75,15,33,73,13] ] - diff --git a/src/Numeric/LinearAlgebra/Devel.hs b/src/Numeric/LinearAlgebra/Devel.hs --- a/src/Numeric/LinearAlgebra/Devel.hs +++ b/src/Numeric/LinearAlgebra/Devel.hs @@ -14,7 +14,7 @@ module Numeric.LinearAlgebra.Devel( -- * FFI tools -- | See @examples/devel@ in the repository. - + createVector, createMatrix, TransArray(..), MatrixOrder(..), orderOf, cmat, fmat, @@ -27,7 +27,7 @@ -- * ST -- | In-place manipulation inside the ST monad. -- See @examples/inplace.hs@ in the repository. - + -- ** Mutable Vectors STVector, newVector, thawVector, freezeVector, runSTVector, readVector, writeVector, modifyVector, liftSTVector, @@ -51,7 +51,7 @@ liftMatrix, liftMatrix2, liftMatrix2Auto, -- * Sparse representation - CSR(..), fromCSR, mkCSR, + CSR(..), fromCSR, mkCSR, impureCSR, GMatrix(..), -- * Misc diff --git a/src/Numeric/LinearAlgebra/Static.hs b/src/Numeric/LinearAlgebra/Static.hs --- a/src/Numeric/LinearAlgebra/Static.hs +++ b/src/Numeric/LinearAlgebra/Static.hs @@ -46,6 +46,13 @@ matrix, -- * Complex ℂ, C, M, Her, her, 𝑖, + toComplex, + fromComplex, + complex, + real, + imag, + sqMagnitude, + magnitude, -- * Products (<>),(#>),(<.>), -- * Linear Systems @@ -74,7 +81,9 @@ (<\>),fromList,takeDiag,svd,eig,eigSH, eigenvalues,eigenvaluesSH,build, qr,size,dot,chol,range,R,C,sym,mTm,unSym, - randomVector,rand,randn,gaussianSample,uniformSample,meanCov) + randomVector,rand,randn,gaussianSample,uniformSample,meanCov, + toComplex, fromComplex, complex, real, magnitude + ) import qualified Numeric.LinearAlgebra as LA import qualified Numeric.LinearAlgebra.Devel as LA import Data.Proxy(Proxy(..)) @@ -243,7 +252,30 @@ where takeDiag x = mkC (LA.takeDiag (extract x)) +-------------------------------------------------------------------------------- + +toComplex :: KnownNat n => (R n, R n) -> C n +toComplex (r,i) = mkC $ LA.toComplex (ud1 r, ud1 i) + +fromComplex :: KnownNat n => C n -> (R n, R n) +fromComplex (C (Dim v)) = let (r,i) = LA.fromComplex v in (mkR r, mkR i) + +complex :: KnownNat n => R n -> C n +complex r = mkC $ LA.toComplex (ud1 r, LA.konst 0 (size r)) + +real :: KnownNat n => C n -> R n +real = fst . fromComplex + +imag :: KnownNat n => C n -> R n +imag = snd . fromComplex + +sqMagnitude :: KnownNat n => C n -> R n +sqMagnitude c = let (r,i) = fromComplex c in r**2 + i**2 + +magnitude :: KnownNat n => C n -> R n +magnitude = sqrt . sqMagnitude + -------------------------------------------------------------------------------- linSolve :: (KnownNat m, KnownNat n) => L m m -> L m n -> Maybe (L m n) @@ -595,9 +627,9 @@ a' = subVector 0 n a b' = subVector 0 n b -mulR (isDiag -> Just (0,a,_)) (extract -> b) = mkL (asColumn a * takeRows (LA.size a) b) +-- mulR (isDiag -> Just (0,a,_)) (extract -> b) = mkL (asColumn a * takeRows (LA.size a) b) -mulR (extract -> a) (isDiag -> Just (0,b,_)) = mkL (takeColumns (LA.size b) a * asRow b) +-- mulR (extract -> a) (isDiag -> Just (0,b,_)) = mkL (takeColumns (LA.size b) a * asRow b) mulR a b = mkL (extract a LA.<> extract b) @@ -658,9 +690,9 @@ a' = subVector 0 n a b' = subVector 0 n b -mulC (isDiagC -> Just (0,a,_)) (extract -> b) = mkM (asColumn a * takeRows (LA.size a) b) +-- mulC (isDiagC -> Just (0,a,_)) (extract -> b) = mkM (asColumn a * takeRows (LA.size a) b) -mulC (extract -> a) (isDiagC -> Just (0,b,_)) = mkM (takeColumns (LA.size b) a * asRow b) +-- mulC (extract -> a) (isDiagC -> Just (0,b,_)) = mkM (takeColumns (LA.size b) a * asRow b) mulC a b = mkM (extract a LA.<> extract b)