{-# LANGUAGE TypeFamilies #-}
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 :: Square sh a -> Vector sh (ComplexOf a)
values =
   Values sh a -> Square sh a -> Vector sh (ComplexOf a)
forall sh a. Values sh a -> Values_ sh a
getValues (Values sh a -> Square sh a -> Vector sh (ComplexOf a))
-> Values sh a -> Square sh a -> Vector sh (ComplexOf 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
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 {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 :: Values_ sh a
valuesAux = (Order -> sh -> ForeignPtr a -> Array sh (Complex ar))
-> Square sh a -> Array sh (Complex ar)
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr a -> Array sh (Complex ar))
 -> Square sh a -> Array sh (Complex ar))
-> (Order -> sh -> ForeignPtr a -> Array sh (Complex ar))
-> Square sh a
-> Array sh (Complex ar)
forall a b. (a -> b) -> a -> b
$ \Order
_order sh
size ForeignPtr a
a ->
      sh -> (Int -> Ptr (Complex ar) -> IO ()) -> Array sh (Complex ar)
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
size ((Int -> Ptr (Complex ar) -> IO ()) -> Array sh (Complex ar))
-> (Int -> Ptr (Complex ar) -> IO ()) -> Array sh (Complex ar)
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr (Complex ar)
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 (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 (Complex ar)
Ptr (Complex (RealOf 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 :: Square sh a -> (Square sh a, Square sh a)
schur =
   Schur sh a -> Square sh a -> (Square sh a, Square sh a)
forall sh a. Schur sh a -> Schur_ sh a
getSchur (Schur sh a -> Square sh a -> (Square sh a, Square sh a))
-> Schur sh a -> Square sh a -> (Square sh a, Square 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
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 {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 :: Schur_ sh a
schurAux = (Order
 -> sh
 -> ForeignPtr a
 -> (Array (Square sh) a, Array (Square sh) a))
-> Schur_ sh a
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order
  -> sh
  -> ForeignPtr a
  -> (Array (Square sh) a, Array (Square sh) a))
 -> Schur_ sh a)
-> (Order
    -> sh
    -> ForeignPtr a
    -> (Array (Square sh) a, Array (Square sh) a))
-> Schur_ 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 (Array (Square sh) a))
-> (Array (Square sh) a, Array (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 (Array (Square sh) a))
 -> (Array (Square sh) a, Array (Square sh) a))
-> (Int -> Ptr a -> IO (Array (Square sh) a))
-> (Array (Square sh) a, Array (Square sh) a)
forall a b. (a -> b) -> a -> b
$ \Int
_ Ptr a
vsPtr ->
      Square sh -> (Ptr a -> IO ()) -> IO (Array (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 (Array (Square sh) a))
-> (Ptr a -> IO ()) -> IO (Array (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 (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 (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 {GEES a -> GEES_ (RealOf a) a
getGEES :: GEES_ (RealOf a) a}

gees :: Class.Floating a => GEES_ (RealOf a) a
gees :: GEES_ (RealOf a) a
gees =
   GEES a -> GEES_ (RealOf a) a
forall a. GEES a -> GEES_ (RealOf a) a
getGEES (GEES a -> GEES_ (RealOf a) a) -> GEES a -> GEES_ (RealOf a) a
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
Class.switchFloating
      (GEES_ (RealOf Float) Float -> GEES Float
forall a. GEES_ (RealOf a) a -> GEES a
GEES 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_ (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_ (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_ (RealOf (Complex Double)) (Complex Double)
forall a. Real a => GEES_ a (Complex a)
geesComplex)

geesReal :: Class.Real a => GEES_ a a
geesReal :: 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 (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 (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 :: 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 (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 :: Square sh a -> System sh ac
decompose =
   Decompose sh a -> Square sh a -> System sh (ComplexOf a)
forall sh a.
Decompose sh a -> Square sh a -> System sh (ComplexOf a)
getDecompose (Decompose sh a -> Square sh a -> System sh (ComplexOf a))
-> Decompose sh a -> Square 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
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 (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 (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 (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 (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 {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 :: Square sh a -> System sh ac
decomposeReal = (Order
 -> sh
 -> ForeignPtr a
 -> (Array (Square sh) (Complex a), Array sh (Complex a),
     Array (Square sh) (Complex a)))
-> Square sh a
-> (Array (Square sh) (Complex a), Array sh (Complex a),
    Array (Square sh) (Complex a))
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order
  -> sh
  -> ForeignPtr a
  -> (Array (Square sh) (Complex a), Array sh (Complex a),
      Array (Square sh) (Complex a)))
 -> Square sh a
 -> (Array (Square sh) (Complex a), Array sh (Complex a),
     Array (Square sh) (Complex a)))
-> (Order
    -> sh
    -> ForeignPtr a
    -> (Array (Square sh) (Complex a), Array sh (Complex a),
        Array (Square sh) (Complex a)))
-> Square sh a
-> (Array (Square sh) (Complex a), Array sh (Complex a),
    Array (Square sh) (Complex a))
forall a b. (a -> b) -> a -> b
$ \Order
order sh
size ForeignPtr a
a ->
   (\(Array sh (Complex a)
w, (Array (Square sh) (Complex a)
vlc,Array (Square sh) (Complex a)
vrc)) -> (Array (Square sh) (Complex a)
vlc, Array sh (Complex a)
w, Array (Square sh) (Complex a) -> Array (Square sh) (Complex a)
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 Array (Square sh) (Complex a)
vrc)) ((Array sh (Complex a),
  (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
 -> (Array (Square sh) (Complex a), Array sh (Complex a),
     Array (Square sh) (Complex a)))
-> (Array sh (Complex a),
    (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
-> (Array (Square sh) (Complex a), Array sh (Complex a),
    Array (Square sh) (Complex a))
forall a b. (a -> b) -> a -> b
$
   sh
-> (Int
    -> Ptr (Complex a)
    -> IO
         (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
-> (Array sh (Complex a),
    (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
forall sh a b.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO b) -> (Array sh a, b)
Array.unsafeCreateWithSizeAndResult sh
size ((Int
  -> Ptr (Complex a)
  -> IO
       (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
 -> (Array sh (Complex a),
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))))
-> (Int
    -> Ptr (Complex a)
    -> IO
         (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
-> (Array sh (Complex a),
    (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr (Complex a)
wPtr ->
   ContT
  (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
  IO
  (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
-> IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT
   (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
   IO
   (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
 -> IO
      (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
-> IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
forall a b. (a -> b) -> a -> b
$ do
      Ptr CChar
jobvlPtr <- Char
-> FortranIO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'V'
      Ptr CChar
jobvrPtr <- Char
-> FortranIO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'V'
      Ptr CInt
nPtr <- Int
-> FortranIO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
aPtr <- Int
-> ForeignPtr a
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     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
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr a
wrPtr <- Int
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     IO
     (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      Ptr a
wiPtr <- Int
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     IO
     (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      Ptr a
vlPtr <- Int
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     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
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr a
vrPtr <- Int
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     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
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      IO ()
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     IO
     ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO ()
 -> ContT
      (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
      IO
      ())
-> IO ()
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     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
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     IO
     ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO ()
 -> ContT
      (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
      IO
      ())
-> IO ()
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     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
      IO (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
 -> ContT
      (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
      IO
      (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
-> IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
forall a b. (a -> b) -> a -> b
$ Order
-> Square sh
-> (Ptr (Complex a) -> Ptr (Complex a) -> IO ())
-> IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
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 (Complex a) -> Ptr (Complex a) -> IO ())
 -> IO
      (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
-> (Ptr (Complex a) -> Ptr (Complex a) -> IO ())
-> IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
forall a b. (a -> b) -> a -> b
$
         \Ptr (Complex a)
vlcPtr Ptr (Complex a)
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 (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 (Complex a)
vrcPtr

eigenvectorsToComplex ::
   (Eq a, Class.Real a) =>
   Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
eigenvectorsToComplex :: 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 (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 (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 :: Square sh ac -> System sh ac
decomposeComplex = (Order
 -> sh
 -> ForeignPtr (Complex a)
 -> (Array (Square sh) (Complex a), Array sh (Complex a),
     Array (Square sh) (Complex a)))
-> Array (Square sh) (Complex a)
-> (Array (Square sh) (Complex a), Array sh (Complex a),
    Array (Square sh) (Complex a))
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order
  -> sh
  -> ForeignPtr (Complex a)
  -> (Array (Square sh) (Complex a), Array sh (Complex a),
      Array (Square sh) (Complex a)))
 -> Array (Square sh) (Complex a)
 -> (Array (Square sh) (Complex a), Array sh (Complex a),
     Array (Square sh) (Complex a)))
-> (Order
    -> sh
    -> ForeignPtr (Complex a)
    -> (Array (Square sh) (Complex a), Array sh (Complex a),
        Array (Square sh) (Complex a)))
-> Array (Square sh) (Complex a)
-> (Array (Square sh) (Complex a), Array sh (Complex a),
    Array (Square sh) (Complex a))
forall a b. (a -> b) -> a -> b
$ \Order
order sh
size ForeignPtr (Complex a)
a ->
   (\(Array sh (Complex a)
w, (Array (Square sh) (Complex a)
vlc,Array (Square sh) (Complex a)
vrc)) -> (Array (Square sh) (Complex a)
vlc, Array sh (Complex a)
w, Array (Square sh) (Complex a) -> Array (Square sh) (Complex a)
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 Array (Square sh) (Complex a)
vrc)) ((Array sh (Complex a),
  (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
 -> (Array (Square sh) (Complex a), Array sh (Complex a),
     Array (Square sh) (Complex a)))
-> (Array sh (Complex a),
    (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
-> (Array (Square sh) (Complex a), Array sh (Complex a),
    Array (Square sh) (Complex a))
forall a b. (a -> b) -> a -> b
$
   sh
-> (Int
    -> Ptr (Complex a)
    -> IO
         (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
-> (Array sh (Complex a),
    (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
forall sh a b.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO b) -> (Array sh a, b)
Array.unsafeCreateWithSizeAndResult sh
size ((Int
  -> Ptr (Complex a)
  -> IO
       (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
 -> (Array sh (Complex a),
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))))
-> (Int
    -> Ptr (Complex a)
    -> IO
         (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
-> (Array sh (Complex a),
    (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr (Complex a)
wPtr ->
   ContT
  (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
  IO
  (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
-> IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT
   (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
   IO
   (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
 -> IO
      (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
-> IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
forall a b. (a -> b) -> a -> b
$ do
      Ptr CChar
jobvlPtr <- Char
-> FortranIO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'V'
      Ptr CChar
jobvrPtr <- Char
-> FortranIO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'V'
      Ptr CInt
nPtr <- Int
-> FortranIO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr (Complex a)
aPtr <- Int
-> ForeignPtr (Complex a)
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     IO
     (Ptr (Complex 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 (Complex a)
a
      Ptr CInt
ldaPtr <- Int
-> FortranIO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr CInt
ldvlPtr <- Int
-> FortranIO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr CInt
ldvrPtr <- Int
-> FortranIO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr a
rworkPtr <- Int
-> FortranIO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     (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 (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
 -> ContT
      (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
      IO
      (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
-> IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
-> ContT
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
     IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
forall a b. (a -> b) -> a -> b
$ Order
-> Square sh
-> (Ptr (Complex a) -> Ptr (Complex a) -> IO ())
-> IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
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 (Complex a) -> Ptr (Complex a) -> IO ())
 -> IO
      (Array (Square sh) (Complex a), Array (Square sh) (Complex a)))
-> (Ptr (Complex a) -> Ptr (Complex a) -> IO ())
-> IO
     (Array (Square sh) (Complex a), Array (Square sh) (Complex a))
forall a b. (a -> b) -> a -> b
$
         \Ptr (Complex a)
vlPtr Ptr (Complex a)
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 (Complex a)
aPtr Ptr CInt
ldaPtr
            Ptr (Complex a)
wPtr Ptr (Complex a)
vlPtr Ptr CInt
ldvlPtr 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 :: 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 (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 (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 :: 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 (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"