diff --git a/THANKS.md b/THANKS.md --- a/THANKS.md +++ b/THANKS.md @@ -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. diff --git a/hmatrix.cabal b/hmatrix.cabal --- a/hmatrix.cabal +++ b/hmatrix.cabal @@ -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 diff --git a/src/Internal/Algorithms.hs b/src/Internal/Algorithms.hs --- a/src/Internal/Algorithms.hs +++ b/src/Internal/Algorithms.hs @@ -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 diff --git a/src/Internal/C/lapack-aux.c b/src/Internal/C/lapack-aux.c --- a/src/Internal/C/lapack-aux.c +++ b/src/Internal/C/lapack-aux.c @@ -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 * diff --git a/src/Internal/C/vector-aux.c b/src/Internal/C/vector-aux.c --- a/src/Internal/C/vector-aux.c +++ b/src/Internal/C/vector-aux.c @@ -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 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 diff --git a/src/Internal/Matrix.hs b/src/Internal/Matrix.hs --- a/src/Internal/Matrix.hs +++ b/src/Internal/Matrix.hs @@ -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 -------------------------------------------------------------------------------- diff --git a/src/Internal/Numeric.hs b/src/Internal/Numeric.hs --- a/src/Internal/Numeric.hs +++ b/src/Internal/Numeric.hs @@ -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 diff --git a/src/Internal/Util.hs b/src/Internal/Util.hs --- a/src/Internal/Util.hs +++ b/src/Internal/Util.hs @@ -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 diff --git a/src/Numeric/LinearAlgebra.hs b/src/Numeric/LinearAlgebra.hs --- a/src/Numeric/LinearAlgebra.hs +++ b/src/Numeric/LinearAlgebra.hs @@ -94,6 +94,11 @@ ldlSolve, ldlPacked, -- ** Positive definite cholSolve, + -- ** Triangular + UpLo(..), + triSolve, + -- ** Tridiagonal + triDiagSolve, -- ** Sparse cgSolve, cgSolve', diff --git a/src/Numeric/LinearAlgebra/Devel.hs b/src/Numeric/LinearAlgebra/Devel.hs --- a/src/Numeric/LinearAlgebra/Devel.hs +++ b/src/Numeric/LinearAlgebra/Devel.hs @@ -55,7 +55,7 @@ GMatrix(..), -- * Misc - toByteString, fromByteString, showInternal + toByteString, fromByteString, showInternal, reorderVector ) where