{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
module Numeric.LAPACK.Matrix.Hermitian.Basic (
   Hermitian,
   Transposition(..),
   fromList,
   autoFromList,
   recheck,
   identity,
   diagonal,
   takeDiagonal,
   forceOrder,

   stack,
   takeTopLeft,
   takeTopRight,
   takeBottomRight,

   multiplyVector,
   multiplyFull,
   square, power,
   outer,
   sumRank1,
   sumRank2,

   toSquare,
   gramian,              gramianAdjoint,
   congruenceDiagonal,   congruenceDiagonalAdjoint,
   congruence,           congruenceAdjoint,
   scaledAnticommutator, scaledAnticommutatorAdjoint,
   addAdjoint,

   takeUpper,
   ) where

import qualified Numeric.LAPACK.Matrix.Symmetric.Private as Symmetric
import qualified Numeric.LAPACK.Matrix.Triangular.Private as Triangular
import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Matrix.Basic as Basic
import qualified Numeric.LAPACK.Split as Split
import Numeric.LAPACK.Matrix.Hermitian.Private (Diagonal(..), TakeDiagonal(..))
import Numeric.LAPACK.Matrix.Triangular.Private
         (forPointers, pack, unpack, unpackToTemp,
          diagonalPointers, diagonalPointerPairs,
          rowMajorPointers, columnMajorPointers)
import Numeric.LAPACK.Matrix.Shape.Private
         (Order(RowMajor,ColumnMajor), flipOrder, sideSwapFromOrder,
          uploFromOrder)
import Numeric.LAPACK.Matrix.Modifier
         (Transposition(NonTransposed, Transposed), transposeOrder,
          Conjugation(Conjugated), conjugatedOnRowMajor)
import Numeric.LAPACK.Matrix.Private
         (Full, General, Square, argSquare, ShapeInt, shapeInt)
import Numeric.LAPACK.Vector (Vector)
import Numeric.LAPACK.Scalar (RealOf, zero, one)
import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))
import Numeric.LAPACK.Private
         (fill, lacgv, realPtr,
          copyConjugate, condConjugate, conjugateToTemp, condConjugateToTemp)

import qualified Numeric.BLAS.FFI.Generic as BlasGen
import qualified Numeric.BLAS.FFI.Complex as BlasComplex
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 as Array
import qualified Data.Array.Comfort.Storable as CheckedArray
import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Storable.Unchecked (Array(Array))
import Data.Array.Comfort.Shape ((:+:)((:+:)))

import Foreign.C.Types (CInt, CChar)
import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
import Foreign.Ptr (Ptr)
import Foreign.Storable (Storable, poke, peek)

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

import Data.Foldable (forM_)

import Data.Function.HT (powerAssociative)


type Hermitian sh = Array (MatrixShape.Hermitian sh)


fromList :: (Shape.C sh, Storable a) => Order -> sh -> [a] -> Hermitian sh a
fromList :: Order -> sh -> [a] -> Hermitian sh a
fromList Order
order sh
sh =
   Hermitian sh -> [a] -> Hermitian sh a
forall sh a. (C sh, Storable a) => sh -> [a] -> Array sh a
CheckedArray.fromList (Order -> sh -> Hermitian sh
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order sh
sh)

autoFromList :: (Storable a) => Order -> [a] -> Hermitian ShapeInt a
autoFromList :: Order -> [a] -> Hermitian ShapeInt a
autoFromList Order
order [a]
xs =
   Order -> ShapeInt -> [a] -> Hermitian ShapeInt a
forall sh a.
(C sh, Storable a) =>
Order -> sh -> [a] -> Hermitian sh a
fromList Order
order
      (Int -> ShapeInt
shapeInt (Int -> ShapeInt) -> Int -> ShapeInt
forall a b. (a -> b) -> a -> b
$ String -> Int -> Int
MatrixShape.triangleExtent String
"Hermitian.autoFromList" (Int -> Int) -> Int -> Int
forall a b. (a -> b) -> a -> b
$
       [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs)
      [a]
xs

uncheck :: Hermitian sh a -> Hermitian (Unchecked sh) a
uncheck :: Hermitian sh a -> Hermitian (Unchecked sh) a
uncheck =
   (Hermitian sh -> Hermitian (Unchecked sh))
-> Hermitian sh a -> Hermitian (Unchecked sh) a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape ((Hermitian sh -> Hermitian (Unchecked sh))
 -> Hermitian sh a -> Hermitian (Unchecked sh) a)
-> (Hermitian sh -> Hermitian (Unchecked sh))
-> Hermitian sh a
-> Hermitian (Unchecked sh) a
forall a b. (a -> b) -> a -> b
$
      \(MatrixShape.Hermitian Order
order sh
sh) ->
         Order -> Unchecked sh -> Hermitian (Unchecked sh)
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order (sh -> Unchecked sh
forall sh. sh -> Unchecked sh
Unchecked sh
sh)

recheck :: Hermitian (Unchecked sh) a -> Hermitian sh a
recheck :: Hermitian (Unchecked sh) a -> Hermitian sh a
recheck =
   (Hermitian (Unchecked sh) -> Hermitian sh)
-> Hermitian (Unchecked sh) a -> Hermitian sh a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape ((Hermitian (Unchecked sh) -> Hermitian sh)
 -> Hermitian (Unchecked sh) a -> Hermitian sh a)
-> (Hermitian (Unchecked sh) -> Hermitian sh)
-> Hermitian (Unchecked sh) a
-> Hermitian sh a
forall a b. (a -> b) -> a -> b
$
      \(MatrixShape.Hermitian Order
order (Unchecked sh
sh)) ->
         Order -> sh -> Hermitian sh
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order sh
sh


identity :: (Shape.C sh, Class.Floating a) => Order -> sh -> Hermitian sh a
identity :: Order -> sh -> Hermitian sh a
identity Order
order sh
sh =
   Hermitian sh -> (Int -> Ptr a -> IO ()) -> Hermitian sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (Order -> sh -> Hermitian sh
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order sh
sh) ((Int -> Ptr a -> IO ()) -> Hermitian sh a)
-> (Int -> Ptr a -> IO ()) -> Hermitian sh a
forall a b. (a -> b) -> a -> b
$
      \Int
triSize Ptr a
aPtr -> do
   a -> Int -> Ptr a -> IO ()
forall a. Floating a => a -> Int -> Ptr a -> IO ()
fill a
forall a. Floating a => a
zero Int
triSize Ptr a
aPtr
   (Ptr a -> IO ()) -> [Ptr a] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ ((Ptr a -> a -> IO ()) -> a -> Ptr a -> IO ()
forall a b c. (a -> b -> c) -> b -> a -> c
flip Ptr a -> a -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke a
forall a. Floating a => a
one) ([Ptr a] -> IO ()) -> [Ptr a] -> IO ()
forall a b. (a -> b) -> a -> b
$ Order -> Int -> Ptr a -> [Ptr a]
forall a. Storable a => Order -> Int -> Ptr a -> [Ptr a]
diagonalPointers Order
order (sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh) Ptr a
aPtr

diagonal ::
   (Shape.C sh, Class.Floating a) =>
   Order -> Vector sh (RealOf a) -> Hermitian sh a
diagonal :: Order -> Vector sh (RealOf a) -> Hermitian sh a
diagonal Order
order =
   Diagonal (Array (Hermitian sh)) (Array sh) a
-> Vector sh (RealOf a) -> Hermitian sh a
forall (f :: * -> *) (g :: * -> *) a.
Diagonal f g a -> g (RealOf a) -> f a
runDiagonal (Diagonal (Array (Hermitian sh)) (Array sh) a
 -> Vector sh (RealOf a) -> Hermitian sh a)
-> Diagonal (Array (Hermitian sh)) (Array sh) a
-> Vector sh (RealOf a)
-> Hermitian sh a
forall a b. (a -> b) -> a -> b
$
   Diagonal (Array (Hermitian sh)) (Array sh) Float
-> Diagonal (Array (Hermitian sh)) (Array sh) Double
-> Diagonal (Array (Hermitian sh)) (Array sh) (Complex Float)
-> Diagonal (Array (Hermitian sh)) (Array sh) (Complex Double)
-> Diagonal (Array (Hermitian sh)) (Array sh) a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      ((Array sh (RealOf Float) -> Array (Hermitian sh) Float)
-> Diagonal (Array (Hermitian sh)) (Array sh) Float
forall (f :: * -> *) (g :: * -> *) a.
(g (RealOf a) -> f a) -> Diagonal f g a
Diagonal ((Array sh (RealOf Float) -> Array (Hermitian sh) Float)
 -> Diagonal (Array (Hermitian sh)) (Array sh) Float)
-> (Array sh (RealOf Float) -> Array (Hermitian sh) Float)
-> Diagonal (Array (Hermitian sh)) (Array sh) Float
forall a b. (a -> b) -> a -> b
$ Order -> Vector sh Float -> Array (Hermitian sh) Float
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Storable ar) =>
Order -> Vector sh ar -> Hermitian sh a
diagonalAux Order
order) ((Array sh (RealOf Double) -> Array (Hermitian sh) Double)
-> Diagonal (Array (Hermitian sh)) (Array sh) Double
forall (f :: * -> *) (g :: * -> *) a.
(g (RealOf a) -> f a) -> Diagonal f g a
Diagonal ((Array sh (RealOf Double) -> Array (Hermitian sh) Double)
 -> Diagonal (Array (Hermitian sh)) (Array sh) Double)
-> (Array sh (RealOf Double) -> Array (Hermitian sh) Double)
-> Diagonal (Array (Hermitian sh)) (Array sh) Double
forall a b. (a -> b) -> a -> b
$ Order -> Vector sh Double -> Array (Hermitian sh) Double
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Storable ar) =>
Order -> Vector sh ar -> Hermitian sh a
diagonalAux Order
order)
      ((Array sh (RealOf (Complex Float))
 -> Array (Hermitian sh) (Complex Float))
-> Diagonal (Array (Hermitian sh)) (Array sh) (Complex Float)
forall (f :: * -> *) (g :: * -> *) a.
(g (RealOf a) -> f a) -> Diagonal f g a
Diagonal ((Array sh (RealOf (Complex Float))
  -> Array (Hermitian sh) (Complex Float))
 -> Diagonal (Array (Hermitian sh)) (Array sh) (Complex Float))
-> (Array sh (RealOf (Complex Float))
    -> Array (Hermitian sh) (Complex Float))
-> Diagonal (Array (Hermitian sh)) (Array sh) (Complex Float)
forall a b. (a -> b) -> a -> b
$ Order -> Vector sh Float -> Array (Hermitian sh) (Complex Float)
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Storable ar) =>
Order -> Vector sh ar -> Hermitian sh a
diagonalAux Order
order) ((Array sh (RealOf (Complex Double))
 -> Array (Hermitian sh) (Complex Double))
-> Diagonal (Array (Hermitian sh)) (Array sh) (Complex Double)
forall (f :: * -> *) (g :: * -> *) a.
(g (RealOf a) -> f a) -> Diagonal f g a
Diagonal ((Array sh (RealOf (Complex Double))
  -> Array (Hermitian sh) (Complex Double))
 -> Diagonal (Array (Hermitian sh)) (Array sh) (Complex Double))
-> (Array sh (RealOf (Complex Double))
    -> Array (Hermitian sh) (Complex Double))
-> Diagonal (Array (Hermitian sh)) (Array sh) (Complex Double)
forall a b. (a -> b) -> a -> b
$ Order -> Vector sh Double -> Array (Hermitian sh) (Complex Double)
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Storable ar) =>
Order -> Vector sh ar -> Hermitian sh a
diagonalAux Order
order)

diagonalAux ::
   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>
   Order -> Vector sh ar -> Hermitian sh a
diagonalAux :: Order -> Vector sh ar -> Hermitian sh a
diagonalAux Order
order (Array sh
sh ForeignPtr ar
x) =
   Hermitian sh -> (Int -> Ptr a -> IO ()) -> Hermitian sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (Order -> sh -> Hermitian sh
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order sh
sh) ((Int -> Ptr a -> IO ()) -> Hermitian sh a)
-> (Int -> Ptr a -> IO ()) -> Hermitian sh a
forall a b. (a -> b) -> a -> b
$
      \Int
triSize Ptr a
aPtr -> do
   a -> Int -> Ptr a -> IO ()
forall a. Floating a => a -> Int -> Ptr a -> IO ()
fill a
forall a. Floating a => a
zero Int
triSize Ptr a
aPtr
   ForeignPtr ar -> (Ptr ar -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr ar
x ((Ptr ar -> IO ()) -> IO ()) -> (Ptr ar -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr ar
xPtr ->
      [(Ptr ar, Ptr a)] -> ((Ptr ar, Ptr a) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ (Order -> Int -> Ptr ar -> Ptr a -> [(Ptr ar, Ptr a)]
forall a b.
(Storable a, Storable b) =>
Order -> Int -> Ptr a -> Ptr b -> [(Ptr a, Ptr b)]
diagonalPointerPairs Order
order (sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh) Ptr ar
xPtr Ptr a
aPtr) (((Ptr ar, Ptr a) -> IO ()) -> IO ())
-> ((Ptr ar, Ptr a) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         \(Ptr ar
srcPtr,Ptr a
dstPtr) -> Ptr ar -> ar -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke (Ptr a -> Ptr (RealOf a)
forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr a
dstPtr) (ar -> IO ()) -> IO ar -> IO ()
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Ptr ar -> IO ar
forall a. Storable a => Ptr a -> IO a
peek Ptr ar
srcPtr


takeDiagonal ::
   (Shape.C sh, Class.Floating a) =>
   Hermitian sh a -> Vector sh (RealOf a)
takeDiagonal :: Hermitian sh a -> Vector sh (RealOf a)
takeDiagonal =
   TakeDiagonal (Array (Hermitian sh)) (Array sh) a
-> Hermitian sh a -> Vector sh (RealOf a)
forall (f :: * -> *) (g :: * -> *) a.
TakeDiagonal f g a -> f a -> g (RealOf a)
runTakeDiagonal (TakeDiagonal (Array (Hermitian sh)) (Array sh) a
 -> Hermitian sh a -> Vector sh (RealOf a))
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) a
-> Hermitian sh a
-> Vector sh (RealOf a)
forall a b. (a -> b) -> a -> b
$
   TakeDiagonal (Array (Hermitian sh)) (Array sh) Float
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) Double
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) (Complex Float)
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) (Complex Double)
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      ((Array (Hermitian sh) Float -> Array sh (RealOf Float))
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) Float
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (RealOf a)) -> TakeDiagonal f g a
TakeDiagonal Array (Hermitian sh) Float -> Array sh (RealOf Float)
forall sh a ar.
(C sh, Storable a, RealOf a ~ ar, Storable ar) =>
Hermitian sh a -> Vector sh ar
takeDiagonalAux) ((Array (Hermitian sh) Double -> Array sh (RealOf Double))
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) Double
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (RealOf a)) -> TakeDiagonal f g a
TakeDiagonal Array (Hermitian sh) Double -> Array sh (RealOf Double)
forall sh a ar.
(C sh, Storable a, RealOf a ~ ar, Storable ar) =>
Hermitian sh a -> Vector sh ar
takeDiagonalAux)
      ((Array (Hermitian sh) (Complex Float)
 -> Array sh (RealOf (Complex Float)))
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) (Complex Float)
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (RealOf a)) -> TakeDiagonal f g a
TakeDiagonal Array (Hermitian sh) (Complex Float)
-> Array sh (RealOf (Complex Float))
forall sh a ar.
(C sh, Storable a, RealOf a ~ ar, Storable ar) =>
Hermitian sh a -> Vector sh ar
takeDiagonalAux) ((Array (Hermitian sh) (Complex Double)
 -> Array sh (RealOf (Complex Double)))
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) (Complex Double)
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (RealOf a)) -> TakeDiagonal f g a
TakeDiagonal Array (Hermitian sh) (Complex Double)
-> Array sh (RealOf (Complex Double))
forall sh a ar.
(C sh, Storable a, RealOf a ~ ar, Storable ar) =>
Hermitian sh a -> Vector sh ar
takeDiagonalAux)

takeDiagonalAux ::
   (Shape.C sh, Storable a, RealOf a ~ ar, Storable ar) =>
   Hermitian sh a -> Vector sh ar
takeDiagonalAux :: Hermitian sh a -> Vector sh ar
takeDiagonalAux (Array (MatrixShape.Hermitian Order
order sh
sh) ForeignPtr a
a) =
   sh -> (Int -> Ptr ar -> IO ()) -> Vector sh ar
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
sh ((Int -> Ptr ar -> IO ()) -> Vector sh ar)
-> (Int -> Ptr ar -> IO ()) -> Vector sh ar
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr ar
xPtr ->
   ForeignPtr a -> (Ptr a -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
a ((Ptr a -> IO ()) -> IO ()) -> (Ptr a -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr a
aPtr ->
      [(Ptr ar, Ptr a)] -> ((Ptr ar, Ptr a) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ (Order -> Int -> Ptr ar -> Ptr a -> [(Ptr ar, Ptr a)]
forall a b.
(Storable a, Storable b) =>
Order -> Int -> Ptr a -> Ptr b -> [(Ptr a, Ptr b)]
diagonalPointerPairs Order
order Int
n Ptr ar
xPtr Ptr a
aPtr) (((Ptr ar, Ptr a) -> IO ()) -> IO ())
-> ((Ptr ar, Ptr a) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         \(Ptr ar
dstPtr,Ptr a
srcPtr) -> Ptr ar -> ar -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr ar
dstPtr (ar -> IO ()) -> IO ar -> IO ()
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Ptr ar -> IO ar
forall a. Storable a => Ptr a -> IO a
peek (Ptr a -> Ptr (RealOf a)
forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr a
srcPtr)


{-
This is not maximally efficient.
It fills up a whole square.
This wastes memory but enables more regular memory access patterns.
Additionally, it fills unused parts of the square with mirrored values.
-}
forceOrder ::
   (Shape.C sh, Class.Floating a) =>
   Order -> Hermitian sh a -> Hermitian sh a
forceOrder :: Order -> Hermitian sh a -> Hermitian sh a
forceOrder Order
newOrder Hermitian sh a
a =
   if Hermitian sh -> Order
forall size. Hermitian size -> Order
MatrixShape.hermitianOrder (Hermitian sh a -> Hermitian sh
forall sh a. Array sh a -> sh
Array.shape Hermitian sh a
a) Order -> Order -> Bool
forall a. Eq a => a -> a -> Bool
== Order
newOrder
      then Hermitian sh a
a
      else Full Small Small sh sh a -> Hermitian sh a
forall vert height width a.
(C vert, C height, C width, Floating a) =>
Full vert Small height width a -> Hermitian width a
fromUpperPart (Full Small Small sh sh a -> Hermitian sh a)
-> Full Small Small sh sh a -> Hermitian sh a
forall a b. (a -> b) -> a -> b
$ Order -> Full Small Small sh sh a -> Full Small Small sh sh a
forall vert horiz height width a.
(C vert, C horiz, C height, C width, Floating a) =>
Order
-> Full vert horiz height width a -> Full vert horiz height width a
Basic.forceOrder Order
newOrder (Full Small Small sh sh a -> Full Small Small sh sh a)
-> Full Small Small sh sh a -> Full Small Small sh sh a
forall a b. (a -> b) -> a -> b
$ Hermitian sh a -> Full Small Small sh sh a
forall sh a. (C sh, Floating a) => Hermitian sh a -> Square sh a
toSquare Hermitian sh a
a

fromUpperPart ::
   (Extent.C vert, Shape.C height, Shape.C width, Class.Floating a) =>
   Full vert Extent.Small height width a -> Hermitian width a
fromUpperPart :: Full vert Small height width a -> Hermitian width a
fromUpperPart = (Order -> width -> Hermitian width)
-> Full vert Small height width a -> Hermitian width a
forall vert height width shape a.
(C vert, C height, C width, C shape, Floating a) =>
(Order -> width -> shape)
-> Full vert Small height width a -> Array shape a
Triangular.fromUpperPart Order -> width -> Hermitian width
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian

{-
Naming is inconsistent to Triangular.takeUpper,
because here Hermitian is the input
and in Triangular.takeUpper, Triangular is the output.
-}
takeUpper ::
   (Shape.C sh, Class.Floating a) =>
   Hermitian sh a ->
   Array (MatrixShape.UpperTriangular MatrixShape.NonUnit sh) a
takeUpper :: Hermitian sh a -> Array (UpperTriangular NonUnit sh) a
takeUpper =
   (Hermitian sh -> UpperTriangular NonUnit sh)
-> Hermitian sh a -> Array (UpperTriangular NonUnit sh) a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape
      (\(MatrixShape.Hermitian Order
order sh
sh) ->
         NonUnit
-> (Empty, Filled) -> Order -> sh -> UpperTriangular NonUnit sh
forall lo diag up size.
diag -> (lo, up) -> Order -> size -> Triangular lo diag up size
MatrixShape.Triangular NonUnit
MatrixShape.NonUnit (Empty, Filled)
MatrixShape.upper Order
order sh
sh)


stack ::
   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>
   Hermitian sh0 a -> General sh0 sh1 a -> Hermitian sh1 a ->
   Hermitian (sh0:+:sh1) a
stack :: Hermitian sh0 a
-> General sh0 sh1 a
-> Hermitian sh1 a
-> Hermitian (sh0 :+: sh1) a
stack Hermitian sh0 a
a General sh0 sh1 a
b Hermitian sh1 a
c =
   let order :: Order
order = Full Big Big sh0 sh1 -> Order
forall vert horiz height width.
Full vert horiz height width -> Order
MatrixShape.fullOrder (Full Big Big sh0 sh1 -> Order) -> Full Big Big sh0 sh1 -> Order
forall a b. (a -> b) -> a -> b
$ General sh0 sh1 a -> Full Big Big sh0 sh1
forall sh a. Array sh a -> sh
Array.shape General sh0 sh1 a
b
   in String
-> ((sh0 :+: sh1) -> Hermitian (sh0 :+: sh1))
-> Hermitian sh0 a
-> General sh0 sh1 a
-> Hermitian sh1 a
-> Hermitian (sh0 :+: sh1) a
forall sh0 height sh1 width sh2 a.
(Box sh0, HeightOf sh0 ~ height, C height, Eq height, Box sh1,
 WidthOf sh1 ~ width, C width, Eq width, C sh2, Floating a) =>
String
-> ((height :+: width) -> sh2)
-> Array sh0 a
-> General height width a
-> Array sh1 a
-> Array sh2 a
Triangular.stack String
"Hermitian" (Order -> (sh0 :+: sh1) -> Hermitian (sh0 :+: sh1)
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order)
         (Order -> Hermitian sh0 a -> Hermitian sh0 a
forall sh a.
(C sh, Floating a) =>
Order -> Hermitian sh a -> Hermitian sh a
forceOrder Order
order Hermitian sh0 a
a) General sh0 sh1 a
b (Order -> Hermitian sh1 a -> Hermitian sh1 a
forall sh a.
(C sh, Floating a) =>
Order -> Hermitian sh a -> Hermitian sh a
forceOrder Order
order Hermitian sh1 a
c)

takeTopLeft ::
   (Shape.C sh0, Shape.C sh1, Class.Floating a) =>
   Hermitian (sh0:+:sh1) a -> Hermitian sh0 a
takeTopLeft :: Hermitian (sh0 :+: sh1) a -> Hermitian sh0 a
takeTopLeft =
   (Hermitian (sh0 :+: sh1) -> (Hermitian sh0, (Order, sh0 :+: sh1)))
-> Hermitian (sh0 :+: sh1) a -> Hermitian sh0 a
forall sh sha height width a.
(C sh, C sha, C height, C width, Floating a) =>
(sh -> (sha, (Order, height :+: width)))
-> Array sh a -> Array sha a
Triangular.takeTopLeft
      (\(MatrixShape.Hermitian Order
order sh :: sh0 :+: sh1
sh@(sh0
sh0:+:sh1
_sh1)) ->
         (Order -> sh0 -> Hermitian sh0
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order sh0
sh0, (Order
order,sh0 :+: sh1
sh)))

takeTopRight ::
   (Shape.C sh0, Shape.C sh1, Class.Floating a) =>
   Hermitian (sh0:+:sh1) a -> General sh0 sh1 a
takeTopRight :: Hermitian (sh0 :+: sh1) a -> General sh0 sh1 a
takeTopRight =
   (Hermitian (sh0 :+: sh1) -> (Order, sh0 :+: sh1))
-> Hermitian (sh0 :+: sh1) a -> General sh0 sh1 a
forall sh height width a.
(C sh, C height, C width, Floating a) =>
(sh -> (Order, height :+: width))
-> Array sh a -> General height width a
Triangular.takeTopRight (\(MatrixShape.Hermitian Order
order sh0 :+: sh1
sh) -> (Order
order,sh0 :+: sh1
sh))

takeBottomRight ::
   (Shape.C sh0, Shape.C sh1, Class.Floating a) =>
   Hermitian (sh0:+:sh1) a -> Hermitian sh1 a
takeBottomRight :: Hermitian (sh0 :+: sh1) a -> Hermitian sh1 a
takeBottomRight =
   (Hermitian (sh0 :+: sh1) -> (Hermitian sh1, (Order, sh0 :+: sh1)))
-> Hermitian (sh0 :+: sh1) a -> Hermitian sh1 a
forall sh sha height width a.
(C sh, C sha, C height, C width, Floating a) =>
(sh -> (sha, (Order, height :+: width)))
-> Array sh a -> Array sha a
Triangular.takeBottomRight
      (\(MatrixShape.Hermitian Order
order sh :: sh0 :+: sh1
sh@(sh0
_sh0:+:sh1
sh1)) ->
         (Order -> sh1 -> Hermitian sh1
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order sh1
sh1, (Order
order,sh0 :+: sh1
sh)))


multiplyVector ::
   (Shape.C sh, Eq sh, Class.Floating a) =>
   Transposition -> Hermitian sh a -> Vector sh a -> Vector sh a
multiplyVector :: Transposition -> Hermitian sh a -> Vector sh a -> Vector sh a
multiplyVector Transposition
transposed
   (Array (MatrixShape.Hermitian Order
order sh
shA) ForeignPtr a
a) (Array sh
shX ForeignPtr a
x) =
      sh -> (Int -> Ptr a -> IO ()) -> Vector sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
shX ((Int -> Ptr a -> IO ()) -> Vector sh a)
-> (Int -> Ptr a -> IO ()) -> Vector sh a
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
yPtr -> do
   String -> Bool -> IO ()
Call.assert String
"Hermitian.multiplyVector: width shapes mismatch" (sh
shA sh -> sh -> Bool
forall a. Eq a => a -> a -> Bool
== sh
shX)
   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 conj :: Conjugation
conj = Order -> Conjugation
conjugatedOnRowMajor (Order -> Conjugation) -> Order -> Conjugation
forall a b. (a -> b) -> a -> b
$ Transposition -> Order -> Order
transposeOrder Transposition
transposed Order
order
      Ptr CChar
uploPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char (Char -> FortranIO () (Ptr CChar))
-> Char -> FortranIO () (Ptr CChar)
forall a b. (a -> b) -> a -> b
$ Order -> Char
uploFromOrder Order
order
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
alphaPtr <- a -> FortranIO () (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
one
      Ptr a
aPtr <- ((Ptr a -> IO ()) -> IO ()) -> FortranIO () (Ptr a)
forall k (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO ()) -> IO ()) -> FortranIO () (Ptr a))
-> ((Ptr a -> IO ()) -> IO ()) -> FortranIO () (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
      Ptr a
xPtr <- Conjugation -> Int -> ForeignPtr a -> FortranIO () (Ptr a)
forall a r.
Floating a =>
Conjugation -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
condConjugateToTemp Conjugation
conj Int
n ForeignPtr a
x
      Ptr CInt
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      Ptr a
betaPtr <- a -> FortranIO () (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
zero
      Ptr CInt
incyPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      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
         Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
BlasGen.hpmv
            Ptr CChar
uploPtr Ptr CInt
nPtr Ptr a
alphaPtr Ptr a
aPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
betaPtr Ptr a
yPtr Ptr CInt
incyPtr
         Conjugation -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Conjugation -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
condConjugate Conjugation
conj Ptr CInt
nPtr Ptr a
yPtr Ptr CInt
incyPtr


square :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a
square :: Hermitian sh a -> Hermitian sh a
square (Array shape :: Hermitian sh
shape@(MatrixShape.Hermitian Order
order sh
sh) ForeignPtr a
a) =
   Hermitian sh -> (Ptr a -> IO ()) -> Hermitian sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate Hermitian sh
shape ((Ptr a -> IO ()) -> Hermitian sh a)
-> (Ptr a -> IO ()) -> Hermitian sh a
forall a b. (a -> b) -> a -> b
$
      Conjugation -> Order -> Int -> ForeignPtr a -> Ptr a -> IO ()
forall a.
Floating a =>
Conjugation -> Order -> Int -> ForeignPtr a -> Ptr a -> IO ()
Symmetric.square Conjugation
Conjugated Order
order (sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh) ForeignPtr a
a

{-
Requires frequent unpacking and packing of triangles.
-}
power ::
   (Shape.C sh, Class.Floating a) => Integer -> Hermitian sh a -> Hermitian sh a
power :: Integer -> Hermitian sh a -> Hermitian sh a
power Integer
n a0 :: Hermitian sh a
a0@(Array (MatrixShape.Hermitian Order
order sh
sh) ForeignPtr a
_) =
   Hermitian (Unchecked sh) a -> Hermitian sh a
forall sh a. Hermitian (Unchecked sh) a -> Hermitian sh a
recheck (Hermitian (Unchecked sh) a -> Hermitian sh a)
-> Hermitian (Unchecked sh) a -> Hermitian sh a
forall a b. (a -> b) -> a -> b
$
   (Hermitian (Unchecked sh) a
 -> Hermitian (Unchecked sh) a -> Hermitian (Unchecked sh) a)
-> Hermitian (Unchecked sh) a
-> Hermitian (Unchecked sh) a
-> Integer
-> Hermitian (Unchecked sh) a
forall a. (a -> a -> a) -> a -> a -> Integer -> a
powerAssociative
      (\Hermitian (Unchecked sh) a
a Hermitian (Unchecked sh) a
b -> Full Small Small (Unchecked sh) (Unchecked sh) a
-> Hermitian (Unchecked sh) a
forall vert height width a.
(C vert, C height, C width, Floating a) =>
Full vert Small height width a -> Hermitian width a
fromUpperPart (Full Small Small (Unchecked sh) (Unchecked sh) a
 -> Hermitian (Unchecked sh) a)
-> Full Small Small (Unchecked sh) (Unchecked sh) a
-> Hermitian (Unchecked sh) a
forall a b. (a -> b) -> a -> b
$ Transposition
-> Hermitian (Unchecked sh) a
-> Full Small Small (Unchecked sh) (Unchecked sh) a
-> Full Small Small (Unchecked sh) (Unchecked sh) a
forall vert horiz height width a.
(C vert, C horiz, C height, Eq height, C width, Floating a) =>
Transposition
-> Hermitian height a
-> Full vert horiz height width a
-> Full vert horiz height width a
multiplyFull Transposition
NonTransposed Hermitian (Unchecked sh) a
a (Full Small Small (Unchecked sh) (Unchecked sh) a
 -> Full Small Small (Unchecked sh) (Unchecked sh) a)
-> Full Small Small (Unchecked sh) (Unchecked sh) a
-> Full Small Small (Unchecked sh) (Unchecked sh) a
forall a b. (a -> b) -> a -> b
$ Hermitian (Unchecked sh) a
-> Full Small Small (Unchecked sh) (Unchecked sh) a
forall sh a. (C sh, Floating a) => Hermitian sh a -> Square sh a
toSquare Hermitian (Unchecked sh) a
b)
      (Order -> Unchecked sh -> Hermitian (Unchecked sh) a
forall sh a. (C sh, Floating a) => Order -> sh -> Hermitian sh a
identity Order
order (Unchecked sh -> Hermitian (Unchecked sh) a)
-> Unchecked sh -> Hermitian (Unchecked sh) a
forall a b. (a -> b) -> a -> b
$ sh -> Unchecked sh
forall sh. sh -> Unchecked sh
Unchecked sh
sh)
      (Hermitian sh a -> Hermitian (Unchecked sh) a
forall sh a. Hermitian sh a -> Hermitian (Unchecked sh) a
uncheck Hermitian sh a
a0)
      Integer
n


multiplyFull ::
   (Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width,
    Class.Floating a) =>
   Transposition -> Hermitian height a ->
   Full vert horiz height width a ->
   Full vert horiz height width a
multiplyFull :: Transposition
-> Hermitian height a
-> Full vert horiz height width a
-> Full vert horiz height width a
multiplyFull Transposition
transposed
   (Array        (MatrixShape.Hermitian Order
orderA height
shA) ForeignPtr a
a)
   (Array shapeB :: Full vert horiz height width
shapeB@(MatrixShape.Full Order
orderB Extent vert horiz height width
extentB) ForeignPtr a
b) =
      Full vert horiz height width
-> (Ptr a -> IO ()) -> Full vert horiz height width a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate Full vert horiz height width
shapeB ((Ptr a -> IO ()) -> Full vert horiz height width a)
-> (Ptr a -> IO ()) -> Full vert horiz height width a
forall a b. (a -> b) -> a -> b
$ \Ptr a
cPtr -> do
   let (height
height,width
width) = Extent vert horiz height width -> (height, width)
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> (height, width)
Extent.dimensions Extent vert horiz height width
extentB
   String -> Bool -> IO ()
Call.assert String
"Hermitian.multiplyFull: shapes mismatch" (height
shA height -> height -> Bool
forall a. Eq a => a -> a -> Bool
== height
height)
   let m0 :: Int
m0 = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
   let n0 :: Int
n0 = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
width
   let size :: Int
size = Int
m0Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
m0
   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 (Char
side,(Int
m,Int
n)) = Order -> (Int, Int) -> (Char, (Int, Int))
forall a. Order -> (a, a) -> (Char, (a, a))
sideSwapFromOrder Order
orderB (Int
m0,Int
n0)
      Ptr CChar
sidePtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
side
      Ptr CChar
uploPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char (Char -> FortranIO () (Ptr CChar))
-> Char -> FortranIO () (Ptr CChar)
forall a b. (a -> b) -> a -> b
$ Order -> Char
uploFromOrder Order
orderA
      Ptr CInt
mPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
alphaPtr <- a -> FortranIO () (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
one
      Ptr a
aPtr <- (Int -> Ptr a -> Ptr a -> IO ())
-> Int -> ForeignPtr a -> FortranIO () (Ptr a)
forall a r.
Storable a =>
(Int -> Ptr a -> Ptr a -> IO ())
-> Int -> ForeignPtr a -> ContT r IO (Ptr a)
unpackToTemp (Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
unpack Order
orderA) Int
m0 ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
m0
      Ptr CInt
incaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      Ptr CInt
sizePtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
size
      Ptr a
bPtr <- ((Ptr a -> IO ()) -> IO ()) -> FortranIO () (Ptr a)
forall k (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO ()) -> IO ()) -> FortranIO () (Ptr a))
-> ((Ptr a -> IO ()) -> IO ()) -> FortranIO () (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
b
      Ptr CInt
ldbPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
m
      Ptr a
betaPtr <- a -> FortranIO () (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
zero
      Ptr CInt
ldcPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
m
      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
         Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Transposition -> Order -> Order
transposeOrder Transposition
transposed Order
orderA Order -> Order -> Bool
forall a. Eq a => a -> a -> Bool
/= Order
orderB) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$
            Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
lacgv Ptr CInt
sizePtr Ptr a
aPtr Ptr CInt
incaPtr
         Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
BlasGen.hemm Ptr CChar
sidePtr Ptr CChar
uploPtr
            Ptr CInt
mPtr Ptr CInt
nPtr Ptr a
alphaPtr Ptr a
aPtr Ptr CInt
ldaPtr
            Ptr a
bPtr Ptr CInt
ldbPtr Ptr a
betaPtr Ptr a
cPtr Ptr CInt
ldcPtr



withConjBuffer ::
   (Shape.C sh, Class.Floating a) =>
   Order -> sh -> Int -> Ptr a ->
   (Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ()) -> ContT r IO ()
withConjBuffer :: Order
-> sh
-> Int
-> Ptr a
-> (Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ())
-> ContT r IO ()
withConjBuffer Order
order sh
sh Int
triSize Ptr a
aPtr Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ()
act = do
   Ptr CChar
uploPtr <- Char -> FortranIO r (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char (Char -> FortranIO r (Ptr CChar))
-> Char -> FortranIO r (Ptr CChar)
forall a b. (a -> b) -> a -> b
$ Order -> Char
uploFromOrder Order
order
   Ptr CInt
nPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (Int -> FortranIO r (Ptr CInt)) -> Int -> FortranIO r (Ptr CInt)
forall a b. (a -> b) -> a -> b
$ sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
   Ptr CInt
incxPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
   Ptr CInt
sizePtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
triSize
   IO () -> ContT r IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT r IO ()) -> IO () -> ContT r IO ()
forall a b. (a -> b) -> a -> b
$ do
      a -> Int -> Ptr a -> IO ()
forall a. Floating a => a -> Int -> Ptr a -> IO ()
fill a
forall a. Floating a => a
zero Int
triSize Ptr a
aPtr
      Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ()
act Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
incxPtr
      Conjugation -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Conjugation -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
condConjugate (Order -> Conjugation
conjugatedOnRowMajor Order
order) Ptr CInt
sizePtr Ptr a
aPtr Ptr CInt
incxPtr

outer ::
   (Shape.C sh, Class.Floating a) => Order -> Vector sh a -> Hermitian sh a
outer :: Order -> Vector sh a -> Hermitian sh a
outer Order
order (Array sh
sh ForeignPtr a
x) =
   Hermitian sh -> (Int -> Ptr a -> IO ()) -> Hermitian sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (Order -> sh -> Hermitian sh
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order sh
sh) ((Int -> Ptr a -> IO ()) -> Hermitian sh a)
-> (Int -> Ptr a -> IO ()) -> Hermitian sh a
forall a b. (a -> b) -> a -> b
$
      \Int
triSize Ptr a
aPtr ->
   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 (RealOf a)
alphaPtr <- Ptr a -> ContT () IO (Ptr (RealOf a))
forall a (f :: * -> *) r.
Floating a =>
f a -> ContT r IO (Ptr (RealOf a))
realOneArg Ptr a
aPtr
      Ptr a
xPtr <- ((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
x
      Order
-> sh
-> Int
-> Ptr a
-> (Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ())
-> ContT () IO ()
forall sh a r.
(C sh, Floating a) =>
Order
-> sh
-> Int
-> Ptr a
-> (Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ())
-> ContT r IO ()
withConjBuffer Order
order sh
sh Int
triSize Ptr a
aPtr ((Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ()) -> ContT () IO ())
-> (Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ()) -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
incxPtr ->
         HPR_ a
forall a. Floating a => HPR_ a
hpr Ptr CChar
uploPtr Ptr CInt
nPtr Ptr (RealOf a)
alphaPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
aPtr


sumRank1 ::
   (Shape.C sh, Eq sh, Class.Floating a) =>
   Order -> sh -> [(RealOf a, Vector sh a)] -> Hermitian sh a
sumRank1 :: Order -> sh -> [(RealOf a, Vector sh a)] -> Hermitian sh a
sumRank1 =
   SumRank1 sh a
-> Order -> sh -> [(RealOf a, Vector sh a)] -> Hermitian sh a
forall sh a. SumRank1 sh a -> SumRank1_ sh (RealOf a) a
getSumRank1 (SumRank1 sh a
 -> Order -> sh -> [(RealOf a, Vector sh a)] -> Hermitian sh a)
-> SumRank1 sh a
-> Order
-> sh
-> [(RealOf a, Vector sh a)]
-> Hermitian sh a
forall a b. (a -> b) -> a -> b
$
   SumRank1 sh Float
-> SumRank1 sh Double
-> SumRank1 sh (Complex Float)
-> SumRank1 sh (Complex Double)
-> SumRank1 sh a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (SumRank1_ sh (RealOf Float) Float -> SumRank1 sh Float
forall sh a. SumRank1_ sh (RealOf a) a -> SumRank1 sh a
SumRank1 SumRank1_ sh (RealOf Float) Float
forall sh a ar.
(C sh, Eq sh, Floating a, RealOf a ~ ar, Storable ar) =>
SumRank1_ sh ar a
sumRank1Aux) (SumRank1_ sh (RealOf Double) Double -> SumRank1 sh Double
forall sh a. SumRank1_ sh (RealOf a) a -> SumRank1 sh a
SumRank1 SumRank1_ sh (RealOf Double) Double
forall sh a ar.
(C sh, Eq sh, Floating a, RealOf a ~ ar, Storable ar) =>
SumRank1_ sh ar a
sumRank1Aux)
      (SumRank1_ sh (RealOf (Complex Float)) (Complex Float)
-> SumRank1 sh (Complex Float)
forall sh a. SumRank1_ sh (RealOf a) a -> SumRank1 sh a
SumRank1 SumRank1_ sh (RealOf (Complex Float)) (Complex Float)
forall sh a ar.
(C sh, Eq sh, Floating a, RealOf a ~ ar, Storable ar) =>
SumRank1_ sh ar a
sumRank1Aux) (SumRank1_ sh (RealOf (Complex Double)) (Complex Double)
-> SumRank1 sh (Complex Double)
forall sh a. SumRank1_ sh (RealOf a) a -> SumRank1 sh a
SumRank1 SumRank1_ sh (RealOf (Complex Double)) (Complex Double)
forall sh a ar.
(C sh, Eq sh, Floating a, RealOf a ~ ar, Storable ar) =>
SumRank1_ sh ar a
sumRank1Aux)

type SumRank1_ sh ar a = Order -> sh -> [(ar, Vector sh a)] -> Hermitian sh a

newtype SumRank1 sh a = SumRank1 {SumRank1 sh a -> SumRank1_ sh (RealOf a) a
getSumRank1 :: SumRank1_ sh (RealOf a) a}

sumRank1Aux ::
   (Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>
   SumRank1_ sh ar a
sumRank1Aux :: SumRank1_ sh ar a
sumRank1Aux Order
order sh
sh [(ar, Vector sh a)]
xs =
   Hermitian sh -> (Int -> Ptr a -> IO ()) -> Array (Hermitian sh) a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (Order -> sh -> Hermitian sh
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order sh
sh) ((Int -> Ptr a -> IO ()) -> Array (Hermitian sh) a)
-> (Int -> Ptr a -> IO ()) -> Array (Hermitian sh) a
forall a b. (a -> b) -> a -> b
$
      \Int
triSize Ptr a
aPtr ->
   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 ar
alphaPtr <- FortranIO () (Ptr ar)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      Order
-> sh
-> Int
-> Ptr a
-> (Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ())
-> ContT () IO ()
forall sh a r.
(C sh, Floating a) =>
Order
-> sh
-> Int
-> Ptr a
-> (Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ())
-> ContT r IO ()
withConjBuffer Order
order sh
sh Int
triSize Ptr a
aPtr ((Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ()) -> ContT () IO ())
-> (Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ()) -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
incxPtr ->
         [(ar, Vector sh a)] -> ((ar, Vector sh a) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [(ar, Vector sh a)]
xs (((ar, Vector sh a) -> IO ()) -> IO ())
-> ((ar, Vector sh a) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \(ar
alpha, Array sh
shX ForeignPtr a
x) ->
            ForeignPtr a -> (Ptr a -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
x ((Ptr a -> IO ()) -> IO ()) -> (Ptr a -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> do
               String -> Bool -> IO ()
Call.assert
                  String
"Hermitian.sumRank1: non-matching vector size" (sh
shsh -> sh -> Bool
forall a. Eq a => a -> a -> Bool
==sh
shX)
               Ptr ar -> ar -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr ar
alphaPtr ar
alpha
               HPR_ a
forall a. Floating a => HPR_ a
hpr Ptr CChar
uploPtr Ptr CInt
nPtr Ptr ar
Ptr (RealOf a)
alphaPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
aPtr


type HPR_ a =
   Ptr CChar -> Ptr CInt ->
   Ptr (RealOf a) -> Ptr a -> Ptr CInt -> Ptr a -> IO ()

newtype HPR a = HPR {HPR a -> HPR_ a
getHPR :: HPR_ a}

hpr :: Class.Floating a => HPR_ a
hpr :: HPR_ a
hpr =
   HPR a -> HPR_ a
forall a. HPR a -> HPR_ a
getHPR (HPR a -> HPR_ a) -> HPR a -> HPR_ a
forall a b. (a -> b) -> a -> b
$
   HPR Float
-> HPR Double
-> HPR (Complex Float)
-> HPR (Complex Double)
-> HPR a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (HPR_ Float -> HPR Float
forall a. HPR_ a -> HPR a
HPR HPR_ Float
forall a.
Real a =>
Ptr CChar
-> Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> IO ()
BlasReal.spr) (HPR_ Double -> HPR Double
forall a. HPR_ a -> HPR a
HPR HPR_ Double
forall a.
Real a =>
Ptr CChar
-> Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> IO ()
BlasReal.spr)
      (HPR_ (Complex Float) -> HPR (Complex Float)
forall a. HPR_ a -> HPR a
HPR HPR_ (Complex Float)
forall a.
Real a =>
Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> IO ()
BlasComplex.hpr) (HPR_ (Complex Double) -> HPR (Complex Double)
forall a. HPR_ a -> HPR a
HPR HPR_ (Complex Double)
forall a.
Real a =>
Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> IO ()
BlasComplex.hpr)


sumRank2 ::
   (Shape.C sh, Eq sh, Class.Floating a) =>
   Order -> sh -> [(a, (Vector sh a, Vector sh a))] -> Hermitian sh a
sumRank2 :: Order -> sh -> [(a, (Vector sh a, Vector sh a))] -> Hermitian sh a
sumRank2 Order
order sh
sh [(a, (Vector sh a, Vector sh a))]
xys =
   Hermitian sh -> (Int -> Ptr a -> IO ()) -> Hermitian sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (Order -> sh -> Hermitian sh
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order sh
sh) ((Int -> Ptr a -> IO ()) -> Hermitian sh a)
-> (Int -> Ptr a -> IO ()) -> Hermitian sh a
forall a b. (a -> b) -> a -> b
$
      \Int
triSize Ptr a
aPtr ->
   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 a
alphaPtr <- FortranIO () (Ptr a)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      Order
-> sh
-> Int
-> Ptr a
-> (Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ())
-> ContT () IO ()
forall sh a r.
(C sh, Floating a) =>
Order
-> sh
-> Int
-> Ptr a
-> (Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ())
-> ContT r IO ()
withConjBuffer Order
order sh
sh Int
triSize Ptr a
aPtr ((Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ()) -> ContT () IO ())
-> (Ptr CChar -> Ptr CInt -> Ptr CInt -> IO ()) -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
incPtr ->
         [(a, (Vector sh a, Vector sh a))]
-> ((a, (Vector sh a, Vector sh a)) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [(a, (Vector sh a, Vector sh a))]
xys (((a, (Vector sh a, Vector sh a)) -> IO ()) -> IO ())
-> ((a, (Vector sh a, Vector sh a)) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \(a
alpha, (Array sh
shX ForeignPtr a
x, Array sh
shY ForeignPtr a
y)) ->
            ForeignPtr a -> (Ptr a -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
x ((Ptr a -> IO ()) -> IO ()) -> (Ptr a -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr ->
            ForeignPtr a -> (Ptr a -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
y ((Ptr a -> IO ()) -> IO ()) -> (Ptr a -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr a
yPtr -> do
               String -> Bool -> IO ()
Call.assert
                  String
"Hermitian.sumRank2: non-matching x vector size" (sh
shsh -> sh -> Bool
forall a. Eq a => a -> a -> Bool
==sh
shX)
               String -> Bool -> IO ()
Call.assert
                  String
"Hermitian.sumRank2: non-matching y vector size" (sh
shsh -> sh -> Bool
forall a. Eq a => a -> a -> Bool
==sh
shY)
               Ptr a -> a -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr a
alphaPtr a
alpha
               Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> IO ()
BlasGen.hpr2 Ptr CChar
uploPtr Ptr CInt
nPtr Ptr a
alphaPtr Ptr a
xPtr Ptr CInt
incPtr Ptr a
yPtr Ptr CInt
incPtr Ptr a
aPtr


{-
It is not strictly necessary to keep the 'order'.
It would be neither more complicated nor less efficient
to change the order via the conversion.
-}
toSquare, _toSquare ::
   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Square sh a
_toSquare :: Hermitian sh a -> Square sh a
_toSquare (Array (MatrixShape.Hermitian Order
order sh
sh) ForeignPtr a
a) =
      Square sh -> (Ptr a -> IO ()) -> Square sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
MatrixShape.square Order
order sh
sh) ((Ptr a -> IO ()) -> Square sh a)
-> (Ptr a -> IO ()) -> Square sh a
forall a b. (a -> b) -> a -> b
$ \Ptr a
bPtr ->
   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 n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
      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
      Ptr a
conjPtr <- Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r. Floating a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
conjugateToTemp (Int -> Int
Shape.triangleSize Int
n) 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
$ do
         Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
unpack (Order -> Order
flipOrder Order
order) Int
n Ptr a
conjPtr Ptr a
bPtr -- wrong
         Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
unpack Order
order Int
n Ptr a
aPtr Ptr a
bPtr

toSquare :: Hermitian sh a -> Square sh a
toSquare (Array (MatrixShape.Hermitian Order
order sh
sh) ForeignPtr a
a) =
      Square sh -> (Ptr a -> IO ()) -> Square sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
MatrixShape.square Order
order sh
sh) ((Ptr a -> IO ()) -> Square sh a)
-> (Ptr a -> IO ()) -> Square sh a
forall a b. (a -> b) -> a -> b
$ \Ptr a
bPtr ->
   ForeignPtr a -> (Ptr a -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
a ((Ptr a -> IO ()) -> IO ()) -> (Ptr a -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr a
aPtr ->
      Conjugation -> Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
Conjugation -> Order -> Int -> Ptr a -> Ptr a -> IO ()
Symmetric.unpack Conjugation
Conjugated Order
order (sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh) Ptr a
aPtr Ptr a
bPtr


gramian ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   General height width a -> Hermitian width a
gramian :: General height width a -> Hermitian width a
gramian (Array (MatrixShape.Full Order
order Extent Big Big height width
extent) ForeignPtr a
a) =
   Hermitian width -> (Ptr a -> IO ()) -> Hermitian width a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> width -> Hermitian width
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order (width -> Hermitian width) -> width -> Hermitian width
forall a b. (a -> b) -> a -> b
$ Extent Big Big height width -> width
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> width
Extent.width Extent Big Big height width
extent) ((Ptr a -> IO ()) -> Hermitian width a)
-> (Ptr a -> IO ()) -> Hermitian width a
forall a b. (a -> b) -> a -> b
$
   \Ptr a
bPtr -> Order
-> ForeignPtr a
-> Ptr a
-> ((Int, Int), (Char, Char, Int))
-> IO ()
forall a.
Floating a =>
Order
-> ForeignPtr a
-> Ptr a
-> ((Int, Int), (Char, Char, Int))
-> IO ()
gramianIO Order
order ForeignPtr a
a Ptr a
bPtr (((Int, Int), (Char, Char, Int)) -> IO ())
-> ((Int, Int), (Char, Char, Int)) -> IO ()
forall a b. (a -> b) -> a -> b
$ Order
-> Extent Big Big height width -> ((Int, Int), (Char, Char, Int))
forall vert horiz height width.
(C vert, C horiz, C height, C width) =>
Order
-> Extent vert horiz height width
-> ((Int, Int), (Char, Char, Int))
gramianParameters Order
order Extent Big Big height width
extent

gramianParameters ::
   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>
   Order ->
   Extent.Extent vert horiz height width ->
   ((Int, Int), (Char, Char, Int))
gramianParameters :: Order
-> Extent vert horiz height width
-> ((Int, Int), (Char, Char, Int))
gramianParameters Order
order Extent vert horiz height width
extent =
   let (height
height, width
width) = Extent vert horiz height width -> (height, width)
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> (height, width)
Extent.dimensions Extent vert horiz height width
extent
       n :: Int
n = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
width
       k :: Int
k = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
    in ((Int
n,Int
k),
         case Order
order of
            Order
ColumnMajor -> (Char
'U', Char
'C', Int
k)
            Order
RowMajor -> (Char
'L', Char
'N', Int
n))


gramianAdjoint ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   General height width a -> Hermitian height a
gramianAdjoint :: General height width a -> Hermitian height a
gramianAdjoint (Array (MatrixShape.Full Order
order Extent Big Big height width
extent) ForeignPtr a
a) =
   Hermitian height -> (Ptr a -> IO ()) -> Hermitian height a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> height -> Hermitian height
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order (height -> Hermitian height) -> height -> Hermitian height
forall a b. (a -> b) -> a -> b
$ Extent Big Big height width -> height
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> height
Extent.height Extent Big Big height width
extent) ((Ptr a -> IO ()) -> Hermitian height a)
-> (Ptr a -> IO ()) -> Hermitian height a
forall a b. (a -> b) -> a -> b
$
   \Ptr a
bPtr -> Order
-> ForeignPtr a
-> Ptr a
-> ((Int, Int), (Char, Char, Int))
-> IO ()
forall a.
Floating a =>
Order
-> ForeignPtr a
-> Ptr a
-> ((Int, Int), (Char, Char, Int))
-> IO ()
gramianIO Order
order ForeignPtr a
a Ptr a
bPtr (((Int, Int), (Char, Char, Int)) -> IO ())
-> ((Int, Int), (Char, Char, Int)) -> IO ()
forall a b. (a -> b) -> a -> b
$ Order
-> Extent Big Big height width -> ((Int, Int), (Char, Char, Int))
forall vert horiz height width.
(C vert, C horiz, C height, C width) =>
Order
-> Extent vert horiz height width
-> ((Int, Int), (Char, Char, Int))
gramianAdjointParameters Order
order Extent Big Big height width
extent

gramianAdjointParameters ::
   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>
   Order ->
   Extent.Extent vert horiz height width ->
   ((Int, Int), (Char, Char, Int))
gramianAdjointParameters :: Order
-> Extent vert horiz height width
-> ((Int, Int), (Char, Char, Int))
gramianAdjointParameters Order
order Extent vert horiz height width
extent =
   let (height
height, width
width) = Extent vert horiz height width -> (height, width)
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> (height, width)
Extent.dimensions Extent vert horiz height width
extent
       n :: Int
n = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
       k :: Int
k = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
width
   in ((Int
n,Int
k),
         case Order
order of
            Order
ColumnMajor -> (Char
'U', Char
'N', Int
n)
            Order
RowMajor -> (Char
'L', Char
'C', Int
k))

{-
Another way to unify 'gramian' and 'gramianAdjoint'
would have been this function:

> gramianConjugation ::
>    Conjugation -> General height width a -> Hermitian width a

with

> gramianAdjoint a = gramianConjugation (transpose a)

but I would like to have

> order (gramianAdjoint a) = order a
-}
gramianIO ::
   (Class.Floating a) =>
   Order ->
   ForeignPtr a -> Ptr a ->
   ((Int, Int), (Char, Char, Int)) -> IO ()
gramianIO :: Order
-> ForeignPtr a
-> Ptr a
-> ((Int, Int), (Char, Char, Int))
-> IO ()
gramianIO Order
order ForeignPtr a
a Ptr a
bPtr ((Int
n,Int
k), (Char
uplo,Char
trans,Int
lda)) =
   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
uploPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
uplo
      Ptr CChar
transPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
trans
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
kPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
k
      Ptr (RealOf a)
alphaPtr <- ForeignPtr a -> ContT () IO (Ptr (RealOf a))
forall a (f :: * -> *) r.
Floating a =>
f a -> ContT r IO (Ptr (RealOf a))
realOneArg ForeignPtr a
a
      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
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
lda
      Ptr (RealOf a)
betaPtr <- ForeignPtr a -> ContT () IO (Ptr (RealOf a))
forall a (f :: * -> *) r.
Floating a =>
f a -> ContT r IO (Ptr (RealOf a))
realZeroArg ForeignPtr a
a
      Ptr a
cPtr <- Int -> ContT () 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
ldcPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim 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
$ do
         HERK_ a
forall a. Floating a => HERK_ a
herk Ptr CChar
uploPtr Ptr CChar
transPtr
            Ptr CInt
nPtr Ptr CInt
kPtr Ptr (RealOf a)
alphaPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr (RealOf a)
betaPtr Ptr a
cPtr Ptr CInt
ldcPtr
         Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
pack Order
order Int
n Ptr a
cPtr Ptr a
bPtr


type HERK_ a =
   Ptr CChar -> Ptr CChar -> Ptr CInt -> Ptr CInt -> Ptr (RealOf a) -> Ptr a ->
   Ptr CInt -> Ptr (RealOf a) -> Ptr a -> Ptr CInt -> IO ()

newtype HERK a = HERK {HERK a -> HERK_ a
getHERK :: HERK_ a}

herk :: Class.Floating a => HERK_ a
herk :: HERK_ a
herk =
   HERK a -> HERK_ a
forall a. HERK a -> HERK_ a
getHERK (HERK a -> HERK_ a) -> HERK a -> HERK_ a
forall a b. (a -> b) -> a -> b
$
   HERK Float
-> HERK Double
-> HERK (Complex Float)
-> HERK (Complex Double)
-> HERK a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (HERK_ Float -> HERK Float
forall a. HERK_ a -> HERK a
HERK HERK_ Float
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
BlasReal.syrk)
      (HERK_ Double -> HERK Double
forall a. HERK_ a -> HERK a
HERK HERK_ Double
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
BlasReal.syrk)
      (HERK_ (Complex Float) -> HERK (Complex Float)
forall a. HERK_ a -> HERK a
HERK HERK_ (Complex Float)
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> IO ()
BlasComplex.herk)
      (HERK_ (Complex Double) -> HERK (Complex Double)
forall a. HERK_ a -> HERK a
HERK HERK_ (Complex Double)
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> IO ()
BlasComplex.herk)


skipCheckCongruence ::
   ((sh -> Unchecked sh) -> matrix0 -> matrix1) ->
   (matrix1 -> Hermitian (Unchecked sh) a) -> matrix0 -> Hermitian sh a
skipCheckCongruence :: ((sh -> Unchecked sh) -> matrix0 -> matrix1)
-> (matrix1 -> Hermitian (Unchecked sh) a)
-> matrix0
-> Hermitian sh a
skipCheckCongruence (sh -> Unchecked sh) -> matrix0 -> matrix1
mapSize matrix1 -> Hermitian (Unchecked sh) a
f matrix0
a =
   Hermitian (Unchecked sh) a -> Hermitian sh a
forall sh a. Hermitian (Unchecked sh) a -> Hermitian sh a
recheck (Hermitian (Unchecked sh) a -> Hermitian sh a)
-> Hermitian (Unchecked sh) a -> Hermitian sh a
forall a b. (a -> b) -> a -> b
$ matrix1 -> Hermitian (Unchecked sh) a
f (matrix1 -> Hermitian (Unchecked sh) a)
-> matrix1 -> Hermitian (Unchecked sh) a
forall a b. (a -> b) -> a -> b
$ (sh -> Unchecked sh) -> matrix0 -> matrix1
mapSize sh -> Unchecked sh
forall sh. sh -> Unchecked sh
Unchecked matrix0
a


congruenceDiagonal ::
   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Vector height (RealOf a) -> General height width a -> Hermitian width a
congruenceDiagonal :: Vector height (RealOf a)
-> General height width a -> Hermitian width a
congruenceDiagonal Vector height (RealOf a)
d =
   ((width -> Unchecked width)
 -> General height width a
 -> Full Big Big height (Unchecked width) a)
-> (Full Big Big height (Unchecked width) a
    -> Hermitian (Unchecked width) a)
-> General height width a
-> Hermitian width a
forall sh matrix0 matrix1 a.
((sh -> Unchecked sh) -> matrix0 -> matrix1)
-> (matrix1 -> Hermitian (Unchecked sh) a)
-> matrix0
-> Hermitian sh a
skipCheckCongruence (width -> Unchecked width)
-> General height width a
-> Full Big Big height (Unchecked width) a
forall vert horiz widthA widthB height a.
(GeneralTallWide vert horiz, GeneralTallWide horiz vert) =>
(widthA -> widthB)
-> Full vert horiz height widthA a
-> Full vert horiz height widthB a
Basic.mapWidth ((Full Big Big height (Unchecked width) a
  -> Hermitian (Unchecked width) a)
 -> General height width a -> Hermitian width a)
-> (Full Big Big height (Unchecked width) a
    -> Hermitian (Unchecked width) a)
-> General height width a
-> Hermitian width a
forall a b. (a -> b) -> a -> b
$ \Full Big Big height (Unchecked width) a
a ->
      a
-> Full Big Big height (Unchecked width) a
-> Full Big Big height (Unchecked width) a
-> Hermitian (Unchecked width) a
forall vert horiz height width a.
(C vert, C horiz, C height, Eq height, C width, Eq width,
 Floating a) =>
a
-> Full vert horiz height width a
-> Full vert horiz height width a
-> Hermitian width a
scaledAnticommutator a
0.5 Full Big Big height (Unchecked width) a
a (Full Big Big height (Unchecked width) a
 -> Hermitian (Unchecked width) a)
-> Full Big Big height (Unchecked width) a
-> Hermitian (Unchecked width) a
forall a b. (a -> b) -> a -> b
$ Vector height (RealOf a)
-> Full Big Big height (Unchecked width) a
-> Full Big Big height (Unchecked width) a
forall vert horiz height width a.
(C vert, C horiz, C height, Eq height, C width, Floating a) =>
Vector height (RealOf a)
-> Full vert horiz height width a -> Full vert horiz height width a
Basic.scaleRowsReal Vector height (RealOf a)
d Full Big Big height (Unchecked width) a
a

congruenceDiagonalAdjoint ::
   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
   General height width a -> Vector width (RealOf a) -> Hermitian height a
congruenceDiagonalAdjoint :: General height width a
-> Vector width (RealOf a) -> Hermitian height a
congruenceDiagonalAdjoint =
   (Vector width (RealOf a)
 -> General height width a -> Hermitian height a)
-> General height width a
-> Vector width (RealOf a)
-> Hermitian height a
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Vector width (RealOf a)
  -> General height width a -> Hermitian height a)
 -> General height width a
 -> Vector width (RealOf a)
 -> Hermitian height a)
-> (Vector width (RealOf a)
    -> General height width a -> Hermitian height a)
-> General height width a
-> Vector width (RealOf a)
-> Hermitian height a
forall a b. (a -> b) -> a -> b
$ \Vector width (RealOf a)
d -> ((height -> Unchecked height)
 -> General height width a
 -> Full Big Big (Unchecked height) width a)
-> (Full Big Big (Unchecked height) width a
    -> Hermitian (Unchecked height) a)
-> General height width a
-> Hermitian height a
forall sh matrix0 matrix1 a.
((sh -> Unchecked sh) -> matrix0 -> matrix1)
-> (matrix1 -> Hermitian (Unchecked sh) a)
-> matrix0
-> Hermitian sh a
skipCheckCongruence (height -> Unchecked height)
-> General height width a
-> Full Big Big (Unchecked height) width a
forall vert horiz heightA heightB width a.
(GeneralTallWide vert horiz, GeneralTallWide horiz vert) =>
(heightA -> heightB)
-> Full vert horiz heightA width a
-> Full vert horiz heightB width a
Basic.mapHeight ((Full Big Big (Unchecked height) width a
  -> Hermitian (Unchecked height) a)
 -> General height width a -> Hermitian height a)
-> (Full Big Big (Unchecked height) width a
    -> Hermitian (Unchecked height) a)
-> General height width a
-> Hermitian height a
forall a b. (a -> b) -> a -> b
$ \Full Big Big (Unchecked height) width a
a ->
      a
-> Full Big Big (Unchecked height) width a
-> Full Big Big (Unchecked height) width a
-> Hermitian (Unchecked height) a
forall vert horiz height width a.
(C vert, C horiz, C height, Eq height, C width, Eq width,
 Floating a) =>
a
-> Full vert horiz height width a
-> Full vert horiz height width a
-> Hermitian height a
scaledAnticommutatorAdjoint a
0.5 Full Big Big (Unchecked height) width a
a (Full Big Big (Unchecked height) width a
 -> Hermitian (Unchecked height) a)
-> Full Big Big (Unchecked height) width a
-> Hermitian (Unchecked height) a
forall a b. (a -> b) -> a -> b
$ Vector width (RealOf a)
-> Full Big Big (Unchecked height) width a
-> Full Big Big (Unchecked height) width a
forall vert horiz height width a.
(C vert, C horiz, C height, C width, Eq width, Floating a) =>
Vector width (RealOf a)
-> Full vert horiz height width a -> Full vert horiz height width a
Basic.scaleColumnsReal Vector width (RealOf a)
d Full Big Big (Unchecked height) width a
a


congruence ::
   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Hermitian height a -> General height width a -> Hermitian width a
congruence :: Hermitian height a -> General height width a -> Hermitian width a
congruence Hermitian height a
b =
   ((width -> Unchecked width)
 -> General height width a
 -> Full Big Big height (Unchecked width) a)
-> (Full Big Big height (Unchecked width) a
    -> Hermitian (Unchecked width) a)
-> General height width a
-> Hermitian width a
forall sh matrix0 matrix1 a.
((sh -> Unchecked sh) -> matrix0 -> matrix1)
-> (matrix1 -> Hermitian (Unchecked sh) a)
-> matrix0
-> Hermitian sh a
skipCheckCongruence (width -> Unchecked width)
-> General height width a
-> Full Big Big height (Unchecked width) a
forall vert horiz widthA widthB height a.
(GeneralTallWide vert horiz, GeneralTallWide horiz vert) =>
(widthA -> widthB)
-> Full vert horiz height widthA a
-> Full vert horiz height widthB a
Basic.mapWidth ((Full Big Big height (Unchecked width) a
  -> Hermitian (Unchecked width) a)
 -> General height width a -> Hermitian width a)
-> (Full Big Big height (Unchecked width) a
    -> Hermitian (Unchecked width) a)
-> General height width a
-> Hermitian width a
forall a b. (a -> b) -> a -> b
$ \Full Big Big height (Unchecked width) a
a ->
      a
-> Full Big Big height (Unchecked width) a
-> Full Big Big height (Unchecked width) a
-> Hermitian (Unchecked width) a
forall vert horiz height width a.
(C vert, C horiz, C height, Eq height, C width, Eq width,
 Floating a) =>
a
-> Full vert horiz height width a
-> Full vert horiz height width a
-> Hermitian width a
scaledAnticommutator a
forall a. Floating a => a
one Full Big Big height (Unchecked width) a
a (Full Big Big height (Unchecked width) a
 -> Hermitian (Unchecked width) a)
-> Full Big Big height (Unchecked width) a
-> Hermitian (Unchecked width) a
forall a b. (a -> b) -> a -> b
$
      Transposition
-> Split Corrupt Small Small height height a
-> Full Big Big height (Unchecked width) a
-> Full Big Big height (Unchecked width) a
forall vertA vert horiz height heightA widthB a lower.
(C vertA, C vert, C horiz, C height, Eq height, C heightA,
 C widthB, Floating a) =>
Transposition
-> Split lower vertA Small heightA height a
-> Full vert horiz height widthB a
-> Full vert horiz height widthB a
Split.tallMultiplyR Transposition
NonTransposed
         ((Hermitian height -> Order)
-> Hermitian height a -> Split Corrupt Small Small height height a
forall symShape sh a.
(Box symShape, HeightOf symShape ~ sh, C sh, Floating a) =>
(symShape -> Order) -> Array symShape a -> Square Corrupt sh a
Split.takeHalf Hermitian height -> Order
forall size. Hermitian size -> Order
MatrixShape.hermitianOrder Hermitian height a
b) Full Big Big height (Unchecked width) a
a

congruenceAdjoint ::
   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
   General height width a -> Hermitian width a -> Hermitian height a
congruenceAdjoint :: General height width a -> Hermitian width a -> Hermitian height a
congruenceAdjoint =
   (Hermitian width a -> General height width a -> Hermitian height a)
-> General height width a
-> Hermitian width a
-> Hermitian height a
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Hermitian width a
  -> General height width a -> Hermitian height a)
 -> General height width a
 -> Hermitian width a
 -> Hermitian height a)
-> (Hermitian width a
    -> General height width a -> Hermitian height a)
-> General height width a
-> Hermitian width a
-> Hermitian height a
forall a b. (a -> b) -> a -> b
$ \Hermitian width a
b -> ((height -> Unchecked height)
 -> General height width a
 -> Full Big Big (Unchecked height) width a)
-> (Full Big Big (Unchecked height) width a
    -> Hermitian (Unchecked height) a)
-> General height width a
-> Hermitian height a
forall sh matrix0 matrix1 a.
((sh -> Unchecked sh) -> matrix0 -> matrix1)
-> (matrix1 -> Hermitian (Unchecked sh) a)
-> matrix0
-> Hermitian sh a
skipCheckCongruence (height -> Unchecked height)
-> General height width a
-> Full Big Big (Unchecked height) width a
forall vert horiz heightA heightB width a.
(GeneralTallWide vert horiz, GeneralTallWide horiz vert) =>
(heightA -> heightB)
-> Full vert horiz heightA width a
-> Full vert horiz heightB width a
Basic.mapHeight ((Full Big Big (Unchecked height) width a
  -> Hermitian (Unchecked height) a)
 -> General height width a -> Hermitian height a)
-> (Full Big Big (Unchecked height) width a
    -> Hermitian (Unchecked height) a)
-> General height width a
-> Hermitian height a
forall a b. (a -> b) -> a -> b
$ \Full Big Big (Unchecked height) width a
a ->
      a
-> Full Big Big (Unchecked height) width a
-> Full Big Big (Unchecked height) width a
-> Hermitian (Unchecked height) a
forall vert horiz height width a.
(C vert, C horiz, C height, Eq height, C width, Eq width,
 Floating a) =>
a
-> Full vert horiz height width a
-> Full vert horiz height width a
-> Hermitian height a
scaledAnticommutatorAdjoint a
forall a. Floating a => a
one Full Big Big (Unchecked height) width a
a (Full Big Big (Unchecked height) width a
 -> Hermitian (Unchecked height) a)
-> Full Big Big (Unchecked height) width a
-> Hermitian (Unchecked height) a
forall a b. (a -> b) -> a -> b
$
      (Split Corrupt Small Small width width a
 -> Full Big Big width (Unchecked height) a
 -> Full Big Big width (Unchecked height) a)
-> Full Big Big (Unchecked height) width a
-> Split Corrupt Small Small width width a
-> Full Big Big (Unchecked height) width a
forall vertA vertB horizA horizB matrix widthA heightA a widthB
       heightB.
(C vertA, C vertB, C horizA, C horizB) =>
(matrix
 -> Full horizA vertA widthA heightA a
 -> Full horizB vertB widthB heightB a)
-> Full vertA horizA heightA widthA a
-> matrix
-> Full vertB horizB heightB widthB a
Basic.swapMultiply (Transposition
-> Split Corrupt Small Small width width a
-> Full Big Big width (Unchecked height) a
-> Full Big Big width (Unchecked height) a
forall vertA vert horiz height heightA widthB a lower.
(C vertA, C vert, C horiz, C height, Eq height, C heightA,
 C widthB, Floating a) =>
Transposition
-> Split lower vertA Small heightA height a
-> Full vert horiz height widthB a
-> Full vert horiz height widthB a
Split.tallMultiplyR Transposition
Transposed)
         Full Big Big (Unchecked height) width a
a ((Hermitian width -> Order)
-> Hermitian width a -> Split Corrupt Small Small width width a
forall symShape sh a.
(Box symShape, HeightOf symShape ~ sh, C sh, Floating a) =>
(symShape -> Order) -> Array symShape a -> Square Corrupt sh a
Split.takeHalf Hermitian width -> Order
forall size. Hermitian size -> Order
MatrixShape.hermitianOrder Hermitian width a
b)


scaledAnticommutator ::
   (Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
   a ->
   Full vert horiz height width a ->
   Full vert horiz height width a -> Hermitian width a
scaledAnticommutator :: a
-> Full vert horiz height width a
-> Full vert horiz height width a
-> Hermitian width a
scaledAnticommutator a
alpha Full vert horiz height width a
arr (Array (MatrixShape.Full Order
order Extent vert horiz height width
extentB) ForeignPtr a
b) = do
   let (Array (MatrixShape.Full Order
_ Extent vert horiz height width
extentA) ForeignPtr a
a) = Order
-> Full vert horiz height width a -> Full vert horiz height width a
forall vert horiz height width a.
(C vert, C horiz, C height, C width, Floating a) =>
Order
-> Full vert horiz height width a -> Full vert horiz height width a
Basic.forceOrder Order
order Full vert horiz height width a
arr
   Hermitian width -> (Ptr a -> IO ()) -> Hermitian width a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> width -> Hermitian width
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order (width -> Hermitian width) -> width -> Hermitian width
forall a b. (a -> b) -> a -> b
$ Extent vert horiz height width -> width
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> width
Extent.width Extent vert horiz height width
extentB) ((Ptr a -> IO ()) -> Hermitian width a)
-> (Ptr a -> IO ()) -> Hermitian width a
forall a b. (a -> b) -> a -> b
$
         \Ptr a
cpPtr -> do
      String -> Bool -> IO ()
Call.assert String
"Hermitian.anticommutator: extents mismatch"
         (Extent vert horiz height width
extentAExtent vert horiz height width
-> Extent vert horiz height width -> Bool
forall a. Eq a => a -> a -> Bool
==Extent vert horiz height width
extentB)
      a
-> Order
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> ((Int, Int), (Char, Char, Int))
-> IO ()
forall a.
Floating a =>
a
-> Order
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> ((Int, Int), (Char, Char, Int))
-> IO ()
scaledAnticommutatorIO a
alpha Order
order ForeignPtr a
a ForeignPtr a
b Ptr a
cpPtr (((Int, Int), (Char, Char, Int)) -> IO ())
-> ((Int, Int), (Char, Char, Int)) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Order
-> Extent vert horiz height width
-> ((Int, Int), (Char, Char, Int))
forall vert horiz height width.
(C vert, C horiz, C height, C width) =>
Order
-> Extent vert horiz height width
-> ((Int, Int), (Char, Char, Int))
gramianParameters Order
order Extent vert horiz height width
extentB

scaledAnticommutatorAdjoint ::
   (Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
   a ->
   Full vert horiz height width a ->
   Full vert horiz height width a -> Hermitian height a
scaledAnticommutatorAdjoint :: a
-> Full vert horiz height width a
-> Full vert horiz height width a
-> Hermitian height a
scaledAnticommutatorAdjoint
      a
alpha Full vert horiz height width a
arr (Array (MatrixShape.Full Order
order Extent vert horiz height width
extentB) ForeignPtr a
b) = do
   let (Array (MatrixShape.Full Order
_ Extent vert horiz height width
extentA) ForeignPtr a
a) = Order
-> Full vert horiz height width a -> Full vert horiz height width a
forall vert horiz height width a.
(C vert, C horiz, C height, C width, Floating a) =>
Order
-> Full vert horiz height width a -> Full vert horiz height width a
Basic.forceOrder Order
order Full vert horiz height width a
arr
   Hermitian height -> (Ptr a -> IO ()) -> Hermitian height a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> height -> Hermitian height
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order (height -> Hermitian height) -> height -> Hermitian height
forall a b. (a -> b) -> a -> b
$ Extent vert horiz height width -> height
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> height
Extent.height Extent vert horiz height width
extentB) ((Ptr a -> IO ()) -> Hermitian height a)
-> (Ptr a -> IO ()) -> Hermitian height a
forall a b. (a -> b) -> a -> b
$
         \Ptr a
cpPtr -> do
      String -> Bool -> IO ()
Call.assert String
"Hermitian.anticommutatorAdjoint: extents mismatch"
         (Extent vert horiz height width
extentAExtent vert horiz height width
-> Extent vert horiz height width -> Bool
forall a. Eq a => a -> a -> Bool
==Extent vert horiz height width
extentB)
      a
-> Order
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> ((Int, Int), (Char, Char, Int))
-> IO ()
forall a.
Floating a =>
a
-> Order
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> ((Int, Int), (Char, Char, Int))
-> IO ()
scaledAnticommutatorIO a
alpha Order
order ForeignPtr a
a ForeignPtr a
b Ptr a
cpPtr (((Int, Int), (Char, Char, Int)) -> IO ())
-> ((Int, Int), (Char, Char, Int)) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Order
-> Extent vert horiz height width
-> ((Int, Int), (Char, Char, Int))
forall vert horiz height width.
(C vert, C horiz, C height, C width) =>
Order
-> Extent vert horiz height width
-> ((Int, Int), (Char, Char, Int))
gramianAdjointParameters Order
order Extent vert horiz height width
extentB

scaledAnticommutatorIO ::
   (Class.Floating a) =>
   a ->
   Order -> ForeignPtr a -> ForeignPtr a -> Ptr a ->
   ((Int, Int), (Char, Char, Int)) -> IO ()
scaledAnticommutatorIO :: a
-> Order
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> ((Int, Int), (Char, Char, Int))
-> IO ()
scaledAnticommutatorIO a
alpha Order
order ForeignPtr a
a ForeignPtr a
b Ptr a
cpPtr ((Int
n,Int
k), (Char
uplo,Char
trans,Int
lda)) =
   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
uploPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
uplo
      Ptr CChar
transPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
trans
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
kPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
k
      Ptr a
alphaPtr <- a -> FortranIO () (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
alpha
      Ptr a
aPtr <- ((Ptr a -> IO ()) -> IO ()) -> FortranIO () (Ptr a)
forall k (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO ()) -> IO ()) -> FortranIO () (Ptr a))
-> ((Ptr a -> IO ()) -> IO ()) -> FortranIO () (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
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
lda
      Ptr a
bPtr <- ((Ptr a -> IO ()) -> IO ()) -> FortranIO () (Ptr a)
forall k (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO ()) -> IO ()) -> FortranIO () (Ptr a))
-> ((Ptr a -> IO ()) -> IO ()) -> FortranIO () (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
b
      let ldbPtr :: Ptr CInt
ldbPtr = Ptr CInt
ldaPtr
      Ptr (RealOf a)
betaPtr <- Ptr a -> ContT () IO (Ptr (RealOf a))
forall a (f :: * -> *) r.
Floating a =>
f a -> ContT r IO (Ptr (RealOf a))
realZeroArg Ptr a
aPtr
      Ptr a
cPtr <- Int -> FortranIO () (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
ldcPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim 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
$ do
         HER2K_ a
forall a. Floating a => HER2K_ a
her2k Ptr CChar
uploPtr Ptr CChar
transPtr Ptr CInt
nPtr Ptr CInt
kPtr Ptr a
alphaPtr
            Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
bPtr Ptr CInt
ldbPtr Ptr (RealOf a)
betaPtr Ptr a
cPtr Ptr CInt
ldcPtr
         Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
pack Order
order Int
n Ptr a
cPtr Ptr a
cpPtr


type HER2K_ a =
   Ptr CChar -> Ptr CChar -> Ptr CInt -> Ptr CInt -> Ptr a ->
   Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->
   Ptr (RealOf a) -> Ptr a -> Ptr CInt -> IO ()

newtype HER2K a = HER2K {HER2K a -> HER2K_ a
getHER2K :: HER2K_ a}

her2k :: Class.Floating a => HER2K_ a
her2k :: HER2K_ a
her2k =
   HER2K a -> HER2K_ a
forall a. HER2K a -> HER2K_ a
getHER2K (HER2K a -> HER2K_ a) -> HER2K a -> HER2K_ a
forall a b. (a -> b) -> a -> b
$
   HER2K Float
-> HER2K Double
-> HER2K (Complex Float)
-> HER2K (Complex Double)
-> HER2K a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (HER2K_ Float -> HER2K Float
forall a. HER2K_ a -> HER2K a
HER2K HER2K_ Float
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
BlasReal.syr2k)
      (HER2K_ Double -> HER2K Double
forall a. HER2K_ a -> HER2K a
HER2K HER2K_ Double
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
BlasReal.syr2k)
      (HER2K_ (Complex Float) -> HER2K (Complex Float)
forall a. HER2K_ a -> HER2K a
HER2K HER2K_ (Complex Float)
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> IO ()
BlasComplex.her2k)
      (HER2K_ (Complex Double) -> HER2K (Complex Double)
forall a. HER2K_ a -> HER2K a
HER2K HER2K_ (Complex Double)
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> IO ()
BlasComplex.her2k)


addAdjoint, _addAdjoint ::
   (Shape.C sh, Class.Floating a) => Square sh a -> Hermitian sh a
_addAdjoint :: Square sh a -> Hermitian sh a
_addAdjoint =
   (Order -> sh -> ForeignPtr a -> Hermitian sh a)
-> Square sh a -> Hermitian sh a
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr a -> Hermitian sh a)
 -> Square sh a -> Hermitian sh a)
-> (Order -> sh -> ForeignPtr a -> Hermitian sh a)
-> Square sh a
-> Hermitian sh a
forall a b. (a -> b) -> a -> b
$ \Order
order sh
sh ForeignPtr a
a ->
      Hermitian sh -> (Int -> Ptr a -> IO ()) -> Hermitian sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (Order -> sh -> Hermitian sh
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order sh
sh) ((Int -> Ptr a -> IO ()) -> Hermitian sh a)
-> (Int -> Ptr a -> IO ()) -> Hermitian sh a
forall a b. (a -> b) -> a -> b
$ \Int
bSize Ptr a
bPtr -> do
   let n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
   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 a
alphaPtr <- a -> FortranIO () (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
one
      Ptr CInt
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      Ptr a
aPtr <- ((Ptr a -> IO ()) -> IO ()) -> FortranIO () (Ptr a)
forall k (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO ()) -> IO ()) -> FortranIO () (Ptr a))
-> ((Ptr a -> IO ()) -> IO ()) -> FortranIO () (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
      Ptr CInt
sizePtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
bSize
      Ptr a
conjPtr <- Int -> FortranIO () (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
bSize
      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
         Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
pack Order
order Int
n Ptr a
aPtr Ptr a
bPtr
         Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
pack (Order -> Order
flipOrder Order
order) Int
n Ptr a
aPtr Ptr a
conjPtr -- wrong
         Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
lacgv Ptr CInt
sizePtr Ptr a
conjPtr Ptr CInt
incxPtr
         Ptr CInt
-> Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt
-> Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.axpy Ptr CInt
sizePtr Ptr a
alphaPtr Ptr a
conjPtr Ptr CInt
incxPtr Ptr a
bPtr Ptr CInt
incxPtr

addAdjoint :: Square sh a -> Hermitian sh a
addAdjoint =
   (Order -> sh -> ForeignPtr a -> Hermitian sh a)
-> Square sh a -> Hermitian sh a
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr a -> Hermitian sh a)
 -> Square sh a -> Hermitian sh a)
-> (Order -> sh -> ForeignPtr a -> Hermitian sh a)
-> Square sh a
-> Hermitian sh a
forall a b. (a -> b) -> a -> b
$ \Order
order sh
sh ForeignPtr a
a ->
      Hermitian sh -> (Ptr a -> IO ()) -> Hermitian sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> sh -> Hermitian sh
forall size. Order -> size -> Hermitian size
MatrixShape.Hermitian Order
order sh
sh) ((Ptr a -> IO ()) -> Hermitian sh a)
-> (Ptr a -> IO ()) -> Hermitian sh a
forall a b. (a -> b) -> a -> b
$ \Ptr a
bPtr -> do
   let n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
   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 a
alphaPtr <- a -> FortranIO () (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
one
      Ptr CInt
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      Ptr CInt
incnPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
aPtr <- ((Ptr a -> IO ()) -> IO ()) -> FortranIO () (Ptr a)
forall k (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO ()) -> IO ()) -> FortranIO () (Ptr a))
-> ((Ptr a -> IO ()) -> IO ()) -> FortranIO () (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
$ case Order
order of
         Order
RowMajor ->
            [(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
aPtr Ptr a
bPtr) ((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
srcPtr,Ptr a
dstPtr) -> do
            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 ()
copyConjugate Ptr CInt
nPtr Ptr a
srcPtr Ptr CInt
incnPtr Ptr a
dstPtr Ptr CInt
incxPtr
            Ptr CInt
-> Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt
-> Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.axpy Ptr CInt
nPtr Ptr a
alphaPtr Ptr a
srcPtr Ptr CInt
incxPtr Ptr a
dstPtr Ptr CInt
incxPtr
         Order
ColumnMajor ->
            [(Int, ((Ptr a, Ptr a), Ptr a))]
-> (Ptr CInt -> ((Ptr a, 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), Ptr a))]
forall a.
Storable a =>
Int -> Ptr a -> Ptr a -> [(Int, ((Ptr a, Ptr a), Ptr a))]
columnMajorPointers Int
n Ptr a
aPtr Ptr a
bPtr) ((Ptr CInt -> ((Ptr a, Ptr a), Ptr a) -> IO ()) -> IO ())
-> (Ptr CInt -> ((Ptr a, Ptr a), Ptr a) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
               \Ptr CInt
nPtr ((Ptr a
srcRowPtr,Ptr a
srcColumnPtr),Ptr a
dstPtr) -> do
            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 ()
copyConjugate Ptr CInt
nPtr Ptr a
srcRowPtr Ptr CInt
incnPtr Ptr a
dstPtr Ptr CInt
incxPtr
            Ptr CInt
-> Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt
-> Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.axpy Ptr CInt
nPtr Ptr a
alphaPtr Ptr a
srcColumnPtr Ptr CInt
incxPtr Ptr a
dstPtr Ptr CInt
incxPtr


_pack :: Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
_pack :: Order -> Int -> Ptr a -> Ptr a -> IO ()
_pack Order
order Int
n Ptr a
fullPtr Ptr a
packedPtr =
   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
      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
$
         case Order
order of
            Order
ColumnMajor ->
               [(Int, ((Ptr a, Ptr a), Ptr a))]
-> (Ptr CInt -> ((Ptr a, 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), Ptr a))]
forall a.
Storable a =>
Int -> Ptr a -> Ptr a -> [(Int, ((Ptr a, Ptr a), Ptr a))]
columnMajorPointers Int
n Ptr a
fullPtr Ptr a
packedPtr) ((Ptr CInt -> ((Ptr a, Ptr a), Ptr a) -> IO ()) -> IO ())
-> (Ptr CInt -> ((Ptr a, Ptr a), Ptr a) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
                  \Ptr CInt
nPtr ((Ptr a
_,Ptr a
srcPtr),Ptr a
dstPtr) ->
                     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
incxPtr
            Order
RowMajor ->
               [(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
srcPtr,Ptr a
dstPtr) ->
                     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
incxPtr


realZeroArg, realOneArg ::
   (Class.Floating a) => f a -> ContT r IO (Ptr (RealOf a))
realZeroArg :: f a -> ContT r IO (Ptr (RealOf a))
realZeroArg =
   RealArg f (ContT r IO) a -> f a -> ContT r IO (Ptr (RealOf a))
forall (f :: * -> *) (g :: * -> *) a.
RealArg f g a -> f a -> g (Ptr (RealOf a))
runRealArg (RealArg f (ContT r IO) a -> f a -> ContT r IO (Ptr (RealOf a)))
-> RealArg f (ContT r IO) a -> f a -> ContT r IO (Ptr (RealOf a))
forall a b. (a -> b) -> a -> b
$
   RealArg f (ContT r IO) Float
-> RealArg f (ContT r IO) Double
-> RealArg f (ContT r IO) (Complex Float)
-> RealArg f (ContT r IO) (Complex Double)
-> RealArg f (ContT r IO) a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      ((f Float -> ContT r IO (Ptr (RealOf Float)))
-> RealArg f (ContT r IO) Float
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (Ptr (RealOf a))) -> RealArg f g a
RealArg ((f Float -> ContT r IO (Ptr (RealOf Float)))
 -> RealArg f (ContT r IO) Float)
-> (f Float -> ContT r IO (Ptr (RealOf Float)))
-> RealArg f (ContT r IO) Float
forall a b. (a -> b) -> a -> b
$ FortranIO r (Ptr Float) -> f Float -> FortranIO r (Ptr Float)
forall a b. a -> b -> a
const (FortranIO r (Ptr Float) -> f Float -> FortranIO r (Ptr Float))
-> FortranIO r (Ptr Float) -> f Float -> FortranIO r (Ptr Float)
forall a b. (a -> b) -> a -> b
$ Float -> FortranIO r (Ptr Float)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number Float
forall a. Floating a => a
zero)
      ((f Double -> ContT r IO (Ptr (RealOf Double)))
-> RealArg f (ContT r IO) Double
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (Ptr (RealOf a))) -> RealArg f g a
RealArg ((f Double -> ContT r IO (Ptr (RealOf Double)))
 -> RealArg f (ContT r IO) Double)
-> (f Double -> ContT r IO (Ptr (RealOf Double)))
-> RealArg f (ContT r IO) Double
forall a b. (a -> b) -> a -> b
$ FortranIO r (Ptr Double) -> f Double -> FortranIO r (Ptr Double)
forall a b. a -> b -> a
const (FortranIO r (Ptr Double) -> f Double -> FortranIO r (Ptr Double))
-> FortranIO r (Ptr Double) -> f Double -> FortranIO r (Ptr Double)
forall a b. (a -> b) -> a -> b
$ Double -> FortranIO r (Ptr Double)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number Double
forall a. Floating a => a
zero)
      ((f (Complex Float) -> ContT r IO (Ptr (RealOf (Complex Float))))
-> RealArg f (ContT r IO) (Complex Float)
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (Ptr (RealOf a))) -> RealArg f g a
RealArg ((f (Complex Float) -> ContT r IO (Ptr (RealOf (Complex Float))))
 -> RealArg f (ContT r IO) (Complex Float))
-> (f (Complex Float) -> ContT r IO (Ptr (RealOf (Complex Float))))
-> RealArg f (ContT r IO) (Complex Float)
forall a b. (a -> b) -> a -> b
$ FortranIO r (Ptr Float)
-> f (Complex Float) -> FortranIO r (Ptr Float)
forall a b. a -> b -> a
const (FortranIO r (Ptr Float)
 -> f (Complex Float) -> FortranIO r (Ptr Float))
-> FortranIO r (Ptr Float)
-> f (Complex Float)
-> FortranIO r (Ptr Float)
forall a b. (a -> b) -> a -> b
$ Float -> FortranIO r (Ptr Float)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number Float
forall a. Floating a => a
zero)
      ((f (Complex Double) -> ContT r IO (Ptr (RealOf (Complex Double))))
-> RealArg f (ContT r IO) (Complex Double)
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (Ptr (RealOf a))) -> RealArg f g a
RealArg ((f (Complex Double) -> ContT r IO (Ptr (RealOf (Complex Double))))
 -> RealArg f (ContT r IO) (Complex Double))
-> (f (Complex Double)
    -> ContT r IO (Ptr (RealOf (Complex Double))))
-> RealArg f (ContT r IO) (Complex Double)
forall a b. (a -> b) -> a -> b
$ FortranIO r (Ptr Double)
-> f (Complex Double) -> FortranIO r (Ptr Double)
forall a b. a -> b -> a
const (FortranIO r (Ptr Double)
 -> f (Complex Double) -> FortranIO r (Ptr Double))
-> FortranIO r (Ptr Double)
-> f (Complex Double)
-> FortranIO r (Ptr Double)
forall a b. (a -> b) -> a -> b
$ Double -> FortranIO r (Ptr Double)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number Double
forall a. Floating a => a
zero)
realOneArg :: f a -> ContT r IO (Ptr (RealOf a))
realOneArg =
   RealArg f (ContT r IO) a -> f a -> ContT r IO (Ptr (RealOf a))
forall (f :: * -> *) (g :: * -> *) a.
RealArg f g a -> f a -> g (Ptr (RealOf a))
runRealArg (RealArg f (ContT r IO) a -> f a -> ContT r IO (Ptr (RealOf a)))
-> RealArg f (ContT r IO) a -> f a -> ContT r IO (Ptr (RealOf a))
forall a b. (a -> b) -> a -> b
$
   RealArg f (ContT r IO) Float
-> RealArg f (ContT r IO) Double
-> RealArg f (ContT r IO) (Complex Float)
-> RealArg f (ContT r IO) (Complex Double)
-> RealArg f (ContT r IO) a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      ((f Float -> ContT r IO (Ptr (RealOf Float)))
-> RealArg f (ContT r IO) Float
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (Ptr (RealOf a))) -> RealArg f g a
RealArg ((f Float -> ContT r IO (Ptr (RealOf Float)))
 -> RealArg f (ContT r IO) Float)
-> (f Float -> ContT r IO (Ptr (RealOf Float)))
-> RealArg f (ContT r IO) Float
forall a b. (a -> b) -> a -> b
$ FortranIO r (Ptr Float) -> f Float -> FortranIO r (Ptr Float)
forall a b. a -> b -> a
const (FortranIO r (Ptr Float) -> f Float -> FortranIO r (Ptr Float))
-> FortranIO r (Ptr Float) -> f Float -> FortranIO r (Ptr Float)
forall a b. (a -> b) -> a -> b
$ Float -> FortranIO r (Ptr Float)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number Float
forall a. Floating a => a
one)
      ((f Double -> ContT r IO (Ptr (RealOf Double)))
-> RealArg f (ContT r IO) Double
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (Ptr (RealOf a))) -> RealArg f g a
RealArg ((f Double -> ContT r IO (Ptr (RealOf Double)))
 -> RealArg f (ContT r IO) Double)
-> (f Double -> ContT r IO (Ptr (RealOf Double)))
-> RealArg f (ContT r IO) Double
forall a b. (a -> b) -> a -> b
$ FortranIO r (Ptr Double) -> f Double -> FortranIO r (Ptr Double)
forall a b. a -> b -> a
const (FortranIO r (Ptr Double) -> f Double -> FortranIO r (Ptr Double))
-> FortranIO r (Ptr Double) -> f Double -> FortranIO r (Ptr Double)
forall a b. (a -> b) -> a -> b
$ Double -> FortranIO r (Ptr Double)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number Double
forall a. Floating a => a
one)
      ((f (Complex Float) -> ContT r IO (Ptr (RealOf (Complex Float))))
-> RealArg f (ContT r IO) (Complex Float)
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (Ptr (RealOf a))) -> RealArg f g a
RealArg ((f (Complex Float) -> ContT r IO (Ptr (RealOf (Complex Float))))
 -> RealArg f (ContT r IO) (Complex Float))
-> (f (Complex Float) -> ContT r IO (Ptr (RealOf (Complex Float))))
-> RealArg f (ContT r IO) (Complex Float)
forall a b. (a -> b) -> a -> b
$ FortranIO r (Ptr Float)
-> f (Complex Float) -> FortranIO r (Ptr Float)
forall a b. a -> b -> a
const (FortranIO r (Ptr Float)
 -> f (Complex Float) -> FortranIO r (Ptr Float))
-> FortranIO r (Ptr Float)
-> f (Complex Float)
-> FortranIO r (Ptr Float)
forall a b. (a -> b) -> a -> b
$ Float -> FortranIO r (Ptr Float)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number Float
forall a. Floating a => a
one)
      ((f (Complex Double) -> ContT r IO (Ptr (RealOf (Complex Double))))
-> RealArg f (ContT r IO) (Complex Double)
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (Ptr (RealOf a))) -> RealArg f g a
RealArg ((f (Complex Double) -> ContT r IO (Ptr (RealOf (Complex Double))))
 -> RealArg f (ContT r IO) (Complex Double))
-> (f (Complex Double)
    -> ContT r IO (Ptr (RealOf (Complex Double))))
-> RealArg f (ContT r IO) (Complex Double)
forall a b. (a -> b) -> a -> b
$ FortranIO r (Ptr Double)
-> f (Complex Double) -> FortranIO r (Ptr Double)
forall a b. a -> b -> a
const (FortranIO r (Ptr Double)
 -> f (Complex Double) -> FortranIO r (Ptr Double))
-> FortranIO r (Ptr Double)
-> f (Complex Double)
-> FortranIO r (Ptr Double)
forall a b. (a -> b) -> a -> b
$ Double -> FortranIO r (Ptr Double)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number Double
forall a. Floating a => a
one)

newtype RealArg f g a = RealArg {RealArg f g a -> f a -> g (Ptr (RealOf a))
runRealArg :: f a -> g (Ptr (RealOf a))}