{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
module Numeric.LAPACK.Matrix.Symmetric.Basic (
   Symmetric,

   gramian,              gramianTransposed,
   congruenceDiagonal,   congruenceDiagonalTransposed,
   congruence,           congruenceTransposed,
   scaledAnticommutator, scaledAnticommutatorTransposed,
   ) where

import qualified Numeric.LAPACK.Matrix.Private as Matrix
import qualified Numeric.LAPACK.Matrix.Basic as Basic
import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Split as Split
import Numeric.LAPACK.Matrix.Triangular.Private (pack, recheck)
import Numeric.LAPACK.Matrix.Shape.Private
         (Order(RowMajor,ColumnMajor), NonUnit(NonUnit))
import Numeric.LAPACK.Matrix.Modifier
         (Transposition(NonTransposed, Transposed))
import Numeric.LAPACK.Matrix.Private (Full)
import Numeric.LAPACK.Vector (Vector)
import Numeric.LAPACK.Scalar (zero, one)
import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))

import qualified Numeric.BLAS.FFI.Generic as BlasGen
import qualified Numeric.Netlib.Utility as Call
import qualified Numeric.Netlib.Class as Class

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

import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
import Foreign.Ptr (Ptr)

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


type Symmetric sh = Array (MatrixShape.FlexSymmetric NonUnit sh)


-- cf. Hermitian.Basic
gramian ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   Matrix.General height width a -> Symmetric width a
gramian :: General height width a -> Symmetric width a
gramian (Array (MatrixShape.Full Order
order Extent Big Big height width
extent) ForeignPtr a
a) =
   Symmetric width -> (Ptr a -> IO ()) -> Symmetric width a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> width -> Symmetric width
forall size. Order -> size -> Symmetric size
symmetricShape Order
order (width -> Symmetric width) -> width -> Symmetric 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 ()) -> Symmetric width a)
-> (Ptr a -> IO ()) -> Symmetric 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
'T', Int
k)
            Order
RowMajor -> (Char
'L', Char
'N', Int
n))


gramianTransposed ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   Matrix.General height width a -> Symmetric height a
gramianTransposed :: General height width a -> Symmetric height a
gramianTransposed (Array (MatrixShape.Full Order
order Extent Big Big height width
extent) ForeignPtr a
a) =
   Symmetric height -> (Ptr a -> IO ()) -> Symmetric height a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> height -> Symmetric height
forall size. Order -> size -> Symmetric size
symmetricShape Order
order (height -> Symmetric height) -> height -> Symmetric 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 ()) -> Symmetric height a)
-> (Ptr a -> IO ()) -> Symmetric 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))
gramianTransposedParameters Order
order Extent Big Big height width
extent

gramianTransposedParameters ::
   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>
   Order ->
   Extent.Extent vert horiz height width ->
   ((Int, Int), (Char, Char, Int))
gramianTransposedParameters :: Order
-> Extent vert horiz height width
-> ((Int, Int), (Char, Char, Int))
gramianTransposedParameters 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
'T', Int
k))

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 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 CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
lda
      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 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
         Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> 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 a
-> Ptr CInt
-> IO ()
BlasGen.syrk Ptr CChar
uploPtr Ptr CChar
transPtr
            Ptr CInt
nPtr Ptr CInt
kPtr Ptr a
alphaPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr 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

skipCheckCongruence ::
   ((sh -> Unchecked sh) -> matrix0 -> matrix1) ->
   (matrix1 -> Symmetric (Unchecked sh) a) -> matrix0 -> Symmetric sh a
skipCheckCongruence :: ((sh -> Unchecked sh) -> matrix0 -> matrix1)
-> (matrix1 -> Symmetric (Unchecked sh) a)
-> matrix0
-> Symmetric sh a
skipCheckCongruence (sh -> Unchecked sh) -> matrix0 -> matrix1
mapSize matrix1 -> Symmetric (Unchecked sh) a
f matrix0
a =
   Symmetric (Unchecked sh) a -> Symmetric sh a
forall lo diag up sh a.
Triangular lo diag up (Unchecked sh) a
-> Triangular lo diag up sh a
recheck (Symmetric (Unchecked sh) a -> Symmetric sh a)
-> Symmetric (Unchecked sh) a -> Symmetric sh a
forall a b. (a -> b) -> a -> b
$ matrix1 -> Symmetric (Unchecked sh) a
f (matrix1 -> Symmetric (Unchecked sh) a)
-> matrix1 -> Symmetric (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 a -> Matrix.General height width a -> Symmetric width a
congruenceDiagonal :: Vector height a -> General height width a -> Symmetric width a
congruenceDiagonal Vector height a
d =
   ((width -> Unchecked width)
 -> General height width a
 -> Full Big Big height (Unchecked width) a)
-> (Full Big Big height (Unchecked width) a
    -> Symmetric (Unchecked width) a)
-> General height width a
-> Symmetric width a
forall sh matrix0 matrix1 a.
((sh -> Unchecked sh) -> matrix0 -> matrix1)
-> (matrix1 -> Symmetric (Unchecked sh) a)
-> matrix0
-> Symmetric 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
  -> Symmetric (Unchecked width) a)
 -> General height width a -> Symmetric width a)
-> (Full Big Big height (Unchecked width) a
    -> Symmetric (Unchecked width) a)
-> General height width a
-> Symmetric 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
-> Symmetric (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
-> Symmetric width a
scaledAnticommutator a
0.5 Full Big Big height (Unchecked width) a
a (Full Big Big height (Unchecked width) a
 -> Symmetric (Unchecked width) a)
-> Full Big Big height (Unchecked width) a
-> Symmetric (Unchecked width) a
forall a b. (a -> b) -> a -> b
$ Vector height 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 a
-> Full vert horiz height width a -> Full vert horiz height width a
Basic.scaleRows Vector height a
d Full Big Big height (Unchecked width) a
a

congruenceDiagonalTransposed ::
   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
   Matrix.General height width a -> Vector width a -> Symmetric height a
congruenceDiagonalTransposed :: General height width a -> Vector width a -> Symmetric height a
congruenceDiagonalTransposed =
   (Vector width a -> General height width a -> Symmetric height a)
-> General height width a -> Vector width a -> Symmetric height a
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Vector width a -> General height width a -> Symmetric height a)
 -> General height width a -> Vector width a -> Symmetric height a)
-> (Vector width a -> General height width a -> Symmetric height a)
-> General height width a
-> Vector width a
-> Symmetric height a
forall a b. (a -> b) -> a -> b
$ \Vector width a
d -> ((height -> Unchecked height)
 -> General height width a
 -> Full Big Big (Unchecked height) width a)
-> (Full Big Big (Unchecked height) width a
    -> Symmetric (Unchecked height) a)
-> General height width a
-> Symmetric height a
forall sh matrix0 matrix1 a.
((sh -> Unchecked sh) -> matrix0 -> matrix1)
-> (matrix1 -> Symmetric (Unchecked sh) a)
-> matrix0
-> Symmetric 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
  -> Symmetric (Unchecked height) a)
 -> General height width a -> Symmetric height a)
-> (Full Big Big (Unchecked height) width a
    -> Symmetric (Unchecked height) a)
-> General height width a
-> Symmetric 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
-> Symmetric (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
-> Symmetric height a
scaledAnticommutatorTransposed a
0.5 Full Big Big (Unchecked height) width a
a (Full Big Big (Unchecked height) width a
 -> Symmetric (Unchecked height) a)
-> Full Big Big (Unchecked height) width a
-> Symmetric (Unchecked height) a
forall a b. (a -> b) -> a -> b
$ Vector width 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 a
-> Full vert horiz height width a -> Full vert horiz height width a
Basic.scaleColumns Vector width a
d Full Big Big (Unchecked height) width a
a


congruence ::
   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Symmetric height a -> Matrix.General height width a -> Symmetric width a
congruence :: Symmetric height a -> General height width a -> Symmetric width a
congruence Symmetric height a
b =
   ((width -> Unchecked width)
 -> General height width a
 -> Full Big Big height (Unchecked width) a)
-> (Full Big Big height (Unchecked width) a
    -> Symmetric (Unchecked width) a)
-> General height width a
-> Symmetric width a
forall sh matrix0 matrix1 a.
((sh -> Unchecked sh) -> matrix0 -> matrix1)
-> (matrix1 -> Symmetric (Unchecked sh) a)
-> matrix0
-> Symmetric 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
  -> Symmetric (Unchecked width) a)
 -> General height width a -> Symmetric width a)
-> (Full Big Big height (Unchecked width) a
    -> Symmetric (Unchecked width) a)
-> General height width a
-> Symmetric 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
-> Symmetric (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
-> Symmetric 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
 -> Symmetric (Unchecked width) a)
-> Full Big Big height (Unchecked width) a
-> Symmetric (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
         ((Triangular Filled NonUnit Filled height -> Order)
-> Symmetric 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 Triangular Filled NonUnit Filled height -> Order
forall lo diag up size. Triangular lo diag up size -> Order
MatrixShape.triangularOrder Symmetric height a
b) Full Big Big height (Unchecked width) a
a

congruenceTransposed ::
   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
   Matrix.General height width a -> Symmetric width a -> Symmetric height a
congruenceTransposed :: General height width a -> Symmetric width a -> Symmetric height a
congruenceTransposed =
   (Symmetric width a -> General height width a -> Symmetric height a)
-> General height width a
-> Symmetric width a
-> Symmetric height a
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Symmetric width a
  -> General height width a -> Symmetric height a)
 -> General height width a
 -> Symmetric width a
 -> Symmetric height a)
-> (Symmetric width a
    -> General height width a -> Symmetric height a)
-> General height width a
-> Symmetric width a
-> Symmetric height a
forall a b. (a -> b) -> a -> b
$ \Symmetric width a
b -> ((height -> Unchecked height)
 -> General height width a
 -> Full Big Big (Unchecked height) width a)
-> (Full Big Big (Unchecked height) width a
    -> Symmetric (Unchecked height) a)
-> General height width a
-> Symmetric height a
forall sh matrix0 matrix1 a.
((sh -> Unchecked sh) -> matrix0 -> matrix1)
-> (matrix1 -> Symmetric (Unchecked sh) a)
-> matrix0
-> Symmetric 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
  -> Symmetric (Unchecked height) a)
 -> General height width a -> Symmetric height a)
-> (Full Big Big (Unchecked height) width a
    -> Symmetric (Unchecked height) a)
-> General height width a
-> Symmetric 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
-> Symmetric (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
-> Symmetric height a
scaledAnticommutatorTransposed a
forall a. Floating a => a
one Full Big Big (Unchecked height) width a
a (Full Big Big (Unchecked height) width a
 -> Symmetric (Unchecked height) a)
-> Full Big Big (Unchecked height) width a
-> Symmetric (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 ((Triangular Filled NonUnit Filled width -> Order)
-> Symmetric 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 Triangular Filled NonUnit Filled width -> Order
forall lo diag up size. Triangular lo diag up size -> Order
MatrixShape.triangularOrder Symmetric 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 -> Symmetric width a
scaledAnticommutator :: a
-> Full vert horiz height width a
-> Full vert horiz height width a
-> Symmetric 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
   Symmetric width -> (Ptr a -> IO ()) -> Symmetric width a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> width -> Symmetric width
forall size. Order -> size -> Symmetric size
symmetricShape Order
order (width -> Symmetric width) -> width -> Symmetric 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 ()) -> Symmetric width a)
-> (Ptr a -> IO ()) -> Symmetric width a
forall a b. (a -> b) -> a -> b
$
         \Ptr a
cpPtr -> do
      String -> Bool -> IO ()
Call.assert String
"Symmetric.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

scaledAnticommutatorTransposed ::
   (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 -> Symmetric height a
scaledAnticommutatorTransposed :: a
-> Full vert horiz height width a
-> Full vert horiz height width a
-> Symmetric height a
scaledAnticommutatorTransposed
      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
   Symmetric height -> (Ptr a -> IO ()) -> Symmetric height a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> height -> Symmetric height
forall size. Order -> size -> Symmetric size
symmetricShape Order
order (height -> Symmetric height) -> height -> Symmetric 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 ()) -> Symmetric height a)
-> (Ptr a -> IO ()) -> Symmetric height a
forall a b. (a -> b) -> a -> b
$
         \Ptr a
cpPtr -> do
      String -> Bool -> IO ()
Call.assert String
"Symmetric.anticommutatorTransposed: 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))
gramianTransposedParameters 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 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 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
         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.syr2k 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 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


symmetricShape :: Order -> size -> MatrixShape.Symmetric size
symmetricShape :: Order -> size -> Symmetric size
symmetricShape = NonUnit -> (Filled, Filled) -> Order -> size -> Symmetric size
forall lo diag up size.
diag -> (lo, up) -> Order -> size -> Triangular lo diag up size
MatrixShape.Triangular NonUnit
NonUnit (Filled, Filled)
forall lo up. (Content lo, Content up) => (lo, up)
MatrixShape.autoUplo