{-# 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
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
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))
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,
Eq height, Shape.C height, Shape.C width, 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