{-# 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