{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
module Numeric.LAPACK.Matrix.Square.Eigen (
   values,
   schur,
   decompose,
   ComplexOf,
   ) where

import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Basic as Basic
import qualified Numeric.LAPACK.Shape as ExtShape
import Numeric.LAPACK.Matrix.Layout.Private (Order(ColumnMajor), swapOnRowMajor)
import Numeric.LAPACK.Matrix.Private (Square, argSquare)
import Numeric.LAPACK.Vector (Vector)
import Numeric.LAPACK.Scalar (ComplexOf, RealOf, zero)
import Numeric.LAPACK.Private
         (copyConjugate, copyToTemp, copyToColumnMajor,
          realPtr, withAutoWorkspaceInfo)

import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
import qualified Numeric.LAPACK.FFI.Real as LapackReal
import qualified Numeric.BLAS.FFI.Real as BlasReal
import qualified Numeric.Netlib.Utility as Call
import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
import qualified Data.Array.Comfort.Storable.Unchecked as Array
import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Storable (Array)

import Foreign.Marshal.Array (advancePtr, peekArray)
import Foreign.C.Types (CInt, CChar)
import Foreign.ForeignPtr (withForeignPtr)
import Foreign.Ptr (Ptr, nullPtr, nullFunPtr)
import Foreign.Storable (Storable)

import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
import Control.Monad.IO.Class (liftIO)

import Data.Complex (Complex)


values ::
   (ExtShape.Permutable sh, Class.Floating a) =>
   Square sh a -> Vector sh (ComplexOf a)
values :: forall sh a.
(Permutable sh, Floating a) =>
Square sh a -> Vector sh (ComplexOf a)
values =
   Values sh a
-> Array (Full Shape Small Small sh sh) a
-> Array sh (Complex (RealOf a))
forall sh a. Values sh a -> Values_ sh a
getValues (Values sh a
 -> Array (Full Shape Small Small sh sh) a
 -> Array sh (Complex (RealOf a)))
-> Values sh a
-> Array (Full Shape Small Small sh sh) a
-> Array sh (Complex (RealOf a))
forall a b. (a -> b) -> a -> b
$
   Values sh Float
-> Values sh Double
-> Values sh (Complex Float)
-> Values sh (Complex Double)
-> Values sh a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
forall (f :: * -> *).
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (Values_ sh Float -> Values sh Float
forall sh a. Values_ sh a -> Values sh a
Values Values_ sh Float
forall sh a ar.
(Permutable sh, Floating a, RealOf a ~ ar, Storable ar) =>
Values_ sh a
valuesAux) (Values_ sh Double -> Values sh Double
forall sh a. Values_ sh a -> Values sh a
Values Values_ sh Double
forall sh a ar.
(Permutable sh, Floating a, RealOf a ~ ar, Storable ar) =>
Values_ sh a
valuesAux)
      (Values_ sh (Complex Float) -> Values sh (Complex Float)
forall sh a. Values_ sh a -> Values sh a
Values Values_ sh (Complex Float)
forall sh a ar.
(Permutable sh, Floating a, RealOf a ~ ar, Storable ar) =>
Values_ sh a
valuesAux) (Values_ sh (Complex Double) -> Values sh (Complex Double)
forall sh a. Values_ sh a -> Values sh a
Values Values_ sh (Complex Double)
forall sh a ar.
(Permutable sh, Floating a, RealOf a ~ ar, Storable ar) =>
Values_ sh a
valuesAux)

type Values_ sh a = Square sh a -> Vector sh (ComplexOf a)

newtype Values sh a = Values {forall sh a. Values sh a -> Values_ sh a
getValues :: Values_ sh a}

valuesAux ::
   (ExtShape.Permutable sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>
   Values_ sh a
valuesAux :: forall sh a ar.
(Permutable sh, Floating a, RealOf a ~ ar, Storable ar) =>
Values_ sh a
valuesAux = (Order -> sh -> ForeignPtr a -> Vector sh (ComplexOf a))
-> Square sh a -> Vector sh (ComplexOf a)
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr a -> Vector sh (ComplexOf a))
 -> Square sh a -> Vector sh (ComplexOf a))
-> (Order -> sh -> ForeignPtr a -> Vector sh (ComplexOf a))
-> Square sh a
-> Vector sh (ComplexOf a)
forall a b. (a -> b) -> a -> b
$ \Order
_order sh
size ForeignPtr a
a ->
      sh
-> (Int -> Ptr (ComplexOf a) -> IO ()) -> Vector sh (ComplexOf a)
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
size ((Int -> Ptr (ComplexOf a) -> IO ()) -> Vector sh (ComplexOf a))
-> (Int -> Ptr (ComplexOf a) -> IO ()) -> Vector sh (ComplexOf a)
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr (ComplexOf a)
wPtr -> do
   let lda :: Int
lda = Int
n
   ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      Ptr CChar
jobvsPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'N'
      Ptr CChar
sortPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'N'
      Ptr a
aPtr <- Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToTemp (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n) ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
lda
      Ptr CInt
sdimPtr <- FortranIO () (Ptr CInt)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      let vsPtr :: Ptr a
vsPtr = Ptr a
forall a. Ptr a
nullPtr
      Ptr CInt
ldvsPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      let bworkPtr :: Ptr a
bworkPtr = Ptr a
forall a. Ptr a
nullPtr
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         String
-> String -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
String
-> String -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo String
eigenMsg String
"gees" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr a
workPtr Ptr CInt
lworkPtr Ptr CInt
infoPtr ->
         GEES_ (RealOf a) a
forall a. Floating a => GEES_ (RealOf a) a
gees
            Ptr CChar
jobvsPtr Ptr CChar
sortPtr Int
n Ptr a
aPtr Ptr CInt
ldaPtr Ptr CInt
sdimPtr
            Ptr (ComplexOf a)
wPtr Ptr a
forall a. Ptr a
vsPtr Ptr CInt
ldvsPtr Ptr a
workPtr Ptr CInt
lworkPtr Ptr Bool
forall a. Ptr a
bworkPtr Ptr CInt
infoPtr


schur ::
   (ExtShape.Permutable sh, Class.Floating a) =>
   Square sh a -> (Square sh a, Square sh a)
schur :: forall sh a.
(Permutable sh, Floating a) =>
Square sh a -> (Square sh a, Square sh a)
schur =
   Schur sh a -> Schur_ sh a
forall sh a. Schur sh a -> Schur_ sh a
getSchur (Schur sh a -> Schur_ sh a) -> Schur sh a -> Schur_ sh a
forall a b. (a -> b) -> a -> b
$
   Schur sh Float
-> Schur sh Double
-> Schur sh (Complex Float)
-> Schur sh (Complex Double)
-> Schur sh a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
forall (f :: * -> *).
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (Schur_ sh Float -> Schur sh Float
forall sh a. Schur_ sh a -> Schur sh a
Schur Schur_ sh Float
forall sh a ar.
(Permutable sh, Floating a, RealOf a ~ ar, Storable ar) =>
Schur_ sh a
schurAux) (Schur_ sh Double -> Schur sh Double
forall sh a. Schur_ sh a -> Schur sh a
Schur Schur_ sh Double
forall sh a ar.
(Permutable sh, Floating a, RealOf a ~ ar, Storable ar) =>
Schur_ sh a
schurAux)
      (Schur_ sh (Complex Float) -> Schur sh (Complex Float)
forall sh a. Schur_ sh a -> Schur sh a
Schur Schur_ sh (Complex Float)
forall sh a ar.
(Permutable sh, Floating a, RealOf a ~ ar, Storable ar) =>
Schur_ sh a
schurAux) (Schur_ sh (Complex Double) -> Schur sh (Complex Double)
forall sh a. Schur_ sh a -> Schur sh a
Schur Schur_ sh (Complex Double)
forall sh a ar.
(Permutable sh, Floating a, RealOf a ~ ar, Storable ar) =>
Schur_ sh a
schurAux)

type Schur_ sh a = Square sh a -> (Square sh a, Square sh a)

newtype Schur sh a = Schur {forall sh a. Schur sh a -> Schur_ sh a
getSchur :: Schur_ sh a}

schurAux ::
   (ExtShape.Permutable sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>
   Schur_ sh a
schurAux :: forall sh a ar.
(Permutable sh, Floating a, RealOf a ~ ar, Storable ar) =>
Schur_ sh a
schurAux = (Order -> sh -> ForeignPtr a -> (Square sh a, Square sh a))
-> Square sh a -> (Square sh a, Square sh a)
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr a -> (Square sh a, Square sh a))
 -> Square sh a -> (Square sh a, Square sh a))
-> (Order -> sh -> ForeignPtr a -> (Square sh a, Square sh a))
-> Square sh a
-> (Square sh a, Square sh a)
forall a b. (a -> b) -> a -> b
$ \Order
order sh
size ForeignPtr a
a ->
   let sh :: Square sh
sh = Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
Layout.square Order
ColumnMajor sh
size
   in Square sh
-> (Int -> Ptr a -> IO (Square sh a)) -> (Square sh a, Square sh a)
forall sh a b.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO b) -> (Array sh a, b)
Array.unsafeCreateWithSizeAndResult Square sh
sh ((Int -> Ptr a -> IO (Square sh a)) -> (Square sh a, Square sh a))
-> (Int -> Ptr a -> IO (Square sh a)) -> (Square sh a, Square sh a)
forall a b. (a -> b) -> a -> b
$ \Int
_ Ptr a
vsPtr ->
      Square sh -> (Ptr a -> IO ()) -> IO (Square sh a)
forall (m :: * -> *) sh a.
(PrimMonad m, C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> m (Array sh a)
ArrayIO.unsafeCreate Square sh
sh ((Ptr a -> IO ()) -> IO (Square sh a))
-> (Ptr a -> IO ()) -> IO (Square sh a)
forall a b. (a -> b) -> a -> b
$ \Ptr a
sPtr -> do

   let n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
size
   let lda :: Int
lda = Int
n
   ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      Ptr CChar
jobvsPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'V'
      Ptr CChar
sortPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'N'
      Ptr a
aPtr <- ((Ptr a -> IO ()) -> IO ()) -> ContT () IO (Ptr a)
forall {k} (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO ()) -> IO ()) -> ContT () IO (Ptr a))
-> ((Ptr a -> IO ()) -> IO ()) -> ContT () IO (Ptr a)
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> (Ptr a -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
a
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyToColumnMajor Order
order Int
n Int
n Ptr a
aPtr Ptr a
sPtr
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
lda
      Ptr CInt
sdimPtr <- FortranIO () (Ptr CInt)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      Ptr (Complex ar)
wPtr <- Int -> FortranIO () (Ptr (Complex ar))
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      Ptr CInt
ldvsPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      let bworkPtr :: Ptr a
bworkPtr = Ptr a
forall a. Ptr a
nullPtr
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         String
-> String -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
String
-> String -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo String
eigenMsg String
"gees" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr a
workPtr Ptr CInt
lworkPtr Ptr CInt
infoPtr ->
         GEES_ (RealOf a) a
forall a. Floating a => GEES_ (RealOf a) a
gees
            Ptr CChar
jobvsPtr Ptr CChar
sortPtr Int
n Ptr a
sPtr Ptr CInt
ldaPtr Ptr CInt
sdimPtr
            Ptr (Complex ar)
Ptr (Complex (RealOf a))
wPtr Ptr a
vsPtr Ptr CInt
ldvsPtr Ptr a
workPtr Ptr CInt
lworkPtr Ptr Bool
forall a. Ptr a
bworkPtr Ptr CInt
infoPtr



type GEES_ ar a =
   Ptr CChar -> Ptr CChar -> Int -> Ptr a -> Ptr CInt ->
   Ptr CInt -> Ptr (Complex ar) -> Ptr a -> Ptr CInt ->
   Ptr a -> Ptr CInt -> Ptr Bool -> Ptr CInt -> IO ()

newtype GEES a = GEES {forall a. GEES a -> GEES_ (RealOf a) a
getGEES :: GEES_ (RealOf a) a}

gees :: Class.Floating a => GEES_ (RealOf a) a
gees :: forall a. Floating a => GEES_ (RealOf a) a
gees =
   GEES a
-> Ptr CChar
-> Ptr CChar
-> Int
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr (Complex (RealOf a))
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr Bool
-> Ptr CInt
-> IO ()
forall a. GEES a -> GEES_ (RealOf a) a
getGEES (GEES a
 -> Ptr CChar
 -> Ptr CChar
 -> Int
 -> Ptr a
 -> Ptr CInt
 -> Ptr CInt
 -> Ptr (Complex (RealOf a))
 -> Ptr a
 -> Ptr CInt
 -> Ptr a
 -> Ptr CInt
 -> Ptr Bool
 -> Ptr CInt
 -> IO ())
-> GEES a
-> Ptr CChar
-> Ptr CChar
-> Int
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr (Complex (RealOf a))
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr Bool
-> Ptr CInt
-> IO ()
forall a b. (a -> b) -> a -> b
$
   GEES Float
-> GEES Double
-> GEES (Complex Float)
-> GEES (Complex Double)
-> GEES a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
forall (f :: * -> *).
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (GEES_ (RealOf Float) Float -> GEES Float
forall a. GEES_ (RealOf a) a -> GEES a
GEES GEES_ Float Float
GEES_ (RealOf Float) Float
forall a. Real a => GEES_ a a
geesReal) (GEES_ (RealOf Double) Double -> GEES Double
forall a. GEES_ (RealOf a) a -> GEES a
GEES GEES_ Double Double
GEES_ (RealOf Double) Double
forall a. Real a => GEES_ a a
geesReal) (GEES_ (RealOf (Complex Float)) (Complex Float)
-> GEES (Complex Float)
forall a. GEES_ (RealOf a) a -> GEES a
GEES GEES_ Float (Complex Float)
GEES_ (RealOf (Complex Float)) (Complex Float)
forall a. Real a => GEES_ a (Complex a)
geesComplex) (GEES_ (RealOf (Complex Double)) (Complex Double)
-> GEES (Complex Double)
forall a. GEES_ (RealOf a) a -> GEES a
GEES GEES_ Double (Complex Double)
GEES_ (RealOf (Complex Double)) (Complex Double)
forall a. Real a => GEES_ a (Complex a)
geesComplex)

geesReal :: Class.Real a => GEES_ a a
geesReal :: forall a. Real a => GEES_ a a
geesReal
      Ptr CChar
jobvsPtr Ptr CChar
sortPtr Int
n Ptr a
aPtr Ptr CInt
ldaPtr Ptr CInt
sdimPtr
      Ptr (Complex a)
wPtr Ptr a
vsPtr Ptr CInt
ldvsPtr Ptr a
workPtr Ptr CInt
lworkPtr Ptr Bool
bworkPtr Ptr CInt
infoPtr =
   ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      let selectPtr :: FunPtr a
selectPtr = FunPtr a
forall a. FunPtr a
nullFunPtr
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
wrPtr <- Int -> FortranIO () (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      Ptr a
wiPtr <- Int -> FortranIO () (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CChar
-> Ptr CChar
-> FunPtr (Ptr a -> Ptr a -> IO Bool)
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr Bool
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> FunPtr (Ptr a -> Ptr a -> IO Bool)
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr Bool
-> Ptr CInt
-> IO ()
LapackReal.gees
            Ptr CChar
jobvsPtr Ptr CChar
sortPtr FunPtr (Ptr a -> Ptr a -> IO Bool)
forall a. FunPtr a
selectPtr Ptr CInt
nPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr CInt
sdimPtr
            Ptr a
wrPtr Ptr a
wiPtr Ptr a
vsPtr Ptr CInt
ldvsPtr Ptr a
workPtr Ptr CInt
lworkPtr Ptr Bool
bworkPtr Ptr CInt
infoPtr
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
forall a.
Real a =>
Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
zipComplex Int
n Ptr a
wrPtr Ptr a
wiPtr Ptr (Complex a)
wPtr

geesComplex :: Class.Real a => GEES_ a (Complex a)
geesComplex :: forall a. Real a => GEES_ a (Complex a)
geesComplex
      Ptr CChar
jobvsPtr Ptr CChar
sortPtr Int
n Ptr (Complex a)
aPtr Ptr CInt
ldaPtr Ptr CInt
sdimPtr
      Ptr (Complex a)
wPtr Ptr (Complex a)
vsPtr Ptr CInt
ldvsPtr Ptr (Complex a)
workPtr Ptr CInt
lworkPtr Ptr Bool
bworkPtr Ptr CInt
infoPtr =
   ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      let selectPtr :: FunPtr a
selectPtr = FunPtr a
forall a. FunPtr a
nullFunPtr
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
rworkPtr <- Int -> FortranIO () (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CChar
-> Ptr CChar
-> FunPtr (Ptr (Complex a) -> IO Bool)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr a
-> Ptr Bool
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> FunPtr (Ptr (Complex a) -> IO Bool)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr a
-> Ptr Bool
-> Ptr CInt
-> IO ()
LapackComplex.gees
            Ptr CChar
jobvsPtr Ptr CChar
sortPtr FunPtr (Ptr (Complex a) -> IO Bool)
forall a. FunPtr a
selectPtr Ptr CInt
nPtr Ptr (Complex a)
aPtr Ptr CInt
ldaPtr Ptr CInt
sdimPtr
            Ptr (Complex a)
wPtr Ptr (Complex a)
vsPtr Ptr CInt
ldvsPtr Ptr (Complex a)
workPtr Ptr CInt
lworkPtr Ptr a
rworkPtr Ptr Bool
bworkPtr Ptr CInt
infoPtr



type System sh a = (Square sh a, Vector sh a, Square sh a)

decompose ::
   (ExtShape.Permutable sh, Class.Floating a, ComplexOf a ~ ac) =>
   Square sh a -> System sh ac
decompose :: forall sh a ac.
(Permutable sh, Floating a, ComplexOf a ~ ac) =>
Square sh a -> System sh ac
decompose =
   Decompose sh a
-> Array (Full Shape Small Small sh sh) a
-> System sh (ComplexOf a)
forall sh a.
Decompose sh a -> Square sh a -> System sh (ComplexOf a)
getDecompose (Decompose sh a
 -> Array (Full Shape Small Small sh sh) a
 -> System sh (ComplexOf a))
-> Decompose sh a
-> Array (Full Shape Small Small sh sh) a
-> System sh (ComplexOf a)
forall a b. (a -> b) -> a -> b
$
   Decompose sh Float
-> Decompose sh Double
-> Decompose sh (Complex Float)
-> Decompose sh (Complex Double)
-> Decompose sh a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
forall (f :: * -> *).
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      ((Square sh Float -> System sh (ComplexOf Float))
-> Decompose sh Float
forall sh a.
(Square sh a -> System sh (ComplexOf a)) -> Decompose sh a
Decompose Square sh Float -> System sh (Complex Float)
Square sh Float -> System sh (ComplexOf Float)
forall sh a ac.
(Permutable sh, Real a, Complex a ~ ac) =>
Square sh a -> System sh ac
decomposeReal)
      ((Square sh Double -> System sh (ComplexOf Double))
-> Decompose sh Double
forall sh a.
(Square sh a -> System sh (ComplexOf a)) -> Decompose sh a
Decompose Square sh Double -> System sh (Complex Double)
Square sh Double -> System sh (ComplexOf Double)
forall sh a ac.
(Permutable sh, Real a, Complex a ~ ac) =>
Square sh a -> System sh ac
decomposeReal)
      ((Square sh (Complex Float)
 -> System sh (ComplexOf (Complex Float)))
-> Decompose sh (Complex Float)
forall sh a.
(Square sh a -> System sh (ComplexOf a)) -> Decompose sh a
Decompose Square sh (Complex Float) -> System sh (Complex Float)
Square sh (Complex Float) -> System sh (ComplexOf (Complex Float))
forall sh a ac.
(Permutable sh, Real a, Complex a ~ ac) =>
Square sh ac -> System sh ac
decomposeComplex)
      ((Square sh (Complex Double)
 -> System sh (ComplexOf (Complex Double)))
-> Decompose sh (Complex Double)
forall sh a.
(Square sh a -> System sh (ComplexOf a)) -> Decompose sh a
Decompose Square sh (Complex Double) -> System sh (Complex Double)
Square sh (Complex Double)
-> System sh (ComplexOf (Complex Double))
forall sh a ac.
(Permutable sh, Real a, Complex a ~ ac) =>
Square sh ac -> System sh ac
decomposeComplex)

newtype Decompose sh a =
   Decompose {forall sh a.
Decompose sh a -> Square sh a -> System sh (ComplexOf a)
getDecompose :: Square sh a -> System sh (ComplexOf a)}

decomposeReal ::
   (ExtShape.Permutable sh, Class.Real a, Complex a ~ ac) =>
   Square sh a -> System sh ac
decomposeReal :: forall sh a ac.
(Permutable sh, Real a, Complex a ~ ac) =>
Square sh a -> System sh ac
decomposeReal = (Order -> sh -> ForeignPtr a -> System sh ac)
-> Square sh a -> System sh ac
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr a -> System sh ac)
 -> Square sh a -> System sh ac)
-> (Order -> sh -> ForeignPtr a -> System sh ac)
-> Square sh a
-> System sh ac
forall a b. (a -> b) -> a -> b
$ \Order
order sh
size ForeignPtr a
a ->
   (\(Vector sh ac
w, (Square sh ac
vlc,Square sh ac
vrc)) -> (Square sh ac
vlc, Vector sh ac
w, Square sh ac -> Square sh ac
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Full meas vert horiz height width a
-> Full meas horiz vert width height a
Basic.adjoint Square sh ac
vrc)) ((Vector sh ac, (Square sh ac, Square sh ac)) -> System sh ac)
-> (Vector sh ac, (Square sh ac, Square sh ac)) -> System sh ac
forall a b. (a -> b) -> a -> b
$
   sh
-> (Int -> Ptr ac -> IO (Square sh ac, Square sh ac))
-> (Vector sh ac, (Square sh ac, Square sh ac))
forall sh a b.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO b) -> (Array sh a, b)
Array.unsafeCreateWithSizeAndResult sh
size ((Int -> Ptr ac -> IO (Square sh ac, Square sh ac))
 -> (Vector sh ac, (Square sh ac, Square sh ac)))
-> (Int -> Ptr ac -> IO (Square sh ac, Square sh ac))
-> (Vector sh ac, (Square sh ac, Square sh ac))
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr ac
wPtr ->
   ContT (Square sh ac, Square sh ac) IO (Square sh ac, Square sh ac)
-> IO (Square sh ac, Square sh ac)
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT (Square sh ac, Square sh ac) IO (Square sh ac, Square sh ac)
 -> IO (Square sh ac, Square sh ac))
-> ContT
     (Square sh ac, Square sh ac) IO (Square sh ac, Square sh ac)
-> IO (Square sh ac, Square sh ac)
forall a b. (a -> b) -> a -> b
$ do
      Ptr CChar
jobvlPtr <- Char -> FortranIO (Square sh ac, Square sh ac) (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'V'
      Ptr CChar
jobvrPtr <- Char -> FortranIO (Square sh ac, Square sh ac) (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'V'
      Ptr CInt
nPtr <- Int -> FortranIO (Square sh ac, Square sh ac) (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
aPtr <- Int
-> ForeignPtr a -> ContT (Square sh ac, Square sh ac) IO (Ptr a)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToTemp (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n) ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO (Square sh ac, Square sh ac) (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr a
wrPtr <- Int -> ContT (Square sh ac, Square sh ac) IO (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      Ptr a
wiPtr <- Int -> ContT (Square sh ac, Square sh ac) IO (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      Ptr a
vlPtr <- Int -> ContT (Square sh ac, Square sh ac) IO (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)
      Ptr CInt
ldvlPtr <- Int -> FortranIO (Square sh ac, Square sh ac) (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr a
vrPtr <- Int -> ContT (Square sh ac, Square sh ac) IO (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)
      Ptr CInt
ldvrPtr <- Int -> FortranIO (Square sh ac, Square sh ac) (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      IO () -> ContT (Square sh ac, Square sh ac) IO ()
forall a. IO a -> ContT (Square sh ac, Square sh ac) IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT (Square sh ac, Square sh ac) IO ())
-> IO () -> ContT (Square sh ac, Square sh ac) IO ()
forall a b. (a -> b) -> a -> b
$ String
-> String -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
String
-> String -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo String
eigenMsg String
"geev" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackReal.geev
            Ptr CChar
jobvlPtr Ptr CChar
jobvrPtr Ptr CInt
nPtr Ptr a
aPtr Ptr CInt
ldaPtr
            Ptr a
wrPtr Ptr a
wiPtr Ptr a
vlPtr Ptr CInt
ldvlPtr Ptr a
vrPtr Ptr CInt
ldvrPtr
      IO () -> ContT (Square sh ac, Square sh ac) IO ()
forall a. IO a -> ContT (Square sh ac, Square sh ac) IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT (Square sh ac, Square sh ac) IO ())
-> IO () -> ContT (Square sh ac, Square sh ac) IO ()
forall a b. (a -> b) -> a -> b
$ Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
forall a.
Real a =>
Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
zipComplex Int
n Ptr a
wrPtr Ptr a
wiPtr Ptr ac
Ptr (Complex a)
wPtr
      IO (Square sh ac, Square sh ac)
-> ContT
     (Square sh ac, Square sh ac) IO (Square sh ac, Square sh ac)
forall a. IO a -> ContT (Square sh ac, Square sh ac) IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (Square sh ac, Square sh ac)
 -> ContT
      (Square sh ac, Square sh ac) IO (Square sh ac, Square sh ac))
-> IO (Square sh ac, Square sh ac)
-> ContT
     (Square sh ac, Square sh ac) IO (Square sh ac, Square sh ac)
forall a b. (a -> b) -> a -> b
$ Order
-> Square sh
-> (Ptr ac -> Ptr ac -> IO ())
-> IO (Square sh ac, Square sh ac)
forall sh a.
(C sh, Storable a) =>
Order
-> sh -> (Ptr a -> Ptr a -> IO ()) -> IO (Array sh a, Array sh a)
createArrayPair Order
order (Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
Layout.square Order
ColumnMajor sh
size) ((Ptr ac -> Ptr ac -> IO ()) -> IO (Square sh ac, Square sh ac))
-> (Ptr ac -> Ptr ac -> IO ()) -> IO (Square sh ac, Square sh ac)
forall a b. (a -> b) -> a -> b
$
         \Ptr ac
vlcPtr Ptr ac
vrcPtr -> do
            Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
forall a.
(Eq a, Real a) =>
Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
eigenvectorsToComplex Int
n Ptr a
wiPtr Ptr a
vlPtr Ptr ac
Ptr (Complex a)
vlcPtr
            Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
forall a.
(Eq a, Real a) =>
Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
eigenvectorsToComplex Int
n Ptr a
wiPtr Ptr a
vrPtr Ptr ac
Ptr (Complex a)
vrcPtr

eigenvectorsToComplex ::
   (Eq a, Class.Real a) =>
   Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
eigenvectorsToComplex :: forall a.
(Eq a, Real a) =>
Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
eigenvectorsToComplex Int
n Ptr a
wiPtr Ptr a
vPtr Ptr (Complex a)
vcPtr = ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
   Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
   Ptr a
zeroPtr <- a -> FortranIO () (Ptr a)
forall a r. Real a => a -> FortranIO r (Ptr a)
Call.real a
forall a. Floating a => a
zero
   Ptr CInt
inc0Ptr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
   Ptr CInt
inc1Ptr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
   Ptr CInt
inc2Ptr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
2
   IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ do
      let go :: Ptr a -> Ptr (Complex a) -> [Bool] -> IO ()
go Ptr a
_ Ptr (Complex a)
_ [] = () -> IO ()
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return ()
          go Ptr a
xPtr Ptr (Complex a)
yPtr (Bool
False:[Bool]
wi) = do
            let yrPtr :: Ptr (RealOf (Complex a))
yrPtr = Ptr (Complex a) -> Ptr (RealOf (Complex a))
forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
yPtr
            let yiPtr :: Ptr a
yiPtr = Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
Ptr (RealOf (Complex a))
yrPtr Int
1
            Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasReal.copy Ptr CInt
nPtr Ptr a
xPtr    Ptr CInt
inc1Ptr Ptr a
Ptr (RealOf (Complex a))
yrPtr Ptr CInt
inc2Ptr
            Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasReal.copy Ptr CInt
nPtr Ptr a
zeroPtr Ptr CInt
inc0Ptr Ptr a
yiPtr Ptr CInt
inc2Ptr
            Ptr a -> Ptr (Complex a) -> [Bool] -> IO ()
go (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
xPtr Int
n) (Ptr (Complex a) -> Int -> Ptr (Complex a)
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (Complex a)
yPtr Int
n) [Bool]
wi
          go Ptr a
xPtr Ptr (Complex a)
yPtr (Bool
True:Bool
True:[Bool]
wi) = do
            let xrPtr :: Ptr a
xrPtr = Ptr a
xPtr
            let xiPtr :: Ptr a
xiPtr = Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
xPtr Int
n
            let yrPtr :: Ptr (RealOf (Complex a))
yrPtr = Ptr (Complex a) -> Ptr (RealOf (Complex a))
forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
yPtr
            let yiPtr :: Ptr a
yiPtr = Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
Ptr (RealOf (Complex a))
yrPtr Int
1
            let y1Ptr :: Ptr (Complex a)
y1Ptr = Ptr (Complex a) -> Int -> Ptr (Complex a)
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (Complex a)
yPtr Int
n
            Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasReal.copy Ptr CInt
nPtr Ptr a
xrPtr Ptr CInt
inc1Ptr Ptr a
Ptr (RealOf (Complex a))
yrPtr Ptr CInt
inc2Ptr
            Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasReal.copy Ptr CInt
nPtr Ptr a
xiPtr Ptr CInt
inc1Ptr Ptr a
yiPtr Ptr CInt
inc2Ptr
            Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
copyConjugate Ptr CInt
nPtr Ptr (Complex a)
yPtr Ptr CInt
inc1Ptr Ptr (Complex a)
y1Ptr Ptr CInt
inc1Ptr
            Ptr a -> Ptr (Complex a) -> [Bool] -> IO ()
go (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
xPtr (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)) (Ptr (Complex a) -> Int -> Ptr (Complex a)
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (Complex a)
yPtr (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)) [Bool]
wi
          go Ptr a
_xPtr Ptr (Complex a)
_yPtr [Bool]
wi =
            String -> IO ()
forall a. HasCallStack => String -> a
error (String -> IO ()) -> String -> IO ()
forall a b. (a -> b) -> a -> b
$ String
"eigenvectorToComplex: invalid non-real pattern " String -> String -> String
forall a. [a] -> [a] -> [a]
++ [Bool] -> String
forall a. Show a => a -> String
show [Bool]
wi
      Ptr a -> Ptr (Complex a) -> [Bool] -> IO ()
go Ptr a
vPtr Ptr (Complex a)
vcPtr ([Bool] -> IO ()) -> ([a] -> [Bool]) -> [a] -> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> Bool) -> [a] -> [Bool]
forall a b. (a -> b) -> [a] -> [b]
map (a
forall a. Floating a => a
zeroa -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=) ([a] -> IO ()) -> IO [a] -> IO ()
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Int -> Ptr a -> IO [a]
forall a. Storable a => Int -> Ptr a -> IO [a]
peekArray Int
n Ptr a
wiPtr

decomposeComplex ::
   (ExtShape.Permutable sh, Class.Real a, Complex a ~ ac) =>
   Square sh ac -> System sh ac
decomposeComplex :: forall sh a ac.
(Permutable sh, Real a, Complex a ~ ac) =>
Square sh ac -> System sh ac
decomposeComplex = (Order -> sh -> ForeignPtr ac -> System sh ac)
-> Square sh ac -> System sh ac
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr ac -> System sh ac)
 -> Square sh ac -> System sh ac)
-> (Order -> sh -> ForeignPtr ac -> System sh ac)
-> Square sh ac
-> System sh ac
forall a b. (a -> b) -> a -> b
$ \Order
order sh
size ForeignPtr ac
a ->
   (\(Vector sh ac
w, (Square sh ac
vlc,Square sh ac
vrc)) -> (Square sh ac
vlc, Vector sh ac
w, Square sh ac -> Square sh ac
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Full meas vert horiz height width a
-> Full meas horiz vert width height a
Basic.adjoint Square sh ac
vrc)) ((Vector sh ac, (Square sh ac, Square sh ac)) -> System sh ac)
-> (Vector sh ac, (Square sh ac, Square sh ac)) -> System sh ac
forall a b. (a -> b) -> a -> b
$
   sh
-> (Int -> Ptr ac -> IO (Square sh ac, Square sh ac))
-> (Vector sh ac, (Square sh ac, Square sh ac))
forall sh a b.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO b) -> (Array sh a, b)
Array.unsafeCreateWithSizeAndResult sh
size ((Int -> Ptr ac -> IO (Square sh ac, Square sh ac))
 -> (Vector sh ac, (Square sh ac, Square sh ac)))
-> (Int -> Ptr ac -> IO (Square sh ac, Square sh ac))
-> (Vector sh ac, (Square sh ac, Square sh ac))
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr ac
wPtr ->
   ContT (Square sh ac, Square sh ac) IO (Square sh ac, Square sh ac)
-> IO (Square sh ac, Square sh ac)
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT (Square sh ac, Square sh ac) IO (Square sh ac, Square sh ac)
 -> IO (Square sh ac, Square sh ac))
-> ContT
     (Square sh ac, Square sh ac) IO (Square sh ac, Square sh ac)
-> IO (Square sh ac, Square sh ac)
forall a b. (a -> b) -> a -> b
$ do
      Ptr CChar
jobvlPtr <- Char -> FortranIO (Square sh ac, Square sh ac) (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'V'
      Ptr CChar
jobvrPtr <- Char -> FortranIO (Square sh ac, Square sh ac) (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'V'
      Ptr CInt
nPtr <- Int -> FortranIO (Square sh ac, Square sh ac) (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr ac
aPtr <- Int
-> ForeignPtr ac -> ContT (Square sh ac, Square sh ac) IO (Ptr ac)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToTemp (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n) ForeignPtr ac
a
      Ptr CInt
ldaPtr <- Int -> FortranIO (Square sh ac, Square sh ac) (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr CInt
ldvlPtr <- Int -> FortranIO (Square sh ac, Square sh ac) (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr CInt
ldvrPtr <- Int -> FortranIO (Square sh ac, Square sh ac) (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr a
rworkPtr <- Int -> FortranIO (Square sh ac, Square sh ac) (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)

      IO (Square sh ac, Square sh ac)
-> ContT
     (Square sh ac, Square sh ac) IO (Square sh ac, Square sh ac)
forall a. IO a -> ContT (Square sh ac, Square sh ac) IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (Square sh ac, Square sh ac)
 -> ContT
      (Square sh ac, Square sh ac) IO (Square sh ac, Square sh ac))
-> IO (Square sh ac, Square sh ac)
-> ContT
     (Square sh ac, Square sh ac) IO (Square sh ac, Square sh ac)
forall a b. (a -> b) -> a -> b
$ Order
-> Square sh
-> (Ptr ac -> Ptr ac -> IO ())
-> IO (Square sh ac, Square sh ac)
forall sh a.
(C sh, Storable a) =>
Order
-> sh -> (Ptr a -> Ptr a -> IO ()) -> IO (Array sh a, Array sh a)
createArrayPair Order
order (Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
Layout.square Order
ColumnMajor sh
size) ((Ptr ac -> Ptr ac -> IO ()) -> IO (Square sh ac, Square sh ac))
-> (Ptr ac -> Ptr ac -> IO ()) -> IO (Square sh ac, Square sh ac)
forall a b. (a -> b) -> a -> b
$
         \Ptr ac
vlPtr Ptr ac
vrPtr ->

         String
-> String
-> (Ptr (Complex a) -> Ptr CInt -> Ptr CInt -> IO ())
-> IO ()
forall a.
Floating a =>
String
-> String -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo String
eigenMsg String
"geev" ((Ptr (Complex a) -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr (Complex a) -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr (Complex a)
workPtr Ptr CInt
lworkPtr Ptr CInt
infoPtr ->
         Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> IO ()
LapackComplex.geev
            Ptr CChar
jobvlPtr Ptr CChar
jobvrPtr Ptr CInt
nPtr Ptr ac
Ptr (Complex a)
aPtr Ptr CInt
ldaPtr
            Ptr ac
Ptr (Complex a)
wPtr Ptr ac
Ptr (Complex a)
vlPtr Ptr CInt
ldvlPtr Ptr ac
Ptr (Complex a)
vrPtr Ptr CInt
ldvrPtr
            Ptr (Complex a)
workPtr Ptr CInt
lworkPtr Ptr a
rworkPtr Ptr CInt
infoPtr


zipComplex ::
   (Class.Real a) => Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
zipComplex :: forall a.
Real a =>
Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
zipComplex Int
n Ptr a
vr Ptr a
vi Ptr (Complex a)
vc =
   ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      Ptr CInt
incyPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
2
      let yPtr :: Ptr (RealOf (Complex a))
yPtr = Ptr (Complex a) -> Ptr (RealOf (Complex a))
forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
vc
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasReal.copy Ptr CInt
nPtr Ptr a
vr Ptr CInt
incxPtr Ptr a
Ptr (RealOf (Complex a))
yPtr Ptr CInt
incyPtr
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasReal.copy Ptr CInt
nPtr Ptr a
vi Ptr CInt
incxPtr (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
Ptr (RealOf (Complex a))
yPtr Int
1) Ptr CInt
incyPtr

createArrayPair ::
   (Shape.C sh, Storable a) =>
   Order -> sh -> (Ptr a -> Ptr a -> IO ()) ->
   IO (Array sh a, Array sh a)
createArrayPair :: forall sh a.
(C sh, Storable a) =>
Order
-> sh -> (Ptr a -> Ptr a -> IO ()) -> IO (Array sh a, Array sh a)
createArrayPair Order
order sh
sh Ptr a -> Ptr a -> IO ()
act =
   ((Array sh a, Array sh a) -> (Array sh a, Array sh a))
-> IO (Array sh a, Array sh a) -> IO (Array sh a, Array sh a)
forall a b. (a -> b) -> IO a -> IO b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Order -> (Array sh a, Array sh a) -> (Array sh a, Array sh a)
forall a. Order -> (a, a) -> (a, a)
swapOnRowMajor Order
order) (IO (Array sh a, Array sh a) -> IO (Array sh a, Array sh a))
-> IO (Array sh a, Array sh a) -> IO (Array sh a, Array sh a)
forall a b. (a -> b) -> a -> b
$
   sh
-> (Int -> Ptr a -> IO (Array sh a)) -> IO (Array sh a, Array sh a)
forall (m :: * -> *) sh a b.
(PrimMonad m, C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO b) -> m (Array sh a, b)
ArrayIO.unsafeCreateWithSizeAndResult sh
sh ((Int -> Ptr a -> IO (Array sh a)) -> IO (Array sh a, Array sh a))
-> (Int -> Ptr a -> IO (Array sh a)) -> IO (Array sh a, Array sh a)
forall a b. (a -> b) -> a -> b
$ \Int
_ Ptr a
vrcPtr ->
   sh -> (Ptr a -> IO ()) -> IO (Array sh a)
forall (m :: * -> *) sh a.
(PrimMonad m, C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> m (Array sh a)
ArrayIO.unsafeCreate sh
sh ((Ptr a -> IO ()) -> IO (Array sh a))
-> (Ptr a -> IO ()) -> IO (Array sh a)
forall a b. (a -> b) -> a -> b
$ \Ptr a
vlcPtr -> Ptr a -> Ptr a -> IO ()
act Ptr a
vlcPtr Ptr a
vrcPtr


eigenMsg :: String
eigenMsg :: String
eigenMsg = String
"only eigenvalues starting with the %d-th one converged"