{-# LANGUAGE GADTs #-}
module Numeric.LAPACK.Matrix.Triangular.Eigen (
   values,
   decompose,
   ) where

import qualified Numeric.LAPACK.Matrix.Mosaic.Basic as Mosaic
import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Vector as Vector
import Numeric.LAPACK.Matrix.Mosaic.Private
         (unpackZero, unpackToTemp, fillTriangle,
          forPointers, rowMajorPointers)
import Numeric.LAPACK.Matrix.Triangular.Basic (TriangularP)
import Numeric.LAPACK.Matrix.Layout.Private
         (Order(ColumnMajor,RowMajor), uploOrder)
import Numeric.LAPACK.Vector (Vector)
import Numeric.LAPACK.Scalar (zero)
import Numeric.LAPACK.Private (copyToColumnMajorTemp, withInfo, errorCodeMsg)

import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
import qualified Numeric.LAPACK.FFI.Real as LapackReal
import qualified Numeric.BLAS.FFI.Generic as BlasGen
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.Unchecked (Array(Array))

import Foreign.C.Types (CInt, CChar)
import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
import Foreign.Ptr (Ptr, nullPtr)

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

import Data.Complex (Complex)
import Data.Tuple.HT (swap, mapPair)


values ::
   (Layout.UpLo uplo, Shape.C sh, Class.Floating a) =>
   TriangularP pack uplo sh a -> Vector sh a
values :: TriangularP pack uplo sh a -> Vector sh a
values = TriangularP pack uplo sh a -> Vector sh a
forall uplo sh a pack mirror.
(UpLo uplo, C sh, Floating a) =>
Mosaic pack mirror uplo sh a -> Vector sh a
Mosaic.takeDiagonal


decompose ::
   (Layout.UpLo uplo, Layout.Packing vpack,
    Shape.C sh, Class.Floating a) =>
   TriangularP pack uplo sh a ->
   (TriangularP vpack uplo sh a, TriangularP vpack uplo sh a)
decompose :: TriangularP pack uplo sh a
-> (TriangularP vpack uplo sh a, TriangularP vpack uplo sh a)
decompose (Array (Layout.Mosaic PackingSingleton pack
packing MirrorSingleton NoMirror
mirror UpLoSingleton uplo
uplo Order
order sh
sh) ForeignPtr a
a) =
   let triShape :: Order -> Mosaic Unpacked NoMirror uplo sh
triShape Order
ord =
         PackingSingleton Unpacked
-> MirrorSingleton NoMirror
-> UpLoSingleton uplo
-> Order
-> sh
-> Mosaic Unpacked NoMirror uplo sh
forall pack mirror uplo size.
PackingSingleton pack
-> MirrorSingleton mirror
-> UpLoSingleton uplo
-> Order
-> size
-> Mosaic pack mirror uplo size
Layout.Mosaic
            PackingSingleton Unpacked
Layout.Unpacked MirrorSingleton NoMirror
mirror UpLoSingleton uplo
uplo (UpLoSingleton uplo -> Order -> Order
forall uplo. UpLoSingleton uplo -> Order -> Order
uploOrder UpLoSingleton uplo
uplo Order
ord) sh
sh
       n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh

   in UpLoSingleton uplo
-> (TriangularP vpack uplo sh a, TriangularP vpack uplo sh a)
-> (TriangularP vpack uplo sh a, TriangularP vpack uplo sh a)
forall uplo a. UpLoSingleton uplo -> (a, a) -> (a, a)
swapUpper UpLoSingleton uplo
uplo ((TriangularP vpack uplo sh a, TriangularP vpack uplo sh a)
 -> (TriangularP vpack uplo sh a, TriangularP vpack uplo sh a))
-> (TriangularP vpack uplo sh a, TriangularP vpack uplo sh a)
-> (TriangularP vpack uplo sh a, TriangularP vpack uplo sh a)
forall a b. (a -> b) -> a -> b
$
      (MosaicUnpacked NoMirror uplo sh a -> TriangularP vpack uplo sh a,
 MosaicUnpacked NoMirror uplo sh a -> TriangularP vpack uplo sh a)
-> (MosaicUnpacked NoMirror uplo sh a,
    MosaicUnpacked NoMirror uplo sh a)
-> (TriangularP vpack uplo sh a, TriangularP vpack uplo sh a)
forall a c b d. (a -> c, b -> d) -> (a, b) -> (c, d)
mapPair (TriangularP vpack uplo sh a -> TriangularP vpack uplo sh a
forall sh a. (C sh, Floating a) => Vector sh a -> Vector sh a
Vector.conjugate (TriangularP vpack uplo sh a -> TriangularP vpack uplo sh a)
-> (MosaicUnpacked NoMirror uplo sh a
    -> TriangularP vpack uplo sh a)
-> MosaicUnpacked NoMirror uplo sh a
-> TriangularP vpack uplo sh a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. MosaicUnpacked NoMirror uplo sh a -> TriangularP vpack uplo sh a
forall pack uplo sh a mirror.
(Packing pack, UpLo uplo, C sh, Floating a) =>
MosaicUnpacked mirror uplo sh a -> Mosaic pack mirror uplo sh a
Mosaic.repack, MosaicUnpacked NoMirror uplo sh a -> TriangularP vpack uplo sh a
forall pack uplo sh a mirror.
(Packing pack, UpLo uplo, C sh, Floating a) =>
MosaicUnpacked mirror uplo sh a -> Mosaic pack mirror uplo sh a
Mosaic.repack) ((MosaicUnpacked NoMirror uplo sh a,
  MosaicUnpacked NoMirror uplo sh a)
 -> (TriangularP vpack uplo sh a, TriangularP vpack uplo sh a))
-> (MosaicUnpacked NoMirror uplo sh a,
    MosaicUnpacked NoMirror uplo sh a)
-> (TriangularP vpack uplo sh a, TriangularP vpack uplo sh a)
forall a b. (a -> b) -> a -> b
$
      Mosaic Unpacked NoMirror uplo sh
-> (Int -> Ptr a -> IO (MosaicUnpacked NoMirror uplo sh a))
-> (MosaicUnpacked NoMirror uplo sh a,
    MosaicUnpacked NoMirror uplo sh a)
forall sh a b.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO b) -> (Array sh a, b)
Array.unsafeCreateWithSizeAndResult (Order -> Mosaic Unpacked NoMirror uplo sh
triShape Order
RowMajor) ((Int -> Ptr a -> IO (MosaicUnpacked NoMirror uplo sh a))
 -> (MosaicUnpacked NoMirror uplo sh a,
     MosaicUnpacked NoMirror uplo sh a))
-> (Int -> Ptr a -> IO (MosaicUnpacked NoMirror uplo sh a))
-> (MosaicUnpacked NoMirror uplo sh a,
    MosaicUnpacked NoMirror uplo sh a)
forall a b. (a -> b) -> a -> b
$ \Int
_ Ptr a
vlPtr ->
      Mosaic Unpacked NoMirror uplo sh
-> (Ptr a -> IO ()) -> IO (MosaicUnpacked NoMirror uplo sh a)
forall (m :: * -> *) sh a.
(PrimMonad m, C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> m (Array sh a)
ArrayIO.unsafeCreate (Order -> Mosaic Unpacked NoMirror uplo sh
triShape Order
ColumnMajor) ((Ptr a -> IO ()) -> IO (MosaicUnpacked NoMirror uplo sh a))
-> (Ptr a -> IO ()) -> IO (MosaicUnpacked NoMirror uplo sh a)
forall a b. (a -> b) -> a -> b
$ \Ptr a
vrPtr ->

   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
sidePtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'B'
      Ptr CChar
howManyPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'A'
      let selectPtr :: Ptr a
selectPtr = Ptr a
forall a. Ptr a
nullPtr
      Ptr a
aPtr <- PackingSingleton pack
-> Order -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a pack.
Floating a =>
PackingSingleton pack
-> Order -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
toColumnMajorTemp PackingSingleton pack
packing (UpLoSingleton uplo -> Order -> Order
forall uplo. UpLoSingleton uplo -> Order -> Order
uploOrder UpLoSingleton uplo
uplo Order
order) Int
n ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr CInt
mmPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
mPtr <- FortranIO () (Ptr CInt)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      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 CInt -> IO ()) -> IO ()
withInfo String
errorCodeMsg String
"trevc" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         TREVC_ a
forall a. Floating a => TREVC_ a
trevc Ptr CChar
sidePtr Ptr CChar
howManyPtr Ptr Bool
forall a. Ptr a
selectPtr Int
n
            Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
vlPtr Ptr CInt
ldaPtr Ptr a
vrPtr Ptr CInt
ldaPtr Ptr CInt
mmPtr Ptr CInt
mPtr

swapUpper :: Layout.UpLoSingleton uplo -> (a,a) -> (a,a)
swapUpper :: UpLoSingleton uplo -> (a, a) -> (a, a)
swapUpper UpLoSingleton uplo
uplo =
   case UpLoSingleton uplo
uplo of
      UpLoSingleton uplo
Layout.Lower -> (a, a) -> (a, a)
forall a. a -> a
id
      UpLoSingleton uplo
Layout.Upper -> (a, a) -> (a, a)
forall a b. (a, b) -> (b, a)
swap

toColumnMajorTemp ::
   (Class.Floating a) =>
   Layout.PackingSingleton pack -> Order ->
   Int -> ForeignPtr a -> ContT () IO (Ptr a)
toColumnMajorTemp :: PackingSingleton pack
-> Order -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
toColumnMajorTemp PackingSingleton pack
packing Order
order Int
n ForeignPtr a
a =
   case PackingSingleton pack
packing of
      PackingSingleton pack
Layout.Packed ->
         let unpk :: Int -> Ptr a -> Ptr a -> IO ()
unpk =
               case Order
order of
                  Order
ColumnMajor -> Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
unpackZero Order
ColumnMajor
                  Order
RowMajor -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
unpackZeroRowMajor
         in (Int -> Ptr a -> Ptr a -> IO ())
-> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Storable a =>
(Int -> Ptr a -> Ptr a -> IO ())
-> Int -> ForeignPtr a -> ContT r IO (Ptr a)
unpackToTemp Int -> Ptr a -> Ptr a -> IO ()
unpk Int
n ForeignPtr a
a
      PackingSingleton pack
Layout.Unpacked ->
         case Order
order of
            Order
ColumnMajor -> ((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
            Order
RowMajor -> Order -> Int -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
order Int
n Int
n ForeignPtr a
a

unpackZeroRowMajor :: Class.Floating a => Int -> Ptr a -> Ptr a -> IO ()
unpackZeroRowMajor :: Int -> Ptr a -> Ptr a -> IO ()
unpackZeroRowMajor Int
n Ptr a
packedPtr Ptr a
fullPtr = do
   a -> Order -> Int -> Ptr a -> IO ()
forall a. Floating a => a -> Order -> Int -> Ptr a -> IO ()
fillTriangle a
forall a. Floating a => a
zero Order
RowMajor Int
n Ptr a
fullPtr
   Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
unpackRowMajor Int
n Ptr a
packedPtr Ptr a
fullPtr

unpackRowMajor :: Class.Floating a => Int -> Ptr a -> Ptr a -> IO ()
unpackRowMajor :: Int -> Ptr a -> Ptr a -> IO ()
unpackRowMajor Int
n Ptr a
packedPtr Ptr a
fullPtr = 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
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
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
$
      [(Int, (Ptr a, Ptr a))]
-> (Ptr CInt -> (Ptr a, Ptr a) -> IO ()) -> IO ()
forall a. [(Int, a)] -> (Ptr CInt -> a -> IO ()) -> IO ()
forPointers (Int -> Ptr a -> Ptr a -> [(Int, (Ptr a, Ptr a))]
forall a.
Storable a =>
Int -> Ptr a -> Ptr a -> [(Int, (Ptr a, Ptr a))]
rowMajorPointers Int
n Ptr a
fullPtr Ptr a
packedPtr) ((Ptr CInt -> (Ptr a, Ptr a) -> IO ()) -> IO ())
-> (Ptr CInt -> (Ptr a, Ptr a) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
            \Ptr CInt
nPtr (Ptr a
dstPtr,Ptr a
srcPtr) ->
         Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.copy Ptr CInt
nPtr Ptr a
srcPtr Ptr CInt
incxPtr Ptr a
dstPtr Ptr CInt
incyPtr


type TREVC_ a =
   Ptr CChar -> Ptr CChar -> Ptr Bool ->
   Int -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->
   Ptr CInt -> Ptr CInt -> Ptr CInt -> IO ()

newtype TREVC a = TREVC {TREVC a -> TREVC_ a
getTREVC :: TREVC_ a}

trevc :: Class.Floating a => TREVC_ a
trevc :: TREVC_ a
trevc =
   TREVC a -> TREVC_ a
forall a. TREVC a -> TREVC_ a
getTREVC (TREVC a -> TREVC_ a) -> TREVC a -> TREVC_ a
forall a b. (a -> b) -> a -> b
$
   TREVC Float
-> TREVC Double
-> TREVC (Complex Float)
-> TREVC (Complex Double)
-> TREVC a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (TREVC_ Float -> TREVC Float
forall a. TREVC_ a -> TREVC a
TREVC TREVC_ Float
forall a. Real a => TREVC_ a
trevcReal) (TREVC_ Double -> TREVC Double
forall a. TREVC_ a -> TREVC a
TREVC TREVC_ Double
forall a. Real a => TREVC_ a
trevcReal)
      (TREVC_ (Complex Float) -> TREVC (Complex Float)
forall a. TREVC_ a -> TREVC a
TREVC TREVC_ (Complex Float)
forall a. Real a => TREVC_ (Complex a)
trevcComplex) (TREVC_ (Complex Double) -> TREVC (Complex Double)
forall a. TREVC_ a -> TREVC a
TREVC TREVC_ (Complex Double)
forall a. Real a => TREVC_ (Complex a)
trevcComplex)

trevcReal :: Class.Real a => TREVC_ a
trevcReal :: TREVC_ a
trevcReal Ptr CChar
sidePtr Ptr CChar
howmnyPtr Ptr Bool
selectPtr Int
n
      Ptr a
tPtr Ptr CInt
ldtPtr Ptr a
vlPtr Ptr CInt
ldvlPtr Ptr a
vrPtr Ptr CInt
ldvrPtr Ptr CInt
mmPtr Ptr CInt
mPtr 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
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
workPtr <- Int -> FortranIO () (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray (Int
3Int -> Int -> Int
forall a. Num a => a -> a -> a
*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
-> Ptr Bool
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr Bool
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> IO ()
LapackReal.trevc Ptr CChar
sidePtr Ptr CChar
howmnyPtr Ptr Bool
selectPtr Ptr CInt
nPtr
            Ptr a
tPtr Ptr CInt
ldtPtr Ptr a
vlPtr Ptr CInt
ldvlPtr Ptr a
vrPtr Ptr CInt
ldvrPtr Ptr CInt
mmPtr Ptr CInt
mPtr Ptr a
workPtr Ptr CInt
infoPtr

trevcComplex :: Class.Real a => TREVC_ (Complex a)
trevcComplex :: TREVC_ (Complex a)
trevcComplex Ptr CChar
sidePtr Ptr CChar
howmnyPtr Ptr Bool
selectPtr Int
n
      Ptr (Complex a)
tPtr Ptr CInt
ldtPtr Ptr (Complex a)
vlPtr Ptr CInt
ldvlPtr Ptr (Complex a)
vrPtr Ptr CInt
ldvrPtr Ptr CInt
mmPtr Ptr CInt
mPtr 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
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr (Complex a)
workPtr <- Int -> FortranIO () (Ptr (Complex 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)
      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
-> Ptr Bool
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr Bool
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr a
-> Ptr CInt
-> IO ()
LapackComplex.trevc Ptr CChar
sidePtr Ptr CChar
howmnyPtr Ptr Bool
selectPtr Ptr CInt
nPtr
            Ptr (Complex a)
tPtr Ptr CInt
ldtPtr Ptr (Complex a)
vlPtr Ptr CInt
ldvlPtr Ptr (Complex a)
vrPtr Ptr CInt
ldvrPtr Ptr CInt
mmPtr Ptr CInt
mPtr
            Ptr (Complex a)
workPtr Ptr a
rworkPtr Ptr CInt
infoPtr