hmatrix 0.20.0.0 → 0.20.1

raw patch · 15 files changed

+668/−321 lines, 15 filesdep +primitivedep ~basePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: primitive

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• hmatrix.cabal +19/−8
• src/Internal/Algorithms.hs +4/−2
• src/Internal/C/lapack-aux.c +461/−234
• src/Internal/C/vector-aux.c +23/−10
• src/Internal/CG.hs +2/−2
• src/Internal/Container.hs +1/−4
• src/Internal/Conversion.hs +0/−1
• src/Internal/Devel.hs +17/−4
• src/Internal/Element.hs +1/−2
• src/Internal/Modular.hs +5/−9
• src/Internal/Sparse.hs +88/−28
• src/Internal/Static.hs +6/−7
• src/Internal/Util.hs +1/−2
• src/Numeric/LinearAlgebra/Devel.hs +3/−3
• src/Numeric/LinearAlgebra/Static.hs +37/−5

@@ -1,5 +1,5 @@ Name:               hmatrix-Version:            0.20.0.0+Version:            0.20.1 License:            BSD3 License-file:       LICENSE Author:             Alberto Ruiz@@ -18,7 +18,7 @@ Category:           Math tested-with:        GHC==8.2, GHC==8.4 -cabal-version:      >=1.8+cabal-version:       >=1.18  build-type:         Simple @@ -29,22 +29,30 @@ flag openblas     description:    Link with OpenBLAS (https://github.com/xianyi/OpenBLAS) optimized libraries.     default:        False-    manual: True+    manual:         True  flag disable-default-paths     description:    When enabled, don't add default hardcoded include/link dirs by default. Needed for hermetic builds like in nix.     default:        False-    manual: True+    manual:         True +flag no-random_r+    description:    When enabled, don't depend on the random_r() C function.+    default:        False+    manual:         True+ library -    Build-Depends:      base >= 4.9 && < 5,+    default-language:   Haskell2010++    Build-Depends:      base >= 4.8 && < 5,                         binary,                         array,                         deepseq,                         random,                         split,                         bytestring,+                        primitive,                         storable-complex,                         semigroups,                         vector >= 0.11@@ -84,7 +92,7 @@                         src/Internal/C/vector-aux.c  -    extensions:         ForeignFunctionInterface+    other-extensions:   ForeignFunctionInterface      ghc-options:        -Wall                         -fno-warn-missing-signatures@@ -99,6 +107,9 @@         cc-options:     -msse2  +    if flag(no-random_r)+        cc-options: -DNO_RANDOM_R+     if os(OSX)         if flag(openblas)             if !flag(disable-default-paths)@@ -128,10 +139,11 @@             extra-lib-dirs: /usr/local/lib             include-dirs: /usr/local/include         extra-libraries: gfortran+        extra-lib-dirs: /usr/local/lib/gcc9 /usr/local/lib/gcc8 /usr/local/lib/gcc7      if os(windows)         if flag(openblas)-            extra-libraries:    libopenblas, libgcc_s_seh-1, libgfortran, libquadmath-0+            extra-libraries:    openblas         else             extra-libraries:    blas lapack @@ -150,4 +162,3 @@ source-repository head     type:     git     location: https://github.com/albertoruiz/hmatrix-

@@ -1,7 +1,6 @@ {-# LANGUAGE FlexibleContexts, FlexibleInstances #-} {-# LANGUAGE CPP #-} {-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE TypeFamilies #-}  {-# OPTIONS_GHC -fno-warn-missing-signatures #-}@@ -27,6 +26,10 @@   UpLo(..) ) where +#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif+ import Internal.Vector import Internal.Matrix import Internal.Element@@ -1159,4 +1162,3 @@ --   for usage in 'chol', 'eigSH', etc. Only a triangular part of the matrix will be used. trustSym :: Matrix t -> Herm t trustSym x = (Herm x)-

@@ -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,9 +470,10 @@              vp,&n,              work, &lwork,              &res);-    CHECK(res,res);+    MARK(res,res);+     free(work);-    OK+    return ret; }  //////////////////// generalized real eigensystem ////////////@@ -425,7 +486,8 @@ 	integer *lwork, integer *info);  int eig_l_G(ODMAT(a), ODMAT(b), CVEC(alpha), DVEC(beta), ODMAT(vl), ODMAT(vr)) {-    integer n = ar;+    integer ret = 0;+    integer n   = ar;     REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE);     REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE);     char jobvl = vlp==NULL?'N':'V';@@ -443,9 +505,12 @@              vlp, &n, vrp, &n,              &ans, &lwork,              &res);+    CHECK(res,res);+     lwork = ceil(ans);     double * work = (double*)malloc(lwork*sizeof(double));     CHECK(!work,MEM);+     dggev_  (&jobvl,&jobvr,              &n,              ap,&n,bp,&n,@@ -453,9 +518,10 @@              vlp, &n, vrp, &n,              work, &lwork,              &res);-    CHECK(res,res);+    MARK(res,res);+     free(work);-    OK+    return ret; }  //////////////////// generalized complex eigensystem ////////////@@ -468,7 +534,8 @@     doublereal *rwork, integer *info);  int eig_l_GC(OCMAT(a), OCMAT(b), CVEC(alpha), CVEC(beta), OCMAT(vl), OCMAT(vr)) {-    integer n = ar;+    integer ret = 0;+    integer n   = ar;     REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE);     REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE);     char jobvl = vlp==NULL?'N':'V';@@ -476,7 +543,8 @@     char jobvr = vrp==NULL?'N':'V';     DEBUGMSG("eig_l_GC");     double *rwork = (double*) malloc(8*n*sizeof(double));-    CHECK(!rwork,MEM);+    UNWIND(!rwork,MEM,cleanup0);+     integer lwork = -1;     integer res;     // ask for optimal lwork@@ -488,9 +556,12 @@              vlp, &n, vrp, &n,              &ans, &lwork,              rwork, &res);+    UNWIND(res,res,cleanup1);+     lwork = ceil(ans.r);     doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));-    CHECK(!work,MEM);+    UNWIND(!work,MEM,cleanup1);+     zggev_  (&jobvl,&jobvr,              &n,              ap,&n,bp,&n,@@ -498,9 +569,14 @@              vlp, &n, vrp, &n,              work, &lwork,              rwork, &res);-    CHECK(res,res);+    UNWIND(res,res,cleanup2);++cleanup2:     free(work);-    OK+cleanup1:+    free(rwork);+cleanup0:+    return ret; }  //////////////////// symmetric real eigensystem ////////////@@ -510,7 +586,8 @@ 	integer *info);  int eig_l_S(int wantV,DVEC(s),ODMAT(v)) {-    integer n = sn;+    integer ret = 0;+    integer n   = sn;     REQUIRES(vr==n && vc==n, BAD_SIZE);     char jobz = wantV?'V':'N';     DEBUGMSG("eig_l_S");@@ -524,17 +601,21 @@              sp,              &ans, &lwork,              &res);+    CHECK(res,res);+     lwork = ceil(ans);     double * work = (double*)malloc(lwork*sizeof(double));     CHECK(!work,MEM);+     dsyev_  (&jobz,&uplo,              &n,vp,&n,              sp,              work, &lwork,              &res);-    CHECK(res,res);+    MARK(res,res);+     free(work);-    OK+    return ret; }  //////////////////// hermitian complex eigensystem ////////////@@ -544,12 +625,15 @@ 	doublereal *rwork, integer *info);  int eig_l_H(int wantV,DVEC(s),OCMAT(v)) {-    integer n = sn;+    integer ret = 0;+    integer n   = sn;+     REQUIRES(vr==n && vc==n, BAD_SIZE);     char jobz = wantV?'V':'N';     DEBUGMSG("eig_l_H");     double *rwork = (double*) malloc((3*n-2)*sizeof(double));-    CHECK(!rwork,MEM);+    UNWIND(!rwork,MEM,cleanup0);+     integer lwork = -1;     char uplo = 'U';     integer res;@@ -561,19 +645,26 @@              &ans, &lwork,              rwork,              &res);+    UNWIND(res,res,cleanup1);+     lwork = ceil(ans.r);     doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));-    CHECK(!work,MEM);+    UNWIND(!work,MEM,cleanup1);+     zheev_  (&jobz,&uplo,              &n,vp,&n,              sp,              work, &lwork,              rwork,              &res);-    CHECK(res,res);+    UNWIND(res,res,cleanup2);++cleanup2:     free(work);+cleanup1:     free(rwork);-    OK+cleanup0:+    return ret; }  //////////////////// general real linear system ////////////@@ -582,23 +673,25 @@ 	*lda, integer *ipiv, doublereal *b, integer *ldb, integer *info);  int linearSolveR_l(ODMAT(a),ODMAT(b)) {-    integer n = ar;+    integer ret  = 0;+    integer n    = ar;     integer nhrs = bc;+     REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);     DEBUGMSG("linearSolveR_l");     integer * ipiv = (integer*)malloc(n*sizeof(integer));+    CHECK(!ipiv,MEM);+     integer res;     dgesv_  (&n,&nhrs,              ap, &n,              ipiv,              bp, &n,              &res);-    if(res>0) {-        return SINGULAR;-    }-    CHECK(res,res);+    CONVERGED(res,SINGULAR);+     free(ipiv);-    OK+    return ret; }  //////////////////// general complex linear system ////////////@@ -608,23 +701,25 @@ 	info);  int linearSolveC_l(OCMAT(a),OCMAT(b)) {-    integer n = ar;+    integer ret  = 0;+    integer n    = ar;     integer nhrs = bc;+     REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);     DEBUGMSG("linearSolveC_l");     integer * ipiv = (integer*)malloc(n*sizeof(integer));+    CHECK(!ipiv,MEM);+     integer res;     zgesv_  (&n,&nhrs,              ap, &n,              ipiv,              bp, &n,              &res);-    if(res>0) {-        return SINGULAR;-    }-    CHECK(res,res);+    CONVERGED(res,SINGULAR);+     free(ipiv);-    OK+    return ret; }  //////// symmetric positive definite real linear system using Cholesky ////////////@@ -768,27 +863,37 @@             integer *info);  int triDiagSolveR_l(DVEC(dl), DVEC(d), DVEC(du), ODMAT(b)) {-    integer n = dn;+    integer ret  = 0;+    integer n    = dn;     integer nhrs = bc;     REQUIRES(n >= 1 && dln == dn - 1 && dun == dn - 1 && br == n, BAD_SIZE);     DEBUGMSG("triDiagSolveR_l");     integer res;     integer* ipiv = (integer*)malloc(n*sizeof(integer));-    double* du2  = (double*)malloc((n - 2)*sizeof(double));+    UNWIND(!ipiv,MEM,cleanup0);++    double* du2 = (double*)malloc((n - 2)*sizeof(double));+    UNWIND(!du2,MEM,cleanup1);+     dgttrf_ (&n,              dlp, dp, dup, du2,              ipiv,              &res);-    CHECK(res,res);+    UNWIND(res,res,cleanup2);+     dgttrs_ ("N",              &n,&nhrs,              dlp, dp, dup, du2,              ipiv, bp, &n,              &res);-    CHECK(res,res);-    free(ipiv);+    UNWIND(res,res,cleanup2);++cleanup2:     free(du2);-    OK+cleanup1:+    free(ipiv);+cleanup0:+    return ret; }  //////// tridiagonal complex linear system ////////////@@ -804,27 +909,37 @@             integer *info);  int triDiagSolveC_l(CVEC(dl), CVEC(d), CVEC(du), OCMAT(b)) {-    integer n = dn;+    integer ret  = 0;+    integer n    = dn;     integer nhrs = bc;     REQUIRES(n >= 1 && dln == dn - 1 && dun == dn - 1 && br == n, BAD_SIZE);     DEBUGMSG("triDiagSolveC_l");     integer res;     integer* ipiv = (integer*)malloc(n*sizeof(integer));+    UNWIND(!ipiv,MEM,cleanup0);+     doublecomplex* du2 = (doublecomplex*)malloc((n - 2)*sizeof(doublecomplex));+    UNWIND(!du2,MEM,cleanup1);+     zgttrf_ (&n,              dlp, dp, dup, du2,              ipiv,              &res);-    CHECK(res,res);+    UNWIND(res,res,cleanup2);+     zgttrs_ ("N",              &n,&nhrs,              dlp, dp, dup, du2,              ipiv, bp, &n,              &res);-    CHECK(res,res);-    free(ipiv);+    UNWIND(res,res,cleanup2);++cleanup2:     free(du2);-    OK+cleanup1:+    free(ipiv);+cleanup0:+    return ret; }  //////////////////// least squares real linear system ////////////@@ -834,10 +949,11 @@ 	doublereal *work, integer *lwork, integer *info);  int linearSolveLSR_l(ODMAT(a),ODMAT(b)) {-    integer m = ar;-    integer n = ac;+    integer ret  = 0;+    integer m    = ar;+    integer n    = ac;     integer nrhs = bc;-    integer ldb = bXc;+    integer ldb  = bXc;     REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE);     DEBUGMSG("linearSolveLSR_l");     integer res;@@ -848,19 +964,21 @@              bp,&ldb,              &ans,&lwork,              &res);+    CHECK(res,res);+     lwork = ceil(ans);     double * work = (double*)malloc(lwork*sizeof(double));+    CHECK(!work,MEM);+     dgels_  ("N",&m,&n,&nrhs,              ap,&m,              bp,&ldb,              work,&lwork,              &res);-    if(res>0) {-        return SINGULAR;-    }-    CHECK(res,res);+    CONVERGED(res,SINGULAR);+     free(work);-    OK+    return ret; }  //////////////////// least squares complex linear system ////////////@@ -870,10 +988,11 @@ 	doublecomplex *work, integer *lwork, integer *info);  int linearSolveLSC_l(OCMAT(a),OCMAT(b)) {-    integer m = ar;-    integer n = ac;+    integer ret  = 0;+    integer m    = ar;+    integer n    = ac;     integer nrhs = bc;-    integer ldb = bXc;+    integer ldb  = bXc;     REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE);     DEBUGMSG("linearSolveLSC_l");     integer res;@@ -884,19 +1003,21 @@              bp,&ldb,              &ans,&lwork,              &res);+    CHECK(res,res);+     lwork = ceil(ans.r);     doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));+    CHECK(!work,MEM);+     zgels_  ("N",&m,&n,&nrhs,              ap,&m,              bp,&ldb,              work,&lwork,              &res);-    if(res>0) {-        return SINGULAR;-    }-    CHECK(res,res);+    CONVERGED(res,SINGULAR);+     free(work);-    OK+    return ret; }  //////////////////// least squares real linear system using SVD ////////////@@ -907,13 +1028,17 @@ 	integer *info);  int linearSolveSVDR_l(double rcond,ODMAT(a),ODMAT(b)) {-    integer m = ar;-    integer n = ac;+    integer ret  = 0;+    integer m    = ar;+    integer n    = ac;     integer nrhs = bc;-    integer ldb = bXc;+    integer ldb  = bXc;     REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE);     DEBUGMSG("linearSolveSVDR_l");-    double*S = (double*)malloc(MIN(m,n)*sizeof(double));++    double * S   = (double*)malloc(MIN(m,n)*sizeof(double));+    UNWIND(!S,MEM,cleanup0);+     integer res;     integer lwork = -1;     integer rank;@@ -925,8 +1050,12 @@              &rcond,&rank,              &ans,&lwork,              &res);+    UNWIND(res,res,cleanup1);+     lwork = ceil(ans);     double * work = (double*)malloc(lwork*sizeof(double));+    UNWIND(!work,MEM,cleanup1);+     dgelss_  (&m,&n,&nrhs,              ap,&m,              bp,&ldb,@@ -934,13 +1063,15 @@              &rcond,&rank,              work,&lwork,              &res);-    if(res>0) {-        return NOCONVER;-    }-    CHECK(res,res);++    CONVERGED(res,NOCONVER);+     free(work);+cleanup1:     free(S);-    OK+cleanup0:+    return ret;+ }  //////////////////// least squares complex linear system using SVD ////////////@@ -952,14 +1083,20 @@     integer *info);  int linearSolveSVDC_l(double rcond, OCMAT(a),OCMAT(b)) {-    integer m = ar;-    integer n = ac;+    integer ret  = 0;+    integer m    = ar;+    integer n    = ac;     integer nrhs = bc;-    integer ldb = bXc;+    integer ldb  = bXc;     REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE);     DEBUGMSG("linearSolveSVDC_l");+     double*S = (double*)malloc(MIN(m,n)*sizeof(double));+    UNWIND(!S,MEM,cleanup0);+     double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double));+    UNWIND(!S,MEM,cleanup1);+     integer res;     integer lwork = -1;     integer rank;@@ -972,8 +1109,12 @@              &ans,&lwork,              RWORK,              &res);+    UNWIND(res,res,cleanup2);+     lwork = ceil(ans.r);     doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));+    UNWIND(!work,MEM,cleanup2);+     zgelss_  (&m,&n,&nrhs,              ap,&m,              bp,&ldb,@@ -982,14 +1123,16 @@              work,&lwork,              RWORK,              &res);-    if(res>0) {-        return NOCONVER;-    }-    CHECK(res,res);+    CONVERGED(res,NOCONVER);+     free(work);+cleanup2:     free(RWORK);+cleanup1:     free(S);-    OK+cleanup0:+    return ret;+ }  //////////////////// Cholesky factorization /////////////////////////@@ -1042,36 +1185,43 @@ 	lda, doublereal *tau, doublereal *work, integer *info);  int qr_l_R(DVEC(tau), ODMAT(r)) {-    integer m = rr;-    integer n = rc;-    integer mn = MIN(m,n);+    integer ret = 0;+    integer m   = rr;+    integer n   = rc;+    integer mn  = MIN(m,n);     REQUIRES(m>=1 && n >=1 && taun == mn, BAD_SIZE);     DEBUGMSG("qr_l_R");     double *WORK = (double*)malloc(n*sizeof(double));     CHECK(!WORK,MEM);+     integer res;     dgeqr2_ (&m,&n,rp,&m,taup,WORK,&res);-    CHECK(res,res);+    MARK(res,res);+     free(WORK);-    OK+    return ret; }  int zgeqr2_(integer *m, integer *n, doublecomplex *a, 	integer *lda, doublecomplex *tau, doublecomplex *work, integer *info);  int qr_l_C(CVEC(tau), OCMAT(r)) {-    integer m = rr;-    integer n = rc;-    integer mn = MIN(m,n);+    integer ret = 0;+    integer m   = rr;+    integer n   = rc;+    integer mn  = MIN(m,n);     REQUIRES(m>=1 && n >=1 && taun == mn, BAD_SIZE);     DEBUGMSG("qr_l_C");+     doublecomplex *WORK = (doublecomplex*)malloc(n*sizeof(doublecomplex));     CHECK(!WORK,MEM);+     integer res;     zgeqr2_ (&m,&n,rp,&m,taup,WORK,&res);-    CHECK(res,res);+    MARK(res,res);+     free(WORK);-    OK+    return ret; }  int dorgqr_(integer *m, integer *n, integer *k, doublereal *@@ -1079,18 +1229,21 @@ 	integer *info);  int c_dorgqr(KDVEC(tau), ODMAT(r)) {-    integer m = rr;-    integer n = MIN(rc,rr);-    integer k = taun;+    integer ret = 0;+    integer m   = rr;+    integer n   = MIN(rc,rr);+    integer k   = taun;     DEBUGMSG("c_dorgqr");     integer lwork = 8*n; // FIXME     double *WORK = (double*)malloc(lwork*sizeof(double));     CHECK(!WORK,MEM);+     integer res;     dorgqr_ (&m,&n,&k,rp,&m,(double*)taup,WORK,&lwork,&res);-    CHECK(res,res);+    MARK(res,res);+     free(WORK);-    OK+    return ret; }  int zungqr_(integer *m, integer *n, integer *k,@@ -1098,18 +1251,21 @@ 	work, integer *lwork, integer *info);  int c_zungqr(KCVEC(tau), OCMAT(r)) {-    integer m = rr;-    integer n = MIN(rc,rr);-    integer k = taun;+    integer ret = 0;+    integer m   = rr;+    integer n   = MIN(rc,rr);+    integer k   = taun;     DEBUGMSG("z_ungqr");     integer lwork = 8*n; // FIXME     doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));     CHECK(!WORK,MEM);+     integer res;     zungqr_ (&m,&n,&k,rp,&m,(doublecomplex*)taup,WORK,&lwork,&res);-    CHECK(res,res);+    MARK(res,res);+     free(WORK);-    OK+    return ret; }  @@ -1120,20 +1276,23 @@ 	integer *lwork, integer *info);  int hess_l_R(DVEC(tau), ODMAT(r)) {-    integer m = rr;-    integer n = rc;-    integer mn = MIN(m,n);+    integer ret = 0;+    integer m   = rr;+    integer n   = rc;+    integer mn  = MIN(m,n);     REQUIRES(m>=1 && n == m && taun == mn-1, BAD_SIZE);     DEBUGMSG("hess_l_R");     integer lwork = 5*n; // FIXME     double *WORK = (double*)malloc(lwork*sizeof(double));     CHECK(!WORK,MEM);+     integer res;     integer one = 1;     dgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res);-    CHECK(res,res);+    MARK(res,res);+     free(WORK);-    OK+    return ret; }  @@ -1142,20 +1301,23 @@ 	work, integer *lwork, integer *info);  int hess_l_C(CVEC(tau), OCMAT(r)) {-    integer m = rr;-    integer n = rc;-    integer mn = MIN(m,n);+    integer ret = 0;+    integer m   = rr;+    integer n   = rc;+    integer mn  = MIN(m,n);     REQUIRES(m>=1 && n == m && taun == mn-1, BAD_SIZE);     DEBUGMSG("hess_l_C");     integer lwork = 5*n; // FIXME     doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));     CHECK(!WORK,MEM);+     integer res;     integer one = 1;     zgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res);-    CHECK(res,res);+    MARK(res,res);+     free(WORK);-    OK+    return ret; }  //////////////////// Schur factorization /////////////////////////@@ -1166,28 +1328,35 @@ 	integer *lwork, logical *bwork, integer *info);  int schur_l_R(ODMAT(u), ODMAT(s)) {-    integer m = sr;-    integer n = sc;+    integer ret = 0;+    integer m   = sr;+    integer n   = sc;     REQUIRES(m>=1 && n==m && ur==n && uc==n, BAD_SIZE);     DEBUGMSG("schur_l_R");     integer lwork = 6*n; // FIXME     double *WORK = (double*)malloc(lwork*sizeof(double));-    double *WR = (double*)malloc(n*sizeof(double));-    double *WI = (double*)malloc(n*sizeof(double));+    UNWIND(!WORK,MEM,cleanup0);+    double *WR   = (double*)malloc(n*sizeof(double));+    UNWIND(!WORK,MEM,cleanup1);+    double *WI   = (double*)malloc(n*sizeof(double));+    UNWIND(!WORK,MEM,cleanup2);     // WR and WI not really required in this call     logical *BWORK = (logical*)malloc(n*sizeof(logical));+    UNWIND(!BWORK,MEM,cleanup3);     integer res;     integer sdim;     dgees_ ("V","N",NULL,&n,sp,&n,&sdim,WR,WI,up,&n,WORK,&lwork,BWORK,&res);-    if(res>0) {-        return NOCONVER;-    }-    CHECK(res,res);-    free(WR);-    free(WI);+    CONVERGED(res,NOCONVER);+     free(BWORK);+cleanup3:+    free(WI);+cleanup2:+    free(WR);+cleanup1:     free(WORK);-    OK+cleanup0:+    return ret; }  @@ -1197,29 +1366,40 @@ 	doublereal *rwork, logical *bwork, integer *info);  int schur_l_C(OCMAT(u), OCMAT(s)) {-    integer m = sr;-    integer n = sc;+    integer ret = 0;+    integer m   = sr;+    integer n   = sc;     REQUIRES(m>=1 && n==m && ur==n && uc==n, BAD_SIZE);     DEBUGMSG("schur_l_C");     integer lwork = 6*n; // FIXME     doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));-    doublecomplex *W = (doublecomplex*)malloc(n*sizeof(doublecomplex));+    UNWIND(!WORK,MEM,cleanup0);++    doublecomplex *W    = (doublecomplex*)malloc(n*sizeof(doublecomplex));+    UNWIND(!W,MEM,cleanup1);+     // W not really required in this call     logical *BWORK = (logical*)malloc(n*sizeof(logical));-    double *RWORK = (double*)malloc(n*sizeof(double));+    UNWIND(!BWORK,MEM,cleanup2);++    double  *RWORK = (double*)malloc(n*sizeof(double));+    UNWIND(!RWORK,MEM,cleanup3);     integer res;     integer sdim;     zgees_ ("V","N",NULL,&n,sp,&n,&sdim,W,                             up,&n,                             WORK,&lwork,RWORK,BWORK,&res);-    if(res>0) {-        return NOCONVER;-    }-    CHECK(res,res);-    free(W);+    CONVERGED(res,NOCONVER);++    free(RWORK);+cleanup3:     free(BWORK);+cleanup2:+    free(W);+cleanup1:     free(WORK);-    OK+cleanup0:+    return ret; }  //////////////////// LU factorization /////////////////////////@@ -1228,24 +1408,30 @@ 	lda, integer *ipiv, integer *info);  int lu_l_R(DVEC(ipiv), ODMAT(r)) {-    integer m = rr;-    integer n = rc;-    integer mn = MIN(m,n);+    integer ret = 0;+    integer m   = rr;+    integer n   = rc;+    integer mn  = MIN(m,n);     REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);     DEBUGMSG("lu_l_R");     integer* auxipiv = (integer*)malloc(mn*sizeof(integer));+    UNWIND(!auxipiv,MEM,cleanup0);+     integer res;     dgetrf_ (&m,&n,rp,&m,auxipiv,&res);     if(res>0) {         res = 0; // FIXME     }-    CHECK(res,res);-    int k;-    for (k=0; k<mn; k++) {+    UNWIND(res,res,cleanup1);++    for (int k=0; k<mn; k++) {         ipivp[k] = auxipiv[k];     }++cleanup1:     free(auxipiv);-    OK+cleanup0:+    return ret; }  @@ -1253,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; }  @@ -1281,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; }  @@ -1306,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; }  @@ -1332,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;@@ -1343,15 +1547,22 @@     dsytrf_ ("L",&n,rp,&lda,auxipiv,&ans,&lwork,&res);     lwork = ceil(ans);     doublereal* work = (doublereal*)malloc(lwork*sizeof(doublereal));+    UNWIND(!work,MEM,cleanup1);+     dsytrf_ ("L",&n,rp,&lda,auxipiv,work,&lwork,&res);-    CHECK(res,res);+    UNWIND(res,res,cleanup2);+     int k;     for (k=0; k<n; k++) {         ipivp[k] = auxipiv[k];     }-    free(auxipiv);++cleanup2:     free(work);-    OK+cleanup1:+    free(auxipiv);+cleanup0:+    return ret; }  @@ -1359,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;@@ -1370,15 +1585,21 @@     zhetrf_ ("L",&n,rp,&lda,auxipiv,&ans,&lwork,&res);     lwork = ceil(ans.r);     doublecomplex* work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));+    UNWIND(!work,MEM,cleanup1);+     zhetrf_ ("L",&n,rp,&lda,auxipiv,work,&lwork,&res);-    CHECK(res,res);+    UNWIND(res,res,cleanup2);     int k;     for (k=0; k<n; k++) {         ipivp[k] = auxipiv[k];     }-    free(auxipiv);++cleanup2:     free(work);-    OK+cleanup1:+    free(auxipiv);+cleanup0:+    return ret;  } @@ -1388,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; }  @@ -1412,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; }  

@@ -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)

@@ -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

@@ -3,10 +3,7 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-}-{-# LANGUAGE UndecidableInstances #-} -{-# OPTIONS_GHC -fno-warn-simplifiable-class-constraints #-}- ----------------------------------------------------------------------------- -- | -- Module      :  Internal.Container@@ -190,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' 

@@ -3,7 +3,6 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-}-{-# LANGUAGE UndecidableInstances #-}  ----------------------------------------------------------------------------- -- |

@@ -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  -------------------------------------------------------------------------------- 

@@ -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)

@@ -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)

@@ -1,11 +1,11 @@-{-# LANGUAGE RecordWildCards #-}-{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE BangPatterns #-} {-# LANGUAGE FlexibleInstances #-}--{-# OPTIONS_GHC -fno-warn-missing-signatures #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE RecordWildCards #-}  module Internal.Sparse(-    GMatrix(..), CSR(..), mkCSR, fromCSR,+    GMatrix(..), CSR(..), mkCSR, fromCSR, impureCSR,     mkSparse, mkDiagR, mkDense,     AssocMatrix,     toDense,@@ -16,10 +16,12 @@ import Internal.Matrix import Internal.Numeric import qualified Data.Vector.Storable as V-import Data.Function(on)+import qualified Data.Vector.Storable.Mutable as M import Control.Arrow((***))-import Control.Monad(when)-import Data.List(groupBy, sort)+import Control.Monad(when, foldM)+import Control.Monad.ST (runST)+import Control.Monad.Primitive (PrimMonad)+import Data.List(sort) import Foreign.C.Types(CInt(..))  import Internal.Devel@@ -27,10 +29,7 @@ import Foreign(Ptr) import Text.Printf(printf) -infixl 0 ~!~-c ~!~ msg = when c (error msg)--type AssocMatrix = [((Int,Int),Double)]+type AssocMatrix = [(IndexOf Matrix, Double)]  data CSR = CSR         { csrVals  :: Vector Double@@ -49,23 +48,79 @@         } deriving Show  +-- | Produce a CSR sparse matrix from a association matrix. mkCSR :: AssocMatrix -> CSR-mkCSR sm' = CSR{..}+mkCSR ms =+  runST $ impureCSR runFold $ sort ms+    where+  runFold next initialise xtract as0 = do+    i0  <- initialise+    acc <- foldM next i0 as0+    xtract acc++-- | Produce a CSR sparse matrix by applying a generic folding function.+--+--   This allows one to build a CSR from an effectful streaming source+--   when combined with libraries like pipes, io-streams, or streaming.+--+--   For example+--+--   > impureCSR Pipes.Prelude.foldM :: PrimMonad m => Producer AssocEntry m () -> m CSR+--   > impureCSR Streaming.Prelude.foldM :: PrimMonad m => Stream (Of AssocEntry) m r -> m (Of CSR r)+--+impureCSR+    :: PrimMonad m+    => (forall x . (x -> (IndexOf Matrix, Double) -> m x) -> m x -> (x -> m CSR) -> r)+    -> r+impureCSR f = f next begin done   where-    sm = sort sm'-    rws = map ((fromList *** fromList)-              . unzip-              . map ((succ.fi.snd) *** id)-              )-        . groupBy ((==) `on` (fst.fst))-        $ sm-    rszs = map (fi . dim . fst) rws-    csrRows = fromList (scanl (+) 1 rszs)-    csrVals = vjoin (map snd rws)-    csrCols = vjoin (map fst rws)-    csrNRows = dim csrRows - 1-    csrNCols = fromIntegral (V.maximum csrCols)+    sfi = succ . fi+    begin = do+      mv <- M.unsafeNew 64+      mr <- M.unsafeNew 64+      mc <- M.unsafeNew 64+      return (mv, mr, mc, 0, 0, 0, -1) +    next (!mv, !mr, !mc, !idxVC, !idxR, !maxC, !curRow) ((r,c),d) = do+      when (r < curRow) $+        error (printf "impureCSR: row %i specified after %i" r curRow)++      let lenVC = M.length mv+          lenR  = M.length mr+          maxC' = max maxC c++      (mv', mc') <-+        if idxVC >= lenVC then do+          mv' <- M.unsafeGrow mv lenVC+          mc' <- M.unsafeGrow mc lenVC+          return (mv', mc')+        else+          return (mv, mc)++      mr' <-+        if idxR >= lenR - 1 then+          M.unsafeGrow mr lenR+        else+          return mr++      M.unsafeWrite mc' idxVC (sfi c)+      M.unsafeWrite mv' idxVC d++      idxR' <-+        foldM+          (\idxR' _ -> idxR' + 1 <$ M.unsafeWrite mr' idxR' (sfi idxVC))+          idxR [1 .. (r-curRow)]++      return (mv', mr', mc', idxVC + 1, idxR', maxC', r)++    done (!mv, !mr, !mc, !idxVC, !idxR, !maxC, !curR) = do+      M.unsafeWrite mr idxR (sfi idxVC)+      vv <- V.unsafeFreeze (M.unsafeTake idxVC mv)+      vc <- V.unsafeFreeze (M.unsafeTake idxVC mc)+      vr <- V.unsafeFreeze (M.unsafeTake (idxR + 1)  mr)+      return $ CSR vv vc vr (succ curR) (succ maxC)++ {- | General matrix with specialized internal representations for      dense, sparse, diagonal, banded, and constant elements. @@ -129,6 +184,7 @@     nCols = csrNCols  +mkDiagR :: Int -> Int -> Vector Double -> GMatrix mkDiagR r c v     | dim v <= min r c = Diag{..}     | otherwise = error $ printf "mkDiagR: incorrect sizes (%d,%d) [%d]" r c (dim v)@@ -144,13 +200,17 @@  gmXv :: GMatrix -> Vector Double -> Vector Double gmXv SparseR { gmCSR = CSR{..}, .. } v = unsafePerformIO $ do-    dim v /= nCols ~!~ printf "gmXv (CSR): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)+    when (dim v /= nCols) $+      error (printf "gmXv (CSR): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v))+     r <- createVector nRows     (csrVals # csrCols # csrRows # v #! r) c_smXv #|"CSRXv"     return r  gmXv SparseC { gmCSC = CSC{..}, .. } v = unsafePerformIO $ do-    dim v /= nCols ~!~ printf "gmXv (CSC): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)+    when (dim v /= nCols) $+      error (printf "gmXv (CSC): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v))+     r <- createVector nRows     (cscVals # cscRows # cscCols # v #! r) c_smTXv #|"CSCXv"     return r

@@ -16,7 +16,6 @@ {-# LANGUAGE DeriveGeneric #-}  {-# OPTIONS_GHC -fno-warn-missing-signatures #-}-{-# OPTIONS_GHC -fno-warn-simplifiable-class-constraints #-}  {- | Module      :  Internal.Static@@ -131,9 +130,9 @@   where     du = fromIntegral . natVal $ (undefined :: Proxy n)     dv = fromIntegral . natVal $ (undefined :: Proxy m)-    u' | du > 1 && LA.size u == 1 = LA.konst (u D.@> 0) du+    u' | du /= 1 && LA.size u == 1 = LA.konst (u D.@> 0) du        | otherwise = u-    v' | dv > 1 && LA.size v == 1 = LA.konst (v D.@> 0) dv+    v' | dv /= 1 && LA.size v == 1 = LA.konst (v D.@> 0) dv        | otherwise = v  @@ -367,7 +366,7 @@     negate = lift1F negate     fromInteger x = Dim (fromInteger x) -instance (Num (Vector t), Num (Matrix t), Fractional t, Numeric t) => Fractional (Dim n (Vector t))+instance (Num (Vector t), Fractional t, Numeric t) => Fractional (Dim n (Vector t))   where     fromRational x = Dim (fromRational x)     (/) = lift2F (/)@@ -392,7 +391,7 @@     pi    = Dim pi  -instance (Num (Matrix t), Numeric t) => Num (Dim m (Dim n (Matrix t)))+instance (Num (Vector t), Numeric t) => Num (Dim m (Dim n (Matrix t)))   where     (+) = (lift2F . lift2F) (+)     (*) = (lift2F . lift2F) (*)@@ -402,12 +401,12 @@     negate = (lift1F . lift1F) negate     fromInteger x = Dim (Dim (fromInteger x)) -instance (Num (Vector t), Num (Matrix t), Fractional t, Numeric t) => Fractional (Dim m (Dim n (Matrix t)))+instance (Num (Vector t), Fractional t, Numeric t) => Fractional (Dim m (Dim n (Matrix t)))   where     fromRational x = Dim (Dim (fromRational x))     (/) = (lift2F.lift2F) (/) -instance (Num (Vector t), Floating (Matrix t), Fractional t, Numeric t) => Floating (Dim m (Dim n (Matrix t))) where+instance (Floating (Vector t), Floating t, Numeric t) => Floating (Dim m (Dim n (Matrix t))) where     sin   = (lift1F . lift1F) sin     cos   = (lift1F . lift1F) cos     tan   = (lift1F . lift1F) tan

@@ -417,7 +417,7 @@  instance Element t => Indexable (Matrix t) (Vector t)   where-    m!j = subVector (j*c) c (flatten m)+    m ! j = subVector (j*c) c (flatten m)       where         c = cols m @@ -912,4 +912,3 @@          , remap r (tr c) p == ep          , tr p ?? (PosCyc (idxs[-5,13]), Pos (idxs[3,7,1])) == (2><3) [35,75,15,33,73,13]          ]-

@@ -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

@@ -46,6 +46,13 @@     matrix,     -- * Complex     ℂ, C, M, Her, her, 𝑖,+    toComplex,+    fromComplex,+    complex,+    real,+    imag,+    sqMagnitude,+    magnitude,     -- * Products     (<>),(#>),(<.>),     -- * Linear Systems@@ -74,7 +81,9 @@     (<\>),fromList,takeDiag,svd,eig,eigSH,     eigenvalues,eigenvaluesSH,build,     qr,size,dot,chol,range,R,C,sym,mTm,unSym,-    randomVector,rand,randn,gaussianSample,uniformSample,meanCov)+    randomVector,rand,randn,gaussianSample,uniformSample,meanCov,+    toComplex, fromComplex, complex, real, magnitude+    ) import qualified Numeric.LinearAlgebra as LA import qualified Numeric.LinearAlgebra.Devel as LA import Data.Proxy(Proxy(..))@@ -243,7 +252,30 @@   where     takeDiag x = mkC (LA.takeDiag (extract x)) +-------------------------------------------------------------------------------- ++toComplex :: KnownNat n => (R n, R n) -> C n+toComplex (r,i) = mkC $ LA.toComplex (ud1 r, ud1 i)++fromComplex :: KnownNat n => C n -> (R n, R n)+fromComplex (C (Dim v)) = let (r,i) = LA.fromComplex v in (mkR r, mkR i)++complex :: KnownNat n => R n -> C n+complex r = mkC $ LA.toComplex (ud1 r, LA.konst 0 (size r))++real :: KnownNat n => C n -> R n+real = fst . fromComplex++imag :: KnownNat n => C n -> R n +imag = snd . fromComplex++sqMagnitude :: KnownNat n => C n -> R n+sqMagnitude c = let (r,i) = fromComplex c in r**2 + i**2++magnitude :: KnownNat n => C n -> R n+magnitude = sqrt . sqMagnitude+ --------------------------------------------------------------------------------  linSolve :: (KnownNat m, KnownNat n) => L m m -> L m n -> Maybe (L m n)@@ -595,9 +627,9 @@     a' = subVector 0 n a     b' = subVector 0 n b -mulR (isDiag -> Just (0,a,_)) (extract -> b) = mkL (asColumn a * takeRows (LA.size a) b)+-- mulR (isDiag -> Just (0,a,_)) (extract -> b) = mkL (asColumn a * takeRows (LA.size a) b) -mulR (extract -> a) (isDiag -> Just (0,b,_)) = mkL (takeColumns (LA.size b) a * asRow b)+-- mulR (extract -> a) (isDiag -> Just (0,b,_)) = mkL (takeColumns (LA.size b) a * asRow b)  mulR a b = mkL (extract a LA.<> extract b) @@ -658,9 +690,9 @@     a' = subVector 0 n a     b' = subVector 0 n b -mulC (isDiagC -> Just (0,a,_)) (extract -> b) = mkM (asColumn a * takeRows (LA.size a) b)+-- mulC (isDiagC -> Just (0,a,_)) (extract -> b) = mkM (asColumn a * takeRows (LA.size a) b) -mulC (extract -> a) (isDiagC -> Just (0,b,_)) = mkM (takeColumns (LA.size b) a * asRow b)+-- mulC (extract -> a) (isDiagC -> Just (0,b,_)) = mkM (takeColumns (LA.size b) a * asRow b)  mulC a b = mkM (extract a LA.<> extract b)