hmatrix 0.17.0.2 → 0.20.2
raw patch · 32 files changed
Files
- CHANGELOG +13/−1
- THANKS.md +34/−4
- hmatrix.cabal +44/−25
- src/Internal/Algorithms.hs +146/−12
- src/Internal/C/lapack-aux.c +684/−214
- src/Internal/C/vector-aux.c +122/−9
- src/Internal/CG.hs +4/−2
- src/Internal/Chain.hs +2/−0
- src/Internal/Container.hs +26/−21
- src/Internal/Conversion.hs +0/−1
- src/Internal/Convolution.hs +8/−2
- src/Internal/Devel.hs +20/−7
- src/Internal/Element.hs +21/−8
- src/Internal/Foreign.hs +0/−102
- src/Internal/IO.hs +13/−3
- src/Internal/LAPACK.hs +144/−26
- src/Internal/Matrix.hs +142/−20
- src/Internal/Modular.hs +22/−15
- src/Internal/Numeric.hs +16/−20
- src/Internal/Random.hs +2/−2
- src/Internal/ST.hs +10/−2
- src/Internal/Sparse.hs +94/−31
- src/Internal/Static.hs +89/−28
- src/Internal/Util.hs +27/−9
- src/Internal/Vector.hs +29/−22
- src/Internal/Vectorized.hs +56/−17
- src/Numeric/LinearAlgebra.hs +20/−6
- src/Numeric/LinearAlgebra/Devel.hs +4/−12
- src/Numeric/LinearAlgebra/HMatrix.hs +6/−0
- src/Numeric/LinearAlgebra/Static.hs +323/−29
- src/Numeric/Matrix.hs +20/−2
- src/Numeric/Vector.hs +12/−2
CHANGELOG view
@@ -1,3 +1,16 @@+0.18.0.0+--------++ * Many new functions and instances in the Static module++ * meanCov and gaussianSample use Herm type++ * thinQR, thinRQ++ * compactSVDTol++ * unitary changed to normalize, also admits Vector (Complex Double)+ 0.17.0.0 -------- @@ -275,4 +288,3 @@ * added NFData instances for Matrix and Vector. * liftVector, liftVector2 replaced by mapVector, zipVector.-
THANKS.md view
@@ -160,7 +160,7 @@ - Denis Laxalde separated the gsl tests from the base ones. - Dominic Steinitz (idontgetoutmuch) reported a bug in the static diagonal creation functions and- added Cholesky to Static.+ added Cholesky to Static. He also added support for tridiagonal matrix solver and fixed several bugs. - Dylan Thurston reported an error in the glpk documentation and ambiguity in the description of linearSolve.@@ -170,7 +170,8 @@ - Ian Ross reported the max/minIndex bug. -- Niklas Hambüchen improved the documentation.+- Niklas Hambüchen improved the documentation and fixed compilation with GHC-8.2+ adding type signatures. Added disable-default-paths flag. - "erdeszt" optimized "conv" using a direct vector reverse. @@ -192,7 +193,8 @@ - Matt Peddie wrote the interfaces to the interpolation and simulated annealing modules. -- "maxc01" solved uninstallability in FreeBSD and improved urandom+- "maxc01" solved uninstallability in FreeBSD, improved urandom, and fixed a Windows+ link error using rand_s. - "ntfrgl" added {take,drop}Last{Rows,Columns} and odeSolveVWith with generalized step control function and fixed link errors related to mod/mod_l.@@ -202,5 +204,33 @@ - Ilan Godik and Douglas McClean helped with Windows support. - Vassil Keremidchiev fixed the cabal options for OpenBlas, fixed several installation- issues, and added support for stack-based build.+ issues, and added support for stack-based build. He also added support for LTS 8.15+ under Windows.++- Greg Nwosu fixed arm compilation++- Patrik Jansson changed meanCov and gaussianSample to use Herm type. Fixed stack.yaml.++- Justin Le added NFData instances for Static types, added mapping and outer product+ methods to Domain, and many other functions to the Static module.++- Sidharth Kapur added Normed and numeric instances for several Static types,+fixed the CPP issue in cabal files, and made many other contributions.++- Matt Renaud improved the documentation.++- Joshua Moerman fixed cabal/stack flags for windows.++- Francesco Mazzoli, Niklas Hambüchen, Patrick Chilton, and Andras Slemmer+ discovered a serious and subtle bug in the wrapper helpers causing memory corruption.+ Andras Slemmer fixed the bug. Thank you all.++- Kevin Slagle implemented thinQR and thinRQ, much faster than the original qr,+ and added compactSVDTol. He also added an optimized reorderVector for hTensor.++- "fedeinthemix" suggested a better name and a more general type for unitary.++- Huw Campbell fixed a bug in equal.++- Hiromi Ishii fixed compilation problems for ghc-8.4
hmatrix.cabal view
@@ -1,11 +1,11 @@ Name: hmatrix-Version: 0.17.0.2+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.0+tested-with: GHC==8.10 -cabal-version: >=1.8+cabal-version: >=1.18 build-type: Simple @@ -29,10 +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++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,@@ -40,8 +52,10 @@ random, split, bytestring,+ primitive, storable-complex,- vector >= 0.8+ semigroups,+ vector >= 0.11 hs-source-dirs: src @@ -55,7 +69,6 @@ Internal.Devel Internal.Vectorized Internal.Matrix- Internal.Foreign Internal.ST Internal.IO Internal.Element@@ -79,13 +92,12 @@ src/Internal/C/vector-aux.c - extensions: ForeignFunctionInterface,- CPP+ other-extensions: ForeignFunctionInterface ghc-options: -Wall -fno-warn-missing-signatures -fno-warn-orphans- -fprof-auto+ -fno-prof-auto cc-options: -O4 -Wall @@ -94,43 +106,51 @@ if arch(i386) cc-options: -msse2 - cpp-options: -DBINARY + if flag(no-random_r)+ cc-options: -DNO_RANDOM_R+ if os(OSX) if flag(openblas)- extra-lib-dirs: /opt/local/lib/openblas/lib+ if !flag(disable-default-paths)+ extra-lib-dirs: /opt/local/lib/openblas/lib extra-libraries: openblas else extra-libraries: blas lapack - extra-lib-dirs: /opt/local/lib/- include-dirs: /opt/local/include/- extra-lib-dirs: /usr/local/lib/- include-dirs: /usr/local/include/+ if !flag(disable-default-paths)+ extra-lib-dirs: /opt/local/lib/+ include-dirs: /opt/local/include/+ extra-lib-dirs: /usr/local/lib/+ include-dirs: /usr/local/include/ if arch(i386) cc-options: -arch i386 frameworks: Accelerate if os(freebsd) if flag(openblas)- extra-lib-dirs: /usr/local/lib/openblas/lib+ if !flag(disable-default-paths)+ extra-lib-dirs: /usr/local/lib/openblas/lib extra-libraries: openblas else extra-libraries: blas lapack - extra-lib-dirs: /usr/local/lib- include-dirs: /usr/local/include- extra-libraries: gfortran+ if !flag(disable-default-paths)+ 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+ extra-libraries: openblas else extra-libraries: blas lapack if os(linux) if flag(openblas)- extra-lib-dirs: /usr/lib/openblas/lib+ if !flag(disable-default-paths)+ extra-lib-dirs: /usr/lib/openblas/lib extra-libraries: openblas else extra-libraries: blas lapack@@ -141,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,9 +1,10 @@ {-# LANGUAGE FlexibleContexts, FlexibleInstances #-} {-# LANGUAGE CPP #-} {-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE TypeFamilies #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-}+ ----------------------------------------------------------------------------- {- | Module : Internal.Algorithms@@ -20,16 +21,24 @@ -} ----------------------------------------------------------------------------- -module Internal.Algorithms where+module Internal.Algorithms (+ module Internal.Algorithms,+ UpLo(..)+) where +#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif+ import Internal.Vector import Internal.Matrix import Internal.Element import Internal.Conversion-import Internal.LAPACK as LAPACK+import Internal.LAPACK import Internal.Numeric import Data.List(foldl1') import qualified Data.Array as A+import qualified Data.Vector.Storable as Vector import Internal.ST import Internal.Vectorized(range) import Control.DeepSeq@@ -57,13 +66,17 @@ mbLinearSolve' :: Matrix t -> Matrix t -> Maybe (Matrix t) linearSolve' :: Matrix t -> Matrix t -> Matrix t cholSolve' :: Matrix t -> Matrix t -> Matrix t+ triSolve' :: UpLo -> Matrix t -> Matrix t -> Matrix t+ triDiagSolve' :: Vector t -> Vector t -> Vector t -> Matrix t -> Matrix t ldlPacked' :: Matrix t -> (Matrix t, [Int]) ldlSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t 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)@@ -82,11 +95,15 @@ linearSolve' = linearSolveR -- (luSolve . luPacked) ?? mbLinearSolve' = mbLinearSolveR cholSolve' = cholSolveR+ triSolve' = triSolveR+ triDiagSolve' = triDiagSolveR linearSolveLS' = linearSolveLSR linearSolveSVD' = linearSolveSVDR Nothing eig' = eigR eigSH'' = eigS+ geig' = eigG eigOnly = eigOnlyR+ geigOnly = eigOnlyG eigOnlySH = eigOnlyS cholSH' = cholS mbCholSH' = mbCholS@@ -111,10 +128,14 @@ linearSolve' = linearSolveC mbLinearSolve' = mbLinearSolveC cholSolve' = cholSolveC+ triSolve' = triSolveC+ triDiagSolve' = triDiagSolveC linearSolveLS' = linearSolveLSC linearSolveSVD' = linearSolveSVDC Nothing eig' = eigC+ geig' = eigGC eigOnly = eigOnlyC+ geigOnly = eigOnlyGC eigSH'' = eigH eigOnlySH = eigOnlyH cholSH' = cholH@@ -158,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@@ -212,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@@ -271,7 +294,8 @@ -0.690 -0.352 >>> s-fromList [35.18264833189422,1.4769076999800903]+[35.18264833189422,1.476907699980091]+it :: Vector Double >>> disp 3 u 5x2@@ -291,9 +315,13 @@ -} compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)-compactSVD m = (u', subVector 0 d s, v') where+compactSVD = compactSVDTol 1++-- | @compactSVDTol r@ is similar to 'compactSVD' (for which @r=1@), but uses tolerance @tol=r*g*eps*(max rows cols)@ to distinguish nonzero singular values, where @g@ is the greatest singular value. If @g<r*eps@, then only one singular value is returned.+compactSVDTol :: Field t => Double -> Matrix t -> (Matrix t, Vector Double, Matrix t)+compactSVDTol r m = (u', subVector 0 d s, v') where (u,s,v) = thinSVD m- d = rankSVD (1*eps) m s `max` 1+ d = rankSVD (r*eps) m s `max` 1 u' = takeColumns d u v' = takeColumns d v @@ -345,6 +373,79 @@ -> Matrix t -- ^ solution cholSolve = {-# SCC "cholSolve" #-} cholSolve' +-- | Solve a triangular linear system. If `Upper` is specified then+-- all elements below the diagonal are ignored; if `Lower` is+-- specified then all elements above the diagonal are ignored.+triSolve+ :: Field t+ => UpLo -- ^ `Lower` or `Upper`+ -> Matrix t -- ^ coefficient matrix+ -> Matrix t -- ^ right hand sides+ -> Matrix t -- ^ solution+triSolve = {-# SCC "triSolve" #-} triSolve'++-- | Solve a tridiagonal linear system. Suppose you wish to solve \(Ax = b\) where+--+-- \[+-- A =+-- \begin{bmatrix}+-- 1.0 & 4.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0+-- \\ 3.0 & 1.0 & 4.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0+-- \\ 0.0 & 3.0 & 1.0 & 4.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0+-- \\ 0.0 & 0.0 & 3.0 & 1.0 & 4.0 & 0.0 & 0.0 & 0.0 & 0.0+-- \\ 0.0 & 0.0 & 0.0 & 3.0 & 1.0 & 4.0 & 0.0 & 0.0 & 0.0+-- \\ 0.0 & 0.0 & 0.0 & 0.0 & 3.0 & 1.0 & 4.0 & 0.0 & 0.0+-- \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 3.0 & 1.0 & 4.0 & 0.0+-- \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 3.0 & 1.0 & 4.0+-- \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 3.0 & 1.0+-- \end{bmatrix}+-- \quad+-- b =+-- \begin{bmatrix}+-- 1.0 & 1.0 & 1.0+-- \\ 1.0 & -1.0 & 2.0+-- \\ 1.0 & 1.0 & 3.0+-- \\ 1.0 & -1.0 & 4.0+-- \\ 1.0 & 1.0 & 5.0+-- \\ 1.0 & -1.0 & 6.0+-- \\ 1.0 & 1.0 & 7.0+-- \\ 1.0 & -1.0 & 8.0+-- \\ 1.0 & 1.0 & 9.0+-- \end{bmatrix}+-- \]+--+-- then+--+-- @+-- dL = fromList [3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0]+-- d = fromList [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]+-- dU = fromList [4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0]+--+-- b = (9><3)+-- [+-- 1.0, 1.0, 1.0,+-- 1.0, -1.0, 2.0,+-- 1.0, 1.0, 3.0,+-- 1.0, -1.0, 4.0,+-- 1.0, 1.0, 5.0,+-- 1.0, -1.0, 6.0,+-- 1.0, 1.0, 7.0,+-- 1.0, -1.0, 8.0,+-- 1.0, 1.0, 9.0+-- ]+--+-- x = triDiagSolve dL d dU b+-- @+--+triDiagSolve+ :: Field t+ => Vector t -- ^ lower diagonal: \(n - 1\) elements+ -> Vector t -- ^ diagonal: \(n\) elements+ -> Vector t -- ^ upper diagonal: \(n - 1\) elements+ -> Matrix t -- ^ right hand sides+ -> Matrix t -- ^ solution+triDiagSolve = {-# SCC "triDiagSolve" #-} triDiagSolve'+ -- | Minimum norm solution of a general linear least squares problem Ax=B using the SVD. Admits rank-deficient systems but it is slower than 'linearSolveLS'. The effective rank of A is determined by treating as zero those singular valures which are less than 'eps' times the largest singular value. linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t linearSolveSVD = {-# SCC "linearSolveSVD" #-} linearSolveSVD'@@ -420,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''@@ -446,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@@ -475,9 +591,14 @@ -- | QR factorization. -- -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.+-- Note: the current implementation is very slow for large matrices. 'thinQR' is much faster. qr :: Field t => Matrix t -> (Matrix t, Matrix t) qr = {-# SCC "qr" #-} unpackQR . qr' +-- | A version of 'qr' which returns only the @min (rows m) (cols m)@ columns of @q@ and rows of @r@.+thinQR :: Field t => Matrix t -> (Matrix t, Matrix t)+thinQR = {-# SCC "thinQR" #-} thinUnpackQR . qr'+ -- | Compute the QR decomposition of a matrix in compact form. qrRaw :: Field t => Matrix t -> QR t qrRaw m = QR x v@@ -494,9 +615,17 @@ -- | RQ factorization. -- -- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.+-- Note: the current implementation is very slow for large matrices. 'thinRQ' is much faster. rq :: Field t => Matrix t -> (Matrix t, Matrix t)-rq m = {-# SCC "rq" #-} (r,q) where- (q',r') = qr $ trans $ rev1 m+rq = {-# SCC "rq" #-} rqFromQR qr++-- | A version of 'rq' which returns only the @min (rows m) (cols m)@ columns of @r@ and rows of @q@.+thinRQ :: Field t => Matrix t -> (Matrix t, Matrix t)+thinRQ = {-# SCC "thinQR" #-} rqFromQR thinQR++rqFromQR :: Field t => (Matrix t -> (Matrix t, Matrix t)) -> Matrix t -> (Matrix t, Matrix t)+rqFromQR qr0 m = (r,q) where+ (q',r') = qr0 $ trans $ rev1 m r = rev2 (trans r') q = rev2 (trans q') rev1 = flipud . fliprl@@ -724,6 +853,12 @@ hs = zipWith haussholder (toList tau) vs q = foldl1' mXm hs +thinUnpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)+thinUnpackQR (pq, tau) = (q, r)+ where mn = uncurry min $ size pq+ q = qrgr mn $ QR pq tau+ r = fromRows $ zipWith (\i v -> Vector.replicate i 0 Vector.++ Vector.drop i v) [0..mn-1] (toRows pq)+ unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) unpackHess hf m | rows m == 1 = ((1><1)[1],m)@@ -1027,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 ////////////@@ -584,6 +766,182 @@ OK } +//////// triangular real linear system ////////////++int dtrtrs_(char *uplo, char *trans, char *diag, integer *n, integer *nrhs,+ doublereal *a, integer *lda, doublereal *b, integer *ldb, integer *+ info);++int triSolveR_l_u(KODMAT(a),ODMAT(b)) {+ integer n = ar;+ integer lda = aXc;+ integer nhrs = bc;+ REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);+ DEBUGMSG("triSolveR_l_u");+ integer res;+ dtrtrs_ ("U",+ "N",+ "N",+ &n,&nhrs,+ (double*)ap, &lda,+ bp, &n,+ &res);+ CHECK(res,res);+ OK+}++int triSolveR_l_l(KODMAT(a),ODMAT(b)) {+ integer n = ar;+ integer lda = aXc;+ integer nhrs = bc;+ REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);+ DEBUGMSG("triSolveR_l_l");+ integer res;+ dtrtrs_ ("L",+ "N",+ "N",+ &n,&nhrs,+ (double*)ap, &lda,+ bp, &n,+ &res);+ CHECK(res,res);+ OK+}++//////// triangular complex linear system ////////////++int ztrtrs_(char *uplo, char *trans, char *diag, integer *n, integer *nrhs,+ doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,+ integer *info);++int triSolveC_l_u(KOCMAT(a),OCMAT(b)) {+ integer n = ar;+ integer lda = aXc;+ integer nhrs = bc;+ REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);+ DEBUGMSG("triSolveC_l_u");+ integer res;+ ztrtrs_ ("U",+ "N",+ "N",+ &n,&nhrs,+ (doublecomplex*)ap, &lda,+ bp, &n,+ &res);+ CHECK(res,res);+ OK+}++int triSolveC_l_l(KOCMAT(a),OCMAT(b)) {+ integer n = ar;+ integer lda = aXc;+ integer nhrs = bc;+ REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);+ DEBUGMSG("triSolveC_l_u");+ integer res;+ ztrtrs_ ("L",+ "N",+ "N",+ &n,&nhrs,+ (doublecomplex*)ap, &lda,+ bp, &n,+ &res);+ CHECK(res,res);+ OK+}++//////// tridiagonal real linear system ////////////++int dgttrf_(integer *n,+ doublereal *dl, doublereal *d, doublereal *du, doublereal *du2,+ integer *ipiv,+ integer *info);++int dgttrs_(char *trans, integer *n, integer *nrhs,+ doublereal *dl, doublereal *d, doublereal *du, doublereal *du2,+ integer *ipiv, doublereal *b, integer *ldb,+ integer *info);++int triDiagSolveR_l(DVEC(dl), DVEC(d), DVEC(du), ODMAT(b)) {+ 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));+ UNWIND(!ipiv,MEM,cleanup0);++ double* du2 = (double*)malloc((n - 2)*sizeof(double));+ UNWIND(!du2,MEM,cleanup1);++ dgttrf_ (&n,+ dlp, dp, dup, du2,+ ipiv,+ &res);+ UNWIND(res,res,cleanup2);++ dgttrs_ ("N",+ &n,&nhrs,+ dlp, dp, dup, du2,+ ipiv, bp, &n,+ &res);+ UNWIND(res,res,cleanup2);++cleanup2:+ free(du2);+cleanup1:+ free(ipiv);+cleanup0:+ return ret;+}++//////// tridiagonal complex linear system ////////////++int zgttrf_(integer *n,+ doublecomplex *dl, doublecomplex *d, doublecomplex *du, doublecomplex *du2,+ integer *ipiv,+ integer *info);++int zgttrs_(char *trans, integer *n, integer *nrhs,+ doublecomplex *dl, doublecomplex *d, doublecomplex *du, doublecomplex *du2,+ integer *ipiv, doublecomplex *b, integer *ldb,+ integer *info);++int triDiagSolveC_l(CVEC(dl), CVEC(d), CVEC(du), OCMAT(b)) {+ 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);+ UNWIND(res,res,cleanup2);++ zgttrs_ ("N",+ &n,&nhrs,+ dlp, dp, dup, du2,+ ipiv, bp, &n,+ &res);+ UNWIND(res,res,cleanup2);++cleanup2:+ free(du2);+cleanup1:+ free(ipiv);+cleanup0:+ return ret;+}+ //////////////////// least squares real linear system //////////// int dgels_(char *trans, integer *m, integer *n, integer *@@ -591,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;@@ -605,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 ////////////@@ -627,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;@@ -641,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 ////////////@@ -664,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;@@ -682,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,@@ -691,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 ////////////@@ -709,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;@@ -729,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,@@ -739,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 /////////////////////////@@ -799,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 *@@ -836,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,@@ -855,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; } @@ -877,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; } @@ -899,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 /////////////////////////@@ -923,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; } @@ -954,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 /////////////////////////@@ -985,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; } @@ -1010,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; } @@ -1038,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; } @@ -1063,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; } @@ -1089,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;@@ -1100,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; } @@ -1116,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;@@ -1127,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; } @@ -1145,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; } @@ -1169,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,19 +932,34 @@ //////////////////////////////////////////////////////////////////////////////// -#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". */+#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) p is matching probalility of random().@@ -959,6 +974,22 @@ return (double)nrand48(state) / (double)max_random; } +#else++#define _CRT_RAND_S+inline double urandom() {+ unsigned int number;+ errno_t err;+ err = rand_s(&number);+ if (err!=0) {+ printf("something wrong\n");+ return -1;+ }+ return (double)number / (double)UINT_MAX;+}++#endif+ double gaussrand(int *phase, double *pV1, double *pV2, double *pS) { double V1=*pV1, V2=*pV2, S=*pS;@@ -985,10 +1016,39 @@ } +#if defined(_WIN32) || defined(WIN32)+ int random_vector(unsigned int seed, int code, DVEC(r)) { int phase = 0; double V1,V2,S; + srand(seed);++ int k;+ switch (code) {+ case 0: { // uniform+ for (k=0; k<rn; k++) {+ rp[k] = urandom();+ }+ OK+ }+ case 1: { // gaussian+ for (k=0; k<rn; k++) {+ rp[k] = gaussrand(&phase,&V1,&V2,&S);+ }+ OK+ }++ default: ERROR(BAD_CODE);+ }+}++#else++int random_vector(unsigned int seed, int code, DVEC(r)) {+ int phase = 0;+ double V1,V2,S;+ srandom(seed); int k;@@ -1010,6 +1070,8 @@ } } +#endif+ #else inline double urandom(struct random_data * buffer) {@@ -1484,3 +1546,54 @@ CHOOSE_IMP } +//////////////////// reorder /////////////////////////++#define REORDER_IMP \+ REQUIRES(kn == stridesn && stridesn == dimsn ,BAD_SIZE); \+ int i,j,l; \+ for (i=1,j=0,l=0;l<kn;++l) { \+ kp[l] = 0; \+ i *= dimsp[l]; \+ j += (dimsp[l]-1) * stridesp[l]; \+ } \+ REQUIRES(i <= vn && j < rn ,BAD_SIZE); \+ for (i=0,j=0;;i++) { \+ rp[i] = vp[j]; \+ for(l=kn-1;;l--) { \+ ++kp[l]; \+ if (kp[l] < dimsp[l]) { \+ j += stridesp[l]; \+ break; \+ } else { \+ if (l == 0) { \+ return 0; \+ } \+ kp[l] = 0; \+ j -= (dimsp[l]-1) * stridesp[l]; \+ } \+ } \+ }++int reorderF(IVEC(k), KIVEC(strides),KIVEC(dims),KFVEC(v),FVEC(r)) {+ REORDER_IMP+}++int reorderD(IVEC(k), KIVEC(strides),KIVEC(dims),KDVEC(v),DVEC(r)) {+ REORDER_IMP+}++int reorderI(IVEC(k), KIVEC(strides),KIVEC(dims),KIVEC(v),IVEC(r)) {+ REORDER_IMP+}++int reorderL(IVEC(k), KIVEC(strides),KIVEC(dims),KLVEC(v),LVEC(r)) {+ REORDER_IMP+}++int reorderC(IVEC(k), KIVEC(strides),KIVEC(dims),KCVEC(v),CVEC(r)) {+ REORDER_IMP+}++int reorderQ(IVEC(k), KIVEC(strides),KIVEC(dims),KQVEC(v),QVEC(r)) {+ REORDER_IMP+}
src/Internal/CG.hs view
@@ -1,6 +1,8 @@ {-# LANGUAGE FlexibleContexts, FlexibleInstances #-} {-# LANGUAGE RecordWildCards #-} +{-# OPTIONS_GHC -fno-warn-orphans #-}+ module Internal.CG( cgSolve, cgSolve', CGState(..), R, V@@ -177,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/Chain.hs view
@@ -1,5 +1,7 @@ {-# LANGUAGE FlexibleContexts #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-}+ ----------------------------------------------------------------------------- -- | -- Module : Internal.Chain
src/Internal/Container.hs view
@@ -1,8 +1,8 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-}-{-# LANGUAGE UndecidableInstances #-} ----------------------------------------------------------------------------- -- |@@ -28,17 +28,21 @@ import Internal.Matrix import Internal.Element import Internal.Numeric-import Internal.Algorithms(Field,linearSolveSVD)-+import Internal.Algorithms(Field,linearSolveSVD,Herm,mTm)+#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif ------------------------------------------------------------------ {- | 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: @@ -82,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 #>@@ -131,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. @@ -162,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: --@@ -179,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' @@ -199,21 +207,20 @@ {- | 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, Matrix Double)+meanCov :: Matrix Double -> (Vector Double, Herm Double) meanCov x = (med,cov) where r = rows x k = 1 / fromIntegral r med = konst k r `vXm` x meds = konst 1 r `outer` med xc = x `sub` meds- cov = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc)+ cov = scale (recip (fromIntegral (r-1))) (mTm xc) -------------------------------------------------------------------------------- @@ -293,5 +300,3 @@ | otherwise = error $ "out of range index in remap" where [i',j'] = conformMs [i,j]- -
src/Internal/Conversion.hs view
@@ -3,7 +3,6 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-}-{-# LANGUAGE UndecidableInstances #-} ----------------------------------------------------------------------------- -- |
src/Internal/Convolution.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} ----------------------------------------------------------------------------- {- |@@ -23,6 +24,9 @@ import Internal.Element import Internal.Conversion import Internal.Container+#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif vectSS :: Element t => Int -> Vector t -> Matrix t@@ -37,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@@ -50,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) @@ -54,14 +55,27 @@ -- | postfix error code check infixl 0 #|+(#|) :: IO CInt -> String -> IO () (#|) = flip check -- | Error capture and conversion to Maybe mbCatch :: IO x -> IO (Maybe x)-mbCatch act = E.catch (Just `fmap` act) f- where f :: SomeException -> IO (Maybe x)- f _ = return Nothing+mbCatch act =+ hush <$>+ Exception.tryJust+ (\e -> if isSyncException e then Just e else Nothing)+ act + where+ hush :: Either a b -> Maybe b+ hush = either (const Nothing) Just++ isSyncException :: SomeException -> Bool+ isSyncException e =+ case Exception.fromException e of+ Just (SomeAsyncException _) -> False+ Nothing -> True+ -------------------------------------------------------------------------------- type CM b r = CInt -> CInt -> Ptr b -> r@@ -80,8 +94,8 @@ where type Trans c b type TransRaw c b- apply :: (Trans c b) -> c -> b- applyRaw :: (TransRaw c b) -> c -> b+ apply :: c -> (b -> IO r) -> (Trans c b) -> IO r+ applyRaw :: c -> (b -> IO r) -> (TransRaw c b) -> IO r infixl 1 `apply`, `applyRaw` instance Storable t => TransArray (Vector t)@@ -92,4 +106,3 @@ {-# INLINE apply #-} applyRaw = avec {-# INLINE applyRaw #-}-
src/Internal/Element.hs view
@@ -1,10 +1,10 @@ {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -fno-warn-orphans #-}+ ----------------------------------------------------------------------------- -- | -- Module : Data.Packed.Matrix@@ -32,14 +32,14 @@ import Foreign.Storable(Storable) import System.IO.Unsafe(unsafePerformIO) import Control.Monad(liftM)+import Foreign.C.Types(CInt) ------------------------------------------------------------------- -#ifdef BINARY 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)@@ -48,7 +48,6 @@ v <- get return (reshape c v) -#endif ------------------------------------------------------------------- @@ -56,8 +55,10 @@ show m | rows m == 0 || cols m == 0 = sizes m ++" []" show m = (sizes m++) . dsp . map (map show) . toLists $ m +sizes :: Matrix t -> [Char] sizes m = "("++show (rows m)++"><"++show (cols m)++")\n" +dsp :: [[[Char]]] -> [Char] dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp where mt = transpose as@@ -76,6 +77,7 @@ rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims +breakAt :: Eq a => a -> [a] -> ([a], [a]) breakAt c l = (a++[c],tail b) where (a,b) = break (==c) l @@ -91,7 +93,8 @@ | Drop Int | DropLast Int deriving Show- ++ppext :: Extractor -> [Char] ppext All = ":" ppext (Range a 1 c) = printf "%d:%d" a c ppext (Range a b c) = printf "%d:%d:%d" a b c@@ -131,10 +134,14 @@ infixl 9 ?? (??) :: Element t => Matrix t -> (Extractor,Extractor) -> Matrix t +minEl :: Vector CInt -> CInt minEl = toScalarI Min+maxEl :: Vector CInt -> CInt maxEl = toScalarI Max+cmodi :: Foreign.C.Types.CInt -> Vector Foreign.C.Types.CInt -> Vector Foreign.C.Types.CInt cmodi = vectorMapValI ModVS +extractError :: Matrix t1 -> (Extractor, Extractor) -> t extractError m (e1,e2)= error $ printf "can't extract (%s,%s) from matrix %dx%d" (ppext e1::String) (ppext e2::String) (rows m) (cols m) m ?? (Range a s b,e) | s /= 1 = m ?? (Pos (idxs [a,a+s .. b]), e)@@ -235,8 +242,10 @@ fromBlocks :: Element t => [[Matrix t]] -> Matrix t fromBlocks = fromBlocksRaw . adaptBlocks +fromBlocksRaw :: Element t => [[Matrix t]] -> Matrix t fromBlocksRaw mms = joinVert . map joinHoriz $ mms +adaptBlocks :: Element t => [[Matrix t]] -> [[Matrix t]] adaptBlocks ms = ms' where bc = case common length ms of Just c -> c@@ -489,6 +498,9 @@ m2' = conformMTo (r,c) m2 -- FIXME do not flatten if equal order+lM :: (Storable t, Element t1, Element t2)+ => (Vector t1 -> Vector t2 -> Vector t)+ -> Matrix t1 -> Matrix t2 -> Matrix t lM f m1 m2 = matrixFromVector RowMajor (max' (rows m1) (rows m2))@@ -507,6 +519,7 @@ ------------------------------------------------------------ +toBlockRows :: Element t => [Int] -> Matrix t -> [Matrix t] toBlockRows [r] m | r == rows m = [m] toBlockRows rs m@@ -516,6 +529,7 @@ szs = map (* cols m) rs g k = (k><0)[] +toBlockCols :: Element t => [Int] -> Matrix t -> [Matrix t] toBlockCols [c] m | c == cols m = [m] toBlockCols cs m = map trans . toBlockRows cs . trans $ m @@ -579,7 +593,7 @@ mapMatrixWithIndexM :: (Element a, Storable b, Monad m) => ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b)-mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m +mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m where c = cols m @@ -601,4 +615,3 @@ mapMatrix :: (Element a, Element b) => (a -> b) -> Matrix a -> Matrix b mapMatrix f = liftMatrix (mapVector f)-
− src/Internal/Foreign.hs
@@ -1,102 +0,0 @@-{-# LANGUAGE MagicHash, UnboxedTuples #-}--- | FFI and hmatrix helpers.------ Sample usage, to upload a perspective matrix to a shader.------ @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3) --- @-----module Internal.Foreign - ( app- , appVector, appVectorLen- , appMatrix, appMatrixLen, appMatrixRaw, appMatrixRawLen- , unsafeMatrixToVector, unsafeMatrixToForeignPtr- ) where--import Foreign.C.Types(CInt)-import Internal.Vector-import Internal.Matrix-import qualified Data.Vector.Storable as S-import Foreign (Ptr, ForeignPtr, Storable)-import GHC.Base (IO(..), realWorld#)--{-# INLINE unsafeInlinePerformIO #-}--- | If we use unsafePerformIO, it may not get inlined, so in a function that returns IO (which are all safe uses of app* in this module), there would be--- unecessary calls to unsafePerformIO or its internals.-unsafeInlinePerformIO :: IO a -> a-unsafeInlinePerformIO (IO f) = case f realWorld# of- (# _, x #) -> x--{-# INLINE app #-}--- | Only useful since it is left associated with a precedence of 1, unlike 'Prelude.$', which is right associative.--- e.g.------ @--- someFunction--- \`appMatrixLen\` m--- \`appVectorLen\` v--- \`app\` other--- \`app\` arguments--- \`app\` go here--- @------ One could also write:------ @--- (someFunction --- \`appMatrixLen\` m--- \`appVectorLen\` v) --- other --- arguments --- (go here)--- @----app :: (a -> b) -> a -> b-app f = f--{-# INLINE appVector #-}-appVector :: Storable a => (Ptr a -> b) -> Vector a -> b-appVector f x = unsafeInlinePerformIO (S.unsafeWith x (return . f))--{-# INLINE appVectorLen #-}-appVectorLen :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b-appVectorLen f x = unsafeInlinePerformIO (S.unsafeWith x (return . f (fromIntegral (S.length x))))--{-# INLINE appMatrix #-}-appMatrix :: Element a => (Ptr a -> b) -> Matrix a -> b-appMatrix f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f))--{-# INLINE appMatrixLen #-}-appMatrixLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b-appMatrixLen f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f r c))- where- r = fromIntegral (rows x)- c = fromIntegral (cols x)--{-# INLINE appMatrixRaw #-}-appMatrixRaw :: Storable a => (Ptr a -> b) -> Matrix a -> b-appMatrixRaw f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f))--{-# INLINE appMatrixRawLen #-}-appMatrixRawLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b-appMatrixRawLen f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f r c))- where- r = fromIntegral (rows x)- c = fromIntegral (cols x)--infixl 1 `app`-infixl 1 `appVector`-infixl 1 `appMatrix`-infixl 1 `appMatrixRaw`--{-# INLINE unsafeMatrixToVector #-}--- | This will disregard the order of the matrix, and simply return it as-is. --- If the order of the matrix is RowMajor, this function is identical to 'flatten'.-unsafeMatrixToVector :: Matrix a -> Vector a-unsafeMatrixToVector = xdat--{-# INLINE unsafeMatrixToForeignPtr #-}-unsafeMatrixToForeignPtr :: Storable a => Matrix a -> (ForeignPtr a, Int)-unsafeMatrixToForeignPtr m = S.unsafeToForeignPtr0 (xdat m)-
src/Internal/IO.hs view
@@ -20,7 +20,7 @@ import Internal.Vector import Internal.Matrix import Internal.Vectorized-import Text.Printf(printf)+import Text.Printf(printf, PrintfArg, PrintfType) import Data.List(intersperse,transpose) import Data.Complex @@ -78,12 +78,18 @@ dispf :: Int -> Matrix Double -> String dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x +sdims :: Matrix t -> [Char] sdims x = show (rows x) ++ "x" ++ show (cols x) +formatFixed :: (Show a, Text.Printf.PrintfArg t, Element t)+ => a -> Matrix t -> String formatFixed d x = format " " (printf ("%."++show d++"f")) $ x +isInt :: Matrix Double -> Bool isInt = all lookslikeInt . toList . flatten +formatScaled :: (Text.Printf.PrintfArg b, RealFrac b, Floating b, Num t, Element b, Show t)+ => t -> Matrix b -> [Char] formatScaled dec t = "E"++show o++"\n" ++ ss where ss = format " " (printf fmt. g) t g x | o >= 0 = x/10^(o::Int)@@ -133,14 +139,18 @@ s2 = if b<0 then "-" else "" s3 = if b<0 then "-" else "+" +shcr :: (Show a, Show t1, Text.Printf.PrintfType t, Text.Printf.PrintfArg t1, RealFrac t1)+ => a -> t1 -> t shcr d a | lookslikeInt a = printf "%.0f" a | otherwise = printf ("%."++show d++"f") a -+lookslikeInt :: (Show a, RealFrac a) => a -> Bool lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx where shx = show x +isZero :: Show a => a -> Bool isZero x = show x `elem` ["0.0","-0.0"]+isOne :: Show a => a -> Bool isOne x = show x `elem` ["1.0","-1.0"] -- | Pretty print a complex matrix with at most n decimal digits.@@ -168,6 +178,6 @@ else return (reshape c v) -+loadMatrix' :: FilePath -> IO (Maybe (Matrix Double)) loadMatrix' name = mbCatch (loadMatrix name)
src/Internal/LAPACK.hs view
@@ -1,6 +1,8 @@ {-# LANGUAGE TypeOperators #-} {-# LANGUAGE ViewPatterns #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-}+ ----------------------------------------------------------------------------- -- | -- Module : Numeric.LinearAlgebra.LAPACK@@ -16,9 +18,11 @@ module Internal.LAPACK where +import Data.Bifunctor (first)+ import Internal.Devel import Internal.Vector-import Internal.Matrix hiding ((#))+import Internal.Matrix hiding ((#), (#!)) import Internal.Conversion import Internal.Element import Foreign.Ptr(nullPtr)@@ -28,10 +32,13 @@ ----------------------------------------------------------------------------------- -infixl 1 #+infixr 1 # a # b = apply a b {-# INLINE (#) #-} +a #! b = a # b # id+{-# INLINE (#!) #-}+ ----------------------------------------------------------------------------------- type TMMM t = t ::> t ::> t ::> Ok@@ -56,7 +63,7 @@ when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++ show (rows a,cols a) ++ " x " ++ show (rows b, cols b) s <- createMatrix ColumnMajor (rows a) (cols b)- f (isT a) (isT b) # (tt a) # (tt b) # s #| st+ ((tt a) # (tt b) #! s) (f (isT a) (isT b)) #| st return s -- | Matrix product based on BLAS's /dgemm/.@@ -80,7 +87,7 @@ when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b s <- createMatrix ColumnMajor (rows a) (cols b)- c_multiplyI m # a # b # s #|"c_multiplyI"+ (a # b #! s) (c_multiplyI m) #|"c_multiplyI" return s multiplyL :: Z -> Matrix Z -> Matrix Z -> Matrix Z@@ -88,7 +95,7 @@ when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b s <- createMatrix ColumnMajor (rows a) (cols b)- c_multiplyL m # a # b # s #|"c_multiplyL"+ (a # b #! s) (c_multiplyL m) #|"c_multiplyL" return s -----------------------------------------------------------------------------@@ -121,7 +128,7 @@ u <- createMatrix ColumnMajor r r s <- createVector (min r c) v <- createMatrix ColumnMajor c c- f # a # u # s # v #| st+ (a # u # s #! v) f #| st return (u,s,v) where r = rows x@@ -149,7 +156,7 @@ u <- createMatrix ColumnMajor r q s <- createVector q v <- createMatrix ColumnMajor q c- f # a # u # s # v #| st+ (a # u # s #! v) f #| st return (u,s,v) where r = rows x@@ -176,7 +183,7 @@ svAux f st x = unsafePerformIO $ do a <- copy ColumnMajor x s <- createVector q- g # a # s #| st+ (a #! s) g #| st return s where r = rows x@@ -197,7 +204,7 @@ a <- copy ColumnMajor x s <- createVector q v <- createMatrix ColumnMajor c c- g # a # s # v #| st+ (a # s #! v) g #| st return (s,v) where r = rows x@@ -218,7 +225,7 @@ a <- copy ColumnMajor x u <- createMatrix ColumnMajor r r s <- createVector q- g # a # u # s #| st+ (a # u #! s) g #| st return (u,s) where r = rows x@@ -229,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 @@ -237,7 +246,7 @@ a <- copy ColumnMajor m l <- createVector r v <- createMatrix ColumnMajor r r- g # a # l # v #| st+ (a # l #! v) g #| st return (l,v) where r = rows m@@ -252,7 +261,7 @@ eigOnlyAux f st m = unsafePerformIO $ do a <- copy ColumnMajor m l <- createVector r- g # a # l #| st+ (a #! l) g #| st return l where r = rows m@@ -277,7 +286,7 @@ a <- copy ColumnMajor m l <- createVector r v <- createMatrix ColumnMajor r r- g # a # l # v #| "eigR"+ (a # l #! v) g #| "eigR" return (l,v) where r = rows m@@ -293,19 +302,75 @@ | 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 l <- createVector r v <- copy ColumnMajor m- f # l # v #| st+ (l #! v) f #| st return (l,v) where r = rows m@@ -356,7 +421,7 @@ | n1==n2 && n1==r = unsafePerformIO . g $ do a' <- copy ColumnMajor a s <- copy ColumnMajor b- f # a' # s #| st+ (a' #! s) f #| st return s | otherwise = error $ st ++ " of nonsquare matrix" where@@ -387,7 +452,7 @@ linearSolveSQAux2 g f st a b | n1==n2 && n1==r = unsafePerformIO . g $ do s <- copy ColumnMajor b- f # a # s #| st+ (a #! s) f #| st return s | otherwise = error $ st ++ " of nonsquare matrix" where@@ -403,6 +468,60 @@ cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) cholSolveC a b = linearSolveSQAux2 id zpotrs "cholSolveC" (fmat a) b +--------------------------------------------------------------------------------+foreign import ccall unsafe "triSolveR_l_u" dtrtrs_u :: R ::> R ::> Ok+foreign import ccall unsafe "triSolveC_l_u" ztrtrs_u :: C ::> C ::> Ok+foreign import ccall unsafe "triSolveR_l_l" dtrtrs_l :: R ::> R ::> Ok+foreign import ccall unsafe "triSolveC_l_l" ztrtrs_l :: C ::> C ::> Ok+++linearSolveTRAux2 g f st a b+ | n1==n2 && n1==r = unsafePerformIO . g $ do+ s <- copy ColumnMajor b+ (a #! s) f #| st+ return s+ | otherwise = error $ st ++ " of nonsquare matrix"+ where+ n1 = rows a+ n2 = cols a+ r = rows b++data UpLo = Lower | Upper++-- | Solves a triangular system of linear equations.+triSolveR :: UpLo -> Matrix Double -> Matrix Double -> Matrix Double+triSolveR Lower a b = linearSolveTRAux2 id dtrtrs_l "triSolveR" (fmat a) b+triSolveR Upper a b = linearSolveTRAux2 id dtrtrs_u "triSolveR" (fmat a) b++-- | Solves a triangular system of linear equations.+triSolveC :: UpLo -> Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)+triSolveC Lower a b = linearSolveTRAux2 id ztrtrs_l "triSolveC" (fmat a) b+triSolveC Upper a b = linearSolveTRAux2 id ztrtrs_u "triSolveC" (fmat a) b++--------------------------------------------------------------------------------+foreign import ccall unsafe "triDiagSolveR_l" dgttrs :: R :> R :> R :> R ::> Ok+foreign import ccall unsafe "triDiagSolveC_l" zgttrs :: C :> C :> C :> C ::> Ok++linearSolveGTAux2 g f st dl d du b+ | 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+ return s+ | otherwise = error $ st ++ " of nonsquare matrix"+ where+ ndl = dim dl+ nd = dim d+ ndu = dim du+ r = rows b++-- | Solves a tridiagonal system of linear equations.+triDiagSolveR dl d du b = linearSolveGTAux2 id dgttrs "triDiagSolveR" dl d du b+triDiagSolveC dl d du b = linearSolveGTAux2 id zgttrs "triDiagSolveC" dl d du b+ ----------------------------------------------------------------------------------- foreign import ccall unsafe "linearSolveLSR_l" dgels :: R ::> R ::> Ok@@ -415,7 +534,7 @@ a' <- copy ColumnMajor a r <- createMatrix ColumnMajor (max m n) nrhs setRect 0 0 b r- f # a' # r #| st+ (a' #! r) f #| st return r | otherwise = error $ "different number of rows in linearSolve ("++st++")" where@@ -458,7 +577,7 @@ cholAux f st a = do r <- copy ColumnMajor a- f # r #| st+ (r # id) f #| st return r -- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.@@ -495,7 +614,7 @@ qrAux f st a = unsafePerformIO $ do r <- copy ColumnMajor a tau <- createVector mn- f # tau # r #| st+ (tau #! r) f #| st return (r,tau) where m = rows a@@ -514,7 +633,7 @@ qrgrAux f st n (a, tau) = unsafePerformIO $ do res <- copy ColumnMajor (subMatrix (0,0) (rows a,n) a)- f # (subVector 0 n tau') # res #| st+ ((subVector 0 n tau') #! res) f #| st return res where tau' = vjoin [tau, constantD 0 n]@@ -534,7 +653,7 @@ hessAux f st a = unsafePerformIO $ do r <- copy ColumnMajor a tau <- createVector (mn-1)- f # tau # r #| st+ (tau #! r) f #| st return (r,tau) where m = rows a@@ -556,7 +675,7 @@ schurAux f st a = unsafePerformIO $ do u <- createMatrix ColumnMajor n n s <- copy ColumnMajor a- f # u # s #| st+ (u #! s) f #| st return (u,s) where n = rows a@@ -576,7 +695,7 @@ luAux f st a = unsafePerformIO $ do lu <- copy ColumnMajor a piv <- createVector (min n m)- f # piv # lu #| st+ (piv #! lu) f #| st return (lu, map (pred.round) (toList piv)) where n = rows a@@ -598,7 +717,7 @@ lusAux f st a piv b | n1==n2 && n2==n =unsafePerformIO $ do x <- copy ColumnMajor b- f # a # piv' # x #| st+ (a # piv' #! x) f #| st return x | otherwise = error st where@@ -622,7 +741,7 @@ ldlAux f st a = unsafePerformIO $ do ldl <- copy ColumnMajor a piv <- createVector (rows a)- f # piv # ldl #| st+ (piv #! ldl) f #| st return (ldl, map (pred.round) (toList piv)) -----------------------------------------------------------------------------------@@ -637,4 +756,3 @@ -- | Solve a complex linear system from a precomputed LDL decomposition ('ldlC'), using LAPACK's /zsytrs/. ldlsC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double) ldlsC a piv b = lusAux zsytrs "ldlsC" (fmat a) piv b-
src/Internal/Matrix.hs view
@@ -5,10 +5,9 @@ {-# LANGUAGE TypeOperators #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE ViewPatterns #-}+{-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE ConstrainedClassMethods #-} -- -- | -- Module : Internal.Matrix -- Copyright : (c) Alberto Ruiz 2007-15@@ -23,7 +22,7 @@ import Internal.Vector import Internal.Devel-import Internal.Vectorized hiding ((#))+import Internal.Vectorized hiding ((#), (#!)) import Foreign.Marshal.Alloc ( free ) import Foreign.Marshal.Array(newArray) import Foreign.Ptr ( Ptr )@@ -58,19 +57,24 @@ cols = icols {-# INLINE cols #-} +size :: Matrix t -> (Int, Int) size m = (irows m, icols m) {-# INLINE size #-} +rowOrder :: Matrix t -> Bool rowOrder m = xCol m == 1 || cols m == 1 {-# INLINE rowOrder #-} +colOrder :: Matrix t -> Bool colOrder m = xRow m == 1 || rows m == 1 {-# INLINE colOrder #-} +is1d :: Matrix t -> Bool is1d (size->(r,c)) = r==1 || c==1 {-# INLINE is1d #-} -- data is not contiguous+isSlice :: Storable t => Matrix t -> Bool isSlice m@(size->(r,c)) = r*c < dim (xdat m) {-# INLINE isSlice #-} @@ -111,15 +115,15 @@ -- C-Haskell matrix adapters {-# INLINE amatr #-}-amatr :: Storable a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b-amatr f x = inlinePerformIO (unsafeWith (xdat x) (return . f r c))+amatr :: Storable a => Matrix a -> (f -> IO r) -> (CInt -> CInt -> Ptr a -> f) -> IO r+amatr x f g = unsafeWith (xdat x) (f . g r c) where r = fi (rows x) c = fi (cols x) {-# INLINE amat #-}-amat :: Storable a => (CInt -> CInt -> CInt -> CInt -> Ptr a -> b) -> Matrix a -> b-amat f x = inlinePerformIO (unsafeWith (xdat x) (return . f r c sr sc))+amat :: Storable a => Matrix a -> (f -> IO r) -> (CInt -> CInt -> CInt -> CInt -> Ptr a -> f) -> IO r+amat x f g = unsafeWith (xdat x) (f . g r c sr sc) where r = fi (rows x) c = fi (cols x)@@ -136,20 +140,28 @@ applyRaw = amatr {-# INLINE applyRaw #-} -infixl 1 #+infixr 1 #+(#) :: TransArray c => c -> (b -> IO r) -> Trans c b -> IO r a # b = apply a b {-# INLINE (#) #-} +(#!) :: (TransArray c, TransArray c1) => c1 -> c -> Trans c1 (Trans c (IO r)) -> IO r+a #! b = a # b # id+{-# INLINE (#!) #-}+ -------------------------------------------------------------------------------- +copy :: Element t => MatrixOrder -> Matrix t -> IO (Matrix t) copy ord m = extractR ord m 0 (idxs[0,rows m-1]) 0 (idxs[0,cols m-1]) +extractAll :: Element t => MatrixOrder -> Matrix t -> Matrix t extractAll ord m = unsafePerformIO (copy ord m) {- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose. >>> flatten (ident 3)-fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]+[1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]+it :: (Num t, Element t) => Vector t -} flatten :: Element t => Matrix t -> Vector t@@ -221,11 +233,13 @@ {-# INLINE (@@>) #-} -- Unsafe matrix access without range checking+atM' :: Storable t => Matrix t -> Int -> Int -> t atM' m i j = xdat m `at'` (i * (xRow m) + j * (xCol m)) {-# INLINE atM' #-} ------------------------------------------------------------------ +matrixFromVector :: Storable t => MatrixOrder -> Int -> Int -> Vector t -> Matrix t matrixFromVector _ 1 _ v@(dim->d) = Matrix { irows = 1, icols = d, xdat = v, xRow = d, xCol = 1 } matrixFromVector _ _ 1 v@(dim->d) = Matrix { irows = d, icols = 1, xdat = v, xRow = 1, xCol = d } matrixFromVector o r c v@@ -283,6 +297,7 @@ remapM :: Matrix CInt -> Matrix CInt -> Matrix a -> Matrix a rowOp :: Int -> a -> Int -> Int -> Int -> Int -> Matrix a -> IO () gemm :: Vector a -> Matrix a -> Matrix a -> Matrix a -> IO ()+ reorderV :: Vector CInt-> Vector CInt-> Vector a -> Vector a -- see reorderVector for documentation instance Element Float where@@ -296,6 +311,7 @@ remapM = remapF rowOp = rowOpAux c_rowOpF gemm = gemmg c_gemmF+ reorderV = reorderAux c_reorderF instance Element Double where constantD = constantAux cconstantR@@ -308,6 +324,7 @@ remapM = remapD rowOp = rowOpAux c_rowOpD gemm = gemmg c_gemmD+ reorderV = reorderAux c_reorderD instance Element (Complex Float) where constantD = constantAux cconstantQ@@ -320,6 +337,7 @@ remapM = remapQ rowOp = rowOpAux c_rowOpQ gemm = gemmg c_gemmQ+ reorderV = reorderAux c_reorderQ instance Element (Complex Double) where constantD = constantAux cconstantC@@ -332,6 +350,7 @@ remapM = remapC rowOp = rowOpAux c_rowOpC gemm = gemmg c_gemmC+ reorderV = reorderAux c_reorderC instance Element (CInt) where constantD = constantAux cconstantI@@ -344,6 +363,7 @@ remapM = remapI rowOp = rowOpAux c_rowOpI gemm = gemmg c_gemmI+ reorderV = reorderAux c_reorderI instance Element Z where constantD = constantAux cconstantL@@ -356,6 +376,7 @@ remapM = remapL rowOp = rowOpAux c_rowOpL gemm = gemmg c_gemmL+ reorderV = reorderAux c_reorderL ------------------------------------------------------------------- @@ -378,18 +399,21 @@ -------------------------------------------------------------------------- +maxZ :: (Num t1, Ord t1, Foldable t) => t t1 -> t1 maxZ xs = if minimum xs == 0 then 0 else maximum xs +conformMs :: Element t => [Matrix t] -> [Matrix t] conformMs ms = map (conformMTo (r,c)) ms where r = maxZ (map rows ms) c = maxZ (map cols ms) -+conformVs :: Element t => [Vector t] -> [Vector t] conformVs vs = map (conformVTo n) vs where n = maxZ (map dim vs) +conformMTo :: Element t => (Int, Int) -> Matrix t -> Matrix t conformMTo (r,c) m | size m == (r,c) = m | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c))@@ -397,18 +421,24 @@ | size m == (1,c) = repRows r m | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to " ++ shDim (r,c) +conformVTo :: Element t => Int -> Vector t -> Vector t conformVTo n v | dim v == n = v | dim v == 1 = constantD (v@>0) n | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n +repRows :: Element t => Int -> Matrix t -> Matrix t repRows n x = fromRows (replicate n (flatten x))+repCols :: Element t => Int -> Matrix t -> Matrix t repCols n x = fromColumns (replicate n (flatten x)) +shSize :: Matrix t -> [Char] shSize = shDim . size +shDim :: (Show a, Show a1) => (a1, a) -> [Char] shDim (r,c) = "(" ++ show r ++"x"++ show c ++")" +emptyM :: Storable t => Int -> Int -> Matrix t emptyM r c = matrixFromVector RowMajor r c (fromList[]) ----------------------------------------------------------------------@@ -423,11 +453,17 @@ --------------------------------------------------------------- +extractAux :: (Eq t3, Eq t2, TransArray c, Storable a, Storable t1,+ Storable t, Num t3, Num t2, Integral t1, Integral t)+ => (t3 -> t2 -> CInt -> Ptr t1 -> CInt -> Ptr t+ -> Trans c (CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt))+ -> MatrixOrder -> c -> t3 -> Vector t1 -> t2 -> Vector t -> IO (Matrix a) extractAux f ord m moder vr modec vc = do let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc r <- createMatrix ord nr nc- f moder modec # vr # vc # m # r #|"extract"+ (vr # vc # m #! r) (f moder modec) #|"extract"+ return r type Extr x = CInt -> CInt -> CIdxs (CIdxs (OM x (OM x (IO CInt))))@@ -441,7 +477,10 @@ --------------------------------------------------------------- -setRectAux f i j m r = f (fi i) (fi j) # m # r #|"setRect"+setRectAux :: (TransArray c1, TransArray c)+ => (CInt -> CInt -> Trans c1 (Trans c (IO CInt)))+ -> Int -> Int -> c1 -> c -> IO ()+setRectAux f i j m r = (m #! r) (f (fi i) (fi j)) #|"setRect" type SetRect x = I -> I -> x ::> x::> Ok @@ -454,19 +493,29 @@ -------------------------------------------------------------------------------- +sortG :: (Storable t, Storable a)+ => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a sortG f v = unsafePerformIO $ do r <- createVector (dim v)- f # v # r #|"sortG"+ (v #! r) f #|"sortG" return r +sortIdxD :: Vector Double -> Vector CInt sortIdxD = sortG c_sort_indexD+sortIdxF :: Vector Float -> Vector CInt sortIdxF = sortG c_sort_indexF+sortIdxI :: Vector CInt -> Vector CInt sortIdxI = sortG c_sort_indexI+sortIdxL :: Vector Z -> Vector I sortIdxL = sortG c_sort_indexL +sortValD :: Vector Double -> Vector Double sortValD = sortG c_sort_valD+sortValF :: Vector Float -> Vector Float sortValF = sortG c_sort_valF+sortValI :: Vector CInt -> Vector CInt sortValI = sortG c_sort_valI+sortValL :: Vector Z -> Vector Z sortValL = sortG c_sort_valL foreign import ccall unsafe "sort_indexD" c_sort_indexD :: CV Double (CV CInt (IO CInt))@@ -481,14 +530,21 @@ -------------------------------------------------------------------------------- +compareG :: (TransArray c, Storable t, Storable a)+ => Trans c (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt)+ -> c -> Vector t -> Vector a compareG f u v = unsafePerformIO $ do r <- createVector (dim v)- f # u # v # r #|"compareG"+ (u # v #! r) f #|"compareG" return r +compareD :: Vector Double -> Vector Double -> Vector CInt compareD = compareG c_compareD+compareF :: Vector Float -> Vector Float -> Vector CInt compareF = compareG c_compareF+compareI :: Vector CInt -> Vector CInt -> Vector CInt compareI = compareG c_compareI+compareL :: Vector Z -> Vector Z -> Vector CInt compareL = compareG c_compareL foreign import ccall unsafe "compareD" c_compareD :: CV Double (CV Double (CV CInt (IO CInt)))@@ -498,16 +554,33 @@ -------------------------------------------------------------------------------- +selectG :: (TransArray c, TransArray c1, TransArray c2, Storable t, Storable a)+ => Trans c2 (Trans c1 (CInt -> Ptr t -> Trans c (CInt -> Ptr a -> IO CInt)))+ -> c2 -> c1 -> Vector t -> c -> Vector a selectG f c u v w = unsafePerformIO $ do r <- createVector (dim v)- f # c # u # v # w # r #|"selectG"+ (c # u # v # w #! r) f #|"selectG" return r +selectD :: Vector CInt -> Vector Double -> Vector Double -> Vector Double -> Vector Double selectD = selectG c_selectD+selectF :: Vector CInt -> Vector Float -> Vector Float -> Vector Float -> Vector Float selectF = selectG c_selectF+selectI :: Vector CInt -> Vector CInt -> Vector CInt -> Vector CInt -> Vector CInt selectI = selectG c_selectI+selectL :: Vector CInt -> Vector Z -> Vector Z -> Vector Z -> Vector Z selectL = selectG c_selectL+selectC :: Vector CInt+ -> Vector (Complex Double)+ -> Vector (Complex Double)+ -> Vector (Complex Double)+ -> Vector (Complex Double) selectC = selectG c_selectC+selectQ :: Vector CInt+ -> Vector (Complex Float)+ -> Vector (Complex Float)+ -> Vector (Complex Float)+ -> Vector (Complex Float) selectQ = selectG c_selectQ type Sel x = CV CInt (CV x (CV x (CV x (CV x (IO CInt)))))@@ -521,16 +594,29 @@ --------------------------------------------------------------------------- +remapG :: (TransArray c, TransArray c1, Storable t, Storable a)+ => (CInt -> CInt -> CInt -> CInt -> Ptr t+ -> Trans c1 (Trans c (CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)))+ -> Matrix t -> c1 -> c -> Matrix a remapG f i j m = unsafePerformIO $ do r <- createMatrix RowMajor (rows i) (cols i)- f # i # j # m # r #|"remapG"+ (i # j # m #! r) f #|"remapG" return r +remapD :: Matrix CInt -> Matrix CInt -> Matrix Double -> Matrix Double remapD = remapG c_remapD+remapF :: Matrix CInt -> Matrix CInt -> Matrix Float -> Matrix Float remapF = remapG c_remapF+remapI :: Matrix CInt -> Matrix CInt -> Matrix CInt -> Matrix CInt remapI = remapG c_remapI+remapL :: Matrix CInt -> Matrix CInt -> Matrix Z -> Matrix Z remapL = remapG c_remapL+remapC :: Matrix CInt+ -> Matrix CInt+ -> Matrix (Complex Double)+ -> Matrix (Complex Double) remapC = remapG c_remapC+remapQ :: Matrix CInt -> Matrix CInt -> Matrix (Complex Float) -> Matrix (Complex Float) remapQ = remapG c_remapQ type Rem x = OM CInt (OM CInt (OM x (OM x (IO CInt))))@@ -544,9 +630,12 @@ -------------------------------------------------------------------------------- +rowOpAux :: (TransArray c, Storable a) =>+ (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Trans c (IO CInt))+ -> Int -> a -> Int -> Int -> Int -> Int -> c -> IO () rowOpAux f c x i1 i2 j1 j2 m = do px <- newArray [x]- f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2) # m #|"rowOp"+ (m # id) (f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2)) #|"rowOp" free px type RowOp x = CInt -> Ptr x -> CInt -> CInt -> CInt -> CInt -> x ::> Ok@@ -562,7 +651,10 @@ -------------------------------------------------------------------------------- -gemmg f v m1 m2 m3 = f # v # m1 # m2 # m3 #|"gemmg"+gemmg :: (TransArray c1, TransArray c, TransArray c2, TransArray c3)+ => Trans c3 (Trans c2 (Trans c1 (Trans c (IO CInt))))+ -> c3 -> c2 -> c1 -> c -> IO ()+gemmg f v m1 m2 m3 = (v # m1 # m2 #! m3) f #|"gemmg" type Tgemm x = x :> x ::> x ::> x ::> Ok @@ -577,6 +669,37 @@ -------------------------------------------------------------------------------- +reorderAux :: (TransArray c, Storable t, Storable a1, Storable t1, Storable a) =>+ (CInt -> Ptr a -> CInt -> Ptr t1+ -> Trans c (CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt))+ -> Vector t1 -> c -> Vector t -> Vector a1+reorderAux f s d v = unsafePerformIO $ do+ k <- createVector (dim s)+ r <- createVector (dim v)+ (k # s # d # v #! r) f #| "reorderV"+ return r++type Reorder x = CV CInt (CV CInt (CV CInt (CV x (CV x (IO CInt)))))++foreign import ccall unsafe "reorderD" c_reorderD :: Reorder Double+foreign import ccall unsafe "reorderF" c_reorderF :: Reorder Float+foreign import ccall unsafe "reorderI" c_reorderI :: Reorder CInt+foreign import ccall unsafe "reorderC" c_reorderC :: Reorder (Complex Double)+foreign import ccall unsafe "reorderQ" c_reorderQ :: Reorder (Complex Float)+foreign import ccall unsafe "reorderL" c_reorderL :: Reorder Z++-- | Transpose an array with dimensions @dims@ by making a copy using @strides@. For example, for an array with 3 indices,+-- @(reorderVector strides dims v) ! ((i * dims ! 1 + j) * dims ! 2 + k) == v ! (i * strides ! 0 + j * strides ! 1 + k * strides ! 2)@+-- This function is intended to be used internally by tensor libraries.+reorderVector :: Element a+ => Vector CInt -- ^ @strides@: array strides+ -> Vector CInt -- ^ @dims@: array dimensions of new array @v@+ -> Vector a -- ^ @v@: flattened input array+ -> Vector a -- ^ @v'@: flattened output array+reorderVector = reorderV++--------------------------------------------------------------------------------+ foreign import ccall unsafe "saveMatrix" c_saveMatrix :: CString -> CString -> Double ::> Ok @@ -590,10 +713,9 @@ saveMatrix name format m = do cname <- newCString name cformat <- newCString format- c_saveMatrix cname cformat # m #|"saveMatrix"+ (m # id) (c_saveMatrix cname cformat) #|"saveMatrix" free cname free cformat return () ---------------------------------------------------------------------------------
src/Internal/Modular.hs view
@@ -1,17 +1,15 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-}-{-# LANGUAGE KindSignatures #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-}-{-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE Rank2Types #-} {-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE UndecidableInstances #-}-{-# LANGUAGE GADTs #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE TypeOperators #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-}+ {- | Module : Internal.Modular Copyright : (c) Alberto Ruiz 2015@@ -38,13 +36,20 @@ gaussElim, gaussElim_1, gaussElim_2, luST, luSolve', luPacked', magnit, invershur) import Internal.ST(mutable)+#if MIN_VERSION_base(4,11,0)+import GHC.TypeLits hiding (Mod)+#else import GHC.TypeLits+#endif import Data.Proxy(Proxy) import Foreign.ForeignPtr(castForeignPtr) import Foreign.Storable import Data.Ratio import Data.Complex import Control.DeepSeq ( NFData(..) )+#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif @@ -72,11 +77,11 @@ where compare a b = compare (unMod a) (unMod b) -instance (Integral 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)@@ -84,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@@ -114,7 +119,7 @@ where show = show . unMod -instance forall n t . (Integral t, KnownNat n) => Num (Mod n t)+instance (Integral t, KnownNat n) => Num (Mod n t) where (+) = l2 (\m a b -> (a + b) `mod` (fromIntegral m)) (*) = l2 (\m a b -> (a * b) `mod` (fromIntegral m))@@ -140,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@@ -157,13 +163,14 @@ 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 forall m . KnownNat m => CTrans (Mod m I)-instance forall m . KnownNat m => CTrans (Mod m Z)+instance KnownNat m => CTrans (Mod m I)+instance KnownNat m => CTrans (Mod m Z) -instance forall m . KnownNat m => Container Vector (Mod m I)+instance KnownNat m => Container Vector (Mod m I) where conj' = id size' = dim@@ -203,7 +210,7 @@ fromZ' = vmod . fromZ' toZ' = toZ' . f2i -instance forall m . KnownNat m => Container Vector (Mod m Z)+instance KnownNat m => Container Vector (Mod m Z) where conj' = id size' = dim@@ -311,7 +318,7 @@ lift2m f a b = liftMatrix vmod (f (f2iM a) (f2iM b)) -instance forall m . KnownNat m => Num (Vector (Mod m I))+instance KnownNat m => Num (Vector (Mod m I)) where (+) = lift2 (+) (*) = lift2 (*)@@ -321,7 +328,7 @@ negate = lift1 negate fromInteger x = fromInt (fromInteger x) -instance forall m . KnownNat m => Num (Vector (Mod m Z))+instance KnownNat m => Num (Vector (Mod m Z)) where (+) = lift2 (+) (*) = lift2 (*)
src/Internal/Numeric.hs view
@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-}@@ -6,6 +5,8 @@ {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE UndecidableInstances #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-}+ ----------------------------------------------------------------------------- -- | -- Module : Data.Packed.Internal.Numeric@@ -26,6 +27,7 @@ import Internal.Vectorized import Internal.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ,multiplyI,multiplyL) import Data.List.Split(chunksOf)+import qualified Data.Vector.Storable as V -------------------------------------------------------------------------------- @@ -103,8 +105,8 @@ add' = vectorZipI Add sub = vectorZipI Sub mul = vectorZipI Mul- equal u v = dim u == dim v && maxElement' (vectorMapI Abs (sub u v)) == 0- scalar' x = fromList [x]+ equal = (==)+ scalar' = V.singleton konst' = constantD build' = buildV cmap' = mapVector@@ -142,8 +144,8 @@ add' = vectorZipL Add sub = vectorZipL Sub mul = vectorZipL Mul- equal u v = dim u == dim v && maxElement' (vectorMapL Abs (sub u v)) == 0- scalar' x = fromList [x]+ equal = (==)+ scalar' = V.singleton konst' = constantD build' = buildV cmap' = mapVector@@ -182,8 +184,8 @@ add' = vectorZipF Add sub = vectorZipF Sub mul = vectorZipF Mul- equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0- scalar' x = fromList [x]+ equal = (==)+ scalar' = V.singleton konst' = constantD build' = buildV cmap' = mapVector@@ -219,8 +221,8 @@ add' = vectorZipR Add sub = vectorZipR Sub mul = vectorZipR Mul- equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0- scalar' x = fromList [x]+ equal = (==)+ scalar' = V.singleton konst' = constantD build' = buildV cmap' = mapVector@@ -256,8 +258,8 @@ add' = vectorZipC Add sub = vectorZipC Sub mul = vectorZipC Mul- equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0- scalar' x = fromList [x]+ equal = (==)+ scalar' = V.singleton konst' = constantD build' = buildV cmap' = mapVector@@ -292,8 +294,8 @@ add' = vectorZipQ Add sub = vectorZipQ Sub mul = vectorZipQ Mul- equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0- scalar' x = fromList [x]+ equal = (==)+ scalar' = V.singleton konst' = constantD build' = buildV cmap' = mapVector@@ -788,13 +790,7 @@ type instance RealOf I = I type instance RealOf Z = Z -type family ComplexOf x--type instance ComplexOf Double = Complex Double-type instance ComplexOf (Complex Double) = Complex Double--type instance ComplexOf Float = Complex Float-type instance ComplexOf (Complex Float) = Complex Float+type ComplexOf x = Complex (RealOf x) type family SingleOf x
src/Internal/Random.hs view
@@ -31,13 +31,13 @@ gaussianSample :: Seed -> Int -- ^ number of rows -> Vector Double -- ^ mean vector- -> Matrix Double -- ^ covariance matrix+ -> Herm Double -- ^ covariance matrix -> Matrix Double -- ^ result gaussianSample seed n med cov = m where c = dim med meds = konst' 1 n `outer` med rs = reshape c $ randomVector seed Gaussian (c * n)- m = rs `mXm` cholSH cov `add` meds+ m = rs `mXm` chol cov `add` meds -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate -- uniform distribution.
src/Internal/ST.hs view
@@ -81,6 +81,8 @@ unsafeFreezeVector (STVector x) = unsafeIOToST . return $ x {-# INLINE safeIndexV #-}+safeIndexV :: Storable t2+ => (STVector s t2 -> Int -> t) -> STVector t1 t2 -> Int -> t safeIndexV f (STVector v) k | k < 0 || k>= dim v = error $ "out of range error in vector (dim=" ++show (dim v)++", pos="++show k++")"@@ -150,9 +152,12 @@ freezeMatrix :: (Element t) => STMatrix s t -> ST s (Matrix t) freezeMatrix m = liftSTMatrix id m +cloneMatrix :: Element t => Matrix t -> IO (Matrix t) cloneMatrix m = copy (orderOf m) m {-# INLINE safeIndexM #-}+safeIndexM :: (STMatrix s t2 -> Int -> Int -> t)+ -> STMatrix t1 t2 -> Int -> Int -> t safeIndexM f (STMatrix m) r c | r<0 || r>=rows m || c<0 || c>=cols m = error $ "out of range error in matrix (size="@@ -184,6 +189,7 @@ | Col Int | FromCol Int +getColRange :: Int -> ColRange -> (Int, Int) getColRange c AllCols = (0,c-1) getColRange c (ColRange a b) = (a `mod` c, b `mod` c) getColRange c (Col a) = (a `mod` c, a `mod` c)@@ -194,6 +200,7 @@ | Row Int | FromRow Int +getRowRange :: Int -> RowRange -> (Int, Int) getRowRange r AllRows = (0,r-1) getRowRange r (RowRange a b) = (a `mod` r, b `mod` r) getRowRange r (Row a) = (a `mod` r, a `mod` r)@@ -223,6 +230,7 @@ i2' = i2 `mod` (rows m) +extractMatrix :: Element a => STMatrix t a -> RowRange -> ColRange -> ST s (Matrix a) extractMatrix (STMatrix m) rr rc = unsafeIOToST (extractR (orderOf m) m 0 (idxs[i1,i2]) 0 (idxs[j1,j2])) where (i1,i2) = getRowRange (rows m) rr@@ -231,6 +239,7 @@ -- | r0 c0 height width data Slice s t = Slice (STMatrix s t) Int Int Int Int +slice :: Element a => Slice t a -> Matrix a slice (Slice (STMatrix m) r0 c0 nr nc) = subMatrix (r0,c0) (nr,nc) m gemmm :: Element t => t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()@@ -238,12 +247,11 @@ where res = unsafeIOToST (gemm v a b r) v = fromList [alpha,beta]- + mutable :: Element t => (forall s . (Int, Int) -> STMatrix s t -> ST s u) -> Matrix t -> (Matrix t,u) mutable f a = runST $ do x <- thawMatrix a info <- f (rows a, cols a) x r <- unsafeFreezeMatrix x return (r,info)-
src/Internal/Sparse.hs view
@@ -1,9 +1,11 @@-{-# LANGUAGE RecordWildCards #-}-{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE BangPatterns #-} {-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE RecordWildCards #-} module Internal.Sparse(- GMatrix(..), CSR(..), mkCSR, fromCSR,+ GMatrix(..), CSR(..), mkCSR, fromCSR, impureCSR, mkSparse, mkDiagR, mkDense, AssocMatrix, toDense,@@ -14,10 +16,12 @@ import Internal.Matrix import Internal.Numeric import qualified Data.Vector.Storable as V-import Data.Function(on)+import qualified Data.Vector.Storable.Mutable as M import Control.Arrow((***))-import Control.Monad(when)-import Data.List(groupBy, sort)+import Control.Monad(when, foldM)+import Control.Monad.ST (runST)+import Control.Monad.Primitive (PrimMonad)+import Data.List(sort) import Foreign.C.Types(CInt(..)) import Internal.Devel@@ -25,10 +29,7 @@ import Foreign(Ptr) import Text.Printf(printf) -infixl 0 ~!~-c ~!~ msg = when c (error msg)--type AssocMatrix = [((Int,Int),Double)]+type AssocMatrix = [(IndexOf Matrix, Double)] data CSR = CSR { csrVals :: Vector Double@@ -47,23 +48,79 @@ } deriving Show +-- | Produce a CSR sparse matrix from a association matrix. mkCSR :: AssocMatrix -> CSR-mkCSR sm' = CSR{..}+mkCSR ms =+ runST $ impureCSR runFold $ sort ms+ where+ runFold next initialise xtract as0 = do+ i0 <- initialise+ acc <- foldM next i0 as0+ xtract acc++-- | Produce a CSR sparse matrix by applying a generic folding function.+--+-- This allows one to build a CSR from an effectful streaming source+-- when combined with libraries like pipes, io-streams, or streaming.+--+-- For example+--+-- > impureCSR Pipes.Prelude.foldM :: PrimMonad m => Producer AssocEntry m () -> m CSR+-- > impureCSR Streaming.Prelude.foldM :: PrimMonad m => Stream (Of AssocEntry) m r -> m (Of CSR r)+--+impureCSR+ :: PrimMonad m+ => (forall x . (x -> (IndexOf Matrix, Double) -> m x) -> m x -> (x -> m CSR) -> r)+ -> r+impureCSR f = f next begin done where- sm = sort sm'- rws = map ((fromList *** fromList)- . unzip- . map ((succ.fi.snd) *** id)- )- . groupBy ((==) `on` (fst.fst))- $ sm- rszs = map (fi . dim . fst) rws- csrRows = fromList (scanl (+) 1 rszs)- csrVals = vjoin (map snd rws)- csrCols = vjoin (map fst rws)- csrNRows = dim csrRows - 1- csrNCols = fromIntegral (V.maximum csrCols)+ sfi = succ . fi+ begin = do+ mv <- M.unsafeNew 64+ mr <- M.unsafeNew 64+ mc <- M.unsafeNew 64+ return (mv, mr, mc, 0, 0, 0, -1) + next (!mv, !mr, !mc, !idxVC, !idxR, !maxC, !curRow) ((r,c),d) = do+ when (r < curRow) $+ error (printf "impureCSR: row %i specified after %i" r curRow)++ let lenVC = M.length mv+ lenR = M.length mr+ maxC' = max maxC c++ (mv', mc') <-+ if idxVC >= lenVC then do+ mv' <- M.unsafeGrow mv lenVC+ mc' <- M.unsafeGrow mc lenVC+ return (mv', mc')+ else+ return (mv, mc)++ mr' <-+ if idxR >= lenR - 1 then+ M.unsafeGrow mr lenR+ else+ return mr++ M.unsafeWrite mc' idxVC (sfi c)+ M.unsafeWrite mv' idxVC d++ idxR' <-+ foldM+ (\idxR' _ -> idxR' + 1 <$ M.unsafeWrite mr' idxR' (sfi idxVC))+ idxR [1 .. (r-curRow)]++ return (mv', mr', mc', idxVC + 1, idxR', maxC', r)++ done (!mv, !mr, !mc, !idxVC, !idxR, !maxC, !curR) = do+ M.unsafeWrite mr idxR (sfi idxVC)+ vv <- V.unsafeFreeze (M.unsafeTake idxVC mv)+ vc <- V.unsafeFreeze (M.unsafeTake idxVC mc)+ vr <- V.unsafeFreeze (M.unsafeTake (idxR + 1) mr)+ return $ CSR vv vc vr (succ curR) (succ maxC)++ {- | General matrix with specialized internal representations for dense, sparse, diagonal, banded, and constant elements. @@ -122,11 +179,12 @@ fromCSR :: CSR -> GMatrix fromCSR csr = SparseR {..} where- gmCSR @ CSR {..} = csr+ gmCSR@CSR {..} = csr nRows = csrNRows nCols = csrNCols +mkDiagR :: Int -> Int -> Vector Double -> GMatrix mkDiagR r c v | dim v <= min r c = Diag{..} | otherwise = error $ printf "mkDiagR: incorrect sizes (%d,%d) [%d]" r c (dim v)@@ -142,15 +200,19 @@ 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- c_smXv # csrVals # csrCols # csrRows # v # r #|"CSRXv"+ (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- c_smTXv # cscVals # cscRows # cscCols # v # r #|"CSCXv"+ (cscVals # cscRows # cscCols # v #! r) c_smTXv #|"CSCXv" return r gmXv Diag{..} v@@ -170,8 +232,10 @@ {- | general matrix - vector product >>> let m = mkSparse [((0,999),1.0),((1,1999),2.0)]+m :: GMatrix >>> m !#> vector [1..2000]-fromList [1000.0,4000.0]+[1000.0,4000.0]+it :: Vector Double -} infixr 8 !#>@@ -211,4 +275,3 @@ tr (Diag v n m) = Diag v m n tr (Dense a n m) = Dense (tr a) m n tr' = tr-
src/Internal/Static.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE CPP #-} #if __GLASGOW_HASKELL__ >= 708 {-# LANGUAGE DataKinds #-}@@ -11,7 +12,11 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE ViewPatterns #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DeriveGeneric #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-}+ {- | Module : Internal.Static Copyright : (c) Alberto Ruiz 2006-14@@ -28,18 +33,35 @@ import Numeric.LinearAlgebra hiding (konst,size,R,C) import Internal.Vector as D hiding (R,C) import Internal.ST+import Control.DeepSeq import Data.Proxy(Proxy) import Foreign.Storable(Storable) import Text.Printf +import Data.Binary+import GHC.Generics (Generic)+import Data.Proxy (Proxy(..))+ -------------------------------------------------------------------------------- type ℝ = Double type ℂ = Complex Double newtype Dim (n :: Nat) t = Dim t- deriving Show+ deriving (Show, Generic) +instance (KnownNat n, Binary a) => Binary (Dim n a) where+ get = do+ k <- get+ let n = natVal (Proxy :: Proxy n)+ if n == k+ then Dim <$> get+ else fail ("Expected dimension " ++ (show n) ++ ", but found dimension " ++ (show k))++ put (Dim x) = do+ put (natVal (Proxy :: Proxy n))+ put x+ lift1F :: (c t -> c t) -> Dim n (c t) -> Dim n (c t)@@ -50,18 +72,22 @@ -> Dim n (c t) -> Dim n (c t) -> Dim n (c t) lift2F f (Dim u) (Dim v) = Dim (f u v) +instance NFData t => NFData (Dim n t) where+ rnf (Dim (force -> !_)) = ()+ -------------------------------------------------------------------------------- newtype R n = R (Dim n (Vector ℝ))- deriving (Num,Fractional,Floating)+ deriving (Num,Fractional,Floating,Generic,Binary) newtype C n = C (Dim n (Vector ℂ))- deriving (Num,Fractional,Floating)+ deriving (Num,Fractional,Floating,Generic) newtype L m n = L (Dim m (Dim n (Matrix ℝ)))--newtype M m n = M (Dim m (Dim n (Matrix ℂ)))+ deriving (Generic, Binary) +newtype M m n = M (Dim m (Dim n (Matrix ℂ)))+ deriving (Generic) mkR :: Vector ℝ -> R n mkR = R . Dim@@ -75,6 +101,18 @@ mkM :: Matrix ℂ -> M m n mkM x = M (Dim (Dim x)) +instance NFData (R n) where+ rnf (R (force -> !_)) = ()++instance NFData (C n) where+ rnf (C (force -> !_)) = ()++instance NFData (L n m) where+ rnf (L (force -> !_)) = ()++instance NFData (M n m) where+ rnf (M (force -> !_)) = ()+ -------------------------------------------------------------------------------- type V n t = Dim n (Vector t)@@ -92,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 @@ -173,7 +211,7 @@ singleM m = rows m == 1 && cols m == 1 -instance forall n. KnownNat n => Sized ℂ (C n) Vector+instance KnownNat n => Sized ℂ (C n) Vector where size _ = fromIntegral . natVal $ (undefined :: Proxy n) konst x = mkC (LA.scalar x)@@ -189,7 +227,7 @@ r = mkC v :: C n -instance forall n. KnownNat n => Sized ℝ (R n) Vector+instance KnownNat n => Sized ℝ (R n) Vector where size _ = fromIntegral . natVal $ (undefined :: Proxy n) konst x = mkR (LA.scalar x)@@ -206,7 +244,7 @@ -instance forall m n . (KnownNat m, KnownNat n) => Sized ℝ (L m n) Matrix+instance (KnownNat m, KnownNat n) => Sized ℝ (L m n) Matrix where size _ = ((fromIntegral . natVal) (undefined :: Proxy m) ,(fromIntegral . natVal) (undefined :: Proxy n))@@ -224,7 +262,7 @@ r = mkL x :: L m n -instance forall m n . (KnownNat m, KnownNat n) => Sized ℂ (M m n) Matrix+instance (KnownNat m, KnownNat n) => Sized ℂ (M m n) Matrix where size _ = ((fromIntegral . natVal) (undefined :: Proxy m) ,(fromIntegral . natVal) (undefined :: Proxy n))@@ -282,43 +320,43 @@ -------------------------------------------------------------------------------- -instance forall n . KnownNat n => Show (R n)+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 forall n . KnownNat n => Show (C n)+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 forall m n . (KnownNat m, KnownNat n) => Show (L m n)+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 forall m n . (KnownNat m, KnownNat n) => Show (M m n)+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 -------------------------------------------------------------------------------- -instance forall n t . (Num (Vector t), Numeric t )=> Num (Dim n (Vector t))+instance (Num (Vector t), Numeric t )=> Num (Dim n (Vector t)) where (+) = lift2F (+) (*) = lift2F (*)@@ -328,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 (/)@@ -353,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) (*)@@ -363,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@@ -482,6 +520,18 @@ (**) = lift2MD (**) pi = M pi +instance Additive (R n) where+ add = (+)++instance Additive (C n) where+ add = (+)++instance (KnownNat m, KnownNat n) => Additive (L m n) where+ add = (+)++instance (KnownNat m, KnownNat n) => Additive (M m n) where+ add = (+)+ -------------------------------------------------------------------------------- @@ -518,6 +568,17 @@ putStr "C " >> putStr (tail . dropWhile (/='x') $ su) --------------------------------------------------------------------------------++overMatL' :: (KnownNat m, KnownNat n)+ => (LA.Matrix ℝ -> LA.Matrix ℝ) -> L m n -> L m n+overMatL' f = mkL . f . unwrap+{-# INLINE overMatL' #-}++overMatM' :: (KnownNat m, KnownNat n)+ => (LA.Matrix ℂ -> LA.Matrix ℂ) -> M m n -> M m n+overMatM' f = mkM . f . unwrap+{-# INLINE overMatM' #-}+ #else
src/Internal/Util.hs view
@@ -1,9 +1,13 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE ViewPatterns #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-}+{-# OPTIONS_GHC -fno-warn-orphans #-} ----------------------------------------------------------------------------- {- |@@ -41,7 +45,7 @@ ℕ,ℤ,ℝ,ℂ,iC, Normed(..), norm_Frob, norm_nuclear, magnit,- unitary,+ normalize, mt, (~!~), pairwiseD2,@@ -79,6 +83,9 @@ import Data.Complex import Data.Function(on) import Internal.ST+#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif type ℝ = Double type ℕ = Int@@ -92,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@@ -255,6 +263,7 @@ -- ^ 2-norm of real vector norm = pnorm PNorm2 +-- | p-norm for vectors, operator norm for matrices class Normed a where norm_0 :: a -> R@@ -319,10 +328,11 @@ norm_2 = norm_2 . double norm_Inf = norm_Inf . double -+-- | Frobenius norm (Schatten p-norm with p=2) norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R norm_Frob = norm_2 . flatten +-- | Sum of singular values (Schatten p-norm with p=1) norm_nuclear :: Field t => Matrix t -> R norm_nuclear = sumElements . singularValues @@ -341,8 +351,8 @@ -- | Obtains a vector in the same direction with 2-norm=1-unitary :: Vector Double -> Vector Double-unitary v = v / scalar (norm v)+normalize :: (Normed (Vector t), Num (Vector t), Field t) => Vector t -> Vector t+normalize v = v / real (scalar (norm_2 v)) -- | trans . inv@@ -369,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@@ -406,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 @@ -611,6 +622,9 @@ s1 = fromRows $ pivotDown (rows x) 0 rs -- interesting s2 = pivotUp (rows x-1) (toRows $ flipud s1) +pivotDown+ :: forall t . (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)+ => Int -> Int -> [Vector t] -> [Vector t] pivotDown t n xs | t == n = [] | otherwise = y : pivotDown t (n+1) ys@@ -618,8 +632,9 @@ 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) | p == 0 = error "gauss: singular!" -- FIXME | otherwise = u : map f zs@@ -630,12 +645,16 @@ redu (_,[]) = [] +pivotUp+ :: forall t . (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)+ => Int -> [Vector t] -> [Vector t] pivotUp n xs | n == -1 = [] | otherwise = y : pivotUp (n-1) ys where y:ys = redu' (n,xs) + redu' :: (Int, [Vector t]) -> [Vector t] redu' (k,x:zs) = u : map f zs where u = x@@ -893,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
@@ -1,6 +1,7 @@-{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-}+{-# LANGUAGE MagicHash, UnboxedTuples, BangPatterns, FlexibleContexts #-} {-# LANGUAGE TypeSynonymInstances #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} -- | -- Module : Internal.Vector@@ -39,12 +40,11 @@ import qualified Data.Vector.Storable as Vector import Data.Vector.Storable(Vector, fromList, unsafeToForeignPtr, unsafeFromForeignPtr, unsafeWith) -#ifdef BINARY import Data.Binary+import Data.Binary.Put import Control.Monad(replicateM) import qualified Data.ByteString.Internal as BS import Data.Vector.Storable.Internal(updPtr)-#endif type I = CInt type Z = Int64@@ -64,13 +64,13 @@ -- | Number of elements dim :: (Storable t) => Vector t -> Int dim = Vector.length+{-# INLINE dim #-} -- C-Haskell vector adapter {-# INLINE avec #-}-avec :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b-avec f v = inlinePerformIO (unsafeWith v (return . f (fromIntegral (Vector.length v))))-infixl 1 `avec`+avec :: Storable a => Vector a -> (f -> IO r) -> ((CInt -> Ptr a -> f) -> IO r)+avec v f g = unsafeWith v $ \ptr -> f (g (fromIntegral (Vector.length v)) ptr) -- allocates memory for a new vector createVector :: Storable a => Int -> IO (Vector a)@@ -94,6 +94,7 @@ -} +safeRead :: Storable a => Vector a -> (Ptr a -> IO c) -> c safeRead v = inlinePerformIO . unsafeWith v {-# INLINE safeRead #-} @@ -115,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@@ -134,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@@ -142,6 +145,7 @@ -> Vector t -- ^ source -> Vector t -- ^ result subVector = Vector.slice+{-# INLINE subVector #-} @@ -167,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@@ -189,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]@@ -201,7 +207,7 @@ --------------------------------------------------------------- --- | transforms a complex vector into a real vector with alternating real and imaginary parts +-- | transforms a complex vector into a real vector with alternating real and imaginary parts asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n) where (fp,i,n) = unsafeToForeignPtr v@@ -246,7 +252,7 @@ {-# INLINE zipVectorWith #-} -- | unzipWith for Vectors-unzipVectorWith :: (Storable (a,b), Storable c, Storable d) +unzipVectorWith :: (Storable (a,b), Storable c, Storable d) => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d) unzipVectorWith f u = unsafePerformIO $ do let n = dim u@@ -257,7 +263,7 @@ unsafeWith w $ \pw -> do let go (-1) = return () go !k = do z <- peekElemOff pu k- let (x,y) = f z + let (x,y) = f z pokeElemOff pv k x pokeElemOff pw k y go (k-1)@@ -284,11 +290,13 @@ go (dim v -1) x {-# INLINE foldVectorWithIndex #-} +foldLoop :: (Int -> t -> t) -> t -> Int -> t foldLoop f s0 d = go (d - 1) s0 where go 0 s = f (0::Int) s go !j !s = go (j - 1) (f j s) +foldVectorG :: Storable t1 => (Int -> (Int -> t1) -> t -> t) -> t -> Vector t1 -> t foldVectorG f s0 v = foldLoop g s0 (dim v) where g !k !s = f k (safeRead v . flip peekElemOff) s {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479)@@ -305,11 +313,11 @@ return w where mapVectorM' w' !k !t | k == t = do- x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k y <- f x return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y | otherwise = do- x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k y <- f x _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y mapVectorM' w' (k+1) t@@ -324,7 +332,7 @@ x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k f x | otherwise = do- x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k _ <- f x mapVectorM' (k+1) t {-# INLINE mapVectorM_ #-}@@ -338,11 +346,11 @@ return w where mapVectorM' w' !k !t | k == t = do- x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k y <- f k x return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y | otherwise = do- x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k y <- f k x _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y mapVectorM' w' (k+1) t@@ -357,7 +365,7 @@ x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k f k x | otherwise = do- x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k _ <- f k x mapVectorM' (k+1) t {-# INLINE mapVectorWithIndexM_ #-}@@ -380,7 +388,6 @@ -------------------------------------------------------------------------------- -#ifdef BINARY -- a 64K cache, with a Double taking 13 bytes in Bytestring, -- implies a chunk size of 5041@@ -392,8 +399,10 @@ m = d `mod` chunk in if m /= 0 then reverse (m:(replicate c chunk)) else (replicate c chunk) +putVector :: (Storable t, Binary t) => Vector t -> Data.Binary.Put.PutM () putVector v = mapM_ put $! toList v +getVector :: (Storable a, Binary a) => Int -> Get (Vector a) getVector d = do xs <- replicateM d get return $! fromList xs@@ -432,7 +441,6 @@ -- get = fmap bs2v get -#endif -------------------------------------------------------------------@@ -458,4 +466,3 @@ unzipVector = unzipVectorWith id --------------------------------------------------------------------
src/Internal/Vectorized.hs view
@@ -27,10 +27,16 @@ import System.IO.Unsafe(unsafePerformIO) import Control.Monad(when) -infixl 1 #+infixr 1 #+(#) :: TransArray c => c -> (b -> IO r) -> TransRaw c b -> IO r a # b = applyRaw a b {-# INLINE (#) #-} +(#!) :: (TransArray c, TransArray c1) => c1 -> c -> TransRaw c1 (TransRaw c (IO r)) -> IO r+a #! b = a # b # id+{-# INLINE (#!) #-}++fromei :: Enum a => a -> CInt fromei x = fromIntegral (fromEnum x) :: CInt data FunCodeV = Sin@@ -97,13 +103,23 @@ sumC :: Vector (Complex Double) -> Complex Double sumC = sumg c_sumC +sumI :: ( TransRaw c (CInt -> Ptr a -> IO CInt) ~ (CInt -> Ptr I -> I :> Ok)+ , TransArray c+ , Storable a+ )+ => I -> c -> a sumI m = sumg (c_sumI m) +sumL :: ( TransRaw c (CInt -> Ptr a -> IO CInt) ~ (CInt -> Ptr Z -> Z :> Ok)+ , TransArray c+ , Storable a+ ) => Z -> c -> a sumL m = sumg (c_sumL m) +sumg :: (TransArray c, Storable a) => TransRaw c (CInt -> Ptr a -> IO CInt) -> c -> a sumg f x = unsafePerformIO $ do r <- createVector 1- f # x # r #| "sum"+ (x #! r) f #| "sum" return $ r @> 0 type TVV t = t :> t :> Ok@@ -137,9 +153,11 @@ prodL :: Z-> Vector Z -> Z prodL = prodg . c_prodL +prodg :: (TransArray c, Storable a)+ => TransRaw c (CInt -> Ptr a -> IO CInt) -> c -> a prodg f x = unsafePerformIO $ do r <- createVector 1- f # x # r #| "prod"+ (x #! r) f #| "prod" return $ r @> 0 @@ -152,26 +170,38 @@ ------------------------------------------------------------------ +toScalarAux :: (Enum a, TransArray c, Storable a1)+ => (CInt -> TransRaw c (CInt -> Ptr a1 -> IO CInt)) -> a -> c -> a1 toScalarAux fun code v = unsafePerformIO $ do r <- createVector 1- fun (fromei code) # v # r #|"toScalarAux"+ (v #! r) (fun (fromei code)) #|"toScalarAux" return (r @> 0) ++vectorMapAux :: (Enum a, Storable t, Storable a1)+ => (CInt -> CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt)+ -> a -> Vector t -> Vector a1 vectorMapAux fun code v = unsafePerformIO $ do r <- createVector (dim v)- fun (fromei code) # v # r #|"vectorMapAux"+ (v #! r) (fun (fromei code)) #|"vectorMapAux" return r +vectorMapValAux :: (Enum a, Storable a2, Storable t, Storable a1)+ => (CInt -> Ptr a2 -> CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt)+ -> a -> a2 -> Vector t -> Vector a1 vectorMapValAux fun code val v = unsafePerformIO $ do r <- createVector (dim v) pval <- newArray [val]- fun (fromei code) pval # v # r #|"vectorMapValAux"+ (v #! r) (fun (fromei code) pval) #|"vectorMapValAux" free pval return r +vectorZipAux :: (Enum a, TransArray c, Storable t, Storable a1)+ => (CInt -> CInt -> Ptr t -> TransRaw c (CInt -> Ptr a1 -> IO CInt))+ -> a -> Vector t -> c -> Vector a1 vectorZipAux fun code u v = unsafePerformIO $ do r <- createVector (dim u)- fun (fromei code) # u # v # r #|"vectorZipAux"+ (u # v #! r) (fun (fromei code)) #|"vectorZipAux" return r ---------------------------------------------------------------------@@ -368,16 +398,17 @@ -> Vector Double randomVector seed dist n = unsafePerformIO $ do r <- createVector n- c_random_vector (fi seed) ((fi.fromEnum) dist) # r #|"randomVector"+ (r # id) (c_random_vector (fi seed) ((fi.fromEnum) dist)) #|"randomVector" return r foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> Double :> Ok -------------------------------------------------------------------------------- +roundVector :: Vector Double -> Vector Double roundVector v = unsafePerformIO $ do r <- createVector (dim v)- c_round_vector # v # r #|"roundVector"+ (v #! r) c_round_vector #|"roundVector" return r foreign import ccall unsafe "round_vector" c_round_vector :: TVV Double@@ -386,12 +417,13 @@ -- | -- >>> range 5--- fromList [0,1,2,3,4]+-- [0,1,2,3,4]+-- it :: Vector I -- range :: Int -> Vector I range n = unsafePerformIO $ do r <- createVector n- c_range_vector # r #|"range"+ (r # id) c_range_vector #|"range" return r foreign import ccall unsafe "range_vector" c_range_vector :: CInt :> Ok@@ -429,9 +461,11 @@ long2intV = tog c_long2int +tog :: (Storable t, Storable a)+ => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a tog f v = unsafePerformIO $ do r <- createVector (dim v)- f # v # r #|"tog"+ (v #! r) f #|"tog" return r foreign import ccall unsafe "float2double" c_float2double :: Float :> Double :> Ok@@ -448,9 +482,11 @@ --------------------------------------------------------------- +stepg :: (Storable t, Storable a)+ => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a stepg f v = unsafePerformIO $ do r <- createVector (dim v)- f # v # r #|"step"+ (v #! r) f #|"step" return r stepD :: Vector Double -> Vector Double@@ -473,9 +509,11 @@ -------------------------------------------------------------------------------- +conjugateAux :: (Storable t, Storable a)+ => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a conjugateAux fun x = unsafePerformIO $ do v <- createVector (dim x)- fun # x # v #|"conjugateAux"+ (x #! v) fun #|"conjugateAux" return v conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)@@ -493,15 +531,17 @@ let n = dim v r <- createVector n let f _ s _ d = copyArray d s n >> return 0- f # v # r #|"cloneVector"+ (v #! r) f #|"cloneVector" return r -------------------------------------------------------------------------------- +constantAux :: (Storable a1, Storable a)+ => (Ptr a1 -> CInt -> Ptr a -> IO CInt) -> a1 -> Int -> Vector a constantAux fun x n = unsafePerformIO $ do v <- createVector n px <- newArray [x]- fun px # v #|"constantAux"+ (v # id) (fun px) #|"constantAux" free px return v @@ -515,4 +555,3 @@ foreign import ccall unsafe "constantL" cconstantL :: TConst Z -----------------------------------------------------------------------
src/Numeric/LinearAlgebra.hs view
@@ -1,5 +1,8 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-}+ ----------------------------------------------------------------------------- {- | Module : Numeric.LinearAlgebra@@ -22,10 +25,12 @@ -- * Numeric classes -- |- -- The standard numeric classes are defined elementwise:+ -- The standard numeric classes are defined elementwise (commonly referred to+ -- 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)@@ -93,6 +98,11 @@ ldlSolve, ldlPacked, -- ** Positive definite cholSolve,+ -- ** Triangular+ UpLo(..),+ triSolve,+ -- ** Tridiagonal+ triDiagSolve, -- ** Sparse cgSolve, cgSolve',@@ -116,16 +126,17 @@ svd, thinSVD, compactSVD,+ compactSVDTol, singularValues, leftSV, rightSV, -- * Eigendecomposition- eig, eigSH,- eigenvalues, eigenvaluesSH,+ eig, geig, eigSH,+ eigenvalues, geigenvalues, eigenvaluesSH, geigSH, -- * QR- qr, rq, qrRaw, qrgr,+ qr, thinQR, rq, thinRQ, qrRaw, qrgr, -- * Cholesky chol, mbChol,@@ -152,7 +163,7 @@ Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample, -- * Misc- meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, magnit,+ meanCov, rowOuters, pairwiseD2, normalize, peps, relativeError, magnit, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, iC, sym, mTm, trustSym, unSym, -- * Auxiliary classes@@ -183,6 +194,9 @@ import Internal.Sparse((!#>)) import Internal.CG import Internal.Conversion+#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif {- | dense matrix product
src/Numeric/LinearAlgebra/Devel.hs view
@@ -12,16 +12,9 @@ -------------------------------------------------------------------------------- module Numeric.LinearAlgebra.Devel(- -- * FFI helpers- -- | Sample usage, to upload a perspective matrix to a shader.- --- -- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3)- -- @- module Internal.Foreign,- -- * FFI tools -- | See @examples/devel@ in the repository.- + createVector, createMatrix, TransArray(..), MatrixOrder(..), orderOf, cmat, fmat,@@ -34,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,@@ -58,15 +51,14 @@ liftMatrix, liftMatrix2, liftMatrix2Auto, -- * Sparse representation- CSR(..), fromCSR, mkCSR,+ CSR(..), fromCSR, mkCSR, impureCSR, GMatrix(..), -- * Misc- toByteString, fromByteString, showInternal+ toByteString, fromByteString, showInternal, reorderVector ) where -import Internal.Foreign import Internal.Devel import Internal.ST import Internal.Vector
src/Numeric/LinearAlgebra/HMatrix.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE CPP #-} -------------------------------------------------------------------------------- {- | Module : Numeric.LinearAlgebra.HMatrix@@ -19,12 +20,17 @@ import Numeric.LinearAlgebra import Internal.Util import Internal.Algorithms(cholSH, mbCholSH, eigSH', eigenvaluesSH', geigSH')+#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif infixr 8 <·> (<·>) :: Numeric t => Vector t -> Vector t -> t (<·>) = dot +app :: Numeric t => Matrix t -> Vector t -> Vector t app m v = m #> v +mul :: Numeric t => Matrix t -> Matrix t -> Matrix t mul a b = a <> b
src/Numeric/LinearAlgebra/Static.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE KindSignatures #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-}@@ -13,6 +14,8 @@ {-# LANGUAGE GADTs #-} {-# LANGUAGE TypeFamilies #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-}+{-# OPTIONS_GHC -fno-warn-orphans #-} {- | Module : Numeric.LinearAlgebra.Static@@ -42,19 +45,31 @@ blockAt, matrix, -- * Complex- C, M, Her, her, 𝑖,+ ℂ, C, M, Her, her, 𝑖,+ toComplex,+ fromComplex,+ complex,+ real,+ imag,+ sqMagnitude,+ magnitude, -- * Products (<>),(#>),(<.>), -- * Linear Systems linSolve, (<\>), -- * Factorizations svd, withCompactSVD, svdTall, svdFlat, Eigen(..),- withNullspace, qr, chol,+ withNullspace, withOrth, qr, chol,+ -- * Norms+ Normed(..),+ -- * Random arrays+ Seed, RandDist(..),+ randomVector, rand, randn, gaussianSample, uniformSample, -- * Misc- mean,+ mean, meanCov, Disp(..), Domain(..),- withVector, withMatrix,- toRows, toColumns,+ withVector, withMatrix, exactLength, exactDims,+ toRows, toColumns, withRows, withColumns, Sized(..), Diag(..), Sym, sym, mTm, unSym, (<·>) ) where @@ -65,18 +80,28 @@ row,col,vector,matrix,linspace,toRows,toColumns, (<\>),fromList,takeDiag,svd,eig,eigSH, eigenvalues,eigenvaluesSH,build,- qr,size,dot,chol,range,R,C,sym,mTm,unSym)+ qr,size,dot,chol,range,R,C,sym,mTm,unSym,+ randomVector,rand,randn,gaussianSample,uniformSample,meanCov,+ toComplex, fromComplex, complex, real, magnitude+ ) import qualified Numeric.LinearAlgebra as LA-import Data.Proxy(Proxy)+import qualified Numeric.LinearAlgebra.Devel as LA+import Data.Proxy(Proxy(..)) import Internal.Static import Control.Arrow((***))+import Text.Printf+import Data.Type.Equality ((:~:)(Refl))+import qualified Data.Bifunctor as BF (first)+#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif ud1 :: R n -> Vector ℝ ud1 (R (Dim v)) = v infixl 4 &-(&) :: forall n . (KnownNat n, 1 <= n)+(&) :: forall n . KnownNat n => R n -> ℝ -> R (n+1) u & x = u # (konst x :: R 1) @@ -204,11 +229,11 @@ infixr 8 <·>-(<·>) :: R n -> R n -> ℝ+(<·>) :: KnownNat n => R n -> R n -> ℝ (<·>) = dotR infixr 8 <.>-(<.>) :: R n -> R n -> ℝ+(<.>) :: KnownNat n => R n -> R n -> ℝ (<.>) = dotR --------------------------------------------------------------------------------@@ -218,21 +243,39 @@ takeDiag :: m -> d -instance forall n . (KnownNat n) => Diag (L n n) (R n)+instance KnownNat n => Diag (L n n) (R n) where- takeDiag m = mkR (LA.takeDiag (extract m))+ takeDiag x = mkR (LA.takeDiag (extract x)) -instance forall m n . (KnownNat m, KnownNat n, m <= n+1) => Diag (L m n) (R m)+instance KnownNat n => Diag (M n n) (C n) where- takeDiag m = mkR (LA.takeDiag (extract m))+ takeDiag x = mkC (LA.takeDiag (extract x)) +-------------------------------------------------------------------------------- -instance forall m n . (KnownNat m, KnownNat n, n <= m+1) => Diag (L m n) (R n)- where- takeDiag m = mkR (LA.takeDiag (extract m)) +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)@@ -289,7 +332,21 @@ her m = Her $ (m + LA.tr m)/2 +instance (KnownNat n) => Disp (Sym n)+ where+ disp n (Sym x) = do+ let a = extract x+ let su = LA.dispf n a+ printf "Sym %d" (cols a) >> putStr (dropWhile (/='\n') $ su) +instance (KnownNat n) => Disp (Her n)+ where+ disp n (Her x) = do+ let a = extract x+ let su = LA.dispcf n a+ printf "Her %d" (cols a) >> putStr (dropWhile (/='\n') $ su)++ instance KnownNat n => Eigen (Sym n) (R n) (L n n) where eigenvalues (Sym (extract -> m)) = mkR . LA.eigenvaluesSH . LA.trustSym $ m@@ -319,6 +376,15 @@ Nothing -> error "static/dynamic mismatch" Just (SomeNat (_ :: Proxy k)) -> f (mkL a :: L n k) +withOrth+ :: forall m n z . (KnownNat m, KnownNat n)+ => L m n+ -> (forall k. (KnownNat k) => L n k -> z)+ -> z+withOrth (LA.orth . extract -> a) f =+ case someNatVal $ fromIntegral $ cols a of+ Nothing -> error "static/dynamic mismatch"+ Just (SomeNat (_ :: Proxy k)) -> f (mkL a :: L n k) withCompactSVD :: forall m n z . (KnownNat m, KnownNat n)@@ -349,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)@@ -365,11 +431,31 @@ toRows :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R n] toRows (LA.toRows . extract -> vs) = map mkR vs +withRows+ :: forall n z . KnownNat n+ => [R n]+ -> (forall m . KnownNat m => L m n -> z)+ -> z+withRows (LA.fromRows . map extract -> m) f =+ case someNatVal $ fromIntegral $ LA.rows m of+ Nothing -> error "static/dynamic mismatch"+ Just (SomeNat (_ :: Proxy m)) -> f (mkL m :: L m n) toColumns :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R m] toColumns (LA.toColumns . extract -> vs) = map mkR vs +withColumns+ :: forall m z . KnownNat m+ => [R m]+ -> (forall n . KnownNat n => L m n -> z)+ -> z+withColumns (LA.fromColumns . map extract -> m) f =+ case someNatVal $ fromIntegral $ LA.cols m of+ Nothing -> error "static/dynamic mismatch"+ Just (SomeNat (_ :: Proxy n)) -> f (mkL m :: L m n) ++ -------------------------------------------------------------------------------- build@@ -392,6 +478,15 @@ Nothing -> error "static/dynamic mismatch" Just (SomeNat (_ :: Proxy m)) -> f (mkR v :: R m) +-- | Useful for constraining two dependently typed vectors to match each+-- other in length when they are unknown at compile-time.+exactLength+ :: forall n m . (KnownNat n, KnownNat m)+ => R m+ -> Maybe (R n)+exactLength v = do+ Refl <- sameNat (Proxy :: Proxy n) (Proxy :: Proxy m)+ return $ mkR (unwrap v) withMatrix :: forall z@@ -407,6 +502,64 @@ Just (SomeNat (_ :: Proxy n)) -> f (mkL a :: L m n) +-- | Useful for constraining two dependently typed matrices to match each+-- other in dimensions when they are unknown at compile-time.+exactDims+ :: forall n m j k . (KnownNat n, KnownNat m, KnownNat j, KnownNat k)+ => L m n+ -> Maybe (L j k)+exactDims m = do+ Refl <- sameNat (Proxy :: Proxy m) (Proxy :: Proxy j)+ Refl <- sameNat (Proxy :: Proxy n) (Proxy :: Proxy k)+ return $ mkL (unwrap m)++randomVector+ :: forall n . KnownNat n+ => Seed+ -> RandDist+ -> R n+randomVector s d = mkR (LA.randomVector s d+ (fromInteger (natVal (Proxy :: Proxy n)))+ )++rand+ :: forall m n . (KnownNat m, KnownNat n)+ => IO (L m n)+rand = mkL <$> LA.rand (fromInteger (natVal (Proxy :: Proxy m)))+ (fromInteger (natVal (Proxy :: Proxy n)))++randn+ :: forall m n . (KnownNat m, KnownNat n)+ => IO (L m n)+randn = mkL <$> LA.randn (fromInteger (natVal (Proxy :: Proxy m)))+ (fromInteger (natVal (Proxy :: Proxy n)))++gaussianSample+ :: forall m n . (KnownNat m, KnownNat n)+ => Seed+ -> R n+ -> Sym n+ -> L m n+gaussianSample s (extract -> mu) (Sym (extract -> sigma)) =+ mkL $ LA.gaussianSample s (fromInteger (natVal (Proxy :: Proxy m)))+ mu (LA.trustSym sigma)++uniformSample+ :: forall m n . (KnownNat m, KnownNat n)+ => Seed+ -> R n -- ^ minimums of each row+ -> R n -- ^ maximums of each row+ -> L m n+uniformSample s (extract -> mins) (extract -> maxs) =+ mkL $ LA.uniformSample s (fromInteger (natVal (Proxy :: Proxy m)))+ (zip (LA.toList mins) (LA.toList maxs))++meanCov+ :: forall m n . (KnownNat m, KnownNat n, 1 <= m)+ => L m n+ -> (R n, Sym n)+meanCov (extract -> vs) = mkR *** (Sym . mkL . LA.unSym) $ LA.meanCov vs+ -------------------------------------------------------------------------------- class Domain field vec mat | mat -> vec field, vec -> mat field, field -> mat vec@@ -416,6 +569,15 @@ dot :: forall n . (KnownNat n) => vec n -> vec n -> field cross :: vec 3 -> vec 3 -> vec 3 diagR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => field -> vec k -> mat m n+ dvmap :: forall n. KnownNat n => (field -> field) -> vec n -> vec n+ dmmap :: forall n m. (KnownNat m, KnownNat n) => (field -> field) -> mat n m -> mat n m+ outer :: forall n m. (KnownNat m, KnownNat n) => vec n -> vec m -> mat n m+ zipWithVector :: forall n. KnownNat n => (field -> field -> field) -> vec n -> vec n -> vec n+ det :: forall n. KnownNat n => mat n n -> field+ invlndet :: forall n. KnownNat n => mat n n -> (mat n n, (field, field))+ expm :: forall n. KnownNat n => mat n n -> mat n n+ sqrtm :: forall n. KnownNat n => mat n n -> mat n n+ inv :: forall n. KnownNat n => mat n n -> mat n n instance Domain ℝ R L@@ -425,6 +587,15 @@ dot = dotR cross = crossR diagR = diagRectR+ dvmap = mapR+ dmmap = mapL+ outer = outerR+ zipWithVector = zipWithR+ det = detL+ invlndet = invlndetL+ expm = expmL+ sqrtm = sqrtmL+ inv = invL instance Domain ℂ C M where@@ -433,6 +604,15 @@ dot = dotC cross = crossC diagR = diagRectC+ dvmap = mapC+ dmmap = mapM'+ outer = outerC+ zipWithVector = zipWithC+ det = detM+ invlndet = invlndetM+ expm = expmM+ sqrtm = sqrtmM+ inv = invM -------------------------------------------------------------------------------- @@ -447,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) @@ -459,10 +639,8 @@ appR m v = mkR (extract m LA.#> extract v) -dotR :: R n -> R n -> ℝ-dotR (ud1 -> u) (ud1 -> v)- | singleV u || singleV v = sumElements (u * v)- | otherwise = udot u v+dotR :: KnownNat n => R n -> R n -> ℝ+dotR (extract -> u) (extract -> v) = LA.dot u v crossR :: R 3 -> R 3 -> R 3@@ -472,6 +650,33 @@ z2 = x!2*y!0-x!0*y!2 z3 = x!0*y!1-x!1*y!0 +outerR :: (KnownNat m, KnownNat n) => R n -> R m -> L n m+outerR (extract -> x) (extract -> y) = mkL (LA.outer x y)++mapR :: KnownNat n => (ℝ -> ℝ) -> R n -> R n+mapR f (unwrap -> v) = mkR (LA.cmap f v)++zipWithR :: KnownNat n => (ℝ -> ℝ -> ℝ) -> R n -> R n -> R n+zipWithR f (extract -> x) (extract -> y) = mkR (LA.zipVectorWith f x y)++mapL :: (KnownNat n, KnownNat m) => (ℝ -> ℝ) -> L n m -> L n m+mapL f = overMatL' (LA.cmap f)++detL :: KnownNat n => Sq n -> ℝ+detL = LA.det . unwrap++invlndetL :: KnownNat n => Sq n -> (L n n, (ℝ, ℝ))+invlndetL = BF.first mkL . LA.invlndet . unwrap++expmL :: KnownNat n => Sq n -> Sq n+expmL = overMatL' LA.expm++sqrtmL :: KnownNat n => Sq n -> Sq n+sqrtmL = overMatL' LA.sqrtm++invL :: KnownNat n => Sq n -> Sq n+invL = overMatL' LA.inv+ -------------------------------------------------------------------------------- mulC :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => M m k -> M k n -> M m n@@ -485,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) @@ -498,9 +703,7 @@ dotC :: KnownNat n => C n -> C n -> ℂ-dotC (unwrap -> u) (unwrap -> v)- | singleV u || singleV v = sumElements (conj u * v)- | otherwise = u LA.<.> v+dotC (extract -> u) (extract -> v) = LA.dot u v crossC :: C 3 -> C 3 -> C 3@@ -510,6 +713,33 @@ z2 = x!2*y!0-x!0*y!2 z3 = x!0*y!1-x!1*y!0 +outerC :: (KnownNat m, KnownNat n) => C n -> C m -> M n m+outerC (extract -> x) (extract -> y) = mkM (LA.outer x y)++mapC :: KnownNat n => (ℂ -> ℂ) -> C n -> C n+mapC f (unwrap -> v) = mkC (LA.cmap f v)++zipWithC :: KnownNat n => (ℂ -> ℂ -> ℂ) -> C n -> C n -> C n+zipWithC f (extract -> x) (extract -> y) = mkC (LA.zipVectorWith f x y)++mapM' :: (KnownNat n, KnownNat m) => (ℂ -> ℂ) -> M n m -> M n m+mapM' f = overMatM' (LA.cmap f)++detM :: KnownNat n => M n n -> ℂ+detM = LA.det . unwrap++invlndetM :: KnownNat n => M n n -> (M n n, (ℂ, ℂ))+invlndetM = BF.first mkM . LA.invlndet . unwrap++expmM :: KnownNat n => M n n -> M n n+expmM = overMatM' LA.expm++sqrtmM :: KnownNat n => M n n -> M n n+sqrtmM = overMatM' LA.sqrtm++invM :: KnownNat n => M n n -> M n n+invM = overMatM' LA.inv+ -------------------------------------------------------------------------------- diagRectR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℝ -> R k -> L m n@@ -616,3 +846,67 @@ where checkT _ = test +--------------------------------------------------------------------------------++instance KnownNat n => Normed (R n)+ where+ norm_0 v = norm_0 (extract v)+ norm_1 v = norm_1 (extract v)+ norm_2 v = norm_2 (extract v)+ norm_Inf v = norm_Inf (extract v)++instance (KnownNat m, KnownNat n) => Normed (L m n)+ where+ norm_0 m = norm_0 (extract m)+ norm_1 m = norm_1 (extract m)+ norm_2 m = norm_2 (extract m)+ norm_Inf m = norm_Inf (extract m)++mkSym f = Sym . f . unSym+mkSym2 f x y = Sym (f (unSym x) (unSym y))++instance KnownNat n => Num (Sym n)+ where+ (+) = mkSym2 (+)+ (*) = mkSym2 (*)+ (-) = mkSym2 (-)+ abs = mkSym abs+ signum = mkSym signum+ negate = mkSym negate+ fromInteger = Sym . fromInteger++instance KnownNat n => Fractional (Sym n)+ where+ fromRational = Sym . fromRational+ (/) = mkSym2 (/)++instance KnownNat n => Floating (Sym n)+ where+ sin = mkSym sin+ cos = mkSym cos+ tan = mkSym tan+ asin = mkSym asin+ acos = mkSym acos+ atan = mkSym atan+ sinh = mkSym sinh+ cosh = mkSym cosh+ tanh = mkSym tanh+ asinh = mkSym asinh+ acosh = mkSym acosh+ atanh = mkSym atanh+ exp = mkSym exp+ log = mkSym log+ sqrt = mkSym sqrt+ (**) = mkSym2 (**)+ pi = Sym pi++instance KnownNat n => Additive (Sym n) where+ add = (+)++instance KnownNat n => Transposable (Sym n) (Sym n) where+ tr = id+ tr' = id++instance KnownNat n => Transposable (Her n) (Her n) where+ tr = id+ tr' (Her m) = Her (tr' m)
src/Numeric/Matrix.hs view
@@ -4,6 +4,8 @@ {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} +{-# OPTIONS_GHC -fno-warn-orphans #-}+ ----------------------------------------------------------------------------- -- | -- Module : Numeric.Matrix@@ -32,8 +34,12 @@ import Internal.Numeric import qualified Data.Monoid as M import Data.List(partition)+import qualified Data.Foldable as F+import qualified Data.Semigroup as S import Internal.Chain+import Foreign.Storable(Storable) + ------------------------------------------------------------------- instance Container Matrix a => Eq (Matrix a) where@@ -77,18 +83,31 @@ -------------------------------------------------------------------------------- +isScalar :: Matrix t -> Bool isScalar m = rows m == 1 && cols m == 1 +adaptScalarM :: (Foreign.Storable.Storable t1, Foreign.Storable.Storable t2)+ => (t1 -> Matrix t2 -> t)+ -> (Matrix t1 -> Matrix t2 -> t)+ -> (Matrix t1 -> t2 -> t)+ -> Matrix t1+ -> Matrix t2+ -> t adaptScalarM f1 f2 f3 x y | isScalar x = f1 (x @@>(0,0) ) y | isScalar y = f3 x (y @@>(0,0) ) | otherwise = f2 x y +instance (Container Vector t, Eq t, Num (Vector t), Product t) => S.Semigroup (Matrix t)+ where+ (<>) = mappend+ sconcat = mconcat . F.toList+ instance (Container Vector t, Eq t, Num (Vector t), Product t) => M.Monoid (Matrix t) where mempty = 1 mappend = adaptScalarM scale mXm (flip scale)- + mconcat xs = work (partition isScalar xs) where work (ss,[]) = product ss@@ -98,4 +117,3 @@ | otherwise = scale x00 m where x00 = x @@> (0,0)-
src/Numeric/Vector.hs view
@@ -3,6 +3,9 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE MultiParamTypeClasses #-}++{-# OPTIONS_GHC -fno-warn-orphans #-}+ ----------------------------------------------------------------------------- -- | -- Module : Numeric.Vector@@ -14,7 +17,7 @@ -- -- Provides instances of standard classes 'Show', 'Read', 'Eq', -- 'Num', 'Fractional', and 'Floating' for 'Vector'.--- +-- ----------------------------------------------------------------------------- module Numeric.Vector () where@@ -23,9 +26,17 @@ import Internal.Vector import Internal.Numeric import Internal.Conversion+import Foreign.Storable(Storable) ------------------------------------------------------------------- +adaptScalar :: (Foreign.Storable.Storable t1, Foreign.Storable.Storable t2)+ => (t1 -> Vector t2 -> t)+ -> (Vector t1 -> Vector t2 -> t)+ -> (Vector t1 -> t2 -> t)+ -> Vector t1+ -> Vector t2+ -> t adaptScalar f1 f2 f3 x y | dim x == 1 = f1 (x@>0) y | dim y == 1 = f3 x (y@>0)@@ -172,4 +183,3 @@ sqrt = vectorMapQ Sqrt (**) = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS)) pi = fromList [pi]-