hmatrix 0.18.0.0 → 0.18.1.0
raw patch · 11 files changed
+416/−22 lines, 11 filesdep ~basePVP ok
version bump matches the API change (PVP)
Dependency ranges changed: base
API changes (from Hackage documentation)
+ Numeric.LinearAlgebra: Lower :: UpLo
+ Numeric.LinearAlgebra: Upper :: UpLo
+ Numeric.LinearAlgebra: data UpLo
+ Numeric.LinearAlgebra: triDiagSolve :: Field t => Vector t -> Vector t -> Vector t -> Matrix t -> Matrix t
+ Numeric.LinearAlgebra: triSolve :: Field t => UpLo -> Matrix t -> Matrix t -> Matrix t
+ Numeric.LinearAlgebra.Devel: reorderVector :: Element a => Vector CInt -> Vector CInt -> Vector a -> Vector a
Files
- THANKS.md +8/−4
- hmatrix.cabal +3/−3
- src/Internal/Algorithms.hs +84/−2
- src/Internal/C/lapack-aux.c +156/−0
- src/Internal/C/vector-aux.c +51/−0
- src/Internal/LAPACK.hs +52/−0
- src/Internal/Matrix.hs +34/−0
- src/Internal/Numeric.hs +13/−12
- src/Internal/Util.hs +9/−0
- src/Numeric/LinearAlgebra.hs +5/−0
- src/Numeric/LinearAlgebra/Devel.hs +1/−1
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. - 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. - "erdeszt" optimized "conv" using a direct vector reverse. @@ -203,7 +204,8 @@ - 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 @@ -224,7 +226,9 @@ Andras Slemmer fixed the bug. Thank you all. - Kevin Slagle implemented thinQR and thinRQ, much faster than the original qr,- and added compactSVDTol.+ 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.
hmatrix.cabal view
@@ -1,5 +1,5 @@ Name: hmatrix-Version: 0.18.0.0+Version: 0.18.1.0 License: BSD3 License-file: LICENSE Author: Alberto Ruiz@@ -83,7 +83,7 @@ ghc-options: -Wall -fno-warn-missing-signatures -fno-warn-orphans- -fprof-auto+ -fno-prof-auto cc-options: -O4 -Wall @@ -121,7 +121,7 @@ if os(windows) if flag(openblas)- extra-libraries: libopenblas, libgcc_s_seh-1, libgfortran-3, libquadmath-0+ extra-libraries: libopenblas, libgcc_s_seh-1, libgfortran, libquadmath-0 else extra-libraries: blas lapack
src/Internal/Algorithms.hs view
@@ -20,13 +20,16 @@ -} ----------------------------------------------------------------------------- -module Internal.Algorithms where+module Internal.Algorithms (+ module Internal.Algorithms,+ UpLo(..)+) where 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@@ -58,6 +61,8 @@ 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@@ -83,6 +88,8 @@ linearSolve' = linearSolveR -- (luSolve . luPacked) ?? mbLinearSolve' = mbLinearSolveR cholSolve' = cholSolveR+ triSolve' = triSolveR+ triDiagSolve' = triDiagSolveR linearSolveLS' = linearSolveLSR linearSolveSVD' = linearSolveSVDR Nothing eig' = eigR@@ -112,6 +119,8 @@ linearSolve' = linearSolveC mbLinearSolve' = mbLinearSolveC cholSolve' = cholSolveC+ triSolve' = triSolveC+ triDiagSolve' = triDiagSolveC linearSolveLS' = linearSolveLSC linearSolveSVD' = linearSolveSVDC Nothing eig' = eigC@@ -349,6 +358,79 @@ -> Matrix t -- ^ right hand sides -> 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
src/Internal/C/lapack-aux.c view
@@ -584,6 +584,162 @@ 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 n = dn;+ integer nhrs = bc;+ REQUIRES(n >= 1 && dln == dn - 1 && dun == dn - 1 && br == n, BAD_SIZE);+ DEBUGMSG("triDiagSolveR_l");+ integer res;+ integer* ipiv = (integer*)malloc(n*sizeof(integer));+ double* du2 = (double*)malloc((n - 2)*sizeof(double));+ dgttrf_ (&n,+ dlp, dp, dup, du2,+ ipiv,+ &res);+ CHECK(res,res);+ dgttrs_ ("N",+ &n,&nhrs,+ dlp, dp, dup, du2,+ ipiv, bp, &n,+ &res);+ CHECK(res,res);+ free(ipiv);+ free(du2);+ OK+}++//////// 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 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));+ doublecomplex* du2 = (doublecomplex*)malloc((n - 2)*sizeof(doublecomplex));+ zgttrf_ (&n,+ dlp, dp, dup, du2,+ ipiv,+ &res);+ CHECK(res,res);+ zgttrs_ ("N",+ &n,&nhrs,+ dlp, dp, dup, du2,+ ipiv, bp, &n,+ &res);+ CHECK(res,res);+ free(ipiv);+ free(du2);+ OK+}+ //////////////////// least squares real linear system //////////// int dgels_(char *trans, integer *m, integer *n, integer *
src/Internal/C/vector-aux.c view
@@ -1533,3 +1533,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/LAPACK.hs view
@@ -406,6 +406,58 @@ 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+ 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
src/Internal/Matrix.hs view
@@ -285,6 +285,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@@ -298,6 +299,7 @@ remapM = remapF rowOp = rowOpAux c_rowOpF gemm = gemmg c_gemmF+ reorderV = reorderAux c_reorderF instance Element Double where constantD = constantAux cconstantR@@ -310,6 +312,7 @@ remapM = remapD rowOp = rowOpAux c_rowOpD gemm = gemmg c_gemmD+ reorderV = reorderAux c_reorderD instance Element (Complex Float) where constantD = constantAux cconstantQ@@ -322,6 +325,7 @@ remapM = remapQ rowOp = rowOpAux c_rowOpQ gemm = gemmg c_gemmQ+ reorderV = reorderAux c_reorderQ instance Element (Complex Double) where constantD = constantAux cconstantC@@ -334,6 +338,7 @@ remapM = remapC rowOp = rowOpAux c_rowOpC gemm = gemmg c_gemmC+ reorderV = reorderAux c_reorderC instance Element (CInt) where constantD = constantAux cconstantI@@ -346,6 +351,7 @@ remapM = remapI rowOp = rowOpAux c_rowOpI gemm = gemmg c_gemmI+ reorderV = reorderAux c_reorderI instance Element Z where constantD = constantAux cconstantL@@ -358,6 +364,7 @@ remapM = remapL rowOp = rowOpAux c_rowOpL gemm = gemmg c_gemmL+ reorderV = reorderAux c_reorderL ------------------------------------------------------------------- @@ -577,6 +584,33 @@ foreign import ccall unsafe "gemm_int64_t" c_gemmL :: Tgemm Z foreign import ccall unsafe "gemm_mod_int32_t" c_gemmMI :: I -> Tgemm I foreign import ccall unsafe "gemm_mod_int64_t" c_gemmML :: Z -> Tgemm Z++--------------------------------------------------------------------------------++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 --------------------------------------------------------------------------------
src/Internal/Numeric.hs view
@@ -25,6 +25,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 -------------------------------------------------------------------------------- @@ -102,8 +103,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@@ -141,8 +142,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@@ -181,8 +182,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@@ -218,8 +219,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@@ -255,8 +256,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@@ -291,8 +292,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
src/Internal/Util.hs view
@@ -2,6 +2,7 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE ViewPatterns #-} @@ -613,6 +614,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@@ -622,6 +626,7 @@ pivot k = (const k &&& id) . 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@@ -632,12 +637,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
src/Numeric/LinearAlgebra.hs view
@@ -94,6 +94,11 @@ ldlSolve, ldlPacked, -- ** Positive definite cholSolve,+ -- ** Triangular+ UpLo(..),+ triSolve,+ -- ** Tridiagonal+ triDiagSolve, -- ** Sparse cgSolve, cgSolve',
src/Numeric/LinearAlgebra/Devel.hs view
@@ -55,7 +55,7 @@ GMatrix(..), -- * Misc- toByteString, fromByteString, showInternal+ toByteString, fromByteString, showInternal, reorderVector ) where