{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE GADTs #-}
module Numeric.LAPACK.Matrix.Symmetric.Unified where

import qualified Numeric.LAPACK.Matrix.Mosaic.Unpacked as Unpacked
import qualified Numeric.LAPACK.Matrix.Mosaic.Packed as Packed
import qualified Numeric.LAPACK.Matrix.Mosaic.Generic as Mosaic
import qualified Numeric.LAPACK.Matrix.Mosaic.Private as MosaicPriv
import qualified Numeric.LAPACK.Matrix.Basic as Basic
import qualified Numeric.LAPACK.Matrix.Private as Matrix
import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Vector as Vector
import qualified Numeric.LAPACK.Scalar as Scalar
import Numeric.LAPACK.Matrix.Mosaic.Private
         (Mosaic, recheck, copyTriangleToTemp, forPointers,
          diagonalPointers, columnMajorPointers, rowMajorPointers,
          withPacking, withPackingLinear, runPacking,
          noLabel, label, applyFuncPair, triArg,
          Labelled(Labelled))
import Numeric.LAPACK.Matrix.Layout.Private
         (Order(RowMajor,ColumnMajor), uploFromOrder, uploOrder)
import Numeric.LAPACK.Matrix.Modifier
         (Conjugation(NonConjugated, Conjugated), conjugatedOnRowMajor,
          Transposition(NonTransposed, Transposed), transposeOrder)
import Numeric.LAPACK.Matrix.Private (Full, General, Square)
import Numeric.LAPACK.Linear.Private (solver, withDeterminantInfo, diagonalMsg)
import Numeric.LAPACK.Vector (Vector)
import Numeric.LAPACK.Scalar (RealOf, zero, one)
import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))
import Numeric.LAPACK.Private
         (fill, copyBlock, copyToTemp,
          condConjugate, copyCondConjugate, condConjugateToTemp,
          withAutoWorkspace, pointerSeq)

import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
import qualified Numeric.LAPACK.FFI.Generic as LapackGen
import qualified Numeric.BLAS.FFI.Complex as BlasComplex
import qualified Numeric.BLAS.FFI.Real as BlasReal
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.Marshal.Array (advancePtr)
import Foreign.C.Types (CInt, CChar)
import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
import Foreign.Ptr (Ptr)
import Foreign.Storable (Storable, peek)

import qualified System.IO.Lazy as LazyIO
import System.IO.Unsafe (unsafePerformIO)

import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
import Control.Monad.IO.Class (liftIO)
import Control.Applicative ((<$>))



data MirrorSingleton mirror where
   SimpleMirror :: MirrorSingleton Layout.SimpleMirror
   ConjugateMirror :: MirrorSingleton Layout.ConjugateMirror

deriving instance Eq (MirrorSingleton mirror)
deriving instance Show (MirrorSingleton mirror)

class (Layout.Mirror mirror) => Mirror mirror where
   autoMirror :: MirrorSingleton mirror
instance Mirror Layout.SimpleMirror where autoMirror = SimpleMirror
instance Mirror Layout.ConjugateMirror where autoMirror = ConjugateMirror

narrowMirror :: (Mirror mirror) => f mirror -> MirrorSingleton mirror
narrowMirror _ = autoMirror

conjugationFromMirror :: (Mirror mirror) => f mirror -> Conjugation
conjugationFromMirror mirror =
   case narrowMirror mirror of
      SimpleMirror -> NonConjugated
      ConjugateMirror -> Conjugated



{- |
> let b = takeHalf a
> ==>
> isTriangular b && a == addTransposed b
-}
takeHalf ::
   (Shape.C sh, Class.Floating a) =>
   Unpacked.Mosaic mirror Shape.Upper sh a ->
   Unpacked.Triangular Shape.Upper sh a
takeHalf (Array shape a) =
   Array.unsafeCreateWithSize
         shape{Layout.mosaicMirror = Layout.NoMirror} $
      \size bPtr -> evalContT $ do
   let n = Shape.size $ Layout.mosaicSize shape
   aPtr <- ContT $ withForeignPtr a
   nPtr <- Call.cint n
   alphaPtr <- Call.number 0.5
   incxPtr <- Call.cint (n+1)
   liftIO $ do
      copyBlock size aPtr bPtr
      BlasGen.scal nPtr alphaPtr bPtr incxPtr


toSquare ::
   (Mirror mirror, Shape.C sh, Class.Floating a) =>
   Packed.Mosaic mirror Shape.Upper sh a -> Square sh a
toSquare (Array (Layout.Mosaic _pack mirror _uplo order sh) a) =
   Array.unsafeCreate (Layout.square order sh) $ \bPtr ->
   withForeignPtr a $ \aPtr -> evalContT $ do

   let conj = conjugationFromMirror mirror
   let n = Shape.size sh
   incxPtr <- Call.cint 1
   incyPtr <- Call.cint n
   liftIO $ case order of
      RowMajor ->
         forPointers (rowMajorPointers n bPtr aPtr) $
               \nPtr (dstPtr,srcPtr) -> do
            copyCondConjugate conj nPtr srcPtr incxPtr dstPtr incyPtr
            BlasGen.copy nPtr srcPtr incxPtr dstPtr incxPtr
      ColumnMajor ->
         forPointers (columnMajorPointers n bPtr aPtr) $
               \nPtr ((dstRowPtr,dstColumnPtr),srcPtr) -> do
            copyCondConjugate conj nPtr srcPtr incxPtr dstRowPtr incyPtr
            BlasGen.copy nPtr srcPtr incxPtr dstColumnPtr incxPtr


addMirrored ::
   (Mirror mirror, Shape.C sh, Class.Floating a) =>
   Square sh a -> Packed.Mosaic mirror Shape.Upper sh a
addMirrored (Array (Layout.Full order extent) a) =
   let sh = Extent.squareSize extent
       mirror = Layout.autoMirror
       shape =
         Layout.Mosaic Layout.Packed mirror Layout.Upper order sh
   in Array.unsafeCreate shape $ \bPtr -> evalContT $ do

   let conj = conjugationFromMirror mirror
   let n = Shape.size sh
   alphaPtr <- Call.number one
   incxPtr <- Call.cint 1
   incnPtr <- Call.cint n
   aPtr <- ContT $ withForeignPtr a
   liftIO $ case order of
      RowMajor ->
         forPointers (rowMajorPointers n aPtr bPtr) $
            \nPtr (srcPtr,dstPtr) -> do
         copyCondConjugate conj nPtr srcPtr incnPtr dstPtr incxPtr
         BlasGen.axpy nPtr alphaPtr srcPtr incxPtr dstPtr incxPtr
      ColumnMajor ->
         forPointers (columnMajorPointers n aPtr bPtr) $
            \nPtr ((srcRowPtr,srcColumnPtr),dstPtr) -> do
         copyCondConjugate conj nPtr srcRowPtr incnPtr dstPtr incxPtr
         BlasGen.axpy nPtr alphaPtr srcColumnPtr incxPtr dstPtr incxPtr


complementUnpacked ::
   (Class.Floating a) =>
   Conjugation -> Order -> Int -> Ptr a -> Ptr a -> IO ()
complementUnpacked conj order n aPtr bPtr = evalContT $ do
   inc1Ptr <- Call.cint 1
   incnPtr <- Call.cint n
   let number = take n . zip [0..]
   liftIO $ case order of
      RowMajor ->
         forPointers (number $ zip (pointerSeq 1 aPtr) (pointerSeq n bPtr)) $
               \nPtr (srcPtr,dstPtr) ->
            copyCondConjugate conj nPtr srcPtr incnPtr dstPtr inc1Ptr
      ColumnMajor ->
         forPointers (number $ zip (pointerSeq n aPtr) (pointerSeq 1 bPtr)) $
               \nPtr (srcPtr,dstPtr) ->
            copyCondConjugate conj nPtr srcPtr inc1Ptr dstPtr incnPtr

complement ::
   (Class.Floating a) =>
   Layout.PackingSingleton pack ->
   Conjugation -> Order -> Int -> Ptr a -> IO ()
complement pack conj order n bPtr = do
   case pack of
      Layout.Packed -> return ()
      Layout.Unpacked -> complementUnpacked conj order n bPtr bPtr


type SPMV ar a =
      Ptr CChar -> Ptr CInt -> Ptr ar -> Ptr a ->
      Ptr a -> Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> IO ()

type SYMV ar a =
      Ptr CChar -> Ptr CInt -> Ptr ar -> Ptr a -> Ptr CInt ->
      Ptr a -> Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> IO ()

multiplyVectorCont ::
   (Class.Floating ar, Class.Floating a) =>
   SPMV ar a ->
   SYMV ar a ->
   Layout.PackingSingleton pack ->
   Layout.UpLoSingleton uplo ->
   Layout.Order ->
   Int ->
   ForeignPtr a -> Ptr a -> Ptr a -> ContT () IO ()
multiplyVectorCont spmv symv pack uplo order n a xPtr yPtr = do
   uploPtr <- Call.char $ uploFromOrder $ uploOrder uplo order
   nPtr <- Call.cint n
   alphaPtr <- Call.number one
   aPtr <- ContT $ withForeignPtr a
   incxPtr <- Call.cint 1
   betaPtr <- Call.number zero
   incyPtr <- Call.cint 1
   withPacking pack $
      applyFuncPair (noLabel spmv) (noLabel symv)
         uploPtr nPtr alphaPtr (triArg aPtr n)
         xPtr incxPtr betaPtr yPtr incyPtr

-- Triangular is a bit different, it has no alpha
multiplyVector ::
   (Mirror mirror, Shape.C sh, Eq sh, Class.Floating a) =>
   Mosaic pack mirror uplo sh a -> Vector sh a -> Vector sh a
multiplyVector
   (Array (Layout.Mosaic pack mirror uplo order shA) a) (Array shX x) =
      Array.unsafeCreateWithSize shX $ \n yPtr -> do
   Call.assert "Symmetric.multiplyVector: width shapes mismatch" (shA == shX)
   evalContT $
      case narrowMirror mirror of
         ConjugateMirror -> do
            let conj = conjugatedOnRowMajor order
            xPtr <- condConjugateToTemp conj n x
            multiplyVectorCont BlasGen.hpmv BlasGen.hemv
               pack uplo order n a xPtr yPtr
            nPtr <- Call.cint n
            incyPtr <- Call.cint 1
            liftIO $ condConjugate conj nPtr yPtr incyPtr
         SimpleMirror -> do
            xPtr <- ContT $ withForeignPtr x
            case Scalar.complexSingletonOfFunctor x of
               Scalar.Real ->
                  multiplyVectorCont BlasReal.spmv BlasReal.symv
                     pack uplo order n a xPtr yPtr
               Scalar.Complex ->
                  multiplyVectorCont LapackComplex.spmv LapackComplex.symv
                     pack uplo order n a xPtr yPtr


multiplyFull ::
   (Layout.UpLo uplo, Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Mosaic pack mirror uplo height a ->
   Full meas vert horiz height width a ->
   Full meas vert horiz height width a
multiplyFull =
   Unpacked.multiplyFull Omni.Arbitrary . Mosaic.unpackDirty


forceUpper ::
   (Layout.Packing pack, Mirror mirror, Shape.C sh, Class.Floating a) =>
   Mosaic pack mirror uplo sh a ->
   Mosaic pack mirror Shape.Upper sh a
forceUpper a =
   case Layout.mosaicUplo $ Array.shape a of
      Layout.Upper -> a
      Layout.Lower -> upperFromLower a

upperFromLower ::
   (Layout.Packing pack, Mirror mirror, Shape.C sh, Class.Floating a) =>
   Mosaic pack mirror Shape.Lower sh a ->
   Mosaic pack mirror Shape.Upper sh a
upperFromLower
      (Array (Layout.Mosaic pack mirror Layout.Lower order sh) a) =

   (case narrowMirror mirror of
      SimpleMirror -> id
      ConjugateMirror -> Vector.conjugate) $
   Array
      (Layout.Mosaic
         pack mirror Layout.Upper (Layout.flipOrder order) sh)
      a


type SYR ar a f =
      Ptr CChar -> Ptr CInt -> Ptr ar -> Ptr a -> Ptr CInt -> Ptr a -> f

withSyrSingleton ::
   (Class.Floating a) => (Scalar.ComplexSingleton a -> SYR a a f) -> SYR a a f
withSyrSingleton f = f Scalar.complexSingleton

spr :: (Class.Floating a) => SYR a a (IO ())
spr = withSyrSingleton $ \sw ->
   case sw of
      Scalar.Real -> BlasReal.spr
      Scalar.Complex -> LapackComplex.spr

syr :: (Class.Floating a) => SYR a a (Ptr CInt -> IO ())
syr = withSyrSingleton $ \sw ->
   case sw of
      Scalar.Real -> BlasReal.syr
      Scalar.Complex -> LapackComplex.syr

outerCont ::
   (Class.Floating a, Class.Floating ar) =>
   SYR ar a (IO ()) ->
   SYR ar a (Ptr CInt -> IO ()) ->
   ForeignPtr a -> Int ->
   Layout.PackingSingleton pack ->
   Layout.UpLoSingleton uplo ->
   Ptr a -> ContT () IO ()
outerCont spr_ syr_ x n pack uplo aPtr = do
   uploPtr <- Call.char $ Layout.uploChar uplo
   nPtr <- Call.cint n
   alphaPtr <- Call.number one
   xPtr <- ContT $ withForeignPtr x
   incxPtr <- Call.cint 1
   withPacking pack $
      applyFuncPair (noLabel spr_) (noLabel syr_)
         uploPtr nPtr alphaPtr xPtr incxPtr (triArg aPtr n)

outer ::
   (Layout.Packing pack, Mirror mirror, Layout.UpLo uplo,
    Shape.C sh, Class.Floating a) =>
   Vector sh a -> Mosaic pack mirror uplo sh a
outer (Array sh x) =
   let n = Shape.size sh
       pack = Layout.autoPacking
       mirror = Layout.autoMirror
       uplo = Layout.autoUplo
       order = Layout.ColumnMajor
   in Array.unsafeCreateWithSize
         (Layout.Mosaic pack mirror uplo order sh) $
      \triSize aPtr -> do
   fill zero triSize aPtr
   evalContT $
      case (Scalar.complexSingletonOfFunctor x, narrowMirror mirror) of
         (Scalar.Complex, ConjugateMirror) ->
            outerCont BlasComplex.hpr BlasComplex.her x n pack uplo aPtr
         _ -> outerCont spr syr x n pack uplo aPtr
   complement pack (conjugationFromMirror mirror) (uploOrder uplo order) n aPtr

outerUpper ::
   (Layout.Packing pack, Mirror mirror, Shape.C sh, Class.Floating a) =>
   Order -> Vector sh a -> Mosaic pack mirror Shape.Upper sh a
outerUpper order x =
   case order of
      Layout.ColumnMajor -> outer x
      Layout.RowMajor -> upperFromLower $ outer x


upperShape ::
   (Layout.Mirror mirror) =>
   Layout.PackingSingleton pack ->
   Layout.MirrorSingleton mirror -> Order -> size ->
   Layout.Mosaic pack mirror Shape.Upper size
upperShape pack mirror = Layout.Mosaic pack mirror Layout.Upper

unpackedShape ::
   (Layout.Mirror mirror) =>
   Layout.MirrorSingleton mirror -> Order -> size ->
   Layout.Mosaic Layout.Unpacked mirror Shape.Upper size
unpackedShape = upperShape Layout.Unpacked

skipCheckCongruence ::
   ((sh -> Unchecked sh) -> matrix0 -> matrix1) ->
   (matrix1 -> Mosaic pack mirror uplo (Unchecked sh) a) ->
   matrix0 -> Mosaic pack mirror uplo sh a
skipCheckCongruence mapSize f  =  recheck . f . mapSize Unchecked


postHook :: (Monad m) => m () -> ContT r m ()
postHook hook = ContT $ \act -> do r <- act (); hook; return r

temporaryUnpacked ::
   (Mirror mirror, Class.Floating a) =>
   Layout.PackingSingleton pack -> MirrorSingleton mirror -> Order ->
   Int -> Ptr a -> ContT () IO (Ptr a)
temporaryUnpacked pack mirror order n cpPtr =
   case pack of
      Layout.Packed -> do
         cPtr <- Call.allocaArray (n*n)
         postHook $ MosaicPriv.pack order n cPtr cpPtr
         return cPtr
      Layout.Unpacked -> do
         postHook $
            complementUnpacked
               (conjugationFromMirror mirror) order n cpPtr cpPtr
         return cpPtr


gramianParameters ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Shape.C width) =>
   MirrorSingleton mirror ->
   Transposition ->
   Order ->
   Extent.Extent meas vert horiz height width ->
   ((Int, Int), (Char, Char, Int))
gramianParameters mirror trans order extent =
   let (height, width) = Extent.dimensions extent
       n = Shape.size width
       k = Shape.size height
       mirrorChar =
         case mirror of
            SimpleMirror -> 'T'
            ConjugateMirror -> 'C'
       transChar =
         case transposeOrder trans order of
            ColumnMajor -> mirrorChar
            RowMajor -> 'N'
       lda =
         case order of
            ColumnMajor -> k
            RowMajor -> n
    in (case trans of NonTransposed -> (n,k); Transposed -> (k,n),
        (uploFromOrder order, transChar, lda))


{-
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 ::
   (Mirror mirror, Class.Floating a) =>
   Layout.PackingSingleton pack -> MirrorSingleton mirror -> Order ->
   ForeignPtr a -> Ptr a ->
   ((Int, Int), (Char, Char, Int)) -> IO ()
gramianIO pack mirror order a cpPtr ((n,k), (uplo,trans,lda)) = evalContT $ do
   uploPtr <- Call.char uplo
   transPtr <- Call.char trans
   nPtr <- Call.cint n
   kPtr <- Call.cint k
   aPtr <- ContT $ withForeignPtr a
   ldaPtr <- Call.leadingDim lda
   cPtr <- temporaryUnpacked pack mirror order n cpPtr
   ldcPtr <- Call.leadingDim n
   case (k, mirror, Scalar.complexSingletonOfFunctor a) of
      (0, _, _) -> liftIO $ fill zero (n*n) cPtr
      (_, ConjugateMirror, Scalar.Complex) -> do
         alphaPtr <- Call.number one
         betaPtr <- Call.number zero
         liftIO $
            BlasComplex.herk uploPtr transPtr
               nPtr kPtr alphaPtr aPtr ldaPtr betaPtr cPtr ldcPtr
      _ -> do
         alphaPtr <- Call.number one
         betaPtr <- Call.number zero
         liftIO $
            BlasGen.syrk uploPtr transPtr
               nPtr kPtr alphaPtr aPtr ldaPtr betaPtr cPtr ldcPtr

gramian ::
   (Layout.Packing pack, Mirror mirror,
    Shape.C height, Shape.C width, Class.Floating a) =>
   Matrix.General height width a ->
   Mosaic pack mirror Shape.Upper width a
gramian (Array (Layout.Full order extent) a) =
   let pack = Layout.autoPacking
       mirror = Layout.autoMirror
   in Array.unsafeCreate (upperShape pack mirror order $ Extent.width extent) $
      \bPtr ->
         gramianIO pack (narrowMirror mirror) order a bPtr $
         gramianParameters (narrowMirror mirror) NonTransposed order extent

gramianTransposed ::
   (Layout.Packing pack, Mirror mirror,
    Shape.C height, Shape.C width, Class.Floating a) =>
   Matrix.General height width a ->
   Mosaic pack mirror Shape.Upper height a
gramianTransposed (Array (Layout.Full order extent) a) =
   let pack = Layout.autoPacking
       mirror = Layout.autoMirror
   in Array.unsafeCreate (upperShape pack mirror order $ Extent.height extent) $
      \bPtr ->
         gramianIO pack (narrowMirror mirror) order a bPtr $
         gramianParameters (narrowMirror mirror) Transposed order extent


scaledAnticommutatorIO ::
   (Mirror mirror, Class.Floating a) =>
   Layout.PackingSingleton pack -> MirrorSingleton mirror -> Order -> a ->
   ForeignPtr a -> ForeignPtr a -> Ptr a ->
   ((Int, Int), (Char, Char, Int)) -> IO ()
scaledAnticommutatorIO
   pack mirror order alpha a b cpPtr ((n,k), (uplo,trans,lda)) =
      evalContT $ do
   uploPtr <- Call.char uplo
   transPtr <- Call.char trans
   nPtr <- Call.cint n
   kPtr <- Call.cint k
   alphaPtr <- Call.number alpha
   aPtr <- ContT $ withForeignPtr a
   ldaPtr <- Call.leadingDim lda
   bPtr <- ContT $ withForeignPtr b
   let ldbPtr = ldaPtr
   cPtr <- temporaryUnpacked pack mirror order n cpPtr
   ldcPtr <- Call.leadingDim n
   case (mirror, Scalar.complexSingletonOfFunctor a) of
      (ConjugateMirror, Scalar.Complex) -> do
         betaPtr <- Call.number zero
         liftIO $
            BlasComplex.her2k uploPtr transPtr nPtr kPtr alphaPtr
               aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr
      _ -> do
         betaPtr <- Call.number zero
         liftIO $
            BlasGen.syr2k uploPtr transPtr nPtr kPtr alphaPtr
               aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr

scaledAnticommutator ::
   (Layout.Packing pack, Mirror mirror,
    Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
   Layout.MirrorSingleton mirror ->
   a ->
   Full meas vert horiz height width a ->
   Full meas vert horiz height width a ->
   Mosaic pack mirror Shape.Upper width a
scaledAnticommutator mirror
      alpha arr (Array (Layout.Full order extentB) b) = do
   let (Array (Layout.Full _ extentA) a) = Basic.forceOrder order arr
       pack = Layout.autoPacking
   Array.unsafeCreate (upperShape pack mirror order $ Extent.width extentB) $
         \cPtr -> do
      Call.assert "Symmetric.anticommutator: extents mismatch"
         (extentA==extentB)
      scaledAnticommutatorIO pack (narrowMirror mirror) order alpha a b cPtr $
         gramianParameters (narrowMirror mirror) NonTransposed order extentB

scaledAnticommutatorTransposed ::
   (Layout.Packing pack, Mirror mirror,
    Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
   Layout.MirrorSingleton mirror ->
   a ->
   Full meas vert horiz height width a ->
   Full meas vert horiz height width a ->
   Mosaic pack mirror Shape.Upper height a
scaledAnticommutatorTransposed mirror
      alpha arr (Array (Layout.Full order extentB) b) = do
   let (Array (Layout.Full _ extentA) a) = Basic.forceOrder order arr
       pack = Layout.autoPacking
   Array.unsafeCreate (upperShape pack mirror order $ Extent.height extentB) $
         \cPtr -> do
      Call.assert "Symmetric.anticommutatorTransposed: extents mismatch"
         (extentA==extentB)
      scaledAnticommutatorIO pack (narrowMirror mirror) order alpha a b cPtr $
         gramianParameters (narrowMirror mirror) Transposed order extentB


congruenceRealDiagonal ::
   (Layout.Packing pack, Mirror mirror,
    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Vector height (RealOf a) ->
   General height width a ->
   Mosaic pack mirror Shape.Upper width a
congruenceRealDiagonal d =
   skipCheckCongruence Basic.mapWidth $ \a ->
      scaledAnticommutator Layout.autoMirror 0.5 a $
         Basic.scaleRowsReal d a

congruenceRealDiagonalTransposed ::
   (Layout.Packing pack, Mirror mirror,
    Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
   General height width a ->
   Vector width (RealOf a) ->
   Mosaic pack mirror Shape.Upper height a
congruenceRealDiagonalTransposed =
   flip $ \d -> skipCheckCongruence Basic.mapHeight $ \a ->
      scaledAnticommutatorTransposed Layout.autoMirror 0.5 a $
         Basic.scaleColumnsReal d a

congruence ::
   (Layout.Packing pack, Mirror mirror,
    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Unpacked.Mosaic mirror Shape.Upper height a ->
   Matrix.General height width a ->
   Mosaic pack mirror Shape.Upper width a
congruence b =
   skipCheckCongruence Basic.mapWidth $ \a ->
      scaledAnticommutator
         (Layout.mosaicMirror $ Array.shape b) one a $
      Unpacked.multiplyFull Omni.Arbitrary (takeHalf b) a

congruenceTransposed ::
   (Layout.Packing pack, Mirror mirror,
    Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
   Matrix.General height width a ->
   Unpacked.Mosaic mirror Shape.Upper width a ->
   Mosaic pack mirror Shape.Upper height a
congruenceTransposed =
   flip $ \b -> skipCheckCongruence Basic.mapHeight $ \a ->
      scaledAnticommutatorTransposed
         (Layout.mosaicMirror $ Array.shape b) one a $
      Basic.swapMultiply
         (Unpacked.multiplyFull Omni.Arbitrary . Mosaic.transpose)
         a (takeHalf b)


solve ::
   (Mirror mirror, Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C width, Shape.C height, Eq height, Class.Floating a) =>
   Mosaic pack mirror Shape.Upper height a ->
   Full meas vert horiz height width a ->
   Full meas vert horiz height width a
solve (Array shape@(Layout.Mosaic pack mirror _upper order sh) a) =
   solver (show mirror) sh $ \n nPtr nrhsPtr xPtr ldxPtr -> do
      uploPtr <- Call.char $ uploFromOrder order
      let conj = conjugationFromMirror mirror
      aPtr <- copyTriangleToTemp conj order (Shape.size shape) a
      ipivPtr <- Call.allocaArray n
      withPackingLinear diagonalMsg pack $ fmap autoWorkspace $
         uncurry applyFuncPair
            (case conj of
               Conjugated ->
                  (label "hpsv" LapackGen.hpsv, label "hesv" LapackGen.hesv)
               NonConjugated ->
                  (label "spsv" LapackGen.spsv, label "sysv" LapackGen.sysv))
            uploPtr nPtr nrhsPtr (triArg aPtr n) ipivPtr xPtr ldxPtr


inverse ::
   (Mirror mirror, Layout.UpLo uplo, Shape.C sh, Class.Floating a) =>
   Mosaic pack mirror uplo sh a ->
   Mosaic pack mirror uplo sh a
inverse (Array shape@(Layout.Mosaic pack mirror uplo order sh) a) =
   Array.unsafeCreateWithSize shape $ \triSize bPtr -> evalContT $ do

   let n = Shape.size sh
   let realOrder = uploOrder uplo order
   uploPtr <- Call.char $ uploFromOrder realOrder
   nPtr <- Call.cint n
   aPtr <- ContT $ withForeignPtr a
   ipivPtr <- Call.allocaArray n
   liftIO $ copyBlock triSize aPtr bPtr
   let conj = conjugationFromMirror mirror
   let (trf,tri) =
         case conj of
            Conjugated ->
               ((label "hptrf" LapackGen.hptrf, label "hetrf" LapackGen.hetrf),
                (label "hptri" LapackGen.hptri, label "hetri" LapackGen.hetri))
            NonConjugated ->
               ((label "sptrf" LapackGen.sptrf, label "sytrf" LapackGen.sytrf),
                (label "sptri" LapackGen.sptri, label "sytri" LapackGen.sytri))
   withPackingLinear diagonalMsg pack $ fmap autoWorkspace $
      uncurry applyFuncPair trf uploPtr nPtr (triArg bPtr n) ipivPtr
   workPtr <- Call.allocaArray n
   withPackingLinear diagonalMsg pack $
      uncurry applyFuncPair tri uploPtr nPtr (triArg bPtr n) ipivPtr workPtr
   liftIO $ complement pack conj realOrder n bPtr


blockDiagonalPointers ::
   (Storable a) =>
   Order -> [(Ptr CInt, Ptr a)] -> LazyIO.T [(Ptr a, Maybe (Ptr a, Ptr a))]
blockDiagonalPointers order =
   let go ((ipiv0Ptr,a0Ptr):ptrs0) = do
         ipiv <- LazyIO.interleave $ peek ipiv0Ptr
         (ext,ptrTuples) <-
            if ipiv >= 0
               then (,) Nothing <$> go ptrs0
               else
                  case ptrs0 of
                     [] -> error "Symmetric.determinant: incomplete 2x2 block"
                     (_ipiv1Ptr,a1Ptr):ptrs1 ->
                        let bPtr =
                              case order of
                                 ColumnMajor -> advancePtr a1Ptr (-1)
                                 RowMajor -> advancePtr a0Ptr 1
                        in (,) (Just (a1Ptr,bPtr)) <$> go ptrs1
         return $ (a0Ptr,ext) : ptrTuples
       go [] = return []
   in go

determinant ::
   (Mirror mirror, Layout.UpLo uplo,
    Shape.C sh, Class.Floating a, Class.Floating ar) =>
   ((Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar) ->
   Mosaic pack mirror uplo sh a -> ar
determinant peekBlockDeterminant
   (Array shape@(Layout.Mosaic pack mirror uplo order sh) a) =
      unsafePerformIO $ evalContT $ do

   let conj = conjugationFromMirror mirror
   let n = Shape.size sh
   uploPtr <- Call.char $ uploFromOrder $ uploOrder uplo order
   nPtr <- Call.cint n
   aPtr <- copyToTemp (Shape.size shape) a
   ipivPtr <- Call.allocaArray n
   let diagPtrs =
         case pack of
            Layout.Packed -> diagonalPointers order n aPtr
            Layout.Unpacked -> take n $ pointerSeq (n+1) aPtr
   let (utrf,ptrf) =
         case conj of
            Conjugated ->
               (label "hptrf" LapackGen.hptrf, label "hetrf" LapackGen.hetrf)
            NonConjugated ->
               (label "sptrf" LapackGen.sptrf, label "sytrf" LapackGen.sytrf)
   let Labelled name makeTrf =
         runPacking pack $ fmap autoWorkspace $
         applyFuncPair utrf ptrf uploPtr nPtr (triArg aPtr n) ipivPtr
   trf <- makeTrf
   liftIO $
      withDeterminantInfo name trf
         ((return $!) =<<
          LazyIO.run
            (fmap product $
             mapM (LazyIO.interleave . peekBlockDeterminant) =<<
             blockDiagonalPointers order (zip (pointerSeq 1 ipivPtr) diagPtrs)))

autoWorkspace ::
   (Class.Floating a) => (Ptr a -> Ptr CInt -> info -> IO ()) -> info -> IO ()
autoWorkspace trf info =
   withAutoWorkspace $ \work lwork -> trf work lwork info