{-# 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 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' = 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'= 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 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' = 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' = 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 = forall t. Matrix t -> Int
rows Matrix t
m forall a. Eq a => a -> a -> Bool
== forall t. Matrix t -> Int
cols Matrix t
m

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

exactHermitian :: Matrix e -> Bool
exactHermitian Matrix e
m = Matrix e
m forall (c :: * -> *) e. Container c e => c e -> c e -> Bool
`equal` 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" #-} forall {m} {c} {a} {b}. Transposable m c => (a, b, m) -> (a, b, c)
g forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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,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" #-} forall {m} {c} {a} {b}. Transposable m c => (a, b, m) -> (a, b, c)
g forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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,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" #-} 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) = forall t.
Field t =>
Matrix t -> (Matrix t, Vector Double, Matrix t)
svd Matrix t
m
    d :: Matrix Double
d = forall t. Storable t => t -> Vector t -> Int -> Int -> Matrix t
diagRect Double
0 Vector Double
s Int
r Int
c
    r :: Int
r = forall t. Matrix t -> Int
rows Matrix t
m
    c :: Int
c = 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 = 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', 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) = forall t.
Field t =>
Matrix t -> (Matrix t, Vector Double, Matrix t)
thinSVD Matrix t
m
    d :: Int
d = forall t. Element t => Double -> Matrix t -> Vector Double -> Int
rankSVD (Double
rforall a. Num a => a -> a -> a
*Double
eps) Matrix t
m Vector Double
s forall a. Ord a => a -> a -> a
`max` Int
1
    u' :: Matrix t
u' = forall t. Element t => Int -> Matrix t -> Matrix t
takeColumns Int
d Matrix t
u
    v' :: Matrix t
v' = 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 | forall {t}. Matrix t -> Bool
vertical Matrix t
m = let (Matrix t
_,Vector Double
s,Matrix t
v) = 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) = 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  | forall {t}. Matrix t -> Bool
vertical Matrix t
m = let (Matrix t
u,Vector Double
s,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
_) = 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
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
showList :: [LU t] -> ShowS
$cshowList :: forall t. (Show t, Element t) => [LU t] -> ShowS
show :: LU t -> String
$cshow :: forall t. (Show t, Element t) => LU t -> String
showsPrec :: Int -> LU t -> ShowS
$cshowsPrec :: forall t. (Show t, Element t) => Int -> LU t -> ShowS
Show

instance (NFData t, Numeric t) => NFData (LU t)
  where
    rnf :: LU t -> ()
rnf (LU Matrix t
m [Int]
_) = 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" #-} forall t. Matrix t -> [Int] -> LU t
LU Matrix t
m [Int]
p
  where
    (Matrix t
m,[Int]
p) = 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" #-} 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" #-} 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" #-} 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" #-} 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" #-} 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" #-} 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" #-} 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" #-} 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
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
showList :: [LDL t] -> ShowS
$cshowList :: forall t. (Show t, Element t) => [LDL t] -> ShowS
show :: LDL t -> String
$cshow :: forall t. (Show t, Element t) => LDL t -> String
showsPrec :: Int -> LDL t -> ShowS
$cshowsPrec :: forall t. (Show t, Element t) => Int -> LDL t -> ShowS
Show

instance (NFData t, Numeric t) => NFData (LDL t)
  where
    rnf :: LDL t -> ()
rnf (LDL Matrix t
m [Int]
_) = 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" #-} forall t. Matrix t -> [Int] -> LDL t
LDL Matrix t
m [Int]
p
  where
   (Matrix t
m,[Int]
p) = 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) = 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" #-} 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" #-} 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" #-} 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" #-} 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" #-} 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'" #-} 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'" #-} 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) = 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) = 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
_) = 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" #-} forall t. Field t => (Matrix t, Vector t) -> (Matrix t, Matrix t)
unpackQR forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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" #-} forall t. Field t => (Matrix t, Vector t) -> (Matrix t, Matrix t)
thinUnpackQR forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 = forall t. Matrix t -> Vector t -> QR t
QR Matrix t
x Vector t
v
  where
    (Matrix t
x,Vector t
v) = 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)
    | forall t. Storable t => Vector t -> Int
dim Vector t
t forall a. Ord a => a -> a -> Bool
> forall a. Ord a => a -> a -> a
min (forall t. Matrix t -> Int
cols Matrix t
a) (forall t. Matrix t -> Int
rows Matrix t
a) Bool -> Bool -> Bool
|| Int
n forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
n forall a. Ord a => a -> a -> Bool
> forall t. Storable t => Vector t -> Int
dim Vector t
t = forall a. HasCallStack => String -> a
error String
"qrgr expects k <= min(rows,cols)"
    | Bool
otherwise = 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" #-} forall t.
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)
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" #-} forall t.
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)
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 forall a b. (a -> b) -> a -> b
$ forall t. Matrix t -> Matrix t
trans forall a b. (a -> b) -> a -> b
$ Matrix t -> Matrix t
rev1 Matrix t
m
    r :: Matrix t
r = Matrix t -> Matrix t
rev2 (forall t. Matrix t -> Matrix t
trans Matrix t
r')
    q :: Matrix t
q = Matrix t -> Matrix t
rev2 (forall t. Matrix t -> Matrix t
trans Matrix t
q')
    rev1 :: Matrix t -> Matrix t
rev1 = forall t. Element t => Matrix t -> Matrix t
flipud forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> Matrix t
fliprl
    rev2 :: Matrix t -> Matrix t
rev2 = forall t. Element t => Matrix t -> Matrix t
fliprl forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 = 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 = 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" #-} 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 = 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" #-} 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" #-} 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 | forall {t}. Matrix t -> Bool
square Matrix t
m = (Matrix t
im,(t
ladm,t
sdm))
           | Bool
otherwise = forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$ String
"invlndet of nonsquare "forall a. [a] -> [a] -> [a]
++ forall t. Matrix t -> String
shSize Matrix t
m forall a. [a] -> [a] -> [a]
++ String
" matrix"
  where
    lp :: LU t
lp@(LU Matrix t
lup [Int]
perm) = forall t. Field t => Matrix t -> LU t
luPacked Matrix t
m
    s :: t
s = forall {a} {b}. (Eq a, Num b, Num a, Enum a) => a -> [a] -> b
signlp (forall t. Matrix t -> Int
rows Matrix t
m) [Int]
perm
    dg :: [t]
dg = forall a. Storable a => Vector a -> [a]
toList forall a b. (a -> b) -> a -> b
$ forall t. Element t => Matrix t -> Vector t
takeDiag forall a b. (a -> b) -> a -> b
$ Matrix t
lup
    ladm :: t
ladm = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a. Floating a => a -> a
logforall b c a. (b -> c) -> (a -> b) -> a -> c
.forall a. Num a => a -> a
abs) [t]
dg
    sdm :: t
sdm = t
sforall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (forall a b. (a -> b) -> [a] -> [b]
map forall a. Num a => a -> a
signum [t]
dg)
    im :: Matrix t
im = forall t. Field t => LU t -> Matrix t -> Matrix t
luSolve LU t
lp (forall a. (Num a, Element a) => Int -> Matrix a
ident (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 | forall {t}. Matrix t -> Bool
square Matrix t
m = {-# SCC "det" #-} t
s forall a. Num a => a -> a -> a
* (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product forall a b. (a -> b) -> a -> b
$ forall a. Storable a => Vector a -> [a]
toList forall a b. (a -> b) -> a -> b
$ forall t. Element t => Matrix t -> Vector t
takeDiag forall a b. (a -> b) -> a -> b
$ Matrix t
lup)
      | Bool
otherwise = forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$ String
"det of nonsquare "forall a. [a] -> [a] -> [a]
++ forall t. Matrix t -> String
shSize Matrix t
m forall a. [a] -> [a] -> [a]
++ String
" matrix"
    where
      LU Matrix t
lup [Int]
perm = forall t. Field t => Matrix t -> LU t
luPacked Matrix t
m
      s :: t
s = forall {a} {b}. (Eq a, Num b, Num a, Enum a) => a -> [a] -> b
signlp (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 = forall t. Numeric t => LU t -> (Matrix t, Matrix t, Matrix t, t)
luFact forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 | forall {t}. Matrix t -> Bool
square Matrix t
m = Matrix t
m forall t. Field t => Matrix t -> Matrix t -> Matrix t
`linearSolve` forall a. (Num a, Element a) => Int -> Matrix a
ident (forall t. Matrix t -> Int
rows Matrix t
m)
      | Bool
otherwise = forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$ String
"inv of nonsquare "forall a. [a] -> [a] -> [a]
++ forall t. Matrix t -> String
shSize Matrix t
m 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 = 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' forall t. Product t => Matrix t -> Matrix t -> Matrix t
`mXm` forall a. (Num a, Element a) => Vector a -> Matrix a
diag Vector t
s' forall t. Product t => Matrix t -> Matrix t -> Matrix t
`mXm` forall t. CTrans t => Matrix t -> Matrix t
ctrans Matrix t
u' where
    (Matrix t
u,Vector Double
s,Matrix t
v) = forall t.
Field t =>
Matrix t -> (Matrix t, Vector Double, Matrix t)
thinSVD Matrix t
m
    sl :: [Double]
sl@(Double
g:[Double]
_) = forall a. Storable a => Vector a -> [a]
toList Vector Double
s
    s' :: Vector t
s' = forall t (c :: * -> *).
(Convert t, Complexable c) =>
c (RealOf t) -> c t
real forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Storable a => [a] -> Vector a
fromList forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
rec forall a b. (a -> b) -> a -> b
$ [Double]
sl
    rec :: Double -> Double
rec Double
x = if Double
x forall a. Ord a => a -> a -> Bool
<= Double
gforall a. Num a => a -> a -> a
*Double
tol then Double
x else Double
1forall a. Fractional a => a -> a -> a
/Double
x
    tol :: Double
tol = (forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall a. Ord a => a -> a -> a
max Int
r Int
c) forall a. Num a => a -> a -> a
* Double
g forall a. Num a => a -> a -> a
* Double
t forall a. Num a => a -> a -> a
* Double
eps)
    r :: Int
r = forall t. Matrix t -> Int
rows Matrix t
m
    c :: Int
c = forall t. Matrix t -> Int
cols Matrix t
m
    d :: Int
d = forall t. Storable t => Vector t -> Int
dim Vector Double
s
    u' :: Matrix t
u' = forall t. Element t => Int -> Matrix t -> Matrix t
takeColumns Int
d Matrix t
u
    v' :: Matrix t
v' = 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 (forall a. Ord a => a -> a -> a
max (forall t. Matrix t -> Int
rows Matrix t
m) (forall t. Matrix t -> Int
cols Matrix t
m)) (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 = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
s
    tol :: Double
tol = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
maxdim forall a. Num a => a -> a -> a
* Double
g forall a. Num a => a -> a -> a
* Double
teps
    s' :: [Double]
s' = forall a. (a -> Bool) -> [a] -> [a]
filter (forall a. Ord a => a -> a -> Bool
>Double
tol) [Double]
s
    k :: Int
k = if Double
g forall a. Ord a => a -> a -> Bool
> Double
teps then 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 forall a. Floating a => a -> a -> a
** forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
1 forall a. Num a => a -> a -> a
- 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
_       -> forall t. Element t => Double -> Matrix t -> Vector Double -> Int
rankSVD Double
tol Matrix t
a Vector Double
s
    vs :: Matrix t
vs = 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 = forall t. Element t => Matrix t -> [Vector t]
toColumns forall a b. (a -> b) -> a -> b
$ forall t.
Field t =>
Either Double Int
-> Matrix t -> (Vector Double, Matrix t) -> Matrix t
nullspaceSVD (forall a b. a -> Either a b
Left (Double
tforall a. Num a => a -> a -> a
*Double
eps)) Matrix t
m (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 = forall a. [a] -> a
last forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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
_       -> forall t. Element t => Double -> Matrix t -> Vector Double -> Int
rankSVD Double
tol Matrix t
a Vector Double
s
    vs :: Matrix t
vs = 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 = forall a. Int -> [a] -> [a]
take Int
r forall a b. (a -> b) -> a -> b
$ forall t. Element t => Matrix t -> [Vector t]
toColumns Matrix t
u
  where
    (Matrix t
u,Vector Double
s,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 (forall a. Ord a => a -> a -> a
max (forall t. Matrix t -> Int
rows Matrix t
m) (forall t. Matrix t -> Int
cols Matrix t
m)) (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 = forall a. (Num a, Element a) => Int -> Matrix a
ident (forall t. Storable t => Vector t -> Int
dim Vector a
v) forall (c :: * -> *) e. Container c e => c e -> c e -> c e
`sub` (a
tau forall t (c :: * -> *). Linear t c => t -> c t -> c t
`scale` (Matrix a
w forall t. Product t => Matrix t -> Matrix t -> Matrix t
`mXm` forall t. CTrans t => Matrix t -> Matrix t
ctrans Matrix a
w))
    where w :: Matrix a
w = forall a. Storable a => Vector a -> Matrix a
asColumn Vector a
v


zh :: Int -> Vector a -> Vector a
zh Int
k Vector a
v = forall a. Storable a => [a] -> Vector a
fromList forall a b. (a -> b) -> a -> b
$ forall a. Int -> a -> [a]
replicate (Int
kforall a. Num a => a -> a -> a
-Int
1) a
0 forall a. [a] -> [a] -> [a]
++ (a
1forall a. a -> [a] -> [a]
:forall a. Int -> [a] -> [a]
drop Int
k [a]
xs)
              where xs :: [a]
xs = 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 = forall t. Storable t => [Vector t] -> Vector t
vjoin [forall t. Storable t => Int -> Int -> Vector t -> Vector t
subVector Int
0 (forall t. Storable t => Vector t -> Int
dim Vector t
v forall a. Num a => a -> a -> a
- Int
k) Vector t
v, forall (c :: * -> *) e. Container c e => e -> IndexOf c -> c e
konst' t
0 Int
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 = forall t. Element t => Matrix t -> [Vector t]
toColumns Matrix t
pq
          m :: Int
m = forall t. Matrix t -> Int
rows Matrix t
pq
          n :: Int
n = forall t. Matrix t -> Int
cols Matrix t
pq
          mn :: Int
mn = forall a. Ord a => a -> a -> a
min Int
m Int
n
          r :: Matrix t
r = forall t. Element t => [Vector t] -> Matrix t
fromColumns forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall {t}.
(Container Vector t, Num t) =>
Int -> Vector t -> Vector t
zt ([Int
mforall a. Num a => a -> a -> a
-Int
1, Int
mforall a. Num a => a -> a -> a
-Int
2 .. Int
1] forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat Int
0) [Vector t]
cs
          vs :: [Vector t]
vs = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall {a}. (Storable a, Num a) => Int -> Vector a -> Vector a
zh [Int
1..Int
mn] [Vector t]
cs
          hs :: [Matrix t]
hs = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Field a => a -> Vector a -> Matrix a
haussholder (forall a. Storable a => Vector a -> [a]
toList Vector t
tau) [Vector t]
vs
          q :: Matrix t
q = forall a. (a -> a -> a) -> [a] -> a
foldl1' 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 = forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry forall a. Ord a => a -> a -> a
min forall a b. (a -> b) -> a -> b
$ forall t. Matrix t -> (Int, Int)
size Matrix t
pq
        q :: Matrix t
q = forall t. Field t => Int -> QR t -> Matrix t
qrgr Int
mn forall a b. (a -> b) -> a -> b
$ forall t. Matrix t -> Vector t -> QR t
QR Matrix t
pq Vector t
tau
        r :: Matrix t
r = forall t. Element t => [Vector t] -> Matrix t
fromRows forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
i Vector t
v -> forall a. Storable a => Int -> a -> Vector a
Vector.replicate Int
i t
0 forall a. Storable a => Vector a -> Vector a -> Vector a
Vector.++ forall a. Storable a => Int -> Vector a -> Vector a
Vector.drop Int
i Vector t
v) [Int
0..Int
mnforall a. Num a => a -> a -> a
-Int
1] (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
    | forall t. Matrix t -> Int
rows Matrix t
m forall a. Eq a => a -> a -> Bool
== Int
1 = ((Int
1forall a. Storable a => Int -> Int -> [a] -> Matrix a
><Int
1)[t
1],Matrix t
m)
    | Bool
otherwise = (forall t. Field t => (Matrix t, Vector t) -> (Matrix t, Matrix t)
uH 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 = forall t. Element t => Matrix t -> [Vector t]
toColumns Matrix t
pq
          m :: Int
m = forall t. Matrix t -> Int
rows Matrix t
pq
          n :: Int
n = forall t. Matrix t -> Int
cols Matrix t
pq
          mn :: Int
mn = forall a. Ord a => a -> a -> a
min Int
m Int
n
          h :: Matrix t
h = forall t. Element t => [Vector t] -> Matrix t
fromColumns forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall {t}.
(Container Vector t, Num t) =>
Int -> Vector t -> Vector t
zt ([Int
mforall a. Num a => a -> a -> a
-Int
2, Int
mforall a. Num a => a -> a -> a
-Int
3 .. Int
1] forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat Int
0) [Vector t]
cs
          vs :: [Vector t]
vs = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall {a}. (Storable a, Num a) => Int -> Vector a -> Vector a
zh [Int
2..Int
mn] [Vector t]
cs
          hs :: [Matrix t]
hs = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Field a => a -> Vector a -> Matrix a
haussholder (forall a. Storable a => Vector a -> [a]
toList Vector t
tau) [Vector t]
vs
          p :: Matrix t
p = forall a. (a -> a -> a) -> [a] -> a
foldl1' 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 = forall a. [a] -> a
last [Double]
s forall a. Fractional a => a -> a -> a
/ forall a. [a] -> a
head [Double]
s
    where s :: [Double]
s = forall a. Storable a => Vector a -> [a]
toList (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 = forall t. Element t => Double -> Matrix t -> Vector Double -> Int
rankSVD Double
eps Matrix t
m (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 forall t. Field t => Matrix t -> Int
rank Matrix (Complex Double)
v forall a. Eq a => a -> a -> Bool
== Int
n
                    then forall a. a -> Maybe a
Just (Vector (Complex Double)
l,Matrix (Complex Double)
v)
                    else forall a. Maybe a
Nothing
    where n :: Int
n = forall t. Matrix t -> Int
rows Matrix (Complex Double)
m
          (Vector (Complex Double)
l,Matrix (Complex Double)
v) = if 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') = forall t. Field t => Herm t -> (Vector Double, Matrix t)
eigSH (forall t. Matrix t -> Herm t
trustSym Matrix (Complex Double)
m) in (forall t (c :: * -> *).
(Convert t, Complexable c) =>
c (RealOf t) -> c t
real Vector Double
l', Matrix (Complex Double)
v')
                    else 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 forall t. Product t => Matrix t -> Matrix t -> Matrix t
`mXm` forall a. (Num a, Element a) => Vector a -> Matrix a
diag (forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
mapVector Complex Double -> Complex Double
f Vector (Complex Double)
l) forall t. Product t => Matrix t -> Matrix t -> Matrix t
`mXm` forall t. Field t => Matrix t -> Matrix t
inv Matrix (Complex Double)
v
    Maybe (Vector (Complex Double), Matrix (Complex Double))
Nothing -> 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 forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
b forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
c where
    a :: Double
a = Double
2forall a b. (Fractional a, Integral b) => a -> b -> a
^^(Integer
3forall a. Num a => a -> a -> a
-Integer
pforall a. Num a => a -> a -> a
-Integer
q)
    b :: Integer
b = forall {a}. (Num a, Enum a) => a -> a
fact Integer
p forall a. Num a => a -> a -> a
* forall {a}. (Num a, Enum a) => a -> a
fact Integer
q
    c :: Integer
c = forall {a}. (Num a, Enum a) => a -> a
fact (Integer
pforall a. Num a => a -> a -> a
+Integer
q) forall a. Num a => a -> a -> a
* forall {a}. (Num a, Enum a) => a -> a
fact (Integer
pforall a. Num a => a -> a -> a
+Integer
qforall a. Num a => a -> a -> a
+Integer
1)
    fact :: a -> a
fact a
n = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [a
1..a
n]

epslist :: [(Int,Double)]
epslist :: [(Int, Double)]
epslist = [ (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 = forall a. [a] -> a
head [ Int
k | (Int
k,Double
g) <- [(Int, Double)]
epslist, Double
gforall 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 = 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 = forall a. (a -> a) -> a -> [a]
iterate Matrix t -> Matrix t
msq Matrix t
f forall a. [a] -> Int -> a
!! Int
j
    where j :: Int
j = forall a. Ord a => a -> a -> a
max Int
0 forall a b. (a -> b) -> a -> b
$ forall a b. (RealFrac a, Integral b) => a -> b
floor forall a b. (a -> b) -> a -> b
$ forall a. Floating a => a -> a -> a
logBase Double
2 forall a b. (a -> b) -> a -> b
$ forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
Infinity Matrix t
m
          a :: Matrix t
a = Matrix t
m forall {t} {c :: * -> *}.
(Linear t c, Fractional t) =>
c t -> t -> c t
*/ forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int
2::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 = forall a. (Num a, Element a) => Int -> Matrix a
ident (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
kforall a. Num a => a -> a -> a
+Int
1
                    c' :: t
c' = t
c forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
qforall a. Num a => a -> a -> a
-Int
kforall a. Num a => a -> a -> a
+Int
1) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int
2forall a. Num a => a -> a -> a
*Int
qforall a. Num a => a -> a -> a
-Int
kforall a. Num a => a -> a -> a
+Int
1)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)forall a b. (Num a, Integral b) => a -> b -> a
^Int
k 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) = 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) forall a. [a] -> Int -> a
!! Int
q
          f :: Matrix t
f = 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
(<>) = forall t. Product t => Matrix t -> Matrix t -> Matrix t
multiply
          c t
v */ :: c t -> t -> c t
*/ t
x = forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale (forall a. Fractional a => a -> a
recip t
x) c t
v
          .* :: t -> Matrix t -> Matrix t
(.*) = forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale
          |+| :: 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 = forall t. Field t => Matrix t -> Matrix t
sqrtmInv

sqrtmInv :: Matrix t -> Matrix t
sqrtmInv Matrix t
x = forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ [(Matrix t, Matrix t)] -> (Matrix t, Matrix t)
fixedPoint forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate (Matrix t, Matrix t) -> (Matrix t, Matrix t)
f (Matrix t
x, forall a. (Num a, Element a) => Int -> Matrix a
ident (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) | forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1 (forall a b. (a, b) -> a
fst (Matrix t, Matrix t)
a Matrix t -> Matrix t -> Matrix t
|-| forall a b. (a, b) -> a
fst (Matrix t, Matrix t)
b) forall a. Ord a => a -> a -> Bool
< 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)
bforall a. a -> [a] -> [a]
:[(Matrix t, Matrix t)]
rest)
          fixedPoint [(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
|+| forall t. Field t => Matrix t -> Matrix t
inv Matrix t
z),
                     t
0.5 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
(.*) = forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale
          |+| :: Matrix t -> Matrix t -> Matrix t
(|+|) = forall c. Additive c => c -> c -> c
add
          |-| :: 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 = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl forall {a} {a}. (Eq a, Num a) => a -> (a, a) -> a
f b
1 (forall a b. [a] -> [b] -> [(a, b)]
zip [a
0..a
rforall 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 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 = (forall t. Element t => [Vector t] -> Matrix t
fromColumns forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> [e]
A.elems Array Int (Vector t)
res, b
sign)
  where
    v :: [Int]
v = [Int
0..Int
rforall a. Num a => a -> a -> a
-Int
1]
    t :: [Vector t]
t = forall t. Element t => Matrix t -> [Vector t]
toColumns (forall a. (Num a, Element a) => Int -> Matrix a
ident Int
r)
    (Array Int (Vector t)
res,b
sign) = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl forall {i} {b} {e}.
(Ix i, Num b) =>
(Array i e, b) -> (i, i) -> (Array i e, b)
swap (forall i e. Ix i => (i, i) -> [e] -> Array i e
A.listArray (Int
0,Int
rforall a. Num a => a -> a -> a
-Int
1) [Vector t]
t, b
1) (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 forall a. Eq a => a -> a -> Bool
/= i
b    = (Array i e
arr forall i e. Ix i => Array i e -> [(i, e)] -> Array i e
A.// [(i
a, Array i e
arr forall i e. Ix i => Array i e -> i -> e
A.! i
b),(i
b,Array i e
arr 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 = forall {b}. (Matrix I, b) -> Vector I
res forall a b. (a -> b) -> a -> b
$ forall t u.
Element t =>
(forall s. (Int, Int) -> STMatrix s t -> ST s u)
-> Matrix t -> (Matrix t, u)
mutable forall t s.
(Num t, Element t) =>
(Int, Int) -> STMatrix s t -> ST s ()
f Matrix I
s0
  where
    s0 :: Matrix I
s0 = forall t. Storable t => Int -> Vector t -> Matrix t
reshape Int
1 (Int -> Vector I
range (forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
s))
    res :: (Matrix I, b) -> Vector I
res = forall t. Element t => Matrix t -> Vector t
flatten forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> a
fst
    swap :: STMatrix s t -> Int -> Int -> ST s ()
swap STMatrix s t
m Int
i Int
j = forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (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 = forall (t :: * -> *) (m :: * -> *) a.
(Foldable t, Monad m) =>
t (m a) -> m ()
sequence_ forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (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
rforall a. Storable a => Int -> Int -> [a] -> Matrix a
><Int
c) [Int -> Int -> a
el Int
s Int
t | Int
s<-[Int
0..Int
rforall a. Num a => a -> a -> a
-Int
1], Int
t<-[Int
0..Int
cforall a. Num a => a -> a -> a
-Int
1]]
    where el :: Int -> Int -> a
el Int
p Int
q = if Int
qforall a. Num a => a -> a -> a
-Int
pforall a. Ord a => a -> a -> Bool
>=Int
h then a
v else a
1 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 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 = forall t. Matrix t -> Int
rows Matrix t
l_u
    c :: Int
c = forall t. Matrix t -> Int
cols Matrix t
l_u
    tu :: Matrix t
tu = 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 = 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 = 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
|+| forall t. Storable t => t -> Vector t -> Int -> Int -> Matrix t
diagRect t
0 (forall (c :: * -> *) e. Container c e => e -> IndexOf c -> c e
konst' t
1 Int
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) = forall {t} {b}.
(Element t, Num t, Num b) =>
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
|+| forall t. Storable t => t -> Vector t -> Int -> Int -> Matrix t
diagRect t
0 (forall (c :: * -> *) e. Container c e => e -> IndexOf c -> c e
konst' t
1 Int
c) Int
r Int
c
    u' :: Matrix t
u' = 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
(|+|) = forall c. Additive c => c -> c -> c
add
    |*| :: 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    = forall e. Product e => Vector e -> RealOf e
norm1
    pnorm NormType
PNorm2    = forall e. (Product e, Floating e) => Vector e -> RealOf e
norm2
    pnorm NormType
Infinity  = forall e. Product e => Vector e -> RealOf e
normInf
    pnorm NormType
Frobenius = 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    = forall e. Product e => Vector e -> RealOf e
norm1
    pnorm NormType
PNorm2    = forall e. (Product e, Floating e) => Vector e -> RealOf e
norm2
    pnorm NormType
Infinity  = forall e. Product e => Vector e -> RealOf e
normInf
    pnorm NormType
Frobenius = 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    = forall e. Product e => Vector e -> RealOf e
norm1
    pnorm NormType
PNorm2    = forall e. (Product e, Floating e) => Vector e -> RealOf e
norm2
    pnorm NormType
Infinity  = forall e. Product e => Vector e -> RealOf e
normInf
    pnorm NormType
Frobenius = 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    = forall e. Product e => Vector e -> RealOf e
norm1
    pnorm NormType
PNorm2    = forall e. (Product e, Floating e) => Vector e -> RealOf e
norm2
    pnorm NormType
Infinity  = forall e. Product e => Vector e -> RealOf e
normInf
    pnorm NormType
Frobenius = 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    = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map (forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> [Vector t]
toColumns
    pnorm NormType
PNorm2    = (forall t. Storable t => Vector t -> Int -> t
@>Int
0) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Field t => Matrix t -> Vector Double
singularValues
    pnorm NormType
Infinity  = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Matrix t -> Matrix t
trans
    pnorm NormType
Frobenius = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm2 forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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    = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map (forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> [Vector t]
toColumns
    pnorm NormType
PNorm2    = (forall t. Storable t => Vector t -> Int -> t
@>Int
0) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Field t => Matrix t -> Vector Double
singularValues
    pnorm NormType
Infinity  = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Matrix t -> Matrix t
trans
    pnorm NormType
Frobenius = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm2 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> Vector t
flatten

instance Normed Matrix Float where
    pnorm :: NormType -> Matrix Float -> RealOf Float
pnorm NormType
PNorm1    = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map (forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> [Vector t]
toColumns
    pnorm NormType
PNorm2    = forall a b. (Real a, Fractional b) => a -> b
realToFrac forall b c a. (b -> c) -> (a -> b) -> a -> c
. (forall t. Storable t => Vector t -> Int -> t
@>Int
0) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Field t => Matrix t -> Vector Double
singularValues forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    pnorm NormType
Infinity  = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Matrix t -> Matrix t
trans
    pnorm NormType
Frobenius = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm2 forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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    = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map (forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Element t => Matrix t -> [Vector t]
toColumns
    pnorm NormType
PNorm2    = forall a b. (Real a, Fractional b) => a -> b
realToFrac forall b c a. (b -> c) -> (a -> b) -> a -> c
. (forall t. Storable t => Vector t -> Int -> t
@>Int
0) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Field t => Matrix t -> Vector Double
singularValues forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
double
    pnorm NormType
Infinity  = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm1 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall t. Matrix t -> Matrix t
trans
    pnorm NormType
Frobenius = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
PNorm2 forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 = forall {b} {a}. (Integral b, Real a) => a -> b
dig (c t -> RealOf t
norm (c t
x forall (c :: * -> *) e. Container c e => c e -> c e -> c e
`sub` c t
y) forall a. Fractional a => a -> a -> a
/ c t -> RealOf t
norm c t
x)
    where norm :: c t -> RealOf t
norm = forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
pnorm NormType
Infinity
          dig :: a -> b
dig a
r = forall a b. (RealFrac a, Integral b) => a -> b
round forall a b. (a -> b) -> a -> b
$ -forall a. Floating a => a -> a -> a
logBase Double
10 (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
aforall a. Num a => a -> a -> a
-a
b)
    mx :: Double
mx = forall a. Ord a => a -> a -> a
max Double
na Double
nb
    mn :: Double
mn = forall a. Ord a => a -> a -> a
min Double
na Double
nb
    r :: Double
r = if Double
mn forall a. Ord a => a -> a -> Bool
< forall x. RealFloat x => x
peps
        then Double
mx
        else Double
nabforall 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) = 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 = forall t. Field t => Matrix t -> Matrix t
cholSH Matrix t
b
    iu :: Matrix t
iu = forall t. Field t => Matrix t -> Matrix t
inv Matrix t
u
    c :: Matrix t
c = 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) = 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
(<>) = 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
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
showList :: [Herm t] -> ShowS
$cshowList :: forall t. (Show t, Element t) => [Herm t] -> ShowS
show :: Herm t -> String
$cshow :: forall t. (Show t, Element t) => Herm t -> String
showsPrec :: Int -> Herm t -> ShowS
$cshowsPrec :: forall t. (Show t, Element t) => Int -> Herm t -> ShowS
Show

instance (NFData t, Numeric t) => NFData (Herm t)
  where
    rnf :: Herm t -> ()
rnf (Herm Matrix t
m) = 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 = forall t. Matrix t -> Herm t
Herm (forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale t
0.5 (forall m mt. Transposable m mt => m -> mt
tr Matrix t
x 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 = forall t. Matrix t -> Herm t
Herm (forall m mt. Transposable m mt => m -> mt
tr Matrix t
x 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) = forall t. Matrix t -> Herm t
Herm (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) = forall t. Matrix t -> Herm t
Herm (Matrix t
a 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 = (forall t. Matrix t -> Herm t
Herm Matrix t
x)