diff --git a/CHANGELOG b/CHANGELOG --- a/CHANGELOG +++ b/CHANGELOG @@ -2,15 +2,15 @@ -------- * Many new functions and instances in the Static module - + * meanCov and gaussianSample use Herm type * thinQR, thinRQ - + * compactSVDTol * unitary changed to normalize, also admits Vector (Complex Double) - + 0.17.0.0 -------- @@ -288,4 +288,3 @@ * added NFData instances for Matrix and Vector. * liftVector, liftVector2 replaced by mapVector, zipVector. - diff --git a/hmatrix.cabal b/hmatrix.cabal --- a/hmatrix.cabal +++ b/hmatrix.cabal @@ -1,11 +1,11 @@ Name: hmatrix -Version: 0.18.2.0 +Version: 0.20.2 License: BSD3 License-file: LICENSE Author: Alberto Ruiz -Maintainer: Alberto Ruiz +Maintainer: Dominic Steinitz Stability: provisional -Homepage: https://github.com/albertoruiz/hmatrix +Homepage: https://github.com/haskell-numerics/hmatrix Synopsis: Numeric Linear Algebra Description: Linear systems, matrix decompositions, and other numerical computations based on BLAS and LAPACK. . @@ -16,9 +16,9 @@ Code examples: Category: Math -tested-with: GHC==8.2, GHC==8.4 +tested-with: GHC==8.10 -cabal-version: >=1.8 +cabal-version: >=1.18 build-type: Simple @@ -29,15 +29,22 @@ 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 + default-language: Haskell2010 + Build-Depends: base >= 4.8 && < 5, binary, array, @@ -45,9 +52,10 @@ random, split, bytestring, + primitive, storable-complex, semigroups, - vector >= 0.8 + vector >= 0.11 hs-source-dirs: src @@ -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 @@ -149,5 +161,4 @@ source-repository head type: git - location: https://github.com/albertoruiz/hmatrix - + location: https://github.com/haskell-numerics/hmatrix diff --git a/src/Internal/Algorithms.hs b/src/Internal/Algorithms.hs --- a/src/Internal/Algorithms.hs +++ b/src/Internal/Algorithms.hs @@ -1,9 +1,10 @@ {-# LANGUAGE FlexibleContexts, FlexibleInstances #-} {-# LANGUAGE CPP #-} {-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE TypeFamilies #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-} + ----------------------------------------------------------------------------- {- | Module : Internal.Algorithms @@ -25,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 @@ -68,8 +73,10 @@ linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t linearSolveLS' :: Matrix t -> Matrix t -> Matrix t eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) + geig' :: Matrix t -> Matrix t -> (Vector (Complex Double), Vector t, Matrix (Complex Double)) eigSH'' :: Matrix t -> (Vector Double, Matrix t) eigOnly :: Matrix t -> Vector (Complex Double) + geigOnly :: Matrix t -> Matrix t -> (Vector (Complex Double), Vector t) eigOnlySH :: Matrix t -> Vector Double cholSH' :: Matrix t -> Matrix t mbCholSH' :: Matrix t -> Maybe (Matrix t) @@ -94,7 +101,9 @@ linearSolveSVD' = linearSolveSVDR Nothing eig' = eigR eigSH'' = eigS + geig' = eigG eigOnly = eigOnlyR + geigOnly = eigOnlyG eigOnlySH = eigOnlyS cholSH' = cholS mbCholSH' = mbCholS @@ -124,7 +133,9 @@ linearSolveLS' = linearSolveLSC linearSolveSVD' = linearSolveSVDC Nothing eig' = eigC + geig' = eigGC eigOnly = eigOnlyC + geigOnly = eigOnlyGC eigSH'' = eigH eigOnlySH = eigOnlyH cholSH' = cholH @@ -168,7 +179,8 @@ -0.690 -0.352 0.433 -0.233 0.398 >>> s -fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15] +[35.18264833189422,1.4769076999800903,1.089145439970417e-15] +it :: Vector Double >>> disp 3 v 3x3 @@ -222,7 +234,8 @@ -0.690 -0.352 0.433 >>> s -fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15] +[35.18264833189422,1.4769076999800903,1.089145439970417e-15] +it :: Vector Double >>> disp 3 v 3x3 @@ -281,7 +294,8 @@ -0.690 -0.352 >>> s -fromList [35.18264833189422,1.4769076999800903] +[35.18264833189422,1.476907699980091] +it :: Vector Double >>> disp 3 u 5x2 @@ -507,10 +521,25 @@ eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) eig = {-# SCC "eig" #-} eig' +-- | Generalized eigenvalues (not ordered) and eigenvectors (as columns) of a pair of nonsymmetric matrices. +-- Eigenvalues are represented as pairs of alpha, beta, where eigenvalue = alpha / beta. Alpha is always +-- complex, but betas has the same type as the input matrix. +-- +-- If @(alphas, betas, v) = geig a b@, then @a \<> v == b \<> v \<> diag (alphas / betas)@ +-- +-- Note that beta can be 0 and that has reasonable interpretation. +geig :: Field t => Matrix t -> Matrix t -> (Vector (Complex Double), Vector t, Matrix (Complex Double)) +geig = {-# SCC "geig" #-} geig' + -- | Eigenvalues (not ordered) of a general square matrix. eigenvalues :: Field t => Matrix t -> Vector (Complex Double) eigenvalues = {-# SCC "eigenvalues" #-} eigOnly +-- | Generalized eigenvalues of a pair of matrices. Represented as pairs of alpha, beta, +-- where eigenvalue is alpha / beta as in 'geig'. +geigenvalues :: Field t => Matrix t -> Matrix t -> (Vector (Complex Double), Vector t) +geigenvalues = {-# SCC "geigenvalues" #-} geigOnly + -- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part. eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t) eigSH' = {-# SCC "eigSH'" #-} eigSH'' @@ -533,7 +562,7 @@ >>> let (l, v) = eigSH a >>> l -fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575] +[11.344814282762075,0.17091518882717918,-0.5157294715892575] >>> disp 3 $ v <> diag l <> tr v 3x3 @@ -1133,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,12 +470,115 @@ vp,&n, work, &lwork, &res); + MARK(res,res); + + free(work); + return ret; +} + +//////////////////// generalized real eigensystem //////////// + +int dggev_(char *jobvl, char *jobvr, integer *n, + doublereal *a, integer *lda, doublereal *b, integer *ldb, + doublereal *alphar, doublereal *alphai, doublereal *beta, + doublereal *vl, integer *ldvl, doublereal *vr, integer *ldvr, + doublereal *work, + integer *lwork, integer *info); + +int eig_l_G(ODMAT(a), ODMAT(b), CVEC(alpha), DVEC(beta), ODMAT(vl), ODMAT(vr)) { + 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'; + REQUIRES(vrp==NULL || (vrr==n && vrc==n), BAD_SIZE); + char jobvr = vrp==NULL?'N':'V'; + DEBUGMSG("eig_l_G"); + integer lwork = -1; + integer res; + // ask for optimal lwork + double ans; + dggev_ (&jobvl,&jobvr, + &n, + ap,&n,bp,&n, + (double*)alphap, (double*)alphap+n, betap, + 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, + (double*)alphap, (double*)alphap+n, betap, + vlp, &n, vrp, &n, + work, &lwork, + &res); + MARK(res,res); + free(work); - OK + return ret; } +//////////////////// generalized complex eigensystem //////////// +int zggev_(char *jobvl, char *jobvr, integer *n, + doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, + doublecomplex *alphar, doublecomplex *beta, + doublecomplex *vl, integer *ldvl, doublecomplex *vr, integer *ldvr, + doublecomplex *work, integer *lwork, + doublereal *rwork, integer *info); + +int eig_l_GC(OCMAT(a), OCMAT(b), CVEC(alpha), CVEC(beta), OCMAT(vl), OCMAT(vr)) { + 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'; + REQUIRES(vrp==NULL || (vrr==n && vrc==n), BAD_SIZE); + char jobvr = vrp==NULL?'N':'V'; + DEBUGMSG("eig_l_GC"); + double *rwork = (double*) malloc(8*n*sizeof(double)); + UNWIND(!rwork,MEM,cleanup0); + + integer lwork = -1; + integer res; + // ask for optimal lwork + doublecomplex ans; + zggev_ (&jobvl,&jobvr, + &n, + ap,&n,bp,&n, + alphap, betap, + vlp, &n, vrp, &n, + &ans, &lwork, + rwork, &res); + UNWIND(res,res,cleanup1); + + lwork = ceil(ans.r); + doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); + UNWIND(!work,MEM,cleanup1); + + zggev_ (&jobvl,&jobvr, + &n, + ap,&n,bp,&n, + alphap, betap, + vlp, &n, vrp, &n, + work, &lwork, + rwork, &res); + UNWIND(res,res,cleanup2); + +cleanup2: + free(work); +cleanup1: + free(rwork); +cleanup0: + return ret; +} + //////////////////// symmetric real eigensystem //////////// int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a, @@ -423,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"); @@ -437,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 //////////// @@ -457,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; @@ -474,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 //////////// @@ -495,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 //////////// @@ -521,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 //////////// @@ -681,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 //////////// @@ -717,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 //////////// @@ -747,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; @@ -761,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 //////////// @@ -783,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; @@ -797,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 //////////// @@ -820,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; @@ -838,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, @@ -847,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 //////////// @@ -865,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; @@ -885,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, @@ -895,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 ///////////////////////// @@ -955,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 * @@ -992,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, @@ -1011,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; } @@ -1033,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; } @@ -1055,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 ///////////////////////// @@ -1079,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; } @@ -1110,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 ///////////////////////// @@ -1141,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; @@ -1256,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; @@ -1283,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>> linspace 5 (-3,7::Double) -fromList [-3.0,-0.5,2.0,4.5,7.0]@ +[-3.0,-0.5,2.0,4.5,7.0] +it :: Vector Double ->>> linspace 5 (8,2+i) :: Vector (Complex Double) -fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0] +>>> linspace 5 (8,3:+2) :: Vector (Complex Double) +[8.0 :+ 0.0,6.75 :+ 0.5,5.5 :+ 1.0,4.25 :+ 1.5,3.0 :+ 2.0] +it :: Vector (Complex Double) Logarithmic spacing can be defined as follows: @@ -85,7 +86,8 @@ >>> let v = vector [10,20,30] >>> m #> v -fromList [140.0,320.0] +[140.0,320.0] +it :: Vector Numeric.LinearAlgebra.Data.R -} infixr 8 #> @@ -134,10 +136,12 @@ >>> let x = a <\> v >>> x -fromList [3.0799999999999996,5.159999999999999] +[3.0799999999999996,5.159999999999999] +it :: Vector Numeric.LinearAlgebra.Data.R >>> a #> x -fromList [13.399999999999999,26.799999999999997,1.0] +[13.399999999999999,26.799999999999997,0.9999999999999991] +it :: Vector Numeric.LinearAlgebra.Data.R It also admits multiple right-hand sides stored as columns in a matrix. @@ -165,7 +169,8 @@ where -- | -- >>> build 5 (**2) :: Vector Double - -- fromList [0.0,1.0,4.0,9.0,16.0] + -- [0.0,1.0,4.0,9.0,16.0] + -- it :: Vector Double -- -- Hilbert matrix of order N: -- @@ -182,7 +187,7 @@ where build = build' -instance Container Matrix e => 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' @@ -202,12 +207,11 @@ {- | Compute mean vector and covariance matrix of the rows of a matrix. ->>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3]) -(fromList [4.010341078059521,5.0197204699640405], -(2><2) - [ 1.9862461923890056, -1.0127225830525157e-2 - , -1.0127225830525157e-2, 3.0373954915729318 ]) - +>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (trustSym $ diagl [2,3]) +([3.9933155655086696,5.061409102770331],Herm (2><2) + [ 1.9963242906624408, -4.227815571404954e-2 + , -4.227815571404954e-2, 3.2003833097832857 ]) +it :: (Vector Double, Herm Double) -} meanCov :: Matrix Double -> (Vector Double, Herm Double) meanCov x = (med,cov) where 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/Convolution.hs b/src/Internal/Convolution.hs --- a/src/Internal/Convolution.hs +++ b/src/Internal/Convolution.hs @@ -41,7 +41,8 @@ {- ^ correlation >>> corr (fromList[1,2,3]) (fromList [1..10]) -fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0] +[14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0] +it :: (Enum t, Product t, Container Vector t) => Vector t -} corr ker v @@ -54,7 +55,8 @@ {- ^ convolution ('corr' with reversed kernel and padded input, equivalent to polynomial product) >>> conv (fromList[1,1]) (fromList [-1,1]) -fromList [-1.0,0.0,1.0] +[-1.0,0.0,1.0] +it :: (Product t, Container Vector t) => Vector t -} conv ker v 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) @@ -54,13 +55,26 @@ -- | postfix error code check infixl 0 #| +(#|) :: IO CInt -> String -> IO () (#|) = flip check -- | 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,9 +1,10 @@ {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} + ----------------------------------------------------------------------------- -- | -- Module : Data.Packed.Matrix @@ -31,13 +32,14 @@ import Foreign.Storable(Storable) import System.IO.Unsafe(unsafePerformIO) import Control.Monad(liftM) +import Foreign.C.Types(CInt) ------------------------------------------------------------------- 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) @@ -53,8 +55,10 @@ show m | rows m == 0 || cols m == 0 = sizes m ++" []" show m = (sizes m++) . dsp . map (map show) . toLists $ m +sizes :: Matrix t -> [Char] sizes m = "("++show (rows m)++"><"++show (cols m)++")\n" +dsp :: [[[Char]]] -> [Char] dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp where mt = transpose as @@ -73,6 +77,7 @@ rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims +breakAt :: Eq a => a -> [a] -> ([a], [a]) breakAt c l = (a++[c],tail b) where (a,b) = break (==c) l @@ -88,7 +93,8 @@ | Drop Int | DropLast Int deriving Show - + +ppext :: Extractor -> [Char] ppext All = ":" ppext (Range a 1 c) = printf "%d:%d" a c ppext (Range a b c) = printf "%d:%d:%d" a b c @@ -128,10 +134,14 @@ infixl 9 ?? (??) :: Element t => Matrix t -> (Extractor,Extractor) -> Matrix t +minEl :: Vector CInt -> CInt minEl = toScalarI Min +maxEl :: Vector CInt -> CInt maxEl = toScalarI Max +cmodi :: Foreign.C.Types.CInt -> Vector Foreign.C.Types.CInt -> Vector Foreign.C.Types.CInt cmodi = vectorMapValI ModVS +extractError :: Matrix t1 -> (Extractor, Extractor) -> t extractError m (e1,e2)= error $ printf "can't extract (%s,%s) from matrix %dx%d" (ppext e1::String) (ppext e2::String) (rows m) (cols m) m ?? (Range a s b,e) | s /= 1 = m ?? (Pos (idxs [a,a+s .. b]), e) @@ -232,8 +242,10 @@ fromBlocks :: Element t => [[Matrix t]] -> Matrix t fromBlocks = fromBlocksRaw . adaptBlocks +fromBlocksRaw :: Element t => [[Matrix t]] -> Matrix t fromBlocksRaw mms = joinVert . map joinHoriz $ mms +adaptBlocks :: Element t => [[Matrix t]] -> [[Matrix t]] adaptBlocks ms = ms' where bc = case common length ms of Just c -> c @@ -486,6 +498,9 @@ m2' = conformMTo (r,c) m2 -- FIXME do not flatten if equal order +lM :: (Storable t, Element t1, Element t2) + => (Vector t1 -> Vector t2 -> Vector t) + -> Matrix t1 -> Matrix t2 -> Matrix t lM f m1 m2 = matrixFromVector RowMajor (max' (rows m1) (rows m2)) @@ -504,6 +519,7 @@ ------------------------------------------------------------ +toBlockRows :: Element t => [Int] -> Matrix t -> [Matrix t] toBlockRows [r] m | r == rows m = [m] toBlockRows rs m @@ -513,6 +529,7 @@ szs = map (* cols m) rs g k = (k><0)[] +toBlockCols :: Element t => [Int] -> Matrix t -> [Matrix t] toBlockCols [c] m | c == cols m = [m] toBlockCols cs m = map trans . toBlockRows cs . trans $ m @@ -576,7 +593,7 @@ mapMatrixWithIndexM :: (Element a, Storable b, Monad m) => ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b) -mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m +mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m where c = cols m @@ -598,4 +615,3 @@ mapMatrix :: (Element a, Element b) => (a -> b) -> Matrix a -> Matrix b mapMatrix f = liftMatrix (mapVector f) - diff --git a/src/Internal/IO.hs b/src/Internal/IO.hs --- a/src/Internal/IO.hs +++ b/src/Internal/IO.hs @@ -20,7 +20,7 @@ import Internal.Vector import Internal.Matrix import Internal.Vectorized -import Text.Printf(printf) +import Text.Printf(printf, PrintfArg, PrintfType) import Data.List(intersperse,transpose) import Data.Complex @@ -78,12 +78,18 @@ dispf :: Int -> Matrix Double -> String dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x +sdims :: Matrix t -> [Char] sdims x = show (rows x) ++ "x" ++ show (cols x) +formatFixed :: (Show a, Text.Printf.PrintfArg t, Element t) + => a -> Matrix t -> String formatFixed d x = format " " (printf ("%."++show d++"f")) $ x +isInt :: Matrix Double -> Bool isInt = all lookslikeInt . toList . flatten +formatScaled :: (Text.Printf.PrintfArg b, RealFrac b, Floating b, Num t, Element b, Show t) + => t -> Matrix b -> [Char] formatScaled dec t = "E"++show o++"\n" ++ ss where ss = format " " (printf fmt. g) t g x | o >= 0 = x/10^(o::Int) @@ -133,14 +139,18 @@ s2 = if b<0 then "-" else "" s3 = if b<0 then "-" else "+" +shcr :: (Show a, Show t1, Text.Printf.PrintfType t, Text.Printf.PrintfArg t1, RealFrac t1) + => a -> t1 -> t shcr d a | lookslikeInt a = printf "%.0f" a | otherwise = printf ("%."++show d++"f") a - +lookslikeInt :: (Show a, RealFrac a) => a -> Bool lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx where shx = show x +isZero :: Show a => a -> Bool isZero x = show x `elem` ["0.0","-0.0"] +isOne :: Show a => a -> Bool isOne x = show x `elem` ["1.0","-1.0"] -- | Pretty print a complex matrix with at most n decimal digits. @@ -168,6 +178,6 @@ else return (reshape c v) - +loadMatrix' :: FilePath -> IO (Maybe (Matrix Double)) loadMatrix' name = mbCatch (loadMatrix name) diff --git a/src/Internal/LAPACK.hs b/src/Internal/LAPACK.hs --- a/src/Internal/LAPACK.hs +++ b/src/Internal/LAPACK.hs @@ -1,6 +1,8 @@ {-# LANGUAGE TypeOperators #-} {-# LANGUAGE ViewPatterns #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-} + ----------------------------------------------------------------------------- -- | -- Module : Numeric.LinearAlgebra.LAPACK @@ -16,6 +18,8 @@ module Internal.LAPACK where +import Data.Bifunctor (first) + import Internal.Devel import Internal.Vector import Internal.Matrix hiding ((#), (#!)) @@ -232,7 +236,9 @@ ----------------------------------------------------------------------------- foreign import ccall unsafe "eig_l_R" dgeev :: R ::> R ::> C :> R ::> Ok +foreign import ccall unsafe "eig_l_G" dggev :: R ::> R ::> C :> R :> R ::> R ::> Ok foreign import ccall unsafe "eig_l_C" zgeev :: C ::> C ::> C :> C ::> Ok +foreign import ccall unsafe "eig_l_GC" zggev :: C ::> C ::> C :> C :> C ::> C ::> Ok foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> R :> R ::> Ok foreign import ccall unsafe "eig_l_H" zheev :: CInt -> R :> C ::> Ok @@ -296,13 +302,69 @@ | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) fixeig _ _ = error "fixeig with impossible inputs" +-- For dggev alpha(i) / beta(i), alpha(i+1) / beta(i+1) form a complex conjugate pair when Im alpha(i) != 0. +-- However, this does not lead to Re alpha(i) == Re alpha(i+1), since beta(i) and beta(i+1) +-- can be different. Therefore old 'fixeig' would fail for 'eigG'. +fixeigG [] _ = [] +fixeigG [_] [v] = [comp' v] +fixeigG ((_:+ai1) : an : as) (v1:v2:vs) + | abs ai1 > 1e-13 = toComplex' (v1, v2) : toComplex' (v1, mapVector negate v2) : fixeigG as vs + | otherwise = comp' v1 : fixeigG (an:as) (v2:vs) +fixeigG _ _ = error "fixeigG with impossible inputs" -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. -- The eigenvalues are not sorted. eigOnlyR :: Matrix Double -> Vector (Complex Double) eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" +-- | Generalized eigenvalues and right eigenvectors of a pair of real matrices, using LAPACK's /dggev/. +-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted. +eigG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double, Matrix (Complex Double)) +eigG a b = (alpha', beta, v'') + where + (alpha, beta, v) = eigGaux dggev a b "eigG" + alpha' = fixeig1 alpha + v' = toRows $ trans v + v'' = fromColumns $ fixeigG (toList alpha') v' +eigGaux f ma mb st = unsafePerformIO $ do + a <- copy ColumnMajor ma + b <- copy ColumnMajor mb + alpha <- createVector r + beta <- createVector r + vr <- createMatrix ColumnMajor r r + + (a # b # alpha # beta #! vr) g #| st + + return (alpha, beta, vr) + where + r = rows ma + g ar ac xra xca pa br bc xrb xcb pb alphan palpha betan pbeta = f ar ac xra xca pa br bc xrb xcb pb alphan palpha betan pbeta 0 0 0 0 nullPtr + +eigGOnlyAux f ma mb st = unsafePerformIO $ do + a <- copy ColumnMajor ma + b <- copy ColumnMajor mb + alpha <- createVector r + beta <- createVector r + + (a # b # alpha #! beta) g #| st + + return (alpha, beta) + where + r = rows ma + g ar ac xra xca pa br bc xrb xcb pb alphan palpha betan pbeta = f ar ac xra xca pa br bc xrb xcb pb alphan palpha betan pbeta 0 0 0 0 nullPtr 0 0 0 0 nullPtr + +-- | Generalized eigenvalues and right eigenvectors of a pair of complex matrices, using LAPACK's /zggev/. +-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted. +eigGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double), Matrix (Complex Double)) +eigGC a b = eigGaux zggev a b "eigGC" + +eigOnlyG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double) +eigOnlyG a b = first fixeig1 $ eigGOnlyAux dggev a b "eigOnlyG" + +eigOnlyGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double)) +eigOnlyGC a b = eigGOnlyAux zggev a b "eigOnlyGC" + ----------------------------------------------------------------------------- eigSHAux f st m = unsafePerformIO $ do @@ -444,8 +506,10 @@ | ndl == nd - 1 && ndu == nd - 1 && nd == r = unsafePerformIO . g $ do + dl' <- head . toRows <$> copy ColumnMajor (fromRows [dl]) + du' <- head . toRows <$> copy ColumnMajor (fromRows [du]) s <- copy ColumnMajor b - (dl # d # du #! s) f #| st + (dl' # d # du' #! s) f #| st return s | otherwise = error $ st ++ " of nonsquare matrix" where diff --git a/src/Internal/Matrix.hs b/src/Internal/Matrix.hs --- a/src/Internal/Matrix.hs +++ b/src/Internal/Matrix.hs @@ -57,19 +57,24 @@ cols = icols {-# INLINE cols #-} +size :: Matrix t -> (Int, Int) size m = (irows m, icols m) {-# INLINE size #-} +rowOrder :: Matrix t -> Bool rowOrder m = xCol m == 1 || cols m == 1 {-# INLINE rowOrder #-} +colOrder :: Matrix t -> Bool colOrder m = xRow m == 1 || rows m == 1 {-# INLINE colOrder #-} +is1d :: Matrix t -> Bool is1d (size->(r,c)) = r==1 || c==1 {-# INLINE is1d #-} -- data is not contiguous +isSlice :: Storable t => Matrix t -> Bool isSlice m@(size->(r,c)) = r*c < dim (xdat m) {-# INLINE isSlice #-} @@ -136,22 +141,27 @@ {-# INLINE applyRaw #-} infixr 1 # +(#) :: TransArray c => c -> (b -> IO r) -> Trans c b -> IO r a # b = apply a b {-# INLINE (#) #-} +(#!) :: (TransArray c, TransArray c1) => c1 -> c -> Trans c1 (Trans c (IO r)) -> IO r a #! b = a # b # id {-# INLINE (#!) #-} -------------------------------------------------------------------------------- +copy :: Element t => MatrixOrder -> Matrix t -> IO (Matrix t) copy ord m = extractR ord m 0 (idxs[0,rows m-1]) 0 (idxs[0,cols m-1]) +extractAll :: Element t => MatrixOrder -> Matrix t -> Matrix t extractAll ord m = unsafePerformIO (copy ord m) {- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose. >>> flatten (ident 3) -fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0] +[1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0] +it :: (Num t, Element t) => Vector t -} flatten :: Element t => Matrix t -> Vector t @@ -223,11 +233,13 @@ {-# INLINE (@@>) #-} -- Unsafe matrix access without range checking +atM' :: Storable t => Matrix t -> Int -> Int -> t atM' m i j = xdat m `at'` (i * (xRow m) + j * (xCol m)) {-# INLINE atM' #-} ------------------------------------------------------------------ +matrixFromVector :: Storable t => MatrixOrder -> Int -> Int -> Vector t -> Matrix t matrixFromVector _ 1 _ v@(dim->d) = Matrix { irows = 1, icols = d, xdat = v, xRow = d, xCol = 1 } matrixFromVector _ _ 1 v@(dim->d) = Matrix { irows = d, icols = 1, xdat = v, xRow = 1, xCol = d } matrixFromVector o r c v @@ -387,18 +399,21 @@ -------------------------------------------------------------------------- +maxZ :: (Num t1, Ord t1, Foldable t) => t t1 -> t1 maxZ xs = if minimum xs == 0 then 0 else maximum xs +conformMs :: Element t => [Matrix t] -> [Matrix t] conformMs ms = map (conformMTo (r,c)) ms where r = maxZ (map rows ms) c = maxZ (map cols ms) - +conformVs :: Element t => [Vector t] -> [Vector t] conformVs vs = map (conformVTo n) vs where n = maxZ (map dim vs) +conformMTo :: Element t => (Int, Int) -> Matrix t -> Matrix t conformMTo (r,c) m | size m == (r,c) = m | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c)) @@ -406,18 +421,24 @@ | size m == (1,c) = repRows r m | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to " ++ shDim (r,c) +conformVTo :: Element t => Int -> Vector t -> Vector t conformVTo n v | dim v == n = v | dim v == 1 = constantD (v@>0) n | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n +repRows :: Element t => Int -> Matrix t -> Matrix t repRows n x = fromRows (replicate n (flatten x)) +repCols :: Element t => Int -> Matrix t -> Matrix t repCols n x = fromColumns (replicate n (flatten x)) +shSize :: Matrix t -> [Char] shSize = shDim . size +shDim :: (Show a, Show a1) => (a1, a) -> [Char] shDim (r,c) = "(" ++ show r ++"x"++ show c ++")" +emptyM :: Storable t => Int -> Int -> Matrix t emptyM r c = matrixFromVector RowMajor r c (fromList[]) ---------------------------------------------------------------------- @@ -432,6 +453,11 @@ --------------------------------------------------------------- +extractAux :: (Eq t3, Eq t2, TransArray c, Storable a, Storable t1, + Storable t, Num t3, Num t2, Integral t1, Integral t) + => (t3 -> t2 -> CInt -> Ptr t1 -> CInt -> Ptr t + -> Trans c (CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)) + -> MatrixOrder -> c -> t3 -> Vector t1 -> t2 -> Vector t -> IO (Matrix a) extractAux f ord m moder vr modec vc = do let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc @@ -451,6 +477,9 @@ --------------------------------------------------------------- +setRectAux :: (TransArray c1, TransArray c) + => (CInt -> CInt -> Trans c1 (Trans c (IO CInt))) + -> Int -> Int -> c1 -> c -> IO () setRectAux f i j m r = (m #! r) (f (fi i) (fi j)) #|"setRect" type SetRect x = I -> I -> x ::> x::> Ok @@ -464,19 +493,29 @@ -------------------------------------------------------------------------------- +sortG :: (Storable t, Storable a) + => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a sortG f v = unsafePerformIO $ do r <- createVector (dim v) (v #! r) f #|"sortG" return r +sortIdxD :: Vector Double -> Vector CInt sortIdxD = sortG c_sort_indexD +sortIdxF :: Vector Float -> Vector CInt sortIdxF = sortG c_sort_indexF +sortIdxI :: Vector CInt -> Vector CInt sortIdxI = sortG c_sort_indexI +sortIdxL :: Vector Z -> Vector I sortIdxL = sortG c_sort_indexL +sortValD :: Vector Double -> Vector Double sortValD = sortG c_sort_valD +sortValF :: Vector Float -> Vector Float sortValF = sortG c_sort_valF +sortValI :: Vector CInt -> Vector CInt sortValI = sortG c_sort_valI +sortValL :: Vector Z -> Vector Z sortValL = sortG c_sort_valL foreign import ccall unsafe "sort_indexD" c_sort_indexD :: CV Double (CV CInt (IO CInt)) @@ -491,14 +530,21 @@ -------------------------------------------------------------------------------- +compareG :: (TransArray c, Storable t, Storable a) + => Trans c (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) + -> c -> Vector t -> Vector a compareG f u v = unsafePerformIO $ do r <- createVector (dim v) (u # v #! r) f #|"compareG" return r +compareD :: Vector Double -> Vector Double -> Vector CInt compareD = compareG c_compareD +compareF :: Vector Float -> Vector Float -> Vector CInt compareF = compareG c_compareF +compareI :: Vector CInt -> Vector CInt -> Vector CInt compareI = compareG c_compareI +compareL :: Vector Z -> Vector Z -> Vector CInt compareL = compareG c_compareL foreign import ccall unsafe "compareD" c_compareD :: CV Double (CV Double (CV CInt (IO CInt))) @@ -508,16 +554,33 @@ -------------------------------------------------------------------------------- +selectG :: (TransArray c, TransArray c1, TransArray c2, Storable t, Storable a) + => Trans c2 (Trans c1 (CInt -> Ptr t -> Trans c (CInt -> Ptr a -> IO CInt))) + -> c2 -> c1 -> Vector t -> c -> Vector a selectG f c u v w = unsafePerformIO $ do r <- createVector (dim v) (c # u # v # w #! r) f #|"selectG" return r +selectD :: Vector CInt -> Vector Double -> Vector Double -> Vector Double -> Vector Double selectD = selectG c_selectD +selectF :: Vector CInt -> Vector Float -> Vector Float -> Vector Float -> Vector Float selectF = selectG c_selectF +selectI :: Vector CInt -> Vector CInt -> Vector CInt -> Vector CInt -> Vector CInt selectI = selectG c_selectI +selectL :: Vector CInt -> Vector Z -> Vector Z -> Vector Z -> Vector Z selectL = selectG c_selectL +selectC :: Vector CInt + -> Vector (Complex Double) + -> Vector (Complex Double) + -> Vector (Complex Double) + -> Vector (Complex Double) selectC = selectG c_selectC +selectQ :: Vector CInt + -> Vector (Complex Float) + -> Vector (Complex Float) + -> Vector (Complex Float) + -> Vector (Complex Float) selectQ = selectG c_selectQ type Sel x = CV CInt (CV x (CV x (CV x (CV x (IO CInt))))) @@ -531,16 +594,29 @@ --------------------------------------------------------------------------- +remapG :: (TransArray c, TransArray c1, Storable t, Storable a) + => (CInt -> CInt -> CInt -> CInt -> Ptr t + -> Trans c1 (Trans c (CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt))) + -> Matrix t -> c1 -> c -> Matrix a remapG f i j m = unsafePerformIO $ do r <- createMatrix RowMajor (rows i) (cols i) (i # j # m #! r) f #|"remapG" return r +remapD :: Matrix CInt -> Matrix CInt -> Matrix Double -> Matrix Double remapD = remapG c_remapD +remapF :: Matrix CInt -> Matrix CInt -> Matrix Float -> Matrix Float remapF = remapG c_remapF +remapI :: Matrix CInt -> Matrix CInt -> Matrix CInt -> Matrix CInt remapI = remapG c_remapI +remapL :: Matrix CInt -> Matrix CInt -> Matrix Z -> Matrix Z remapL = remapG c_remapL +remapC :: Matrix CInt + -> Matrix CInt + -> Matrix (Complex Double) + -> Matrix (Complex Double) remapC = remapG c_remapC +remapQ :: Matrix CInt -> Matrix CInt -> Matrix (Complex Float) -> Matrix (Complex Float) remapQ = remapG c_remapQ type Rem x = OM CInt (OM CInt (OM x (OM x (IO CInt)))) @@ -554,6 +630,9 @@ -------------------------------------------------------------------------------- +rowOpAux :: (TransArray c, Storable a) => + (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Trans c (IO CInt)) + -> Int -> a -> Int -> Int -> Int -> Int -> c -> IO () rowOpAux f c x i1 i2 j1 j2 m = do px <- newArray [x] (m # id) (f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2)) #|"rowOp" @@ -572,6 +651,9 @@ -------------------------------------------------------------------------------- +gemmg :: (TransArray c1, TransArray c, TransArray c2, TransArray c3) + => Trans c3 (Trans c2 (Trans c1 (Trans c (IO CInt)))) + -> c3 -> c2 -> c1 -> c -> IO () gemmg f v m1 m2 m3 = (v # m1 # m2 #! m3) f #|"gemmg" type Tgemm x = x :> x ::> x ::> x ::> Ok @@ -587,6 +669,10 @@ -------------------------------------------------------------------------------- +reorderAux :: (TransArray c, Storable t, Storable a1, Storable t1, Storable a) => + (CInt -> Ptr a -> CInt -> Ptr t1 + -> Trans c (CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt)) + -> Vector t1 -> c -> Vector t -> Vector a1 reorderAux f s d v = unsafePerformIO $ do k <- createVector (dim s) r <- createVector (dim v) diff --git a/src/Internal/Modular.hs b/src/Internal/Modular.hs --- a/src/Internal/Modular.hs +++ b/src/Internal/Modular.hs @@ -1,18 +1,15 @@ {-# 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 #-} + {- | Module : Internal.Modular Copyright : (c) Alberto Ruiz 2015 @@ -80,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) @@ -92,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 @@ -148,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 @@ -165,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/Numeric.hs b/src/Internal/Numeric.hs --- a/src/Internal/Numeric.hs +++ b/src/Internal/Numeric.hs @@ -5,6 +5,8 @@ {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE UndecidableInstances #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-} + ----------------------------------------------------------------------------- -- | -- Module : Data.Packed.Internal.Numeric @@ -788,13 +790,7 @@ type instance RealOf I = I type instance RealOf Z = Z -type family ComplexOf x - -type instance ComplexOf Double = Complex Double -type instance ComplexOf (Complex Double) = Complex Double - -type instance ComplexOf Float = Complex Float -type instance ComplexOf (Complex Float) = Complex Float +type ComplexOf x = Complex (RealOf x) type family SingleOf x diff --git a/src/Internal/ST.hs b/src/Internal/ST.hs --- a/src/Internal/ST.hs +++ b/src/Internal/ST.hs @@ -81,6 +81,8 @@ unsafeFreezeVector (STVector x) = unsafeIOToST . return $ x {-# INLINE safeIndexV #-} +safeIndexV :: Storable t2 + => (STVector s t2 -> Int -> t) -> STVector t1 t2 -> Int -> t safeIndexV f (STVector v) k | k < 0 || k>= dim v = error $ "out of range error in vector (dim=" ++show (dim v)++", pos="++show k++")" @@ -150,9 +152,12 @@ freezeMatrix :: (Element t) => STMatrix s t -> ST s (Matrix t) freezeMatrix m = liftSTMatrix id m +cloneMatrix :: Element t => Matrix t -> IO (Matrix t) cloneMatrix m = copy (orderOf m) m {-# INLINE safeIndexM #-} +safeIndexM :: (STMatrix s t2 -> Int -> Int -> t) + -> STMatrix t1 t2 -> Int -> Int -> t safeIndexM f (STMatrix m) r c | r<0 || r>=rows m || c<0 || c>=cols m = error $ "out of range error in matrix (size=" @@ -184,6 +189,7 @@ | Col Int | FromCol Int +getColRange :: Int -> ColRange -> (Int, Int) 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) @@ -194,6 +200,7 @@ | Row Int | FromRow Int +getRowRange :: Int -> RowRange -> (Int, 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) @@ -223,6 +230,7 @@ i2' = i2 `mod` (rows m) +extractMatrix :: Element a => STMatrix t a -> RowRange -> ColRange -> ST s (Matrix a) extractMatrix (STMatrix m) rr rc = unsafeIOToST (extractR (orderOf m) m 0 (idxs[i1,i2]) 0 (idxs[j1,j2])) where (i1,i2) = getRowRange (rows m) rr @@ -231,6 +239,7 @@ -- | r0 c0 height width data Slice s t = Slice (STMatrix s t) Int Int Int Int +slice :: Element a => Slice t a -> Matrix a slice (Slice (STMatrix m) r0 c0 nr nc) = subMatrix (r0,c0) (nr,nc) m gemmm :: Element t => t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s () @@ -238,12 +247,11 @@ where res = unsafeIOToST (gemm v a b r) v = fromList [alpha,beta] - + mutable :: Element t => (forall s . (Int, Int) -> STMatrix s t -> ST s u) -> Matrix t -> (Matrix t,u) mutable f a = runST $ do x <- thawMatrix a info <- f (rows a, cols a) x r <- unsafeFreezeMatrix x return (r,info) - diff --git a/src/Internal/Sparse.hs b/src/Internal/Sparse.hs --- a/src/Internal/Sparse.hs +++ b/src/Internal/Sparse.hs @@ -1,9 +1,11 @@ -{-# LANGUAGE RecordWildCards #-} -{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE BangPatterns #-} {-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE RecordWildCards #-} module Internal.Sparse( - GMatrix(..), CSR(..), mkCSR, fromCSR, + GMatrix(..), CSR(..), mkCSR, fromCSR, impureCSR, mkSparse, mkDiagR, mkDense, AssocMatrix, toDense, @@ -14,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 @@ -25,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 @@ -47,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. @@ -122,11 +179,12 @@ fromCSR :: CSR -> GMatrix fromCSR csr = SparseR {..} where - gmCSR @ CSR {..} = csr + gmCSR@CSR {..} = csr nRows = csrNRows 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) @@ -142,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 @@ -170,8 +232,10 @@ {- | general matrix - vector product >>> let m = mkSparse [((0,999),1.0),((1,1999),2.0)] +m :: GMatrix >>> m !#> vector [1..2000] -fromList [1000.0,4000.0] +[1000.0,4000.0] +it :: Vector Double -} infixr 8 !#> diff --git a/src/Internal/Static.hs b/src/Internal/Static.hs --- a/src/Internal/Static.hs +++ b/src/Internal/Static.hs @@ -15,6 +15,8 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE DeriveGeneric #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-} + {- | Module : Internal.Static Copyright : (c) Alberto Ruiz 2006-14 @@ -128,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 @@ -321,34 +323,34 @@ instance KnownNat n => Show (R n) where show s@(R (Dim v)) - | singleV v = "("++show (v!0)++" :: R "++show d++")" - | otherwise = "(vector"++ drop 8 (show v)++" :: R "++show d++")" + | singleV v = "(" ++ show (v!0) ++ " :: R " ++ show d ++ ")" + | otherwise = "(vector " ++ show v ++ " :: R " ++ show d ++")" where d = size s instance KnownNat n => Show (C n) where show s@(C (Dim v)) - | singleV v = "("++show (v!0)++" :: C "++show d++")" - | otherwise = "(vector"++ drop 8 (show v)++" :: C "++show d++")" + | singleV v = "(" ++ show (v!0) ++ " :: C " ++ show d ++ ")" + | otherwise = "(vector " ++ show v ++ " :: C " ++ show d ++")" where d = size s instance (KnownNat m, KnownNat n) => Show (L m n) where - show (isDiag -> Just (z,y,(m',n'))) = printf "(diag %s %s :: L %d %d)" (show z) (drop 9 $ show y) m' n' + show (isDiag -> Just (z,y,(m',n'))) = printf "(diag %s %s :: L %d %d)" (show z) (show y) m' n' show s@(L (Dim (Dim x))) | singleM x = printf "(%s :: L %d %d)" (show (x `atIndex` (0,0))) m' n' - | otherwise = "(matrix"++ dropWhile (/='\n') (show x)++" :: L "++show m'++" "++show n'++")" + | otherwise = "(matrix" ++ dropWhile (/='\n') (show x) ++ " :: L " ++ show m' ++ " " ++ show n' ++ ")" where (m',n') = size s instance (KnownNat m, KnownNat n) => Show (M m n) where - show (isDiagC -> Just (z,y,(m',n'))) = printf "(diag %s %s :: M %d %d)" (show z) (drop 9 $ show y) m' n' + show (isDiagC -> Just (z,y,(m',n'))) = printf "(diag %s %s :: M %d %d)" (show z) (show y) m' n' show s@(M (Dim (Dim x))) | singleM x = printf "(%s :: M %d %d)" (show (x `atIndex` (0,0))) m' n' - | otherwise = "(matrix"++ dropWhile (/='\n') (show x)++" :: M "++show m'++" "++show n'++")" + | otherwise = "(matrix" ++ dropWhile (/='\n') (show x) ++ " :: M " ++ show m' ++ " " ++ show n' ++ ")" where (m',n') = size s @@ -364,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 (/) @@ -389,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) (*) @@ -399,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 @@ -6,6 +6,8 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE ViewPatterns #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} ----------------------------------------------------------------------------- {- | @@ -97,7 +99,8 @@ {- | Create a real vector. >>> vector [1..5] -fromList [1.0,2.0,3.0,4.0,5.0] +[1.0,2.0,3.0,4.0,5.0] +it :: Vector R -} vector :: [R] -> Vector R @@ -376,7 +379,8 @@ On a matrix it gets the k-th row as a vector: >>> matrix 5 [1..15] ! 1 -fromList [6.0,7.0,8.0,9.0,10.0] +[6.0,7.0,8.0,9.0,10.0] +it :: Vector Double >>> matrix 5 [1..15] ! 1 ! 3 9.0 @@ -413,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 @@ -628,7 +632,7 @@ y:ys = redu (pivot n xs) pivot k = (const k &&& id) - . sortBy (flip compare `on` (abs. (!k))) + . sortBy (flip compare `on` (abs. (! k))) redu :: (Int, [Vector t]) -> [Vector t] redu (k,x:zs) @@ -908,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/Internal/Vector.hs b/src/Internal/Vector.hs --- a/src/Internal/Vector.hs +++ b/src/Internal/Vector.hs @@ -1,6 +1,7 @@ {-# LANGUAGE MagicHash, UnboxedTuples, BangPatterns, FlexibleContexts #-} {-# LANGUAGE TypeSynonymInstances #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} -- | -- Module : Internal.Vector @@ -40,6 +41,7 @@ import Data.Vector.Storable(Vector, fromList, unsafeToForeignPtr, unsafeFromForeignPtr, unsafeWith) import Data.Binary +import Data.Binary.Put import Control.Monad(replicateM) import qualified Data.ByteString.Internal as BS import Data.Vector.Storable.Internal(updPtr) @@ -92,6 +94,7 @@ -} +safeRead :: Storable a => Vector a -> (Ptr a -> IO c) -> c safeRead v = inlinePerformIO . unsafeWith v {-# INLINE safeRead #-} @@ -113,7 +116,8 @@ be used, for instance, with infinite lists. >>> 5 |> [1..] -fromList [1.0,2.0,3.0,4.0,5.0] +[1.0,2.0,3.0,4.0,5.0] +it :: (Enum a, Num a, Foreign.Storable.Storable a) => Vector a -} (|>) :: (Storable a) => Int -> [a] -> Vector a @@ -132,7 +136,8 @@ {- | takes a number of consecutive elements from a Vector >>> subVector 2 3 (fromList [1..10]) -fromList [3.0,4.0,5.0] +[3.0,4.0,5.0] +it :: (Enum t, Num t, Foreign.Storable.Storable t) => Vector t -} subVector :: Storable t => Int -- ^ index of the starting element @@ -166,7 +171,8 @@ {- | concatenate a list of vectors >>> vjoin [fromList [1..5::Double], konst 1 3] -fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0] +[1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0] +it :: Vector Double -} vjoin :: Storable t => [Vector t] -> Vector t @@ -188,7 +194,8 @@ {- | Extract consecutive subvectors of the given sizes. >>> takesV [3,4] (linspace 10 (1,10::Double)) -[fromList [1.0,2.0,3.0],fromList [4.0,5.0,6.0,7.0]] +[[1.0,2.0,3.0],[4.0,5.0,6.0,7.0]] +it :: [Vector Double] -} takesV :: Storable t => [Int] -> Vector t -> [Vector t] @@ -283,11 +290,13 @@ go (dim v -1) x {-# INLINE foldVectorWithIndex #-} +foldLoop :: (Int -> t -> t) -> t -> Int -> t foldLoop f s0 d = go (d - 1) s0 where go 0 s = f (0::Int) s go !j !s = go (j - 1) (f j s) +foldVectorG :: Storable t1 => (Int -> (Int -> t1) -> t -> t) -> t -> Vector t1 -> t foldVectorG f s0 v = foldLoop g s0 (dim v) where g !k !s = f k (safeRead v . flip peekElemOff) s {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479) @@ -390,8 +399,10 @@ m = d `mod` chunk in if m /= 0 then reverse (m:(replicate c chunk)) else (replicate c chunk) +putVector :: (Storable t, Binary t) => Vector t -> Data.Binary.Put.PutM () putVector v = mapM_ put $! toList v +getVector :: (Storable a, Binary a) => Int -> Get (Vector a) getVector d = do xs <- replicateM d get return $! fromList xs diff --git a/src/Internal/Vectorized.hs b/src/Internal/Vectorized.hs --- a/src/Internal/Vectorized.hs +++ b/src/Internal/Vectorized.hs @@ -28,12 +28,15 @@ import Control.Monad(when) infixr 1 # +(#) :: TransArray c => c -> (b -> IO r) -> TransRaw c b -> IO r a # b = applyRaw a b {-# INLINE (#) #-} +(#!) :: (TransArray c, TransArray c1) => c1 -> c -> TransRaw c1 (TransRaw c (IO r)) -> IO r a #! b = a # b # id {-# INLINE (#!) #-} +fromei :: Enum a => a -> CInt fromei x = fromIntegral (fromEnum x) :: CInt data FunCodeV = Sin @@ -100,10 +103,20 @@ sumC :: Vector (Complex Double) -> Complex Double sumC = sumg c_sumC +sumI :: ( TransRaw c (CInt -> Ptr a -> IO CInt) ~ (CInt -> Ptr I -> I :> Ok) + , TransArray c + , Storable a + ) + => I -> c -> a sumI m = sumg (c_sumI m) +sumL :: ( TransRaw c (CInt -> Ptr a -> IO CInt) ~ (CInt -> Ptr Z -> Z :> Ok) + , TransArray c + , Storable a + ) => Z -> c -> a sumL m = sumg (c_sumL m) +sumg :: (TransArray c, Storable a) => TransRaw c (CInt -> Ptr a -> IO CInt) -> c -> a sumg f x = unsafePerformIO $ do r <- createVector 1 (x #! r) f #| "sum" @@ -140,6 +153,8 @@ prodL :: Z-> Vector Z -> Z prodL = prodg . c_prodL +prodg :: (TransArray c, Storable a) + => TransRaw c (CInt -> Ptr a -> IO CInt) -> c -> a prodg f x = unsafePerformIO $ do r <- createVector 1 (x #! r) f #| "prod" @@ -155,16 +170,25 @@ ------------------------------------------------------------------ +toScalarAux :: (Enum a, TransArray c, Storable a1) + => (CInt -> TransRaw c (CInt -> Ptr a1 -> IO CInt)) -> a -> c -> a1 toScalarAux fun code v = unsafePerformIO $ do r <- createVector 1 (v #! r) (fun (fromei code)) #|"toScalarAux" return (r @> 0) + +vectorMapAux :: (Enum a, Storable t, Storable a1) + => (CInt -> CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt) + -> a -> Vector t -> Vector a1 vectorMapAux fun code v = unsafePerformIO $ do r <- createVector (dim v) (v #! r) (fun (fromei code)) #|"vectorMapAux" return r +vectorMapValAux :: (Enum a, Storable a2, Storable t, Storable a1) + => (CInt -> Ptr a2 -> CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt) + -> a -> a2 -> Vector t -> Vector a1 vectorMapValAux fun code val v = unsafePerformIO $ do r <- createVector (dim v) pval <- newArray [val] @@ -172,6 +196,9 @@ free pval return r +vectorZipAux :: (Enum a, TransArray c, Storable t, Storable a1) + => (CInt -> CInt -> Ptr t -> TransRaw c (CInt -> Ptr a1 -> IO CInt)) + -> a -> Vector t -> c -> Vector a1 vectorZipAux fun code u v = unsafePerformIO $ do r <- createVector (dim u) (u # v #! r) (fun (fromei code)) #|"vectorZipAux" @@ -378,6 +405,7 @@ -------------------------------------------------------------------------------- +roundVector :: Vector Double -> Vector Double roundVector v = unsafePerformIO $ do r <- createVector (dim v) (v #! r) c_round_vector #|"roundVector" @@ -389,7 +417,8 @@ -- | -- >>> range 5 --- fromList [0,1,2,3,4] +-- [0,1,2,3,4] +-- it :: Vector I -- range :: Int -> Vector I range n = unsafePerformIO $ do @@ -432,6 +461,8 @@ long2intV = tog c_long2int +tog :: (Storable t, Storable a) + => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a tog f v = unsafePerformIO $ do r <- createVector (dim v) (v #! r) f #|"tog" @@ -451,6 +482,8 @@ --------------------------------------------------------------- +stepg :: (Storable t, Storable a) + => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a stepg f v = unsafePerformIO $ do r <- createVector (dim v) (v #! r) f #|"step" @@ -476,6 +509,8 @@ -------------------------------------------------------------------------------- +conjugateAux :: (Storable t, Storable a) + => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a conjugateAux fun x = unsafePerformIO $ do v <- createVector (dim x) (x #! v) fun #|"conjugateAux" @@ -501,6 +536,8 @@ -------------------------------------------------------------------------------- +constantAux :: (Storable a1, Storable a) + => (Ptr a1 -> CInt -> Ptr a -> IO CInt) -> a1 -> Int -> Vector a constantAux fun x n = unsafePerformIO $ do v <- createVector n px <- newArray [x] diff --git a/src/Numeric/LinearAlgebra.hs b/src/Numeric/LinearAlgebra.hs --- a/src/Numeric/LinearAlgebra.hs +++ b/src/Numeric/LinearAlgebra.hs @@ -1,6 +1,8 @@ {-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-} + ----------------------------------------------------------------------------- {- | Module : Numeric.LinearAlgebra @@ -27,7 +29,8 @@ -- as the Hadamard product or the Schur product): -- -- >>> vector [1,2,3] * vector [3,0,-2] - -- fromList [3.0,0.0,-6.0] + -- [3.0,0.0,-6.0] + -- it :: Vector R -- -- >>> matrix 3 [1..9] * ident 3 -- (3><3) @@ -128,8 +131,8 @@ leftSV, rightSV, -- * Eigendecomposition - eig, eigSH, - eigenvalues, eigenvaluesSH, + eig, geig, eigSH, + eigenvalues, geigenvalues, eigenvaluesSH, geigSH, -- * QR 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/HMatrix.hs b/src/Numeric/LinearAlgebra/HMatrix.hs --- a/src/Numeric/LinearAlgebra/HMatrix.hs +++ b/src/Numeric/LinearAlgebra/HMatrix.hs @@ -28,7 +28,9 @@ (<ยท>) :: Numeric t => Vector t -> Vector t -> t (<ยท>) = dot +app :: Numeric t => Matrix t -> Vector t -> Vector t app m v = m #> v +mul :: Numeric t => Matrix t -> Matrix t -> Matrix t mul a b = a <> b 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 @@ -14,6 +14,8 @@ {-# LANGUAGE GADTs #-} {-# LANGUAGE TypeFamilies #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} {- | Module : Numeric.LinearAlgebra.Static @@ -44,6 +46,13 @@ matrix, -- * Complex โ„‚, C, M, Her, her, ๐‘–, + toComplex, + fromComplex, + complex, + real, + imag, + sqMagnitude, + magnitude, -- * Products (<>),(#>),(<.>), -- * Linear Systems @@ -72,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(..)) @@ -90,7 +101,7 @@ infixl 4 & -(&) :: forall n . (KnownNat n, 1 <= n) +(&) :: forall n . KnownNat n => R n -> โ„ -> R (n+1) u & x = u # (konst x :: R 1) @@ -241,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) @@ -381,7 +415,7 @@ headTail :: (KnownNat n, 1<=n) => R n -> (โ„, R (n-1)) -headTail = ((!0) . extract *** id) . split +headTail = ((! 0) . extract *** id) . split splitRows :: forall p m n . (KnownNat p, KnownNat m, KnownNat n, p<=m) => L m n -> (L p n, L (m-p) n) @@ -593,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) @@ -656,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) diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs --- a/src/Numeric/Matrix.hs +++ b/src/Numeric/Matrix.hs @@ -4,6 +4,8 @@ {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} + ----------------------------------------------------------------------------- -- | -- Module : Numeric.Matrix @@ -35,6 +37,7 @@ import qualified Data.Foldable as F import qualified Data.Semigroup as S import Internal.Chain +import Foreign.Storable(Storable) ------------------------------------------------------------------- @@ -80,8 +83,16 @@ -------------------------------------------------------------------------------- +isScalar :: Matrix t -> Bool isScalar m = rows m == 1 && cols m == 1 +adaptScalarM :: (Foreign.Storable.Storable t1, Foreign.Storable.Storable t2) + => (t1 -> Matrix t2 -> t) + -> (Matrix t1 -> Matrix t2 -> t) + -> (Matrix t1 -> t2 -> t) + -> Matrix t1 + -> Matrix t2 + -> t adaptScalarM f1 f2 f3 x y | isScalar x = f1 (x @@>(0,0) ) y | isScalar y = f3 x (y @@>(0,0) ) @@ -96,7 +107,7 @@ where mempty = 1 mappend = adaptScalarM scale mXm (flip scale) - + mconcat xs = work (partition isScalar xs) where work (ss,[]) = product ss @@ -106,4 +117,3 @@ | otherwise = scale x00 m where x00 = x @@> (0,0) - diff --git a/src/Numeric/Vector.hs b/src/Numeric/Vector.hs --- a/src/Numeric/Vector.hs +++ b/src/Numeric/Vector.hs @@ -3,6 +3,9 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} + +{-# OPTIONS_GHC -fno-warn-orphans #-} + ----------------------------------------------------------------------------- -- | -- Module : Numeric.Vector @@ -14,7 +17,7 @@ -- -- Provides instances of standard classes 'Show', 'Read', 'Eq', -- 'Num', 'Fractional', and 'Floating' for 'Vector'. --- +-- ----------------------------------------------------------------------------- module Numeric.Vector () where @@ -23,9 +26,17 @@ import Internal.Vector import Internal.Numeric import Internal.Conversion +import Foreign.Storable(Storable) ------------------------------------------------------------------- +adaptScalar :: (Foreign.Storable.Storable t1, Foreign.Storable.Storable t2) + => (t1 -> Vector t2 -> t) + -> (Vector t1 -> Vector t2 -> t) + -> (Vector t1 -> t2 -> t) + -> Vector t1 + -> Vector t2 + -> t adaptScalar f1 f2 f3 x y | dim x == 1 = f1 (x@>0) y | dim y == 1 = f3 x (y@>0) @@ -172,4 +183,3 @@ sqrt = vectorMapQ Sqrt (**) = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS)) pi = fromList [pi] -