10 files changed, 446 insertions, 25 deletions

diff --git a/packages/base/THANKS.md b/packages/base/THANKS.md
index 7e82d40..45ded81 100644
--- a/packages/base/THANKS.md
+++ b/packages/base/THANKS.md
@@ -160,7 +160,7 @@ module reorganization, monadic mapVectorM, and many other improvements.
 - 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 @@ module reorganization, monadic mapVectorM, and many other improvements.
 
 - 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 @@ module reorganization, monadic mapVectorM, and many other improvements.
 - 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 @@ fixed the CPP issue in cabal files, and made many other contributions.
   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/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index 4d01e57..968a00e 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/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
@@ -130,7 +130,7 @@ library
 
     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/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs
index 70d65d7..99c9e34 100644
--- a/packages/base/src/Internal/Algorithms.hs
+++ b/packages/base/src/Internal/Algorithms.hs
@@ -20,13 +20,16 @@ imported from "Numeric.LinearAlgebra.LAPACK".
 -}
 -----------------------------------------------------------------------------
 
-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 @@ class (Numeric t,
     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 @@ instance Field Double where
     linearSolve' = linearSolveR                 -- (luSolve . luPacked) ??
     mbLinearSolve' = mbLinearSolveR
     cholSolve' = cholSolveR
+    triSolve' = triSolveR
+    triDiagSolve' = triDiagSolveR
     linearSolveLS' = linearSolveLSR
     linearSolveSVD' = linearSolveSVDR Nothing
     eig' = eigR
@@ -112,6 +119,8 @@ instance Field (Complex Double) where
     linearSolve' = linearSolveC
     mbLinearSolve' = mbLinearSolveC
     cholSolve' = cholSolveC
+    triSolve' = triSolveC
+    triDiagSolve' = triDiagSolveC
     linearSolveLS' = linearSolveLSC
     linearSolveSVD' = linearSolveSVDC Nothing
     eig' = eigC
@@ -350,6 +359,79 @@ cholSolve
     -> 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'
diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index ff7ad92..5018a45 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -584,6 +584,162 @@ int cholSolveC_l(KOCMAT(a),OCMAT(b)) {
     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/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index 231109a..e306454 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -406,6 +406,58 @@ cholSolveR a b = linearSolveSQAux2 id dpotrs "cholSolveR" (fmat a) b
 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
diff --git a/packages/base/src/Internal/Numeric.hs b/packages/base/src/Internal/Numeric.hs
index dd74165..c9ef0c5 100644
--- a/packages/base/src/Internal/Numeric.hs
+++ b/packages/base/src/Internal/Numeric.hs
@@ -25,6 +25,7 @@ import Internal.Conversion
 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 @@ instance Container Vector I
     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 @@ instance Container Vector Z
     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 @@ instance Container Vector Float
     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 @@ instance Container Vector Double
     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 @@ instance Container Vector (Complex Double)
     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 @@ instance Container Vector (Complex Float)
     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/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index 17d3e50..ec21fe4 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -2,6 +2,7 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE ViewPatterns #-}
 
 
@@ -613,6 +614,9 @@ gaussElim_1 x y = dropColumns (rows x) (flipud $ fromRows s2)
     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 @@ pivotDown t n xs
     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 @@ pivotDown t n xs
     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/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index badf8f9..fd100e0 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -94,6 +94,11 @@ module Numeric.LinearAlgebra (
     ldlSolve, ldlPacked,
     -- ** Positive definite
     cholSolve,
+    -- ** Triangular
+    UpLo(..),
+    triSolve,
+    -- ** Tridiagonal
+    triDiagSolve,
     -- ** Sparse
     cgSolve,
     cgSolve',
diff --git a/packages/gsl/src/Numeric/GSL/gsl-ode.c b/packages/gsl/src/Numeric/GSL/gsl-ode.c
index a6bdb55..72c8617 100644
--- a/packages/gsl/src/Numeric/GSL/gsl-ode.c
+++ b/packages/gsl/src/Numeric/GSL/gsl-ode.c
@@ -178,11 +178,16 @@ int ode(int method, int control, double h,
            status = gsl_odeiv2_driver_apply (d, &t, ti, y);
      
            if (status != GSL_SUCCESS) {
-                  printf ("error in ode, return value=%d\n", status);
-                  break;
-                }
-
-//           printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
+             int k;
+             printf ("error in ode, return value=%d\n", status);
+             printf("last successful values are:\n");
+             printf("t = %.5e\n", t);
+             for (k=0; k < xin; k++)
+               {
+                 printf("y[%d] = %.5e\n", k, y[k]);
+               }
+             break;
+           }
            
            for(j=0; j<xin; j++) {
                solp[i*xin + j] = y[j];
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 043ebf3..6d54f4d 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -131,6 +131,111 @@ mbCholTest = utest "mbCholTest" (ok1 && ok2) where
     ok1 = mbChol (trustSym m1) == Nothing
     ok2 = mbChol (trustSym m2) == Just (chol $ trustSym m2)
 
+-----------------------------------------------------
+
+triTest = utest "triTest" ok1 where
+
+  a :: Matrix R
+  a = (4><4)
+    [
+       4.30,  0.00,  0.00, 0.00,
+      -3.96, -4.87,  0.00, 0.00,
+       0.40,  0.31, -8.02, 0.00,
+      -0.27,  0.07, -5.95, 0.12
+    ]
+
+  w :: Matrix R
+  w = (4><2)
+    [
+      -12.90, -21.50,
+       16.75,  14.93,
+      -17.55,   6.33,
+      -11.04,   8.09
+    ]
+
+  v :: Matrix R
+  v = triSolve Lower a w
+
+  e :: Matrix R
+  e = (4><2)
+    [
+      -3.0000, -5.0000,
+      -1.0000,  1.0000,
+       2.0000, -1.0000,
+       1.0000,  6.0000
+    ]
+
+  ok1 = (norm_Inf . flatten $ e - v) <= 1e-13
+
+-----------------------------------------------------
+
+triDiagTest = utest "triDiagTest" (ok1 && ok2) where
+
+  dL, d, dU :: Vector Double
+  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 :: Matrix R
+  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
+    ]
+
+  y :: Matrix R
+  y = (9><9)
+    [
+      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
+    ]
+
+  x :: Matrix R
+  x = triDiagSolve dL d dU b
+
+  z :: Matrix C
+  z = (4><4)
+    [
+      1.0 :+ 1.0, 4.0 :+ 4.0, 0.0 :+ 0.0, 0.0 :+ 0.0,
+      3.0 :+ 3.0, 1.0 :+ 1.0, 4.0 :+ 4.0, 0.0 :+ 0.0,
+      0.0 :+ 0.0, 3.0 :+ 3.0, 1.0 :+ 1.0, 4.0 :+ 4.0,
+      0.0 :+ 0.0, 0.0 :+ 0.0, 3.0 :+ 3.0, 1.0 :+ 1.0
+    ]
+
+  zDL, zD, zDu :: Vector C
+  zDL = fromList [3.0 :+ 3.0, 3.0 :+ 3.0, 3.0 :+ 3.0]
+  zD  = fromList [1.0 :+ 1.0, 1.0 :+ 1.0, 1.0 :+ 1.0, 1.0 :+ 1.0]
+  zDu = fromList [4.0 :+ 4.0, 4.0 :+ 4.0, 4.0 :+ 4.0]
+
+  zB :: Matrix C
+  zB = (4><3)
+    [
+      1.0 :+ 1.0,   1.0  :+   1.0,  1.0 :+ (-1.0),
+      1.0 :+ 1.0, (-1.0) :+ (-1.0), 1.0 :+ (-1.0),
+      1.0 :+ 1.0,   1.0  :+   1.0,  1.0 :+ (-1.0),
+      1.0 :+ 1.0, (-1.0) :+ (-1.0), 1.0 :+ (-1.0)
+    ]
+
+  u :: Matrix C
+  u = triDiagSolve zDL zD zDu zB
+
+  ok1 = (maximum $ map abs $ concat $ toLists $ b - (y <> x)) <= 1e-15
+  ok2 = (maximum $ map magnitude $ concat $ toLists $ zB - (z <> u)) <= 1e-15
+
 ---------------------------------------------------------------------
 
 randomTestGaussian = (unSym c) :~3~: unSym (snd (meanCov dat))
@@ -715,6 +820,8 @@ runTests n = do
                        && rank ((2><3)[1,0,0,1,7*peps,0::Double]) == 2
         , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)
         , mbCholTest
+        , triTest
+        , triDiagTest
         , utest "offset" offsetTest
         , normsVTest
         , normsMTest