{-# 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 t -> Bool
exactHermitian Matrix t
m = Matrix t
m Matrix t -> Matrix t -> Bool
forall (c :: * -> *) e. Container c e => c e -> c e -> Bool
`equal` Matrix t -> Matrix t
forall t. CTrans t => Matrix t -> Matrix t
ctrans Matrix t
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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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
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]
_) = 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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
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]
_) = 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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' :: 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' :: 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: (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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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 :: 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 (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 (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 :: 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 (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 :: 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 :: 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 :: 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 :: 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
forall t (c :: * -> *).
(Convert t, 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 :: 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 (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 (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 :: 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 :: 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 :: 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 :: Matrix t -> Vector t
nullVector = [Vector t] -> Vector t
forall a. [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 :: 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 :: 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 :: 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 :: (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. (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 :: (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 :: (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. (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 :: Matrix t -> Double
rcond Matrix t
m = [Double] -> Double
forall a. [a] -> a
last [Double]
s Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
forall a. [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 :: 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 t. (Num t, Container Vector t, CTrans t) => Matrix t -> 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
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 (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. [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 :: 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 :: 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. [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 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. [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 :: 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 (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl b -> (a, a) -> b
forall a p. (Eq a, Num p) => p -> (a, a) -> p
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 :: p -> (a, a) -> p
f p
s (a
a,a
b) | a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
b    = -p
s
                    | Bool
otherwise =  p
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 (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 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 (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 :: (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 :: 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.
(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
|+| 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 (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 (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 (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
forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
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 (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)
forall t (c :: * -> *).
(Convert t, Complexable c) =>
c t -> c (DoubleOf t)
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' :: 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 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 :: (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 :: 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' :: 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
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) = 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 :: 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 :: 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 :: 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 :: Matrix t -> Herm t
trustSym Matrix t
x = (Matrix t -> Herm t
forall t. Matrix t -> Herm t
Herm Matrix t
x)