    {-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
    {-# LANGUAGE CPP #-}
    {-# LANGUAGE MultiParamTypeClasses #-}
    {-# LANGUAGE TypeFamilies #-}

    {-# OPTIONS_GHC -fno-warn-missing-signatures #-}

    -----------------------------------------------------------------------------
    {- |
    Module      :  Internal.Algorithms
    Copyright   :  (c) Alberto Ruiz 2006-14
    License     :  BSD3
    Maintainer  :  Alberto Ruiz
    Stability   :  provisional

    High level generic interface to common matrix computations.

    Specific functions for particular base types can also be explicitly
    imported from "Numeric.LinearAlgebra.LAPACK".

    -}
    -----------------------------------------------------------------------------

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

    {- | Generic linear algebra functions for double precision real and complex matrices.

    (Single precision data can be converted using 'single' and 'double').

    -}
    class (Numeric t,
           Convert t,
           Normed Matrix t,
           Normed Vector t,
           Floating t,
           Linear t Vector,
           Linear t Matrix,
           Additive (Vector t),
           Additive (Matrix t),
           RealOf t ~ Double) => Field t where
        svd'         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
        thinSVD'     :: Matrix t -> (Matrix t, Vector Double, Matrix t)
        sv'          :: Matrix t -> Vector Double
        luPacked'    :: Matrix t -> (Matrix t, [Int])
        luSolve'     :: (Matrix t, [Int]) -> Matrix t -> Matrix 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
        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)
        qr'          :: Matrix t -> (Matrix t, Vector t)
        qrgr'        :: Int -> (Matrix t, Vector t) -> Matrix t
        hess'        :: Matrix t -> (Matrix t, Matrix t)
        schur'       :: Matrix t -> (Matrix t, Matrix t)

    instance Field Double where
        svd' :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
    svd' = Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
    svdRd
        thinSVD' :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
    thinSVD' = Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
    thinSVDRd
        sv' :: Matrix Double -> Vector Double
    sv' = Matrix Double -> Vector Double
    svR
        luPacked' :: Matrix Double -> (Matrix Double, [Int])
    luPacked' = Matrix Double -> (Matrix Double, [Int])
    luR
        luSolve' :: (Matrix Double, [Int]) -> Matrix Double -> Matrix Double
    luSolve' (Matrix Double
    l_u,[Int]
    perm) = Matrix Double -> [Int] -> Matrix Double -> Matrix Double
    lusR Matrix Double
    l_u [Int]
    perm
        linearSolve' :: Matrix Double -> Matrix Double -> Matrix Double
    linearSolve' = Matrix Double -> Matrix Double -> Matrix Double
    linearSolveR                 -- (luSolve . luPacked) ??
        mbLinearSolve' :: Matrix Double -> Matrix Double -> Maybe (Matrix Double)
    mbLinearSolve' = Matrix Double -> Matrix Double -> Maybe (Matrix Double)
    mbLinearSolveR
        cholSolve' :: Matrix Double -> Matrix Double -> Matrix Double
    cholSolve' = Matrix Double -> Matrix Double -> Matrix Double
    cholSolveR
        triSolve' :: UpLo -> Matrix Double -> Matrix Double -> Matrix Double
    triSolve' = UpLo -> Matrix Double -> Matrix Double -> Matrix Double
    triSolveR
        triDiagSolve' :: Vector Double
    -> Vector Double -> Vector Double -> Matrix Double -> Matrix Double
    triDiagSolve' = Vector Double
    -> Vector Double -> Vector Double -> Matrix Double -> Matrix Double
    triDiagSolveR
        linearSolveLS' :: Matrix Double -> Matrix Double -> Matrix Double
    linearSolveLS' = Matrix Double -> Matrix Double -> Matrix Double
    linearSolveLSR
        linearSolveSVD' :: Matrix Double -> Matrix Double -> Matrix Double
    linearSolveSVD' = Maybe Double -> Matrix Double -> Matrix Double -> Matrix Double
    linearSolveSVDR Maybe Double
    forall a. Maybe a
    Nothing
        eig' :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
    eig' = Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
    eigR
        eigSH'' :: Matrix Double -> (Vector Double, Matrix Double)
    eigSH'' = Matrix Double -> (Vector Double, Matrix Double)
    eigS
        geig' :: Matrix Double
    -> Matrix Double
    -> (Vector (Complex Double), Vector Double,
        Matrix (Complex Double))
    geig' = Matrix Double
    -> Matrix Double
    -> (Vector (Complex Double), Vector Double,
        Matrix (Complex Double))
    eigG
        eigOnly :: Matrix Double -> Vector (Complex Double)
    eigOnly = Matrix Double -> Vector (Complex Double)
    eigOnlyR
        geigOnly :: Matrix Double
    -> Matrix Double -> (Vector (Complex Double), Vector Double)
    geigOnly = Matrix Double
    -> Matrix Double -> (Vector (Complex Double), Vector Double)
    eigOnlyG
        eigOnlySH :: Matrix Double -> Vector Double
    eigOnlySH = Matrix Double -> Vector Double
    eigOnlyS
        cholSH' :: Matrix Double -> Matrix Double
    cholSH' = Matrix Double -> Matrix Double
    cholS
        mbCholSH' :: Matrix Double -> Maybe (Matrix Double)
    mbCholSH' = Matrix Double -> Maybe (Matrix Double)
    mbCholS
        qr' :: Matrix Double -> (Matrix Double, Vector Double)
    qr' = Matrix Double -> (Matrix Double, Vector Double)
    qrR
        qrgr' :: Int -> (Matrix Double, Vector Double) -> Matrix Double
    qrgr' = Int -> (Matrix Double, Vector Double) -> Matrix Double
    qrgrR
        hess' :: Matrix Double -> (Matrix Double, Matrix Double)
    hess' = (Matrix Double -> (Matrix Double, Vector Double))
    -> Matrix Double -> (Matrix Double, Matrix Double)
    forall t.
    Field t =>
    (Matrix t -> (Matrix t, Vector t))
    -> Matrix t -> (Matrix t, Matrix t)
    unpackHess Matrix Double -> (Matrix Double, Vector Double)
    hessR
        schur' :: Matrix Double -> (Matrix Double, Matrix Double)
    schur' = Matrix Double -> (Matrix Double, Matrix Double)
    schurR
        ldlPacked' :: Matrix Double -> (Matrix Double, [Int])
    ldlPacked' = Matrix Double -> (Matrix Double, [Int])
    ldlR
        ldlSolve' :: (Matrix Double, [Int]) -> Matrix Double -> Matrix Double
    ldlSolve'= (Matrix Double -> [Int] -> Matrix Double -> Matrix Double)
    -> (Matrix Double, [Int]) -> Matrix Double -> Matrix Double
    forall a b c. (a -> b -> c) -> (a, b) -> c
    uncurry Matrix Double -> [Int] -> Matrix Double -> Matrix Double
    ldlsR

    instance Field (Complex Double) where
    #ifdef NOZGESDD
        svd' = svdC
        thinSVD' = thinSVDC
    #else
        svd' :: Matrix (Complex Double)
    -> (Matrix (Complex Double), Vector Double,
        Matrix (Complex Double))
    svd' = Matrix (Complex Double)
    -> (Matrix (Complex Double), Vector Double,
        Matrix (Complex Double))
    svdCd
        thinSVD' :: Matrix (Complex Double)
    -> (Matrix (Complex Double), Vector Double,
        Matrix (Complex Double))
    thinSVD' = Matrix (Complex Double)
    -> (Matrix (Complex Double), Vector Double,
        Matrix (Complex Double))
    thinSVDCd
    #endif
        sv' :: Matrix (Complex Double) -> Vector Double
    sv' = Matrix (Complex Double) -> Vector Double
    svC
        luPacked' :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
    luPacked' = Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
    luC
        luSolve' :: (Matrix (Complex Double), [Int])
    -> Matrix (Complex Double) -> Matrix (Complex Double)
    luSolve' (Matrix (Complex Double)
    l_u,[Int]
    perm) = Matrix (Complex Double)
    -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
    lusC Matrix (Complex Double)
    l_u [Int]
    perm
        linearSolve' :: Matrix (Complex Double)
    -> Matrix (Complex Double) -> Matrix (Complex Double)
    linearSolve' = Matrix (Complex Double)
    -> Matrix (Complex Double) -> Matrix (Complex Double)
    linearSolveC
        mbLinearSolve' :: Matrix (Complex Double)
    -> Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
    mbLinearSolve' = Matrix (Complex Double)
    -> Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
    mbLinearSolveC
        cholSolve' :: Matrix (Complex Double)
    -> Matrix (Complex Double) -> Matrix (Complex Double)
    cholSolve' = Matrix (Complex Double)
    -> Matrix (Complex Double) -> Matrix (Complex Double)
    cholSolveC
        triSolve' :: UpLo
    -> Matrix (Complex Double)
    -> Matrix (Complex Double)
    -> Matrix (Complex Double)
    triSolve' = UpLo
    -> Matrix (Complex Double)
    -> Matrix (Complex Double)
    -> Matrix (Complex Double)
    triSolveC
        triDiagSolve' :: Vector (Complex Double)
    -> Vector (Complex Double)
    -> Vector (Complex Double)
    -> Matrix (Complex Double)
    -> Matrix (Complex Double)
    triDiagSolve' = Vector (Complex Double)
    -> Vector (Complex Double)
    -> Vector (Complex Double)
    -> Matrix (Complex Double)
    -> Matrix (Complex Double)
    triDiagSolveC
        linearSolveLS' :: Matrix (Complex Double)
    -> Matrix (Complex Double) -> Matrix (Complex Double)
    linearSolveLS' = Matrix (Complex Double)
    -> Matrix (Complex Double) -> Matrix (Complex Double)
    linearSolveLSC
        linearSolveSVD' :: Matrix (Complex Double)
    -> Matrix (Complex Double) -> Matrix (Complex Double)
    linearSolveSVD' = Maybe Double
    -> Matrix (Complex Double)
    -> Matrix (Complex Double)
    -> Matrix (Complex Double)
    linearSolveSVDC Maybe Double
    forall a. Maybe a
    Nothing
        eig' :: Matrix (Complex Double)
    -> (Vector (Complex Double), Matrix (Complex Double))
    eig' = Matrix (Complex Double)
    -> (Vector (Complex Double), Matrix (Complex Double))
    eigC
        geig' :: Matrix (Complex Double)
    -> Matrix (Complex Double)
    -> (Vector (Complex Double), Vector (Complex Double),
        Matrix (Complex Double))
    geig' = Matrix (Complex Double)
    -> Matrix (Complex Double)
    -> (Vector (Complex Double), Vector (Complex Double),
        Matrix (Complex Double))
    eigGC
        eigOnly :: Matrix (Complex Double) -> Vector (Complex Double)
    eigOnly = Matrix (Complex Double) -> Vector (Complex Double)
    eigOnlyC
        geigOnly :: Matrix (Complex Double)
    -> Matrix (Complex Double)
    -> (Vector (Complex Double), Vector (Complex Double))
    geigOnly = Matrix (Complex Double)
    -> Matrix (Complex Double)
    -> (Vector (Complex Double), Vector (Complex Double))
    eigOnlyGC
        eigSH'' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
    eigSH'' = Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
    eigH
        eigOnlySH :: Matrix (Complex Double) -> Vector Double
    eigOnlySH = Matrix (Complex Double) -> Vector Double
    eigOnlyH
        cholSH' :: Matrix (Complex Double) -> Matrix (Complex Double)
    cholSH' = Matrix (Complex Double) -> Matrix (Complex Double)
    cholH
        mbCholSH' :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
    mbCholSH' = Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
    mbCholH
        qr' :: Matrix (Complex Double)
    -> (Matrix (Complex Double), Vector (Complex Double))
    qr' = Matrix (Complex Double)
    -> (Matrix (Complex Double), Vector (Complex Double))
    qrC
        qrgr' :: Int
    -> (Matrix (Complex Double), Vector (Complex Double))
    -> Matrix (Complex Double)
    qrgr' = Int
    -> (Matrix (Complex Double), Vector (Complex Double))
    -> Matrix (Complex Double)
    qrgrC
        hess' :: Matrix (Complex Double)
    -> (Matrix (Complex Double), Matrix (Complex Double))
    hess' = (Matrix (Complex Double)
     -> (Matrix (Complex Double), Vector (Complex Double)))
    -> Matrix (Complex Double)
    -> (Matrix (Complex Double), Matrix (Complex Double))
    forall t.
    Field t =>
    (Matrix t -> (Matrix t, Vector t))
    -> Matrix t -> (Matrix t, Matrix t)
    unpackHess Matrix (Complex Double)
    -> (Matrix (Complex Double), Vector (Complex Double))
    hessC
        schur' :: Matrix (Complex Double)
    -> (Matrix (Complex Double), Matrix (Complex Double))
    schur' = Matrix (Complex Double)
    -> (Matrix (Complex Double), Matrix (Complex Double))
    schurC
        ldlPacked' :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
    ldlPacked' = Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
    ldlC
        ldlSolve' :: (Matrix (Complex Double), [Int])
    -> Matrix (Complex Double) -> Matrix (Complex Double)
    ldlSolve' = (Matrix (Complex Double)
     -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double))
    -> (Matrix (Complex Double), [Int])
    -> Matrix (Complex Double)
    -> Matrix (Complex Double)
    forall a b c. (a -> b -> c) -> (a, b) -> c
    uncurry Matrix (Complex Double)
    -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
    ldlsC

    --------------------------------------------------------------

    square :: Matrix t -> Bool
    square Matrix t
    m = Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    == Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    m

    vertical :: Matrix t -> Bool
    vertical Matrix t
    m = Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    >= Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    m

    exactHermitian :: Matrix e -> Bool
    exactHermitian Matrix e
    m = Matrix e
    m Matrix e -> Matrix e -> Bool
    forall (c :: * -> *) e. Container c e => c e -> c e -> Bool
    `equal` Matrix e -> Matrix e
    forall t. CTrans t => Matrix t -> Matrix t
    ctrans Matrix e
    m

    --------------------------------------------------------------

    {- | Full singular value decomposition.

    @
    a = (5><3)
     [  1.0,  2.0,  3.0
     ,  4.0,  5.0,  6.0
     ,  7.0,  8.0,  9.0
     , 10.0, 11.0, 12.0
     , 13.0, 14.0, 15.0 ] :: Matrix Double
    @

    >>> let (u,s,v) = svd a

    >>> disp 3 u
    5x5
    -0.101   0.768   0.614   0.028  -0.149
    -0.249   0.488  -0.503   0.172   0.646
    -0.396   0.208  -0.405  -0.660  -0.449
    -0.543  -0.072  -0.140   0.693  -0.447
    -0.690  -0.352   0.433  -0.233   0.398

    >>> s
    [35.18264833189422,1.4769076999800903,1.089145439970417e-15]
    it :: Vector Double

    >>> disp 3 v
    3x3
    -0.519  -0.751   0.408
    -0.576  -0.046  -0.816
    -0.632   0.659   0.408

    >>> let d = diagRect 0 s 5 3
    >>> disp 3 d
    5x3
    35.183  0.000  0.000
     0.000  1.477  0.000
     0.000  0.000  0.000
     0.000  0.000  0.000

    >>> disp 3 $ u <> d <> tr v
    5x3
     1.000   2.000   3.000
     4.000   5.000   6.000
     7.000   8.000   9.000
    10.000  11.000  12.000
    13.000  14.000  15.000

    -}
    svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
    svd :: forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    svd = {-# SCC "svd" #-} (Matrix t, Vector Double, Matrix t)
    -> (Matrix t, Vector Double, Matrix t)
    forall {m} {c} {a} {b}. Transposable m c => (a, b, m) -> (a, b, c)
    g ((Matrix t, Vector Double, Matrix t)
     -> (Matrix t, Vector Double, Matrix t))
    -> (Matrix t -> (Matrix t, Vector Double, Matrix t))
    -> Matrix t
    -> (Matrix t, Vector Double, Matrix t)
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> (Matrix t, Vector Double, Matrix t)
    forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    svd'
      where
        g :: (a, b, m) -> (a, b, c)
    g (a
    u,b
    s,m
    v) = (a
    u,b
    s,m -> c
    forall m mt. Transposable m mt => m -> mt
    tr m
    v)

    {- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@.

    If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> tr v@.

    @
    a = (5><3)
     [  1.0,  2.0,  3.0
     ,  4.0,  5.0,  6.0
     ,  7.0,  8.0,  9.0
     , 10.0, 11.0, 12.0
     , 13.0, 14.0, 15.0 ] :: Matrix Double
    @

    >>> let (u,s,v) = thinSVD a

    >>> disp 3 u
    5x3
    -0.101   0.768   0.614
    -0.249   0.488  -0.503
    -0.396   0.208  -0.405
    -0.543  -0.072  -0.140
    -0.690  -0.352   0.433

    >>> s
    [35.18264833189422,1.4769076999800903,1.089145439970417e-15]
    it :: Vector Double

    >>> disp 3 v
    3x3
    -0.519  -0.751   0.408
    -0.576  -0.046  -0.816
    -0.632   0.659   0.408

    >>> disp 3 $ u <> diag s <> tr v
    5x3
     1.000   2.000   3.000
     4.000   5.000   6.000
     7.000   8.000   9.000
    10.000  11.000  12.000
    13.000  14.000  15.000

    -}
    thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
    thinSVD :: forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    thinSVD = {-# SCC "thinSVD" #-} (Matrix t, Vector Double, Matrix t)
    -> (Matrix t, Vector Double, Matrix t)
    forall {m} {c} {a} {b}. Transposable m c => (a, b, m) -> (a, b, c)
    g ((Matrix t, Vector Double, Matrix t)
     -> (Matrix t, Vector Double, Matrix t))
    -> (Matrix t -> (Matrix t, Vector Double, Matrix t))
    -> Matrix t
    -> (Matrix t, Vector Double, Matrix t)
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> (Matrix t, Vector Double, Matrix t)
    forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    thinSVD'
      where
        g :: (a, b, m) -> (a, b, c)
    g (a
    u,b
    s,m
    v) = (a
    u,b
    s,m -> c
    forall m mt. Transposable m mt => m -> mt
    tr m
    v)

    -- | Singular values only.
    singularValues :: Field t => Matrix t -> Vector Double
    singularValues :: forall t. Field t => Matrix t -> Vector Double
    singularValues = {-# SCC "singularValues" #-} Matrix t -> Vector Double
    forall t. Field t => Matrix t -> Vector Double
    sv'

    -- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
    --
    -- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> tr v@.
    fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
    fullSVD :: forall t.
    Field t =>
    Matrix t -> (Matrix t, Matrix Double, Matrix t)
    fullSVD Matrix t
    m = (Matrix t
    u,Matrix Double
    d,Matrix t
    v) where
        (Matrix t
    u,Vector Double
    s,Matrix t
    v) = Matrix t -> (Matrix t, Vector Double, Matrix t)
    forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    svd Matrix t
    m
        d :: Matrix Double
    d = Double -> Vector Double -> Int -> Int -> Matrix Double
    forall t. Storable t => t -> Vector t -> Int -> Int -> Matrix t
    diagRect Double
    0 Vector Double
    s Int
    r Int
    c
        r :: Int
    r = Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m
        c :: Int
    c = Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    m

    {- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors.

    @
    a = (5><3)
     [  1.0,  2.0,  3.0
     ,  4.0,  5.0,  6.0
     ,  7.0,  8.0,  9.0
     , 10.0, 11.0, 12.0
     , 13.0, 14.0, 15.0 ] :: Matrix Double
    @

    >>> let (u,s,v) = compactSVD a

    >>> disp 3 u
    5x2
    -0.101   0.768
    -0.249   0.488
    -0.396   0.208
    -0.543  -0.072
    -0.690  -0.352

    >>> s
    [35.18264833189422,1.476907699980091]
    it :: Vector Double

    >>> disp 3 u
    5x2
    -0.101   0.768
    -0.249   0.488
    -0.396   0.208
    -0.543  -0.072
    -0.690  -0.352

    >>> disp 3 $ u <> diag s <> tr v
    5x3
     1.000   2.000   3.000
     4.000   5.000   6.000
     7.000   8.000   9.000
    10.000  11.000  12.000
    13.000  14.000  15.000

    -}
    compactSVD :: Field t  => Matrix t -> (Matrix t, Vector Double, Matrix t)
    compactSVD :: forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    compactSVD = Double -> Matrix t -> (Matrix t, Vector Double, Matrix t)
    forall t.
    Field t =>
    Double -> Matrix t -> (Matrix t, Vector Double, Matrix t)
    compactSVDTol Double
    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 :: forall t.
    Field t =>
    Double -> Matrix t -> (Matrix t, Vector Double, Matrix t)
    compactSVDTol Double
    r Matrix t
    m = (Matrix t
    u', Int -> Int -> Vector Double -> Vector Double
    forall t. Storable t => Int -> Int -> Vector t -> Vector t
    subVector Int
    0 Int
    d Vector Double
    s, Matrix t
    v') where
        (Matrix t
    u,Vector Double
    s,Matrix t
    v) = Matrix t -> (Matrix t, Vector Double, Matrix t)
    forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    thinSVD Matrix t
    m
        d :: Int
    d = Double -> Matrix t -> Vector Double -> Int
    forall t. Element t => Double -> Matrix t -> Vector Double -> Int
    rankSVD (Double
    rDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    eps) Matrix t
    m Vector Double
    s Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    `max` Int
    1
        u' :: Matrix t
    u' = Int -> Matrix t -> Matrix t
    forall t. Element t => Int -> Matrix t -> Matrix t
    takeColumns Int
    d Matrix t
    u
        v' :: Matrix t
    v' = Int -> Matrix t -> Matrix t
    forall t. Element t => Int -> Matrix t -> Matrix t
    takeColumns Int
    d Matrix t
    v

    -- | Singular values and all right singular vectors (as columns).
    rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
    rightSV :: forall t. Field t => Matrix t -> (Vector Double, Matrix t)
    rightSV Matrix t
    m | Matrix t -> Bool
    forall {t}. Matrix t -> Bool
    vertical Matrix t
    m = let (Matrix t
    _,Vector Double
    s,Matrix t
    v) = Matrix t -> (Matrix t, Vector Double, Matrix t)
    forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    thinSVD Matrix t
    m in (Vector Double
    s,Matrix t
    v)
              | Bool
    otherwise  = let (Matrix t
    _,Vector Double
    s,Matrix t
    v) = Matrix t -> (Matrix t, Vector Double, Matrix t)
    forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    svd Matrix t
    m     in (Vector Double
    s,Matrix t
    v)

    -- | Singular values and all left singular vectors (as columns).
    leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
    leftSV :: forall t. Field t => Matrix t -> (Matrix t, Vector Double)
    leftSV Matrix t
    m  | Matrix t -> Bool
    forall {t}. Matrix t -> Bool
    vertical Matrix t
    m = let (Matrix t
    u,Vector Double
    s,Matrix t
    _) = Matrix t -> (Matrix t, Vector Double, Matrix t)
    forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    svd Matrix t
    m     in (Matrix t
    u,Vector Double
    s)
              | Bool
    otherwise  = let (Matrix t
    u,Vector Double
    s,Matrix t
    _) = Matrix t -> (Matrix t, Vector Double, Matrix t)
    forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    thinSVD Matrix t
    m in (Matrix t
    u,Vector Double
    s)

    --------------------------------------------------------------

    -- | LU decomposition of a matrix in a compact format.
    data LU t = LU (Matrix t) [Int] deriving Int -> LU t -> ShowS
    [LU t] -> ShowS
    LU t -> String
    (Int -> LU t -> ShowS)
    -> (LU t -> String) -> ([LU t] -> ShowS) -> Show (LU t)
    forall t. (Show t, Element t) => Int -> LU t -> ShowS
    forall t. (Show t, Element t) => [LU t] -> ShowS
    forall t. (Show t, Element t) => LU t -> String
    forall a.
    (Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
    $cshowsPrec :: forall t. (Show t, Element t) => Int -> LU t -> ShowS
    showsPrec :: Int -> LU t -> ShowS
    $cshow :: forall t. (Show t, Element t) => LU t -> String
    show :: LU t -> String
    $cshowList :: forall t. (Show t, Element t) => [LU t] -> ShowS
    showList :: [LU t] -> ShowS
    Show

    instance (NFData t, Numeric t) => NFData (LU t)
      where
        rnf :: LU t -> ()
    rnf (LU Matrix t
    m [Int]
    _) = Matrix t -> ()
    forall a. NFData a => a -> ()
    rnf Matrix t
    m

    -- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
    luPacked :: Field t => Matrix t -> LU t
    luPacked :: forall t. Field t => Matrix t -> LU t
    luPacked Matrix t
    x = {-# SCC "luPacked" #-} Matrix t -> [Int] -> LU t
    forall t. Matrix t -> [Int] -> LU t
    LU Matrix t
    m [Int]
    p
      where
        (Matrix t
    m,[Int]
    p) = Matrix t -> (Matrix t, [Int])
    forall t. Field t => Matrix t -> (Matrix t, [Int])
    luPacked' Matrix t
    x

    -- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'.
    luSolve :: Field t => LU t -> Matrix t -> Matrix t
    luSolve :: forall t. Field t => LU t -> Matrix t -> Matrix t
    luSolve (LU Matrix t
    m [Int]
    p) = {-# SCC "luSolve" #-} (Matrix t, [Int]) -> Matrix t -> Matrix t
    forall t. Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
    luSolve' (Matrix t
    m,[Int]
    p)

    -- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'.
    -- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.
    linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t
    linearSolve :: forall t. Field t => Matrix t -> Matrix t -> Matrix t
    linearSolve = {-# SCC "linearSolve" #-} Matrix t -> Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t -> Matrix t
    linearSolve'

    -- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'.
    mbLinearSolve :: Field t => Matrix t -> Matrix t -> Maybe (Matrix t)
    mbLinearSolve :: forall t. Field t => Matrix t -> Matrix t -> Maybe (Matrix t)
    mbLinearSolve = {-# SCC "linearSolve" #-} Matrix t -> Matrix t -> Maybe (Matrix t)
    forall t. Field t => Matrix t -> Matrix t -> Maybe (Matrix t)
    mbLinearSolve'

    -- | Solve a symmetric or Hermitian positive definite linear system using a precomputed Cholesky decomposition obtained by 'chol'.
    cholSolve
        :: Field t
        => Matrix t -- ^ Cholesky decomposition of the coefficient matrix
        -> Matrix t -- ^ right hand sides
        -> Matrix t -- ^ solution
    cholSolve :: forall t. Field t => Matrix t -> Matrix t -> Matrix t
    cholSolve = {-# SCC "cholSolve" #-} Matrix t -> Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t -> Matrix t
    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 :: forall t. Field t => UpLo -> Matrix t -> Matrix t -> Matrix t
    triSolve = {-# SCC "triSolve" #-} UpLo -> Matrix t -> Matrix t -> Matrix t
    forall t. Field t => UpLo -> Matrix t -> Matrix t -> Matrix t
    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 :: forall t.
    Field t =>
    Vector t -> Vector t -> Vector t -> Matrix t -> Matrix t
    triDiagSolve = {-# SCC "triDiagSolve" #-} Vector t -> Vector t -> Vector t -> Matrix t -> Matrix t
    forall t.
    Field t =>
    Vector t -> Vector t -> Vector t -> Matrix t -> Matrix t
    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 :: forall t. Field t => Matrix t -> Matrix t -> Matrix t
    linearSolveSVD = {-# SCC "linearSolveSVD" #-} Matrix t -> Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t -> Matrix t
    linearSolveSVD'

    -- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'.
    linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t
    linearSolveLS :: forall t. Field t => Matrix t -> Matrix t -> Matrix t
    linearSolveLS = {-# SCC "linearSolveLS" #-} Matrix t -> Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t -> Matrix t
    linearSolveLS'

    --------------------------------------------------------------------------------

    -- | LDL decomposition of a complex Hermitian or real symmetric matrix in a compact format.
    data LDL t = LDL (Matrix t) [Int] deriving Int -> LDL t -> ShowS
    [LDL t] -> ShowS
    LDL t -> String
    (Int -> LDL t -> ShowS)
    -> (LDL t -> String) -> ([LDL t] -> ShowS) -> Show (LDL t)
    forall t. (Show t, Element t) => Int -> LDL t -> ShowS
    forall t. (Show t, Element t) => [LDL t] -> ShowS
    forall t. (Show t, Element t) => LDL t -> String
    forall a.
    (Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
    $cshowsPrec :: forall t. (Show t, Element t) => Int -> LDL t -> ShowS
    showsPrec :: Int -> LDL t -> ShowS
    $cshow :: forall t. (Show t, Element t) => LDL t -> String
    show :: LDL t -> String
    $cshowList :: forall t. (Show t, Element t) => [LDL t] -> ShowS
    showList :: [LDL t] -> ShowS
    Show

    instance (NFData t, Numeric t) => NFData (LDL t)
      where
        rnf :: LDL t -> ()
    rnf (LDL Matrix t
    m [Int]
    _) = Matrix t -> ()
    forall a. NFData a => a -> ()
    rnf Matrix t
    m

    -- | Similar to 'ldlPacked', without checking that the input matrix is hermitian or symmetric. It works with the lower triangular part.
    ldlPackedSH :: Field t => Matrix t -> LDL t
    ldlPackedSH :: forall t. Field t => Matrix t -> LDL t
    ldlPackedSH Matrix t
    x = {-# SCC "ldlPacked" #-} Matrix t -> [Int] -> LDL t
    forall t. Matrix t -> [Int] -> LDL t
    LDL Matrix t
    m [Int]
    p
      where
       (Matrix t
    m,[Int]
    p) = Matrix t -> (Matrix t, [Int])
    forall t. Field t => Matrix t -> (Matrix t, [Int])
    ldlPacked' Matrix t
    x

    -- | Obtains the LDL decomposition of a matrix in a compact data structure suitable for 'ldlSolve'.
    ldlPacked :: Field t => Herm t -> LDL t
    ldlPacked :: forall t. Field t => Herm t -> LDL t
    ldlPacked (Herm Matrix t
    m) = Matrix t -> LDL t
    forall t. Field t => Matrix t -> LDL t
    ldlPackedSH Matrix t
    m

    -- | Solution of a linear system (for several right hand sides) from a precomputed LDL factorization obtained by 'ldlPacked'.
    --
    --   Note: this can be slower than the general solver based on the LU decomposition.
    ldlSolve :: Field t => LDL t -> Matrix t -> Matrix t
    ldlSolve :: forall t. Field t => LDL t -> Matrix t -> Matrix t
    ldlSolve (LDL Matrix t
    m [Int]
    p) = {-# SCC "ldlSolve" #-} (Matrix t, [Int]) -> Matrix t -> Matrix t
    forall t. Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
    ldlSolve' (Matrix t
    m,[Int]
    p)

    --------------------------------------------------------------

    {- | Eigenvalues (not ordered) and eigenvectors (as columns) of a general square matrix.

    If @(s,v) = eig m@ then @m \<> v == v \<> diag s@

    @
    a = (3><3)
     [ 3, 0, -2
     , 4, 5, -1
     , 3, 1,  0 ] :: Matrix Double
    @

    >>> let (l, v) = eig a

    >>> putStr . dispcf 3 . asRow $ l
    1x3
    1.925+1.523i  1.925-1.523i  4.151

    >>> putStr . dispcf 3 $ v
    3x3
    -0.455+0.365i  -0.455-0.365i   0.181
            0.603          0.603  -0.978
     0.033+0.543i   0.033-0.543i  -0.104

    >>> putStr . dispcf 3 $ complex a <> v
    3x3
    -1.432+0.010i  -1.432-0.010i   0.753
     1.160+0.918i   1.160-0.918i  -4.059
    -0.763+1.096i  -0.763-1.096i  -0.433

    >>> putStr . dispcf 3 $ v <> diag l
    3x3
    -1.432+0.010i  -1.432-0.010i   0.753
     1.160+0.918i   1.160-0.918i  -4.059
    -0.763+1.096i  -0.763-1.096i  -0.433

    -}
    eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
    eig :: forall t.
    Field t =>
    Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
    eig = {-# SCC "eig" #-} Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
    forall t.
    Field t =>
    Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
    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 :: forall t.
    Field t =>
    Matrix t
    -> Matrix t
    -> (Vector (Complex Double), Vector t, Matrix (Complex Double))
    geig = {-# SCC "geig" #-} Matrix t
    -> Matrix t
    -> (Vector (Complex Double), Vector t, Matrix (Complex Double))
    forall t.
    Field t =>
    Matrix t
    -> Matrix t
    -> (Vector (Complex Double), Vector t, Matrix (Complex Double))
    geig'

    -- | Eigenvalues (not ordered) of a general square matrix.
    eigenvalues :: Field t => Matrix t -> Vector (Complex Double)
    eigenvalues :: forall t. Field t => Matrix t -> Vector (Complex Double)
    eigenvalues = {-# SCC "eigenvalues" #-} Matrix t -> Vector (Complex Double)
    forall t. Field t => Matrix t -> Vector (Complex Double)
    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 :: forall t.
    Field t =>
    Matrix t -> Matrix t -> (Vector (Complex Double), Vector t)
    geigenvalues = {-# SCC "geigenvalues" #-} Matrix t -> Matrix t -> (Vector (Complex Double), Vector t)
    forall t.
    Field t =>
    Matrix t -> Matrix t -> (Vector (Complex Double), Vector t)
    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' :: forall t. Field t => Matrix t -> (Vector Double, Matrix t)
    eigSH' = {-# SCC "eigSH'" #-} Matrix t -> (Vector Double, Matrix t)
    forall t. Field t => Matrix t -> (Vector Double, Matrix t)
    eigSH''

    -- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
    eigenvaluesSH' :: Field t => Matrix t -> Vector Double
    eigenvaluesSH' :: forall t. Field t => Matrix t -> Vector Double
    eigenvaluesSH' = {-# SCC "eigenvaluesSH'" #-} Matrix t -> Vector Double
    forall t. Field t => Matrix t -> Vector Double
    eigOnlySH

    {- | Eigenvalues and eigenvectors (as columns) of a complex hermitian or real symmetric matrix, in descending order.

    If @(s,v) = eigSH m@ then @m == v \<> diag s \<> tr v@

    @
    a = (3><3)
     [ 1.0, 2.0, 3.0
     , 2.0, 4.0, 5.0
     , 3.0, 5.0, 6.0 ]
    @

    >>> let (l, v) = eigSH a

    >>> l
    [11.344814282762075,0.17091518882717918,-0.5157294715892575]

    >>> disp 3 $ v <> diag l <> tr v
    3x3
    1.000  2.000  3.000
    2.000  4.000  5.000
    3.000  5.000  6.000

    -}
    eigSH :: Field t => Herm t -> (Vector Double, Matrix t)
    eigSH :: forall t. Field t => Herm t -> (Vector Double, Matrix t)
    eigSH (Herm Matrix t
    m) = Matrix t -> (Vector Double, Matrix t)
    forall t. Field t => Matrix t -> (Vector Double, Matrix t)
    eigSH' Matrix t
    m

    -- | Eigenvalues (in descending order) of a complex hermitian or real symmetric matrix.
    eigenvaluesSH :: Field t => Herm t -> Vector Double
    eigenvaluesSH :: forall t. Field t => Herm t -> Vector Double
    eigenvaluesSH (Herm Matrix t
    m) = Matrix t -> Vector Double
    forall t. Field t => Matrix t -> Vector Double
    eigenvaluesSH' Matrix t
    m

    --------------------------------------------------------------

    -- | QR decomposition of a matrix in compact form. (The orthogonal matrix is not explicitly formed.)
    data QR t = QR (Matrix t) (Vector t)

    instance (NFData t, Numeric t) => NFData (QR t)
      where
        rnf :: QR t -> ()
    rnf (QR Matrix t
    m Vector t
    _) = Matrix t -> ()
    forall a. NFData a => a -> ()
    rnf Matrix t
    m

    -- | 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 :: forall t. Field t => Matrix t -> (Matrix t, Matrix t)
    qr = {-# SCC "qr" #-} (Matrix t, Vector t) -> (Matrix t, Matrix t)
    forall t. Field t => (Matrix t, Vector t) -> (Matrix t, Matrix t)
    unpackQR ((Matrix t, Vector t) -> (Matrix t, Matrix t))
    -> (Matrix t -> (Matrix t, Vector t))
    -> Matrix t
    -> (Matrix t, Matrix t)
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> (Matrix t, Vector t)
    forall t. Field t => Matrix t -> (Matrix t, Vector t)
    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 :: forall t. Field t => Matrix t -> (Matrix t, Matrix t)
    thinQR = {-# SCC "thinQR" #-} (Matrix t, Vector t) -> (Matrix t, Matrix t)
    forall t. Field t => (Matrix t, Vector t) -> (Matrix t, Matrix t)
    thinUnpackQR ((Matrix t, Vector t) -> (Matrix t, Matrix t))
    -> (Matrix t -> (Matrix t, Vector t))
    -> Matrix t
    -> (Matrix t, Matrix t)
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> (Matrix t, Vector t)
    forall t. Field t => Matrix t -> (Matrix t, Vector t)
    qr'

    -- | Compute the QR decomposition of a matrix in compact form.
    qrRaw :: Field t => Matrix t -> QR t
    qrRaw :: forall t. Field t => Matrix t -> QR t
    qrRaw Matrix t
    m = Matrix t -> Vector t -> QR t
    forall t. Matrix t -> Vector t -> QR t
    QR Matrix t
    x Vector t
    v
      where
        (Matrix t
    x,Vector t
    v) = Matrix t -> (Matrix t, Vector t)
    forall t. Field t => Matrix t -> (Matrix t, Vector t)
    qr' Matrix t
    m

    -- | generate a matrix with k orthogonal columns from the compact QR decomposition obtained by 'qrRaw'.
    --
    qrgr :: Field t => Int -> QR t -> Matrix t
    qrgr :: forall t. Field t => Int -> QR t -> Matrix t
    qrgr Int
    n (QR Matrix t
    a Vector t
    t)
        | Vector t -> Int
    forall t. Storable t => Vector t -> Int
    dim Vector t
    t Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    > Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    min (Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    a) (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    a) Bool -> Bool -> Bool
    || Int
    n Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    < Int
    0 Bool -> Bool -> Bool
    || Int
    n Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    > Vector t -> Int
    forall t. Storable t => Vector t -> Int
    dim Vector t
    t = String -> Matrix t
    forall a. HasCallStack => String -> a
    error String
    "qrgr expects k <= min(rows,cols)"
        | Bool
    otherwise = Int -> (Matrix t, Vector t) -> Matrix t
    forall t. Field t => Int -> (Matrix t, Vector t) -> Matrix t
    qrgr' Int
    n (Matrix t
    a,Vector t
    t)

    -- | 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 :: forall t. Field t => Matrix t -> (Matrix t, Matrix t)
    rq = {-# SCC "rq" #-} (Matrix t -> (Matrix t, Matrix t))
    -> Matrix t -> (Matrix t, Matrix t)
    forall t.
    Field t =>
    (Matrix t -> (Matrix t, Matrix t))
    -> Matrix t -> (Matrix t, Matrix t)
    rqFromQR Matrix t -> (Matrix t, Matrix t)
    forall t. Field t => Matrix t -> (Matrix t, Matrix t)
    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 :: forall t. Field t => Matrix t -> (Matrix t, Matrix t)
    thinRQ = {-# SCC "thinQR" #-} (Matrix t -> (Matrix t, Matrix t))
    -> Matrix t -> (Matrix t, Matrix t)
    forall t.
    Field t =>
    (Matrix t -> (Matrix t, Matrix t))
    -> Matrix t -> (Matrix t, Matrix t)
    rqFromQR Matrix t -> (Matrix t, Matrix t)
    forall t. Field t => Matrix t -> (Matrix t, Matrix t)
    thinQR

    rqFromQR :: Field t => (Matrix t -> (Matrix t, Matrix t)) -> Matrix t -> (Matrix t, Matrix t)
    rqFromQR :: forall t.
    Field t =>
    (Matrix t -> (Matrix t, Matrix t))
    -> Matrix t -> (Matrix t, Matrix t)
    rqFromQR Matrix t -> (Matrix t, Matrix t)
    qr0 Matrix t
    m = (Matrix t
    r,Matrix t
    q) where
        (Matrix t
    q',Matrix t
    r') = Matrix t -> (Matrix t, Matrix t)
    qr0 (Matrix t -> (Matrix t, Matrix t))
    -> Matrix t -> (Matrix t, Matrix t)
    forall a b. (a -> b) -> a -> b
    $ Matrix t -> Matrix t
    forall t. Matrix t -> Matrix t
    trans (Matrix t -> Matrix t) -> Matrix t -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ Matrix t -> Matrix t
    rev1 Matrix t
    m
        r :: Matrix t
    r = Matrix t -> Matrix t
    rev2 (Matrix t -> Matrix t
    forall t. Matrix t -> Matrix t
    trans Matrix t
    r')
        q :: Matrix t
    q = Matrix t -> Matrix t
    rev2 (Matrix t -> Matrix t
    forall t. Matrix t -> Matrix t
    trans Matrix t
    q')
        rev1 :: Matrix t -> Matrix t
    rev1 = Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t
    flipud (Matrix t -> Matrix t)
    -> (Matrix t -> Matrix t) -> Matrix t -> Matrix t
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t
    fliprl
        rev2 :: Matrix t -> Matrix t
    rev2 = Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t
    fliprl (Matrix t -> Matrix t)
    -> (Matrix t -> Matrix t) -> Matrix t -> Matrix t
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t
    flipud

    -- | Hessenberg factorization.
    --
    -- If @(p,h) = hess m@ then @m == p \<> h \<> tr p@, where p is unitary
    -- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal).
    hess        :: Field t => Matrix t -> (Matrix t, Matrix t)
    hess :: forall t. Field t => Matrix t -> (Matrix t, Matrix t)
    hess = Matrix t -> (Matrix t, Matrix t)
    forall t. Field t => Matrix t -> (Matrix t, Matrix t)
    hess'

    -- | Schur factorization.
    --
    -- If @(u,s) = schur m@ then @m == u \<> s \<> tr u@, where u is unitary
    -- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is
    -- upper triangular in 2x2 blocks.
    --
    -- \"Anything that the Jordan decomposition can do, the Schur decomposition
    -- can do better!\" (Van Loan)
    schur       :: Field t => Matrix t -> (Matrix t, Matrix t)
    schur :: forall t. Field t => Matrix t -> (Matrix t, Matrix t)
    schur = Matrix t -> (Matrix t, Matrix t)
    forall t. Field t => Matrix t -> (Matrix t, Matrix t)
    schur'

    -- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'.
    mbCholSH :: Field t => Matrix t -> Maybe (Matrix t)
    mbCholSH :: forall t. Field t => Matrix t -> Maybe (Matrix t)
    mbCholSH = {-# SCC "mbCholSH" #-} Matrix t -> Maybe (Matrix t)
    forall t. Field t => Matrix t -> Maybe (Matrix t)
    mbCholSH'

    -- | Similar to 'chol', without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
    cholSH      :: Field t => Matrix t -> Matrix t
    cholSH :: forall t. Field t => Matrix t -> Matrix t
    cholSH = Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t
    cholSH'

    -- | Cholesky factorization of a positive definite hermitian or symmetric matrix.
    --
    -- If @c = chol m@ then @c@ is upper triangular and @m == tr c \<> c@.
    chol :: Field t => Herm t ->  Matrix t
    chol :: forall t. Field t => Herm t -> Matrix t
    chol (Herm Matrix t
    m) = {-# SCC "chol" #-} Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t
    cholSH' Matrix t
    m

    -- | Similar to 'chol', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'.
    mbChol :: Field t => Herm t -> Maybe (Matrix t)
    mbChol :: forall t. Field t => Herm t -> Maybe (Matrix t)
    mbChol (Herm Matrix t
    m) = {-# SCC "mbChol" #-} Matrix t -> Maybe (Matrix t)
    forall t. Field t => Matrix t -> Maybe (Matrix t)
    mbCholSH' Matrix t
    m

    -- | Joint computation of inverse and logarithm of determinant of a square matrix.
    invlndet :: Field t
             => Matrix t
             -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det))
    invlndet :: forall t. Field t => Matrix t -> (Matrix t, (t, t))
    invlndet Matrix t
    m | Matrix t -> Bool
    forall {t}. Matrix t -> Bool
    square Matrix t
    m = (Matrix t
    im,(t
    ladm,t
    sdm))
               | Bool
    otherwise = String -> (Matrix t, (t, t))
    forall a. HasCallStack => String -> a
    error (String -> (Matrix t, (t, t))) -> String -> (Matrix t, (t, t))
    forall a b. (a -> b) -> a -> b
    $ String
    "invlndet of nonsquare "String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ Matrix t -> String
    forall t. Matrix t -> String
    shSize Matrix t
    m String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ String
    " matrix"
      where
        lp :: LU t
    lp@(LU Matrix t
    lup [Int]
    perm) = Matrix t -> LU t
    forall t. Field t => Matrix t -> LU t
    luPacked Matrix t
    m
        s :: t
    s = Int -> [Int] -> t
    forall {a} {b}. (Eq a, Num b, Num a, Enum a) => a -> [a] -> b
    signlp (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m) [Int]
    perm
        dg :: [t]
    dg = Vector t -> [t]
    forall a. Storable a => Vector a -> [a]
    toList (Vector t -> [t]) -> Vector t -> [t]
    forall a b. (a -> b) -> a -> b
    $ Matrix t -> Vector t
    forall t. Element t => Matrix t -> Vector t
    takeDiag (Matrix t -> Vector t) -> Matrix t -> Vector t
    forall a b. (a -> b) -> a -> b
    $ Matrix t
    lup
        ladm :: t
    ladm = [t] -> t
    forall a. Num a => [a] -> a
    forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
    sum ([t] -> t) -> [t] -> t
    forall a b. (a -> b) -> a -> b
    $ (t -> t) -> [t] -> [t]
    forall a b. (a -> b) -> [a] -> [b]
    map (t -> t
    forall a. Floating a => a -> a
    log(t -> t) -> (t -> t) -> t -> t
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    .t -> t
    forall a. Num a => a -> a
    abs) [t]
    dg
        sdm :: t
    sdm = t
    st -> t -> t
    forall a. Num a => a -> a -> a
    * [t] -> t
    forall a. Num a => [a] -> a
    forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
    product ((t -> t) -> [t] -> [t]
    forall a b. (a -> b) -> [a] -> [b]
    map t -> t
    forall a. Num a => a -> a
    signum [t]
    dg)
        im :: Matrix t
    im = LU t -> Matrix t -> Matrix t
    forall t. Field t => LU t -> Matrix t -> Matrix t
    luSolve LU t
    lp (Int -> Matrix t
    forall a. (Num a, Element a) => Int -> Matrix a
    ident (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m))

    -- | Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'.
    det :: Field t => Matrix t -> t
    det :: forall t. Field t => Matrix t -> t
    det Matrix t
    m | Matrix t -> Bool
    forall {t}. Matrix t -> Bool
    square Matrix t
    m = {-# SCC "det" #-} t
    s t -> t -> t
    forall a. Num a => a -> a -> a
    * ([t] -> t
    forall a. Num a => [a] -> a
    forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
    product ([t] -> t) -> [t] -> t
    forall a b. (a -> b) -> a -> b
    $ Vector t -> [t]
    forall a. Storable a => Vector a -> [a]
    toList (Vector t -> [t]) -> Vector t -> [t]
    forall a b. (a -> b) -> a -> b
    $ Matrix t -> Vector t
    forall t. Element t => Matrix t -> Vector t
    takeDiag (Matrix t -> Vector t) -> Matrix t -> Vector t
    forall a b. (a -> b) -> a -> b
    $ Matrix t
    lup)
          | Bool
    otherwise = String -> t
    forall a. HasCallStack => String -> a
    error (String -> t) -> String -> t
    forall a b. (a -> b) -> a -> b
    $ String
    "det of nonsquare "String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ Matrix t -> String
    forall t. Matrix t -> String
    shSize Matrix t
    m String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ String
    " matrix"
        where
          LU Matrix t
    lup [Int]
    perm = Matrix t -> LU t
    forall t. Field t => Matrix t -> LU t
    luPacked Matrix t
    m
          s :: t
    s = Int -> [Int] -> t
    forall {a} {b}. (Eq a, Num b, Num a, Enum a) => a -> [a] -> b
    signlp (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m) [Int]
    perm

    -- | Explicit LU factorization of a general matrix.
    --
    -- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular,
    -- u is upper triangular, p is a permutation matrix and s is the signature of the permutation.
    lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
    lu :: forall t. Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
    lu = LU t -> (Matrix t, Matrix t, Matrix t, t)
    forall t. Numeric t => LU t -> (Matrix t, Matrix t, Matrix t, t)
    luFact (LU t -> (Matrix t, Matrix t, Matrix t, t))
    -> (Matrix t -> LU t)
    -> Matrix t
    -> (Matrix t, Matrix t, Matrix t, t)
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> LU t
    forall t. Field t => Matrix t -> LU t
    luPacked

    -- | Inverse of a square matrix. See also 'invlndet'.
    inv :: Field t => Matrix t -> Matrix t
    inv :: forall t. Field t => Matrix t -> Matrix t
    inv Matrix t
    m | Matrix t -> Bool
    forall {t}. Matrix t -> Bool
    square Matrix t
    m = Matrix t
    m Matrix t -> Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t -> Matrix t
    `linearSolve` Int -> Matrix t
    forall a. (Num a, Element a) => Int -> Matrix a
    ident (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m)
          | Bool
    otherwise = String -> Matrix t
    forall a. HasCallStack => String -> a
    error (String -> Matrix t) -> String -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ String
    "inv of nonsquare "String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ Matrix t -> String
    forall t. Matrix t -> String
    shSize Matrix t
    m String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ String
    " matrix"

    -- | Pseudoinverse of a general matrix with default tolerance ('pinvTol' 1, similar to GNU-Octave).
    pinv :: Field t => Matrix t -> Matrix t
    pinv :: forall t. Field t => Matrix t -> Matrix t
    pinv = Double -> Matrix t -> Matrix t
    forall t. Field t => Double -> Matrix t -> Matrix t
    pinvTol Double
    1

    {- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value.

    @
    m = (3><3) [ 1, 0,    0
               , 0, 1,    0
               , 0, 0, 1e-10] :: Matrix Double
    @

    >>> pinv m
    1. 0.           0.
    0. 1.           0.
    0. 0. 10000000000.

    >>> pinvTol 1E8 m
    1. 0. 0.
    0. 1. 0.
    0. 0. 1.

    -}

    pinvTol :: Field t => Double -> Matrix t -> Matrix t
    pinvTol :: forall t. Field t => Double -> Matrix t -> Matrix t
    pinvTol Double
    t Matrix t
    m = Matrix t
    v' Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    `mXm` Vector t -> Matrix t
    forall a. (Num a, Element a) => Vector a -> Matrix a
    diag Vector t
    s' Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    `mXm` Matrix t -> Matrix t
    forall t. CTrans t => Matrix t -> Matrix t
    ctrans Matrix t
    u' where
        (Matrix t
    u,Vector Double
    s,Matrix t
    v) = Matrix t -> (Matrix t, Vector Double, Matrix t)
    forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    thinSVD Matrix t
    m
        sl :: [Double]
    sl@(Double
    g:[Double]
    _) = Vector Double -> [Double]
    forall a. Storable a => Vector a -> [a]
    toList Vector Double
    s
        s' :: Vector t
    s' = Vector Double -> Vector t
    Vector (RealOf t) -> Vector t
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c (RealOf t) -> c t
    forall (c :: * -> *). Complexable c => c (RealOf t) -> c t
    real (Vector Double -> Vector t)
    -> ([Double] -> Vector Double) -> [Double] -> Vector t
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [Double] -> Vector Double
    forall a. Storable a => [a] -> Vector a
    fromList ([Double] -> Vector Double)
    -> ([Double] -> [Double]) -> [Double] -> Vector Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Double -> Double) -> [Double] -> [Double]
    forall a b. (a -> b) -> [a] -> [b]
    map Double -> Double
    rec ([Double] -> Vector t) -> [Double] -> Vector t
    forall a b. (a -> b) -> a -> b
    $ [Double]
    sl
        rec :: Double -> Double
    rec Double
    x = if Double
    x Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Double
    gDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    tol then Double
    x else Double
    1Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    /Double
    x
        tol :: Double
    tol = (Int -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    max Int
    r Int
    c) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    g Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    t Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    eps)
        r :: Int
    r = Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m
        c :: Int
    c = Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    m
        d :: Int
    d = Vector Double -> Int
    forall t. Storable t => Vector t -> Int
    dim Vector Double
    s
        u' :: Matrix t
    u' = Int -> Matrix t -> Matrix t
    forall t. Element t => Int -> Matrix t -> Matrix t
    takeColumns Int
    d Matrix t
    u
        v' :: Matrix t
    v' = Int -> Matrix t -> Matrix t
    forall t. Element t => Int -> Matrix t -> Matrix t
    takeColumns Int
    d Matrix t
    v

    -- | Numeric rank of a matrix from the SVD decomposition.
    rankSVD :: Element t
            => Double   -- ^ numeric zero (e.g. 1*'eps')
            -> Matrix t -- ^ input matrix m
            -> Vector Double -- ^ 'sv' of m
            -> Int      -- ^ rank of m
    rankSVD :: forall t. Element t => Double -> Matrix t -> Vector Double -> Int
    rankSVD Double
    teps Matrix t
    m Vector Double
    s = Double -> Int -> [Double] -> Int
    ranksv Double
    teps (Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    max (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m) (Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    m)) (Vector Double -> [Double]
    forall a. Storable a => Vector a -> [a]
    toList Vector Double
    s)

    -- | Numeric rank of a matrix from its singular values.
    ranksv ::  Double   -- ^ numeric zero (e.g. 1*'eps')
            -> Int      -- ^ maximum dimension of the matrix
            -> [Double] -- ^ singular values
            -> Int      -- ^ rank of m
    ranksv :: Double -> Int -> [Double] -> Int
    ranksv Double
    teps Int
    maxdim [Double]
    s = Int
    k where
        g :: Double
    g = [Double] -> Double
    forall a. Ord a => [a] -> a
    forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
    maximum [Double]
    s
        tol :: Double
    tol = Int -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral Int
    maxdim Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    g Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    teps
        s' :: [Double]
    s' = (Double -> Bool) -> [Double] -> [Double]
    forall a. (a -> Bool) -> [a] -> [a]
    filter (Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    >Double
    tol) [Double]
    s
        k :: Int
    k = if Double
    g Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    > Double
    teps then [Double] -> Int
    forall a. [a] -> Int
    forall (t :: * -> *) a. Foldable t => t a -> Int
    length [Double]
    s' else Int
    0

    -- | The machine precision of a Double: @eps = 2.22044604925031e-16@ (the value used by GNU-Octave).
    eps :: Double
    eps :: Double
    eps =  Double
    2.22044604925031e-16

    -- | 1 + 0.5*peps == 1,  1 + 0.6*peps /= 1
    peps :: RealFloat x => x
    peps :: forall x. RealFloat x => x
    peps = x
    x where x :: x
    x = x
    2.0 x -> x -> x
    forall a. Floating a => a -> a -> a
    ** Int -> x
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Int
    1 Int -> Int -> Int
    forall a. Num a => a -> a -> a
    - x -> Int
    forall a. RealFloat a => a -> Int
    floatDigits x
    x)

    -----------------------------------------------------------------------

    -- | The nullspace of a matrix from its precomputed SVD decomposition.
    nullspaceSVD :: Field t
                 => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
                                      --   or Right \"theoretical\" matrix rank.
                 -> Matrix t          -- ^ input matrix m
                 -> (Vector Double, Matrix t) -- ^ 'rightSV' of m
                 -> Matrix t          -- ^ nullspace
    nullspaceSVD :: forall t.
    Field t =>
    Either Double Int
    -> Matrix t -> (Vector Double, Matrix t) -> Matrix t
    nullspaceSVD Either Double Int
    hint Matrix t
    a (Vector Double
    s,Matrix t
    v) = Matrix t
    vs where
        tol :: Double
    tol = case Either Double Int
    hint of
            Left Double
    t -> Double
    t
            Either Double Int
    _      -> Double
    eps
        k :: Int
    k = case Either Double Int
    hint of
            Right Int
    t -> Int
    t
            Either Double Int
    _       -> Double -> Matrix t -> Vector Double -> Int
    forall t. Element t => Double -> Matrix t -> Vector Double -> Int
    rankSVD Double
    tol Matrix t
    a Vector Double
    s
        vs :: Matrix t
    vs = Int -> Matrix t -> Matrix t
    forall t. Element t => Int -> Matrix t -> Matrix t
    dropColumns Int
    k Matrix t
    v

    -- | The nullspace of a matrix. See also 'nullspaceSVD'.
    nullspacePrec :: Field t
                  => Double     -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps')
                  -> Matrix t   -- ^ input matrix
                  -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
    nullspacePrec :: forall t. Field t => Double -> Matrix t -> [Vector t]
    nullspacePrec Double
    t Matrix t
    m = Matrix t -> [Vector t]
    forall t. Element t => Matrix t -> [Vector t]
    toColumns (Matrix t -> [Vector t]) -> Matrix t -> [Vector t]
    forall a b. (a -> b) -> a -> b
    $ Either Double Int
    -> Matrix t -> (Vector Double, Matrix t) -> Matrix t
    forall t.
    Field t =>
    Either Double Int
    -> Matrix t -> (Vector Double, Matrix t) -> Matrix t
    nullspaceSVD (Double -> Either Double Int
    forall a b. a -> Either a b
    Left (Double
    tDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    eps)) Matrix t
    m (Matrix t -> (Vector Double, Matrix t)
    forall t. Field t => Matrix t -> (Vector Double, Matrix t)
    rightSV Matrix t
    m)

    -- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision.
    nullVector :: Field t => Matrix t -> Vector t
    nullVector :: forall t. Field t => Matrix t -> Vector t
    nullVector = [Vector t] -> Vector t
    forall a. HasCallStack => [a] -> a
    last ([Vector t] -> Vector t)
    -> (Matrix t -> [Vector t]) -> Matrix t -> Vector t
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Double -> Matrix t -> [Vector t]
    forall t. Field t => Double -> Matrix t -> [Vector t]
    nullspacePrec Double
    1

    -- | The range space a matrix from its precomputed SVD decomposition.
    orthSVD :: Field t
                 => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
                                      --   or Right \"theoretical\" matrix rank.
                 -> Matrix t          -- ^ input matrix m
                 -> (Matrix t, Vector Double) -- ^ 'leftSV' of m
                 -> Matrix t          -- ^ orth
    orthSVD :: forall t.
    Field t =>
    Either Double Int
    -> Matrix t -> (Matrix t, Vector Double) -> Matrix t
    orthSVD Either Double Int
    hint Matrix t
    a (Matrix t
    v,Vector Double
    s) = Matrix t
    vs where
        tol :: Double
    tol = case Either Double Int
    hint of
            Left Double
    t -> Double
    t
            Either Double Int
    _      -> Double
    eps
        k :: Int
    k = case Either Double Int
    hint of
            Right Int
    t -> Int
    t
            Either Double Int
    _       -> Double -> Matrix t -> Vector Double -> Int
    forall t. Element t => Double -> Matrix t -> Vector Double -> Int
    rankSVD Double
    tol Matrix t
    a Vector Double
    s
        vs :: Matrix t
    vs = Int -> Matrix t -> Matrix t
    forall t. Element t => Int -> Matrix t -> Matrix t
    takeColumns Int
    k Matrix t
    v

    orth :: Field t => Matrix t -> [Vector t]
    -- ^ Return an orthonormal basis of the range space of a matrix
    orth :: forall t. Field t => Matrix t -> [Vector t]
    orth Matrix t
    m = Int -> [Vector t] -> [Vector t]
    forall a. Int -> [a] -> [a]
    take Int
    r ([Vector t] -> [Vector t]) -> [Vector t] -> [Vector t]
    forall a b. (a -> b) -> a -> b
    $ Matrix t -> [Vector t]
    forall t. Element t => Matrix t -> [Vector t]
    toColumns Matrix t
    u
      where
        (Matrix t
    u,Vector Double
    s,Matrix t
    _) = Matrix t -> (Matrix t, Vector Double, Matrix t)
    forall t.
    Field t =>
    Matrix t -> (Matrix t, Vector Double, Matrix t)
    compactSVD Matrix t
    m
        r :: Int
    r = Double -> Int -> [Double] -> Int
    ranksv Double
    eps (Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    max (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m) (Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    m)) (Vector Double -> [Double]
    forall a. Storable a => Vector a -> [a]
    toList Vector Double
    s)

    ------------------------------------------------------------------------

    -- many thanks, quickcheck!

    haussholder :: (Field a) => a -> Vector a -> Matrix a
    haussholder :: forall a. Field a => a -> Vector a -> Matrix a
    haussholder a
    tau Vector a
    v = Int -> Matrix a
    forall a. (Num a, Element a) => Int -> Matrix a
    ident (Vector a -> Int
    forall t. Storable t => Vector t -> Int
    dim Vector a
    v) Matrix a -> Matrix a -> Matrix a
    forall (c :: * -> *) e. Container c e => c e -> c e -> c e
    `sub` (a
    tau a -> Matrix a -> Matrix a
    forall t (c :: * -> *). Linear t c => t -> c t -> c t
    `scale` (Matrix a
    w Matrix a -> Matrix a -> Matrix a
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    `mXm` Matrix a -> Matrix a
    forall t. CTrans t => Matrix t -> Matrix t
    ctrans Matrix a
    w))
        where w :: Matrix a
    w = Vector a -> Matrix a
    forall a. Storable a => Vector a -> Matrix a
    asColumn Vector a
    v

    zh :: Int -> Vector a -> Vector a
    zh Int
    k Vector a
    v = [a] -> Vector a
    forall a. Storable a => [a] -> Vector a
    fromList ([a] -> Vector a) -> [a] -> Vector a
    forall a b. (a -> b) -> a -> b
    $ Int -> a -> [a]
    forall a. Int -> a -> [a]
    replicate (Int
    kInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1) a
    0 [a] -> [a] -> [a]
    forall a. [a] -> [a] -> [a]
    ++ (a
    1a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    :Int -> [a] -> [a]
    forall a. Int -> [a] -> [a]
    drop Int
    k [a]
    xs)
                  where xs :: [a]
    xs = Vector a -> [a]
    forall a. Storable a => Vector a -> [a]
    toList Vector a
    v

    zt :: Int -> Vector t -> Vector t
    zt Int
    0 Vector t
    v = Vector t
    v
    zt Int
    k Vector t
    v = [Vector t] -> Vector t
    forall t. Storable t => [Vector t] -> Vector t
    vjoin [Int -> Int -> Vector t -> Vector t
    forall t. Storable t => Int -> Int -> Vector t -> Vector t
    subVector Int
    0 (Vector t -> Int
    forall t. Storable t => Vector t -> Int
    dim Vector t
    v Int -> Int -> Int
    forall a. Num a => a -> a -> a
    - Int
    k) Vector t
    v, t -> IndexOf Vector -> Vector t
    forall (c :: * -> *) e. Container c e => e -> IndexOf c -> c e
    konst' t
    0 Int
    IndexOf Vector
    k]

    unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
    unpackQR :: forall t. Field t => (Matrix t, Vector t) -> (Matrix t, Matrix t)
    unpackQR (Matrix t
    pq, Vector t
    tau) =  {-# SCC "unpackQR" #-} (Matrix t
    q,Matrix t
    r)
        where cs :: [Vector t]
    cs = Matrix t -> [Vector t]
    forall t. Element t => Matrix t -> [Vector t]
    toColumns Matrix t
    pq
              m :: Int
    m = Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    pq
              n :: Int
    n = Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    pq
              mn :: Int
    mn = Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    min Int
    m Int
    n
              r :: Matrix t
    r = [Vector t] -> Matrix t
    forall t. Element t => [Vector t] -> Matrix t
    fromColumns ([Vector t] -> Matrix t) -> [Vector t] -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ (Int -> Vector t -> Vector t) -> [Int] -> [Vector t] -> [Vector t]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith Int -> Vector t -> Vector t
    forall {t}.
    (Container Vector t, Num t) =>
    Int -> Vector t -> Vector t
    zt ([Int
    mInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1, Int
    mInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    2 .. Int
    1] [Int] -> [Int] -> [Int]
    forall a. [a] -> [a] -> [a]
    ++ Int -> [Int]
    forall a. a -> [a]
    repeat Int
    0) [Vector t]
    cs
              vs :: [Vector t]
    vs = (Int -> Vector t -> Vector t) -> [Int] -> [Vector t] -> [Vector t]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith Int -> Vector t -> Vector t
    forall {a}. (Storable a, Num a) => Int -> Vector a -> Vector a
    zh [Int
    1..Int
    mn] [Vector t]
    cs
              hs :: [Matrix t]
    hs = (t -> Vector t -> Matrix t) -> [t] -> [Vector t] -> [Matrix t]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith t -> Vector t -> Matrix t
    forall a. Field a => a -> Vector a -> Matrix a
    haussholder (Vector t -> [t]
    forall a. Storable a => Vector a -> [a]
    toList Vector t
    tau) [Vector t]
    vs
              q :: Matrix t
    q = (Matrix t -> Matrix t -> Matrix t) -> [Matrix t] -> Matrix t
    forall a. HasCallStack => (a -> a -> a) -> [a] -> a
    foldl1' Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    mXm [Matrix t]
    hs

    thinUnpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
    thinUnpackQR :: forall t. Field t => (Matrix t, Vector t) -> (Matrix t, Matrix t)
    thinUnpackQR (Matrix t
    pq, Vector t
    tau) = (Matrix t
    q, Matrix t
    r)
      where mn :: Int
    mn = (Int -> Int -> Int) -> (Int, Int) -> Int
    forall a b c. (a -> b -> c) -> (a, b) -> c
    uncurry Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    min ((Int, Int) -> Int) -> (Int, Int) -> Int
    forall a b. (a -> b) -> a -> b
    $ Matrix t -> (Int, Int)
    forall t. Matrix t -> (Int, Int)
    size Matrix t
    pq
            q :: Matrix t
    q = Int -> QR t -> Matrix t
    forall t. Field t => Int -> QR t -> Matrix t
    qrgr Int
    mn (QR t -> Matrix t) -> QR t -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ Matrix t -> Vector t -> QR t
    forall t. Matrix t -> Vector t -> QR t
    QR Matrix t
    pq Vector t
    tau
            r :: Matrix t
    r = [Vector t] -> Matrix t
    forall t. Element t => [Vector t] -> Matrix t
    fromRows ([Vector t] -> Matrix t) -> [Vector t] -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ (Int -> Vector t -> Vector t) -> [Int] -> [Vector t] -> [Vector t]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith (\Int
    i Vector t
    v -> Int -> t -> Vector t
    forall a. Storable a => Int -> a -> Vector a
    Vector.replicate Int
    i t
    0 Vector t -> Vector t -> Vector t
    forall a. Storable a => Vector a -> Vector a -> Vector a
    Vector.++ Int -> Vector t -> Vector t
    forall a. Storable a => Int -> Vector a -> Vector a
    Vector.drop Int
    i Vector t
    v) [Int
    0..Int
    mnInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1] (Matrix t -> [Vector t]
    forall t. Element t => Matrix t -> [Vector t]
    toRows Matrix t
    pq)

    unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t)
    unpackHess :: forall t.
    Field t =>
    (Matrix t -> (Matrix t, Vector t))
    -> Matrix t -> (Matrix t, Matrix t)
    unpackHess Matrix t -> (Matrix t, Vector t)
    hf Matrix t
    m
        | Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    == Int
    1 = ((Int
    1Int -> Int -> [t] -> Matrix t
    forall a. Storable a => Int -> Int -> [a] -> Matrix a
    ><Int
    1)[t
    1],Matrix t
    m)
        | Bool
    otherwise = ((Matrix t, Vector t) -> (Matrix t, Matrix t)
    forall t. Field t => (Matrix t, Vector t) -> (Matrix t, Matrix t)
    uH ((Matrix t, Vector t) -> (Matrix t, Matrix t))
    -> (Matrix t -> (Matrix t, Vector t))
    -> Matrix t
    -> (Matrix t, Matrix t)
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> (Matrix t, Vector t)
    hf) Matrix t
    m

    uH :: (Matrix t, Vector t) -> (Matrix t, Matrix t)
    uH (Matrix t
    pq, Vector t
    tau) = (Matrix t
    p,Matrix t
    h)
        where cs :: [Vector t]
    cs = Matrix t -> [Vector t]
    forall t. Element t => Matrix t -> [Vector t]
    toColumns Matrix t
    pq
              m :: Int
    m = Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    pq
              n :: Int
    n = Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    pq
              mn :: Int
    mn = Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    min Int
    m Int
    n
              h :: Matrix t
    h = [Vector t] -> Matrix t
    forall t. Element t => [Vector t] -> Matrix t
    fromColumns ([Vector t] -> Matrix t) -> [Vector t] -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ (Int -> Vector t -> Vector t) -> [Int] -> [Vector t] -> [Vector t]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith Int -> Vector t -> Vector t
    forall {t}.
    (Container Vector t, Num t) =>
    Int -> Vector t -> Vector t
    zt ([Int
    mInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    2, Int
    mInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    3 .. Int
    1] [Int] -> [Int] -> [Int]
    forall a. [a] -> [a] -> [a]
    ++ Int -> [Int]
    forall a. a -> [a]
    repeat Int
    0) [Vector t]
    cs
              vs :: [Vector t]
    vs = (Int -> Vector t -> Vector t) -> [Int] -> [Vector t] -> [Vector t]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith Int -> Vector t -> Vector t
    forall {a}. (Storable a, Num a) => Int -> Vector a -> Vector a
    zh [Int
    2..Int
    mn] [Vector t]
    cs
              hs :: [Matrix t]
    hs = (t -> Vector t -> Matrix t) -> [t] -> [Vector t] -> [Matrix t]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith t -> Vector t -> Matrix t
    forall a. Field a => a -> Vector a -> Matrix a
    haussholder (Vector t -> [t]
    forall a. Storable a => Vector a -> [a]
    toList Vector t
    tau) [Vector t]
    vs
              p :: Matrix t
    p = (Matrix t -> Matrix t -> Matrix t) -> [Matrix t] -> Matrix t
    forall a. HasCallStack => (a -> a -> a) -> [a] -> a
    foldl1' Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    mXm [Matrix t]
    hs

    --------------------------------------------------------------------------

    -- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.
    rcond :: Field t => Matrix t -> Double
    rcond :: forall t. Field t => Matrix t -> Double
    rcond Matrix t
    m = [Double] -> Double
    forall a. HasCallStack => [a] -> a
    last [Double]
    s Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / [Double] -> Double
    forall a. HasCallStack => [a] -> a
    head [Double]
    s
        where s :: [Double]
    s = Vector Double -> [Double]
    forall a. Storable a => Vector a -> [a]
    toList (Matrix t -> Vector Double
    forall t. Field t => Matrix t -> Vector Double
    singularValues Matrix t
    m)

    -- | Number of linearly independent rows or columns. See also 'ranksv'
    rank :: Field t => Matrix t -> Int
    rank :: forall t. Field t => Matrix t -> Int
    rank Matrix t
    m = Double -> Matrix t -> Vector Double -> Int
    forall t. Element t => Double -> Matrix t -> Vector Double -> Int
    rankSVD Double
    eps Matrix t
    m (Matrix t -> Vector Double
    forall t. Field t => Matrix t -> Vector Double
    singularValues Matrix t
    m)

    {-
    expm' m = case diagonalize (complex m) of
        Just (l,v) -> v `mXm` diag (exp l) `mXm` inv v
        Nothing -> error "Sorry, expm not yet implemented for non-diagonalizable matrices"
      where exp = vectorMapC Exp
    -}

    diagonalize :: Matrix (Complex Double)
    -> Maybe (Vector (Complex Double), Matrix (Complex Double))
    diagonalize Matrix (Complex Double)
    m = if Matrix (Complex Double) -> Int
    forall t. Field t => Matrix t -> Int
    rank Matrix (Complex Double)
    v Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    == Int
    n
                        then (Vector (Complex Double), Matrix (Complex Double))
    -> Maybe (Vector (Complex Double), Matrix (Complex Double))
    forall a. a -> Maybe a
    Just (Vector (Complex Double)
    l,Matrix (Complex Double)
    v)
                        else Maybe (Vector (Complex Double), Matrix (Complex Double))
    forall a. Maybe a
    Nothing
        where n :: Int
    n = Matrix (Complex Double) -> Int
    forall t. Matrix t -> Int
    rows Matrix (Complex Double)
    m
              (Vector (Complex Double)
    l,Matrix (Complex Double)
    v) = if Matrix (Complex Double) -> Bool
    forall {e}.
    (Num e, Container Vector e, CTrans e) =>
    Matrix e -> Bool
    exactHermitian Matrix (Complex Double)
    m
                        then let (Vector Double
    l',Matrix (Complex Double)
    v') = Herm (Complex Double) -> (Vector Double, Matrix (Complex Double))
    forall t. Field t => Herm t -> (Vector Double, Matrix t)
    eigSH (Matrix (Complex Double) -> Herm (Complex Double)
    forall t. Matrix t -> Herm t
    trustSym Matrix (Complex Double)
    m) in (Vector (RealOf (Complex Double)) -> Vector (Complex Double)
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c (RealOf t) -> c t
    forall (c :: * -> *).
    Complexable c =>
    c (RealOf (Complex Double)) -> c (Complex Double)
    real Vector Double
    Vector (RealOf (Complex Double))
    l', Matrix (Complex Double)
    v')
                        else Matrix (Complex Double)
    -> (Vector (Complex Double), Matrix (Complex Double))
    forall t.
    Field t =>
    Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
    eig Matrix (Complex Double)
    m

    -- | Generic matrix functions for diagonalizable matrices. For instance:
    --
    -- @logm = matFunc log@
    --
    matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
    matFunc :: (Complex Double -> Complex Double)
    -> Matrix (Complex Double) -> Matrix (Complex Double)
    matFunc Complex Double -> Complex Double
    f Matrix (Complex Double)
    m = case Matrix (Complex Double)
    -> Maybe (Vector (Complex Double), Matrix (Complex Double))
    diagonalize Matrix (Complex Double)
    m of
        Just (Vector (Complex Double)
    l,Matrix (Complex Double)
    v) -> Matrix (Complex Double)
    v Matrix (Complex Double)
    -> Matrix (Complex Double) -> Matrix (Complex Double)
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    `mXm` Vector (Complex Double) -> Matrix (Complex Double)
    forall a. (Num a, Element a) => Vector a -> Matrix a
    diag ((Complex Double -> Complex Double)
    -> Vector (Complex Double) -> Vector (Complex Double)
    forall a b.
    (Storable a, Storable b) =>
    (a -> b) -> Vector a -> Vector b
    mapVector Complex Double -> Complex Double
    f Vector (Complex Double)
    l) Matrix (Complex Double)
    -> Matrix (Complex Double) -> Matrix (Complex Double)
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    `mXm` Matrix (Complex Double) -> Matrix (Complex Double)
    forall t. Field t => Matrix t -> Matrix t
    inv Matrix (Complex Double)
    v
        Maybe (Vector (Complex Double), Matrix (Complex Double))
    Nothing -> String -> Matrix (Complex Double)
    forall a. HasCallStack => String -> a
    error String
    "Sorry, matFunc requires a diagonalizable matrix"

    --------------------------------------------------------------

    golubeps :: Integer -> Integer -> Double
    golubeps :: Integer -> Integer -> Double
    golubeps Integer
    p Integer
    q = Double
    a Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Integer -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral Integer
    b Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Integer -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral Integer
    c where
        a :: Double
    a = Double
    2Double -> Integer -> Double
    forall a b. (Fractional a, Integral b) => a -> b -> a
    ^^(Integer
    3Integer -> Integer -> Integer
    forall a. Num a => a -> a -> a
    -Integer
    pInteger -> Integer -> Integer
    forall a. Num a => a -> a -> a
    -Integer
    q)
        b :: Integer
    b = Integer -> Integer
    forall {a}. (Num a, Enum a) => a -> a
    fact Integer
    p Integer -> Integer -> Integer
    forall a. Num a => a -> a -> a
    * Integer -> Integer
    forall {a}. (Num a, Enum a) => a -> a
    fact Integer
    q
        c :: Integer
    c = Integer -> Integer
    forall {a}. (Num a, Enum a) => a -> a
    fact (Integer
    pInteger -> Integer -> Integer
    forall a. Num a => a -> a -> a
    +Integer
    q) Integer -> Integer -> Integer
    forall a. Num a => a -> a -> a
    * Integer -> Integer
    forall {a}. (Num a, Enum a) => a -> a
    fact (Integer
    pInteger -> Integer -> Integer
    forall a. Num a => a -> a -> a
    +Integer
    qInteger -> Integer -> Integer
    forall a. Num a => a -> a -> a
    +Integer
    1)
        fact :: a -> a
    fact a
    n = [a] -> a
    forall a. Num a => [a] -> a
    forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
    product [a
    1..a
    n]

    epslist :: [(Int,Double)]
    epslist :: [(Int, Double)]
    epslist = [ (Integer -> Int
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral Integer
    k, Integer -> Integer -> Double
    golubeps Integer
    k Integer
    k) | Integer
    k <- [Integer
    1..]]

    geps :: Double -> Int
    geps Double
    delta = [Int] -> Int
    forall a. HasCallStack => [a] -> a
    head [ Int
    k | (Int
    k,Double
    g) <- [(Int, Double)]
    epslist, Double
    gDouble -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    <Double
    delta]

    {- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,
         based on a scaled Pade approximation.
    -}
    expm :: Field t => Matrix t -> Matrix t
    expm :: forall t. Field t => Matrix t -> Matrix t
    expm = Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t
    expGolub

    expGolub :: Field t => Matrix t -> Matrix t
    expGolub :: forall t. Field t => Matrix t -> Matrix t
    expGolub Matrix t
    m = (Matrix t -> Matrix t) -> Matrix t -> [Matrix t]
    forall a. (a -> a) -> a -> [a]
    iterate Matrix t -> Matrix t
    msq Matrix t
    f [Matrix t] -> Int -> Matrix t
    forall a. HasCallStack => [a] -> Int -> a
    !! Int
    j
        where j :: Int
    j = Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    max Int
    0 (Int -> Int) -> Int -> Int
    forall a b. (a -> b) -> a -> b
    $ Double -> Int
    forall b. Integral b => Double -> b
    forall a b. (RealFrac a, Integral b) => a -> b
    floor (Double -> Int) -> Double -> Int
    forall a b. (a -> b) -> a -> b
    $ Double -> Double -> Double
    forall a. Floating a => a -> a -> a
    logBase Double
    2 (Double -> Double) -> Double -> Double
    forall a b. (a -> b) -> a -> b
    $ NormType -> Matrix t -> RealOf t
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    Infinity Matrix t
    m
              a :: Matrix t
    a = Matrix t
    m Matrix t -> t -> Matrix t
    forall {t} {c :: * -> *}.
    (Linear t c, Fractional t) =>
    c t -> t -> c t
    */ Int -> t
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral ((Int
    2::Int)Int -> Int -> Int
    forall a b. (Num a, Integral b) => a -> b -> a
    ^Int
    j)
              q :: Int
    q = Double -> Int
    geps Double
    eps -- 7 steps
              eye :: Matrix t
    eye = Int -> Matrix t
    forall a. (Num a, Element a) => Int -> Matrix a
    ident (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m)
              work :: (Int, t, Matrix t, Matrix t, Matrix t)
    -> (Int, t, Matrix t, Matrix t, Matrix t)
    work (Int
    k,t
    c,Matrix t
    x,Matrix t
    n,Matrix t
    d) = (Int
    k',t
    c',Matrix t
    x',Matrix t
    n',Matrix t
    d')
                  where k' :: Int
    k' = Int
    kInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1
                        c' :: t
    c' = t
    c t -> t -> t
    forall a. Num a => a -> a -> a
    * Int -> t
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Int
    qInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    kInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1) t -> t -> t
    forall a. Fractional a => a -> a -> a
    / Int -> t
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral ((Int
    2Int -> Int -> Int
    forall a. Num a => a -> a -> a
    *Int
    qInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    kInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1)Int -> Int -> Int
    forall a. Num a => a -> a -> a
    *Int
    k)
                        x' :: Matrix t
    x' = Matrix t
    a Matrix t -> Matrix t -> Matrix t
    <> Matrix t
    x
                        n' :: Matrix t
    n' = Matrix t
    n Matrix t -> Matrix t -> Matrix t
    |+| (t
    c' t -> Matrix t -> Matrix t
    .* Matrix t
    x')
                        d' :: Matrix t
    d' = Matrix t
    d Matrix t -> Matrix t -> Matrix t
    |+| (((-t
    1)t -> Int -> t
    forall a b. (Num a, Integral b) => a -> b -> a
    ^Int
    k t -> t -> t
    forall a. Num a => a -> a -> a
    * t
    c') t -> Matrix t -> Matrix t
    .* Matrix t
    x')
              (Int
    _,t
    _,Matrix t
    _,Matrix t
    nf,Matrix t
    df) = ((Int, t, Matrix t, Matrix t, Matrix t)
     -> (Int, t, Matrix t, Matrix t, Matrix t))
    -> (Int, t, Matrix t, Matrix t, Matrix t)
    -> [(Int, t, Matrix t, Matrix t, Matrix t)]
    forall a. (a -> a) -> a -> [a]
    iterate (Int, t, Matrix t, Matrix t, Matrix t)
    -> (Int, t, Matrix t, Matrix t, Matrix t)
    work (Int
    1,t
    1,Matrix t
    eye,Matrix t
    eye,Matrix t
    eye) [(Int, t, Matrix t, Matrix t, Matrix t)]
    -> Int -> (Int, t, Matrix t, Matrix t, Matrix t)
    forall a. HasCallStack => [a] -> Int -> a
    !! Int
    q
              f :: Matrix t
    f = Matrix t -> Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t -> Matrix t
    linearSolve Matrix t
    df Matrix t
    nf
              msq :: Matrix t -> Matrix t
    msq Matrix t
    x = Matrix t
    x Matrix t -> Matrix t -> Matrix t
    <> Matrix t
    x

              <> :: Matrix t -> Matrix t -> Matrix t
    (<>) = Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    multiply
              c t
    v */ :: c t -> t -> c t
    */ t
    x = t -> c t -> c t
    forall t (c :: * -> *). Linear t c => t -> c t -> c t
    scale (t -> t
    forall a. Fractional a => a -> a
    recip t
    x) c t
    v
              .* :: t -> Matrix t -> Matrix t
    (.*) = t -> Matrix t -> Matrix t
    forall t (c :: * -> *). Linear t c => t -> c t -> c t
    scale
              |+| :: Matrix t -> Matrix t -> Matrix t
    (|+|) = Matrix t -> Matrix t -> Matrix t
    forall c. Additive c => c -> c -> c
    add

    --------------------------------------------------------------

    {- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia.
    It only works with invertible matrices that have a real solution.

    @m = (2><2) [4,9
               ,0,4] :: Matrix Double@

    >>> sqrtm m
    (2><2)
     [ 2.0, 2.25
     , 0.0,  2.0 ]

    For diagonalizable matrices you can try 'matFunc' @sqrt@:

    >>> matFunc sqrt ((2><2) [1,0,0,-1])
    (2><2)
     [ 1.0 :+ 0.0, 0.0 :+ 0.0
     , 0.0 :+ 0.0, 0.0 :+ 1.0 ]

    -}
    sqrtm ::  Field t => Matrix t -> Matrix t
    sqrtm :: forall t. Field t => Matrix t -> Matrix t
    sqrtm = Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t
    sqrtmInv

    sqrtmInv :: Matrix t -> Matrix t
    sqrtmInv Matrix t
    x = (Matrix t, Matrix t) -> Matrix t
    forall a b. (a, b) -> a
    fst ((Matrix t, Matrix t) -> Matrix t)
    -> (Matrix t, Matrix t) -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ [(Matrix t, Matrix t)] -> (Matrix t, Matrix t)
    fixedPoint ([(Matrix t, Matrix t)] -> (Matrix t, Matrix t))
    -> [(Matrix t, Matrix t)] -> (Matrix t, Matrix t)
    forall a b. (a -> b) -> a -> b
    $ ((Matrix t, Matrix t) -> (Matrix t, Matrix t))
    -> (Matrix t, Matrix t) -> [(Matrix t, Matrix t)]
    forall a. (a -> a) -> a -> [a]
    iterate (Matrix t, Matrix t) -> (Matrix t, Matrix t)
    f (Matrix t
    x, Int -> Matrix t
    forall a. (Num a, Element a) => Int -> Matrix a
    ident (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    x))
        where fixedPoint :: [(Matrix t, Matrix t)] -> (Matrix t, Matrix t)
    fixedPoint ((Matrix t, Matrix t)
    a:(Matrix t, Matrix t)
    b:[(Matrix t, Matrix t)]
    rest) | NormType -> Matrix t -> RealOf t
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1 ((Matrix t, Matrix t) -> Matrix t
    forall a b. (a, b) -> a
    fst (Matrix t, Matrix t)
    a Matrix t -> Matrix t -> Matrix t
    |-| (Matrix t, Matrix t) -> Matrix t
    forall a b. (a, b) -> a
    fst (Matrix t, Matrix t)
    b) RealOf t -> RealOf t -> Bool
    forall a. Ord a => a -> a -> Bool
    < RealOf t
    forall x. RealFloat x => x
    peps   = (Matrix t, Matrix t)
    a
                                    | Bool
    otherwise = [(Matrix t, Matrix t)] -> (Matrix t, Matrix t)
    fixedPoint ((Matrix t, Matrix t)
    b(Matrix t, Matrix t)
    -> [(Matrix t, Matrix t)] -> [(Matrix t, Matrix t)]
    forall a. a -> [a] -> [a]
    :[(Matrix t, Matrix t)]
    rest)
              fixedPoint [(Matrix t, Matrix t)]
    _ = String -> (Matrix t, Matrix t)
    forall a. HasCallStack => String -> a
    error String
    "fixedpoint with impossible inputs"
              f :: (Matrix t, Matrix t) -> (Matrix t, Matrix t)
    f (Matrix t
    y,Matrix t
    z) = (t
    0.5 t -> Matrix t -> Matrix t
    .* (Matrix t
    y Matrix t -> Matrix t -> Matrix t
    |+| Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t
    inv Matrix t
    z),
                         t
    0.5 t -> Matrix t -> Matrix t
    .* (Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t
    inv Matrix t
    y Matrix t -> Matrix t -> Matrix t
    |+| Matrix t
    z))
              .* :: t -> Matrix t -> Matrix t
    (.*) = t -> Matrix t -> Matrix t
    forall t (c :: * -> *). Linear t c => t -> c t -> c t
    scale
              |+| :: Matrix t -> Matrix t -> Matrix t
    (|+|) = Matrix t -> Matrix t -> Matrix t
    forall c. Additive c => c -> c -> c
    add
              |-| :: Matrix t -> Matrix t -> Matrix t
    (|-|) = Matrix t -> Matrix t -> Matrix t
    forall (c :: * -> *) e. Container c e => c e -> c e -> c e
    sub

    ------------------------------------------------------------------

    signlp :: a -> [a] -> b
    signlp a
    r [a]
    vals = (b -> (a, a) -> b) -> b -> [(a, a)] -> b
    forall b a. (b -> a -> b) -> b -> [a] -> b
    forall (t :: * -> *) b a.
    Foldable t =>
    (b -> a -> b) -> b -> t a -> b
    foldl b -> (a, a) -> b
    forall {a} {a}. (Eq a, Num a) => a -> (a, a) -> a
    f b
    1 ([a] -> [a] -> [(a, a)]
    forall a b. [a] -> [b] -> [(a, b)]
    zip [a
    0..a
    ra -> a -> a
    forall a. Num a => a -> a -> a
    -a
    1] [a]
    vals)
        where f :: a -> (a, a) -> a
    f a
    s (a
    a,a
    b) | a
    a a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    /= a
    b    = -a
    s
                        | Bool
    otherwise =  a
    s

    fixPerm :: Int -> [Int] -> (Matrix t, b)
    fixPerm Int
    r [Int]
    vals = ([Vector t] -> Matrix t
    forall t. Element t => [Vector t] -> Matrix t
    fromColumns ([Vector t] -> Matrix t) -> [Vector t] -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ Array Int (Vector t) -> [Vector t]
    forall i e. Array i e -> [e]
    A.elems Array Int (Vector t)
    res, b
    sign)
      where
        v :: [Int]
    v = [Int
    0..Int
    rInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1]
        t :: [Vector t]
    t = Matrix t -> [Vector t]
    forall t. Element t => Matrix t -> [Vector t]
    toColumns (Int -> Matrix t
    forall a. (Num a, Element a) => Int -> Matrix a
    ident Int
    r)
        (Array Int (Vector t)
    res,b
    sign) = ((Array Int (Vector t), b)
     -> (Int, Int) -> (Array Int (Vector t), b))
    -> (Array Int (Vector t), b)
    -> [(Int, Int)]
    -> (Array Int (Vector t), b)
    forall b a. (b -> a -> b) -> b -> [a] -> b
    forall (t :: * -> *) b a.
    Foldable t =>
    (b -> a -> b) -> b -> t a -> b
    foldl (Array Int (Vector t), b)
    -> (Int, Int) -> (Array Int (Vector t), b)
    forall {i} {b} {e}.
    (Ix i, Num b) =>
    (Array i e, b) -> (i, i) -> (Array i e, b)
    swap ((Int, Int) -> [Vector t] -> Array Int (Vector t)
    forall i e. Ix i => (i, i) -> [e] -> Array i e
    A.listArray (Int
    0,Int
    rInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1) [Vector t]
    t, b
    1) ([Int] -> [Int] -> [(Int, Int)]
    forall a b. [a] -> [b] -> [(a, b)]
    zip [Int]
    v [Int]
    vals)
        swap :: (Array i e, b) -> (i, i) -> (Array i e, b)
    swap (Array i e
    arr,b
    s) (i
    a,i
    b)
          | i
    a i -> i -> Bool
    forall a. Eq a => a -> a -> Bool
    /= i
    b    = (Array i e
    arr Array i e -> [(i, e)] -> Array i e
    forall i e. Ix i => Array i e -> [(i, e)] -> Array i e
    A.// [(i
    a, Array i e
    arr Array i e -> i -> e
    forall i e. Ix i => Array i e -> i -> e
    A.! i
    b),(i
    b,Array i e
    arr Array i e -> i -> e
    forall i e. Ix i => Array i e -> i -> e
    A.! i
    a)],-b
    s)
          | Bool
    otherwise = (Array i e
    arr,b
    s)

    fixPerm' :: [Int] -> Vector I
    fixPerm' :: [Int] -> Vector I
    fixPerm' [Int]
    s = (Matrix I, ()) -> Vector I
    forall {b}. (Matrix I, b) -> Vector I
    res ((Matrix I, ()) -> Vector I) -> (Matrix I, ()) -> Vector I
    forall a b. (a -> b) -> a -> b
    $ (forall s. (Int, Int) -> STMatrix s I -> ST s ())
    -> Matrix I -> (Matrix I, ())
    forall t u.
    Element t =>
    (forall s. (Int, Int) -> STMatrix s t -> ST s u)
    -> Matrix t -> (Matrix t, u)
    mutable (Int, Int) -> STMatrix s I -> ST s ()
    forall s. (Int, Int) -> STMatrix s I -> ST s ()
    forall t s.
    (Num t, Element t) =>
    (Int, Int) -> STMatrix s t -> ST s ()
    f Matrix I
    s0
      where
        s0 :: Matrix I
    s0 = Int -> Vector I -> Matrix I
    forall t. Storable t => Int -> Vector t -> Matrix t
    reshape Int
    1 (Int -> Vector I
    range ([Int] -> Int
    forall a. [a] -> Int
    forall (t :: * -> *) a. Foldable t => t a -> Int
    length [Int]
    s))
        res :: (Matrix I, b) -> Vector I
    res = Matrix I -> Vector I
    forall t. Element t => Matrix t -> Vector t
    flatten (Matrix I -> Vector I)
    -> ((Matrix I, b) -> Matrix I) -> (Matrix I, b) -> Vector I
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Matrix I, b) -> Matrix I
    forall a b. (a, b) -> a
    fst
        swap :: STMatrix s t -> Int -> Int -> ST s ()
    swap STMatrix s t
    m Int
    i Int
    j = RowOper t -> STMatrix s t -> ST s ()
    forall t s.
    (Num t, Element t) =>
    RowOper t -> STMatrix s t -> ST s ()
    rowOper (Int -> Int -> ColRange -> RowOper t
    forall t. Int -> Int -> ColRange -> RowOper t
    SWAP Int
    i Int
    j ColRange
    AllCols) STMatrix s t
    m
        f :: (Num t, Element t) => (Int, Int) -> STMatrix s t -> ST s () -- needed because of TypeFamilies
        f :: forall t s.
    (Num t, Element t) =>
    (Int, Int) -> STMatrix s t -> ST s ()
    f (Int, Int)
    _ STMatrix s t
    p = [ST s ()] -> ST s ()
    forall (t :: * -> *) (m :: * -> *) a.
    (Foldable t, Monad m) =>
    t (m a) -> m ()
    sequence_ ([ST s ()] -> ST s ()) -> [ST s ()] -> ST s ()
    forall a b. (a -> b) -> a -> b
    $ (Int -> Int -> ST s ()) -> [Int] -> [Int] -> [ST s ()]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith (STMatrix s t -> Int -> Int -> ST s ()
    forall {t} {s}.
    (Num t, Element t) =>
    STMatrix s t -> Int -> Int -> ST s ()
    swap STMatrix s t
    p) [Int
    0..] [Int]
    s

    triang :: Int -> Int -> Int -> a -> Matrix a
    triang Int
    r Int
    c Int
    h a
    v = (Int
    rInt -> Int -> [a] -> Matrix a
    forall a. Storable a => Int -> Int -> [a] -> Matrix a
    ><Int
    c) [Int -> Int -> a
    el Int
    s Int
    t | Int
    s<-[Int
    0..Int
    rInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1], Int
    t<-[Int
    0..Int
    cInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1]]
        where el :: Int -> Int -> a
    el Int
    p Int
    q = if Int
    qInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    pInt -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    >=Int
    h then a
    v else a
    1 a -> a -> a
    forall a. Num a => a -> a -> a
    - a
    v

    -- | Compute the explicit LU decomposition from the compact one obtained by 'luPacked'.
    luFact :: Numeric t => LU t -> (Matrix t, Matrix t, Matrix t, t)
    luFact :: forall t. Numeric t => LU t -> (Matrix t, Matrix t, Matrix t, t)
    luFact (LU Matrix t
    l_u [Int]
    perm)
        | Int
    r Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Int
    c    = (Matrix t
    l ,Matrix t
    u ,Matrix t
    p, t
    s)
        | Bool
    otherwise = (Matrix t
    l',Matrix t
    u',Matrix t
    p, t
    s)
      where
        r :: Int
    r = Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    l_u
        c :: Int
    c = Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    l_u
        tu :: Matrix t
    tu = Int -> Int -> Int -> t -> Matrix t
    forall {a}.
    (Num a, Storable a) =>
    Int -> Int -> Int -> a -> Matrix a
    triang Int
    r Int
    c Int
    0 t
    1
        tl :: Matrix t
    tl = Int -> Int -> Int -> t -> Matrix t
    forall {a}.
    (Num a, Storable a) =>
    Int -> Int -> Int -> a -> Matrix a
    triang Int
    r Int
    c Int
    0 t
    0
        l :: Matrix t
    l = Int -> Matrix t -> Matrix t
    forall t. Element t => Int -> Matrix t -> Matrix t
    takeColumns Int
    r (Matrix t
    l_u Matrix t -> Matrix t -> Matrix t
    |*| Matrix t
    tl) Matrix t -> Matrix t -> Matrix t
    |+| t -> Vector t -> Int -> Int -> Matrix t
    forall t. Storable t => t -> Vector t -> Int -> Int -> Matrix t
    diagRect t
    0 (t -> IndexOf Vector -> Vector t
    forall (c :: * -> *) e. Container c e => e -> IndexOf c -> c e
    konst' t
    1 Int
    IndexOf Vector
    r) Int
    r Int
    r
        u :: Matrix t
    u = Matrix t
    l_u Matrix t -> Matrix t -> Matrix t
    |*| Matrix t
    tu
        (Matrix t
    p,t
    s) = Int -> [Int] -> (Matrix t, t)
    forall {t} {b}.
    (Num t, Num b, Element t) =>
    Int -> [Int] -> (Matrix t, b)
    fixPerm Int
    r [Int]
    perm
        l' :: Matrix t
    l' = (Matrix t
    l_u Matrix t -> Matrix t -> Matrix t
    |*| Matrix t
    tl) Matrix t -> Matrix t -> Matrix t
    |+| t -> Vector t -> Int -> Int -> Matrix t
    forall t. Storable t => t -> Vector t -> Int -> Int -> Matrix t
    diagRect t
    0 (t -> IndexOf Vector -> Vector t
    forall (c :: * -> *) e. Container c e => e -> IndexOf c -> c e
    konst' t
    1 Int
    IndexOf Vector
    c) Int
    r Int
    c
        u' :: Matrix t
    u' = Int -> Matrix t -> Matrix t
    forall t. Element t => Int -> Matrix t -> Matrix t
    takeRows Int
    c (Matrix t
    l_u Matrix t -> Matrix t -> Matrix t
    |*| Matrix t
    tu)
        |+| :: Matrix t -> Matrix t -> Matrix t
    (|+|) = Matrix t -> Matrix t -> Matrix t
    forall c. Additive c => c -> c -> c
    add
        |*| :: Matrix t -> Matrix t -> Matrix t
    (|*|) = Matrix t -> Matrix t -> Matrix t
    forall (c :: * -> *) e. Container c e => c e -> c e -> c e
    mul

    ---------------------------------------------------------------------------

    data NormType = Infinity | PNorm1 | PNorm2 | Frobenius

    class (RealFloat (RealOf t)) => Normed c t where
        pnorm :: NormType -> c t -> RealOf t

    instance Normed Vector Double where
        pnorm :: NormType -> Vector Double -> RealOf Double
    pnorm NormType
    PNorm1    = Vector Double -> RealOf Double
    forall e. Product e => Vector e -> RealOf e
    norm1
        pnorm NormType
    PNorm2    = Vector Double -> RealOf Double
    forall e. (Product e, Floating e) => Vector e -> RealOf e
    norm2
        pnorm NormType
    Infinity  = Vector Double -> RealOf Double
    forall e. Product e => Vector e -> RealOf e
    normInf
        pnorm NormType
    Frobenius = Vector Double -> RealOf Double
    forall e. (Product e, Floating e) => Vector e -> RealOf e
    norm2

    instance Normed Vector (Complex Double) where
        pnorm :: NormType -> Vector (Complex Double) -> RealOf (Complex Double)
    pnorm NormType
    PNorm1    = Vector (Complex Double) -> RealOf (Complex Double)
    forall e. Product e => Vector e -> RealOf e
    norm1
        pnorm NormType
    PNorm2    = Vector (Complex Double) -> RealOf (Complex Double)
    forall e. (Product e, Floating e) => Vector e -> RealOf e
    norm2
        pnorm NormType
    Infinity  = Vector (Complex Double) -> RealOf (Complex Double)
    forall e. Product e => Vector e -> RealOf e
    normInf
        pnorm NormType
    Frobenius = NormType -> Vector (Complex Double) -> RealOf (Complex Double)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm2

    instance Normed Vector Float where
        pnorm :: NormType -> Vector Float -> RealOf Float
    pnorm NormType
    PNorm1    = Vector Float -> RealOf Float
    forall e. Product e => Vector e -> RealOf e
    norm1
        pnorm NormType
    PNorm2    = Vector Float -> RealOf Float
    forall e. (Product e, Floating e) => Vector e -> RealOf e
    norm2
        pnorm NormType
    Infinity  = Vector Float -> RealOf Float
    forall e. Product e => Vector e -> RealOf e
    normInf
        pnorm NormType
    Frobenius = NormType -> Vector Float -> RealOf Float
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm2

    instance Normed Vector (Complex Float) where
        pnorm :: NormType -> Vector (Complex Float) -> RealOf (Complex Float)
    pnorm NormType
    PNorm1    = Vector (Complex Float) -> RealOf (Complex Float)
    forall e. Product e => Vector e -> RealOf e
    norm1
        pnorm NormType
    PNorm2    = Vector (Complex Float) -> RealOf (Complex Float)
    forall e. (Product e, Floating e) => Vector e -> RealOf e
    norm2
        pnorm NormType
    Infinity  = Vector (Complex Float) -> RealOf (Complex Float)
    forall e. Product e => Vector e -> RealOf e
    normInf
        pnorm NormType
    Frobenius = NormType -> Vector (Complex Float) -> RealOf (Complex Float)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm2

    instance Normed Matrix Double where
        pnorm :: NormType -> Matrix Double -> RealOf Double
    pnorm NormType
    PNorm1    = [Double] -> Double
    forall a. Ord a => [a] -> a
    forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
    maximum ([Double] -> Double)
    -> (Matrix Double -> [Double]) -> Matrix Double -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Vector Double -> Double) -> [Vector Double] -> [Double]
    forall a b. (a -> b) -> [a] -> [b]
    map (NormType -> Vector Double -> RealOf Double
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1) ([Vector Double] -> [Double])
    -> (Matrix Double -> [Vector Double]) -> Matrix Double -> [Double]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Double -> [Vector Double]
    forall t. Element t => Matrix t -> [Vector t]
    toColumns
        pnorm NormType
    PNorm2    = (Vector Double -> Int -> Double
    forall t. Storable t => Vector t -> Int -> t
    @>Int
    0) (Vector Double -> Double)
    -> (Matrix Double -> Vector Double) -> Matrix Double -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Double -> Vector Double
    forall t. Field t => Matrix t -> Vector Double
    singularValues
        pnorm NormType
    Infinity  = NormType -> Matrix Double -> RealOf Double
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1 (Matrix Double -> Double)
    -> (Matrix Double -> Matrix Double) -> Matrix Double -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Double -> Matrix Double
    forall t. Matrix t -> Matrix t
    trans
        pnorm NormType
    Frobenius = NormType -> Vector Double -> RealOf Double
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm2 (Vector Double -> Double)
    -> (Matrix Double -> Vector Double) -> Matrix Double -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Double -> Vector Double
    forall t. Element t => Matrix t -> Vector t
    flatten

    instance Normed Matrix (Complex Double) where
        pnorm :: NormType -> Matrix (Complex Double) -> RealOf (Complex Double)
    pnorm NormType
    PNorm1    = [Double] -> Double
    forall a. Ord a => [a] -> a
    forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
    maximum ([Double] -> Double)
    -> (Matrix (Complex Double) -> [Double])
    -> Matrix (Complex Double)
    -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Vector (Complex Double) -> Double)
    -> [Vector (Complex Double)] -> [Double]
    forall a b. (a -> b) -> [a] -> [b]
    map (NormType -> Vector (Complex Double) -> RealOf (Complex Double)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1) ([Vector (Complex Double)] -> [Double])
    -> (Matrix (Complex Double) -> [Vector (Complex Double)])
    -> Matrix (Complex Double)
    -> [Double]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix (Complex Double) -> [Vector (Complex Double)]
    forall t. Element t => Matrix t -> [Vector t]
    toColumns
        pnorm NormType
    PNorm2    = (Vector Double -> Int -> Double
    forall t. Storable t => Vector t -> Int -> t
    @>Int
    0) (Vector Double -> Double)
    -> (Matrix (Complex Double) -> Vector Double)
    -> Matrix (Complex Double)
    -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix (Complex Double) -> Vector Double
    forall t. Field t => Matrix t -> Vector Double
    singularValues
        pnorm NormType
    Infinity  = NormType -> Matrix (Complex Double) -> RealOf (Complex Double)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1 (Matrix (Complex Double) -> Double)
    -> (Matrix (Complex Double) -> Matrix (Complex Double))
    -> Matrix (Complex Double)
    -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix (Complex Double) -> Matrix (Complex Double)
    forall t. Matrix t -> Matrix t
    trans
        pnorm NormType
    Frobenius = NormType -> Vector (Complex Double) -> RealOf (Complex Double)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm2 (Vector (Complex Double) -> Double)
    -> (Matrix (Complex Double) -> Vector (Complex Double))
    -> Matrix (Complex Double)
    -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix (Complex Double) -> Vector (Complex Double)
    forall t. Element t => Matrix t -> Vector t
    flatten

    instance Normed Matrix Float where
        pnorm :: NormType -> Matrix Float -> RealOf Float
    pnorm NormType
    PNorm1    = [Float] -> Float
    forall a. Ord a => [a] -> a
    forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
    maximum ([Float] -> Float)
    -> (Matrix Float -> [Float]) -> Matrix Float -> Float
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Vector Float -> Float) -> [Vector Float] -> [Float]
    forall a b. (a -> b) -> [a] -> [b]
    map (NormType -> Vector Float -> RealOf Float
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1) ([Vector Float] -> [Float])
    -> (Matrix Float -> [Vector Float]) -> Matrix Float -> [Float]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Float -> [Vector Float]
    forall t. Element t => Matrix t -> [Vector t]
    toColumns
        pnorm NormType
    PNorm2    = Double -> Float
    forall a b. (Real a, Fractional b) => a -> b
    realToFrac (Double -> Float)
    -> (Matrix Float -> Double) -> Matrix Float -> Float
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Vector Double -> Int -> Double
    forall t. Storable t => Vector t -> Int -> t
    @>Int
    0) (Vector Double -> Double)
    -> (Matrix Float -> Vector Double) -> Matrix Float -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Double -> Vector Double
    forall t. Field t => Matrix t -> Vector Double
    singularValues (Matrix Double -> Vector Double)
    -> (Matrix Float -> Matrix Double) -> Matrix Float -> Vector Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Float -> Matrix Double
    Matrix Float -> Matrix (DoubleOf Float)
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c t -> c (DoubleOf t)
    forall (c :: * -> *).
    Complexable c =>
    c Float -> c (DoubleOf Float)
    double
        pnorm NormType
    Infinity  = NormType -> Matrix Float -> RealOf Float
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1 (Matrix Float -> Float)
    -> (Matrix Float -> Matrix Float) -> Matrix Float -> Float
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Float -> Matrix Float
    forall t. Matrix t -> Matrix t
    trans
        pnorm NormType
    Frobenius = NormType -> Vector Float -> RealOf Float
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm2 (Vector Float -> Float)
    -> (Matrix Float -> Vector Float) -> Matrix Float -> Float
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Float -> Vector Float
    forall t. Element t => Matrix t -> Vector t
    flatten

    instance Normed Matrix (Complex Float) where
        pnorm :: NormType -> Matrix (Complex Float) -> RealOf (Complex Float)
    pnorm NormType
    PNorm1    = [Float] -> Float
    forall a. Ord a => [a] -> a
    forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
    maximum ([Float] -> Float)
    -> (Matrix (Complex Float) -> [Float])
    -> Matrix (Complex Float)
    -> Float
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Vector (Complex Float) -> Float)
    -> [Vector (Complex Float)] -> [Float]
    forall a b. (a -> b) -> [a] -> [b]
    map (NormType -> Vector (Complex Float) -> RealOf (Complex Float)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1) ([Vector (Complex Float)] -> [Float])
    -> (Matrix (Complex Float) -> [Vector (Complex Float)])
    -> Matrix (Complex Float)
    -> [Float]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix (Complex Float) -> [Vector (Complex Float)]
    forall t. Element t => Matrix t -> [Vector t]
    toColumns
        pnorm NormType
    PNorm2    = Double -> Float
    forall a b. (Real a, Fractional b) => a -> b
    realToFrac (Double -> Float)
    -> (Matrix (Complex Float) -> Double)
    -> Matrix (Complex Float)
    -> Float
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Vector Double -> Int -> Double
    forall t. Storable t => Vector t -> Int -> t
    @>Int
    0) (Vector Double -> Double)
    -> (Matrix (Complex Float) -> Vector Double)
    -> Matrix (Complex Float)
    -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix (Complex Double) -> Vector Double
    forall t. Field t => Matrix t -> Vector Double
    singularValues (Matrix (Complex Double) -> Vector Double)
    -> (Matrix (Complex Float) -> Matrix (Complex Double))
    -> Matrix (Complex Float)
    -> Vector Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix (Complex Float) -> Matrix (Complex Double)
    Matrix (Complex Float) -> Matrix (DoubleOf (Complex Float))
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c t -> c (DoubleOf t)
    forall (c :: * -> *).
    Complexable c =>
    c (Complex Float) -> c (DoubleOf (Complex Float))
    double
        pnorm NormType
    Infinity  = NormType -> Matrix (Complex Float) -> RealOf (Complex Float)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1 (Matrix (Complex Float) -> Float)
    -> (Matrix (Complex Float) -> Matrix (Complex Float))
    -> Matrix (Complex Float)
    -> Float
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix (Complex Float) -> Matrix (Complex Float)
    forall t. Matrix t -> Matrix t
    trans
        pnorm NormType
    Frobenius = NormType -> Vector (Complex Float) -> RealOf (Complex Float)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm2 (Vector (Complex Float) -> Float)
    -> (Matrix (Complex Float) -> Vector (Complex Float))
    -> Matrix (Complex Float)
    -> Float
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix (Complex Float) -> Vector (Complex Float)
    forall t. Element t => Matrix t -> Vector t
    flatten

    -- | Approximate number of common digits in the maximum element.
    relativeError' :: (Normed c t, Container c t) => c t -> c t -> Int
    relativeError' :: forall (c :: * -> *) t.
    (Normed c t, Container c t) =>
    c t -> c t -> Int
    relativeError' c t
    x c t
    y = RealOf t -> Int
    forall {b} {a}. (Integral b, Real a) => a -> b
    dig (c t -> RealOf t
    norm (c t
    x c t -> c t -> c t
    forall (c :: * -> *) e. Container c e => c e -> c e -> c e
    `sub` c t
    y) RealOf t -> RealOf t -> RealOf t
    forall a. Fractional a => a -> a -> a
    / c t -> RealOf t
    norm c t
    x)
        where norm :: c t -> RealOf t
    norm = NormType -> c t -> RealOf t
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    Infinity
              dig :: a -> b
    dig a
    r = Double -> b
    forall b. Integral b => Double -> b
    forall a b. (RealFrac a, Integral b) => a -> b
    round (Double -> b) -> Double -> b
    forall a b. (a -> b) -> a -> b
    $ -Double -> Double -> Double
    forall a. Floating a => a -> a -> a
    logBase Double
    10 (a -> Double
    forall a b. (Real a, Fractional b) => a -> b
    realToFrac a
    r :: Double)

    relativeError :: Num a => (a -> Double) -> a -> a -> Double
    relativeError :: forall a. Num a => (a -> Double) -> a -> a -> Double
    relativeError a -> Double
    norm a
    a a
    b = Double
    r
      where
        na :: Double
    na = a -> Double
    norm a
    a
        nb :: Double
    nb = a -> Double
    norm a
    b
        nab :: Double
    nab = a -> Double
    norm (a
    aa -> a -> a
    forall a. Num a => a -> a -> a
    -a
    b)
        mx :: Double
    mx = Double -> Double -> Double
    forall a. Ord a => a -> a -> a
    max Double
    na Double
    nb
        mn :: Double
    mn = Double -> Double -> Double
    forall a. Ord a => a -> a -> a
    min Double
    na Double
    nb
        r :: Double
    r = if Double
    mn Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    < Double
    forall x. RealFloat x => x
    peps
            then Double
    mx
            else Double
    nabDouble -> Double -> Double
    forall a. Fractional a => a -> a -> a
    /Double
    mx

    ----------------------------------------------------------------------

    -- | Generalized symmetric positive definite eigensystem Av = lBv,
    -- for A and B symmetric, B positive definite.
    geigSH :: Field t
            => Herm t -- ^ A
            -> Herm t -- ^ B
            -> (Vector Double, Matrix t)
    geigSH :: forall t. Field t => Herm t -> Herm t -> (Vector Double, Matrix t)
    geigSH (Herm Matrix t
    a) (Herm Matrix t
    b) = Matrix t -> Matrix t -> (Vector Double, Matrix t)
    forall t.
    Field t =>
    Matrix t -> Matrix t -> (Vector Double, Matrix t)
    geigSH' Matrix t
    a Matrix t
    b

    geigSH' :: Field t
            => Matrix t -- ^ A
            -> Matrix t -- ^ B
            -> (Vector Double, Matrix t)
    geigSH' :: forall t.
    Field t =>
    Matrix t -> Matrix t -> (Vector Double, Matrix t)
    geigSH' Matrix t
    a Matrix t
    b = (Vector Double
    l,Matrix t
    v')
      where
        u :: Matrix t
    u = Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t
    cholSH Matrix t
    b
        iu :: Matrix t
    iu = Matrix t -> Matrix t
    forall t. Field t => Matrix t -> Matrix t
    inv Matrix t
    u
        c :: Matrix t
    c = Matrix t -> Matrix t
    forall t. CTrans t => Matrix t -> Matrix t
    ctrans Matrix t
    iu Matrix t -> Matrix t -> Matrix t
    <> Matrix t
    a Matrix t -> Matrix t -> Matrix t
    <> Matrix t
    iu
        (Vector Double
    l,Matrix t
    v) = Matrix t -> (Vector Double, Matrix t)
    forall t. Field t => Matrix t -> (Vector Double, Matrix t)
    eigSH' Matrix t
    c
        v' :: Matrix t
    v' = Matrix t
    iu Matrix t -> Matrix t -> Matrix t
    <> Matrix t
    v
        <> :: Matrix t -> Matrix t -> Matrix t
    (<>) = Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    mXm

    --------------------------------------------------------------------------------

    -- | A matrix that, by construction, it is known to be complex Hermitian or real symmetric.
    --
    --   It can be created using 'sym', 'mTm', or 'trustSym', and the matrix can be extracted using 'unSym'.
    newtype Herm t = Herm (Matrix t) deriving Int -> Herm t -> ShowS
    [Herm t] -> ShowS
    Herm t -> String
    (Int -> Herm t -> ShowS)
    -> (Herm t -> String) -> ([Herm t] -> ShowS) -> Show (Herm t)
    forall t. (Show t, Element t) => Int -> Herm t -> ShowS
    forall t. (Show t, Element t) => [Herm t] -> ShowS
    forall t. (Show t, Element t) => Herm t -> String
    forall a.
    (Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
    $cshowsPrec :: forall t. (Show t, Element t) => Int -> Herm t -> ShowS
    showsPrec :: Int -> Herm t -> ShowS
    $cshow :: forall t. (Show t, Element t) => Herm t -> String
    show :: Herm t -> String
    $cshowList :: forall t. (Show t, Element t) => [Herm t] -> ShowS
    showList :: [Herm t] -> ShowS
    Show

    instance (NFData t, Numeric t) => NFData (Herm t)
      where
        rnf :: Herm t -> ()
    rnf (Herm Matrix t
    m) = Matrix t -> ()
    forall a. NFData a => a -> ()
    rnf Matrix t
    m

    -- | Extract the general matrix from a 'Herm' structure, forgetting its symmetric or Hermitian property.
    unSym :: Herm t -> Matrix t
    unSym :: forall t. Herm t -> Matrix t
    unSym (Herm Matrix t
    x) = Matrix t
    x

    -- | Compute the complex Hermitian or real symmetric part of a square matrix (@(x + tr x)/2@).
    sym :: Field t => Matrix t -> Herm t
    sym :: forall t. Field t => Matrix t -> Herm t
    sym Matrix t
    x = Matrix t -> Herm t
    forall t. Matrix t -> Herm t
    Herm (t -> Matrix t -> Matrix t
    forall t (c :: * -> *). Linear t c => t -> c t -> c t
    scale t
    0.5 (Matrix t -> Matrix t
    forall m mt. Transposable m mt => m -> mt
    tr Matrix t
    x Matrix t -> Matrix t -> Matrix t
    forall c. Additive c => c -> c -> c
    `add` Matrix t
    x))

    -- | Compute the contraction @tr x <> x@ of a general matrix.
    mTm :: Numeric t => Matrix t -> Herm t
    mTm :: forall t. Numeric t => Matrix t -> Herm t
    mTm Matrix t
    x = Matrix t -> Herm t
    forall t. Matrix t -> Herm t
    Herm (Matrix t -> Matrix t
    forall m mt. Transposable m mt => m -> mt
    tr Matrix t
    x Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    `mXm` Matrix t
    x)

    instance Field t => Linear t Herm where
        scale :: t -> Herm t -> Herm t
    scale  t
    x (Herm Matrix t
    m) = Matrix t -> Herm t
    forall t. Matrix t -> Herm t
    Herm (t -> Matrix t -> Matrix t
    forall t (c :: * -> *). Linear t c => t -> c t -> c t
    scale t
    x Matrix t
    m)

    instance Field t => Additive (Herm t) where
        add :: Herm t -> Herm t -> Herm t
    add (Herm Matrix t
    a) (Herm Matrix t
    b) = Matrix t -> Herm t
    forall t. Matrix t -> Herm t
    Herm (Matrix t
    a Matrix t -> Matrix t -> Matrix t
    forall c. Additive c => c -> c -> c
    `add` Matrix t
    b)

    -- | At your own risk, declare that a matrix is complex Hermitian or real symmetric
    --   for usage in 'chol', 'eigSH', etc. Only a triangular part of the matrix will be used.
    trustSym :: Matrix t -> Herm t
    trustSym :: forall t. Matrix t -> Herm t
    trustSym Matrix t
    x = (Matrix t -> Herm t
    forall t. Matrix t -> Herm t
    Herm Matrix t
    x)