hmatrix 0.19.0.0 → 0.20.2
raw patch · 21 files changed
Files
- hmatrix.cabal +23/−12
- src/Internal/Algorithms.hs +32/−6
- src/Internal/C/lapack-aux.c +539/−225
- src/Internal/C/vector-aux.c +23/−10
- src/Internal/CG.hs +2/−2
- src/Internal/Container.hs +19/−17
- src/Internal/Conversion.hs +0/−1
- src/Internal/Convolution.hs +4/−2
- src/Internal/Devel.hs +17/−4
- src/Internal/Element.hs +1/−2
- src/Internal/LAPACK.hs +63/−1
- src/Internal/Matrix.hs +2/−1
- src/Internal/Modular.hs +5/−9
- src/Internal/Sparse.hs +92/−30
- src/Internal/Static.hs +14/−15
- src/Internal/Util.hs +6/−5
- src/Internal/Vector.hs +8/−4
- src/Internal/Vectorized.hs +2/−1
- src/Numeric/LinearAlgebra.hs +4/−3
- src/Numeric/LinearAlgebra/Devel.hs +3/−3
- src/Numeric/LinearAlgebra/Static.hs +39/−7
hmatrix.cabal view
@@ -1,11 +1,11 @@ Name: hmatrix-Version: 0.19.0.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: <http://dis.um.es/~alberto/hmatrix/hmatrix.html> 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
src/Internal/Algorithms.hs view
@@ -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@@ -70,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)@@ -96,7 +101,9 @@ linearSolveSVD' = linearSolveSVDR Nothing eig' = eigR eigSH'' = eigS+ geig' = eigG eigOnly = eigOnlyR+ geigOnly = eigOnlyG eigOnlySH = eigOnlyS cholSH' = cholS mbCholSH' = mbCholS@@ -126,7 +133,9 @@ linearSolveLS' = linearSolveLSC linearSolveSVD' = linearSolveSVDC Nothing eig' = eigC+ geig' = eigGC eigOnly = eigOnlyC+ geigOnly = eigOnlyGC eigSH'' = eigH eigOnlySH = eigOnlyH cholSH' = cholH@@ -170,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@@ -224,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@@ -283,7 +294,8 @@ -0.690 -0.352 >>> s-fromList [35.18264833189422,1.4769076999800903]+[35.18264833189422,1.476907699980091]+it :: Vector Double >>> disp 3 u 5x2@@ -509,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''@@ -535,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@@ -1135,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)-
src/Internal/C/lapack-aux.c view
@@ -45,6 +45,9 @@ for(q=0;q<M##r*M##c;q++) printf("%.1f ",M##p[q]); printf("\n");} #define CHECK(RES,CODE) MACRO(if(RES) return CODE;)+#define MARK(RES,CODE) MACRO(if(RES) { ret = CODE; })+#define CONVERGED(RES,CODE) MACRO(if(RES > 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<mn; k++) {+ UNWIND(res,res,cleanup1);++ for (int k=0; k<mn; k++) { ipivp[k] = auxipiv[k]; }++cleanup1: free(auxipiv);- OK+cleanup0:+ return ret; } @@ -1166,24 +1439,31 @@ integer *lda, integer *ipiv, integer *info); int lu_l_C(DVEC(ipiv), 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 && 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<mn; k++) {+ UNWIND(res,res,cleanup1);++ for (int k=0; k<mn; k++) { ipivp[k] = auxipiv[k]; }++cleanup1: free(auxipiv);- OK+cleanup0:+ return ret; } @@ -1194,23 +1474,26 @@ ldb, integer *info); int luS_l_R(KODMAT(a), KDVEC(ipiv), ODMAT(b)) {- integer m = ar;- integer n = ac;- integer lda = aXc;- integer mrhs = br;- integer nrhs = bc;+ integer ret = 0;+ integer m = ar;+ integer n = ac;+ integer lda = aXc;+ integer mrhs = br;+ integer nrhs = bc; - REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);- integer* auxipiv = (integer*)malloc(n*sizeof(integer));- int k;- for (k=0; k<n; k++) {- auxipiv[k] = (integer)ipivp[k];- }- integer res;- dgetrs_ ("N",&n,&nrhs,(/*no const (!?)*/ double*)ap,&lda,auxipiv,bp,&mrhs,&res);- CHECK(res,res);- free(auxipiv);- OK+ REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);+ integer* auxipiv = (integer*)malloc(n*sizeof(integer));+ CHECK(!auxipiv,MEM);++ for (int k=0; k<n; k++) {+ auxipiv[k] = (integer)ipivp[k];+ }+ integer res;+ dgetrs_ ("N",&n,&nrhs,(/*no const (!?)*/ double*)ap,&lda,auxipiv,bp,&mrhs,&res);+ MARK(res,res);++ free(auxipiv);+ return ret; } @@ -1219,23 +1502,26 @@ integer *ldb, integer *info); int luS_l_C(KOCMAT(a), KDVEC(ipiv), OCMAT(b)) {- integer m = ar;- integer n = ac;- integer lda = aXc;+ integer ret = 0;+ integer m = ar;+ integer n = ac;+ integer lda = aXc; integer mrhs = br; integer nrhs = bc; REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE); integer* auxipiv = (integer*)malloc(n*sizeof(integer));- int k;- for (k=0; k<n; k++) {+ CHECK(!auxipiv,MEM);++ for (int k=0; k<n; k++) { auxipiv[k] = (integer)ipivp[k]; } integer res; zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&lda,auxipiv,bp,&mrhs,&res);- CHECK(res,res);+ MARK(res,res);+ free(auxipiv);- OK+ return ret; } @@ -1245,10 +1531,15 @@ doublereal *work, integer *lwork, integer *info); int ldl_R(DVEC(ipiv), ODMAT(r)) {- integer n = rr;+ integer ret = 0;+ integer n = rr;+ REQUIRES(n>=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<n; k++) { ipivp[k] = auxipiv[k]; }- free(auxipiv);++cleanup2: free(work);- OK+cleanup1:+ free(auxipiv);+cleanup0:+ return ret; } @@ -1272,10 +1570,14 @@ doublecomplex *work, integer *lwork, integer *info); int ldl_C(DVEC(ipiv), OCMAT(r)) {- integer n = rr;+ integer ret = 0;+ integer n = rr;+ REQUIRES(n>=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<n; k++) { ipivp[k] = auxipiv[k]; }- free(auxipiv);++cleanup2: free(work);- OK+cleanup1:+ free(auxipiv);+cleanup0:+ return ret; } @@ -1301,23 +1609,26 @@ integer *ipiv, doublereal *b, integer *ldb, integer *info); int ldl_S_R(KODMAT(a), KDVEC(ipiv), ODMAT(b)) {- integer m = ar;- integer n = ac;- integer lda = aXc;- integer mrhs = br;- integer nrhs = bc;+ integer ret = 0;+ integer m = ar;+ integer n = ac;+ integer lda = aXc;+ integer mrhs = br;+ integer nrhs = bc; - REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);- integer* auxipiv = (integer*)malloc(n*sizeof(integer));- int k;- for (k=0; k<n; k++) {- auxipiv[k] = (integer)ipivp[k];- }- integer res;- dsytrs_ ("L",&n,&nrhs,(/*no const (!?)*/ double*)ap,&lda,auxipiv,bp,&mrhs,&res);- CHECK(res,res);- free(auxipiv);- OK+ REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);+ integer* auxipiv = (integer*)malloc(n*sizeof(integer));+ CHECK(!auxipiv,MEM);++ for (int k=0; k<n; k++) {+ auxipiv[k] = (integer)ipivp[k];+ }+ integer res;+ dsytrs_ ("L",&n,&nrhs,(/*no const (!?)*/ double*)ap,&lda,auxipiv,bp,&mrhs,&res);+ MARK(res,res);++ free(auxipiv);+ return ret; } @@ -1325,23 +1636,26 @@ integer *ipiv, doublecomplex *b, integer *ldb, integer *info); int ldl_S_C(KOCMAT(a), KDVEC(ipiv), OCMAT(b)) {- integer m = ar;- integer n = ac;- integer lda = aXc;+ integer ret = 0;+ integer m = ar;+ integer n = ac;+ integer lda = aXc; integer mrhs = br; integer nrhs = bc; REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE); integer* auxipiv = (integer*)malloc(n*sizeof(integer));- int k;- for (k=0; k<n; k++) {+ CHECK(!auxipiv,MEM);++ for (int k=0; k<n; k++) { auxipiv[k] = (integer)ipivp[k]; } integer res; zhetrs_ ("L",&n,&nrhs,(doublecomplex*)ap,&lda,auxipiv,bp,&mrhs,&res);- CHECK(res,res);+ MARK(res,res);+ free(auxipiv);- OK+ return ret; }
src/Internal/C/vector-aux.c view
@@ -932,20 +932,33 @@ //////////////////////////////////////////////////////////////////////////////// -#if defined (__APPLE__) || (__FreeBSD__)-/* FreeBSD and Mac OS X do not provide random_r(), thread safety cannot be- guaranteed.+#if defined (__APPLE__) || (__FreeBSD__) || defined(NO_RANDOM_R) || defined(_WIN32) || defined(WIN32)+/* Windows use thread-safe random+ See: http://stackoverflow.com/questions/143108/is-windows-rand-s-thread-safe+*/+#if defined (__APPLE__) || (__FreeBSD__) || defined(NO_RANDOM_R)++/* For FreeBSD, Mac OS X, and other libcs (like `musl`) that do not provide+ random_r(), or if the use of random_r() is explicitly disabled, thread safety+ cannot be guaranteed.+ As per current understanding, this should at worst lead to less "random"+ numbers being generated, in particular+ * if another thread somebody calls lcong48() at the same time as nrand48()+ is called+ * in addition to that, for glibc with NO_RANDOM_R enabled when ndrand48()+ is called for the first time by multiple threads in parallel due to the+ initialisation function placed within it+ See: http://www.evanjones.ca/random-thread-safe.html+ For FreeBSD and Mac OS X, nrand48() is much better than random(). See: http://www.evanjones.ca/random-thread-safe.html-*/-#pragma message "randomVector is not thread-safe in OSX and FreeBSD"-#endif -#if defined (__APPLE__) || (__FreeBSD__) || defined(_WIN32) || defined(WIN32)-/* Windows use thread-safe random- See: http://stackoverflow.com/questions/143108/is-windows-rand-s-thread-safe+ TODO: As mentioned in the linked article, this could be fixed:+ "the best solution for truly portable applications is to include+ your own random number generator implementation,+ and not rely on the system's C library". */-#if defined (__APPLE__) || (__FreeBSD__)+#pragma message "randomVector is not thread-safe in OSX and FreeBSD or with NO_RANDOM_R; this likely leads to less random numbers at worst; see http://www.evanjones.ca/random-thread-safe.html" inline double urandom() { /* the probalility of matching will be theoretically p^3(in fact, it is not)
src/Internal/CG.hs view
@@ -179,9 +179,9 @@ where m1 = convomat n k m2 = map (((+n) *** id) *** id) m1- + testb n = vect $ take n $ cycle ([0..10]++[9,8..1])- + denseSolve a = flatten . linearSolveLS a . asColumn -- mkDiag v = mkDiagR (dim v) (dim v) v
src/Internal/Container.hs view
@@ -3,10 +3,7 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-}-{-# LANGUAGE UndecidableInstances #-} -{-# OPTIONS_GHC -fno-warn-simplifiable-class-constraints #-}- ----------------------------------------------------------------------------- -- | -- Module : Internal.Container@@ -40,10 +37,12 @@ {- | Creates a real vector containing a range of values: >>> 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: @@ -87,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 #>@@ -136,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. @@ -167,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: --@@ -184,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' @@ -204,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
src/Internal/Conversion.hs view
@@ -3,7 +3,6 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-}-{-# LANGUAGE UndecidableInstances #-} ----------------------------------------------------------------------------- -- |
src/Internal/Convolution.hs view
@@ -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
src/Internal/Devel.hs view
@@ -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 --------------------------------------------------------------------------------
src/Internal/Element.hs view
@@ -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)
src/Internal/LAPACK.hs view
@@ -18,6 +18,8 @@ module Internal.LAPACK where +import Data.Bifunctor (first)+ import Internal.Devel import Internal.Vector import Internal.Matrix hiding ((#), (#!))@@ -234,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 @@ -298,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@@ -446,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
src/Internal/Matrix.hs view
@@ -160,7 +160,8 @@ {- | 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
src/Internal/Modular.hs view
@@ -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)
src/Internal/Sparse.hs view
@@ -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. @@ -124,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)@@ -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@@ -172,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 !#>
src/Internal/Static.hs view
@@ -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 @@ -324,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 @@ -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
src/Internal/Util.hs view
@@ -99,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@@ -378,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@@ -415,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 @@ -630,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)@@ -910,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] ]-
src/Internal/Vector.hs view
@@ -116,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@@ -135,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@@ -169,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@@ -191,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]
src/Internal/Vectorized.hs view
@@ -417,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
src/Numeric/LinearAlgebra.hs view
@@ -29,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)@@ -130,8 +131,8 @@ leftSV, rightSV, -- * Eigendecomposition- eig, eigSH,- eigenvalues, eigenvaluesSH,+ eig, geig, eigSH,+ eigenvalues, geigenvalues, eigenvaluesSH, geigSH, -- * QR
src/Numeric/LinearAlgebra/Devel.hs view
@@ -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
src/Numeric/LinearAlgebra/Static.hs view
@@ -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(..))@@ -92,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) @@ -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)@@ -383,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)@@ -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)