{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE GADTs #-}
module Numeric.LAPACK.Matrix.BandedHermitian.Basic (
   BandedHermitian,
   StaticVector,
   Transposition(..),
   fromList,
   identity,
   identityFatOrder,
   diagonal,
   takeDiagonal,
   toHermitian,
   toBanded,
   forceOrder,
   takeTopLeft,
   takeBottomRight,
   multiplyVector,
   multiplyFull,
   gramian,
   sumRank1,
   takeUpper,
   ) where

import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Extent.Private as ExtentPriv
import qualified Numeric.LAPACK.Matrix.Extent as Extent
import qualified Numeric.LAPACK.Matrix.Banded.Basic as Banded
import qualified Numeric.LAPACK.Matrix.Mosaic.Private as Mos
import qualified Numeric.LAPACK.Matrix.RowMajor as RowMajor
import qualified Numeric.LAPACK.Matrix.Private as Matrix
import qualified Numeric.LAPACK.Vector.Private as VectorPriv
import qualified Numeric.LAPACK.Vector as Vector
import qualified Data.Array.Comfort.Shape.Static as ShapeStatic
import Numeric.LAPACK.Matrix.Hermitian.Private (TakeDiagonal(..))
import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian)
import Numeric.LAPACK.Matrix.Layout.Private
         (Order(RowMajor,ColumnMajor), uploFromOrder,
          UnaryProxy, natFromProxy)
import Numeric.BLAS.Matrix.Modifier
         (Transposition(NonTransposed, Transposed), transposeOrder,
          Conjugation(NonConjugated), conjugatedOnRowMajor)
import Numeric.LAPACK.Vector (Vector)
import Numeric.LAPACK.Scalar (RealOf, zero, one)
import Numeric.LAPACK.Private
         (fill, lacgv, caseRealComplexFunc, realPtr,
          copyBlock, copyConjugate, condConjugate, condConjugateToTemp,
          pointerSeq, pokeCInt, copySubMatrix)

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 Type.Data.Num.Unary.Literal as TypeNum
import qualified Type.Data.Num.Unary.Proof as Proof
import qualified Type.Data.Num.Unary as Unary
import Type.Data.Num.Unary ((:+:))
import Type.Data.Num (integralFromProxy)
import Type.Base.Proxy (Proxy(Proxy))

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.Marshal.Array (advancePtr)
import Foreign.C.Types (CInt, CChar)
import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
import Foreign.Ptr (Ptr)
import Foreign.Storable (Storable, poke, peek, peekElemOff)

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

import Data.Foldable (for_)
import Data.Tuple.HT (mapPair)

import Data.Complex (Complex, conjugate)


type BandedHermitian offDiag size =
      Array (Layout.BandedHermitian offDiag size)

type Diagonal size = BandedHermitian TypeNum.U0 size


fromList ::
   (Unary.Natural offDiag, Shape.C size, Storable a) =>
   UnaryProxy offDiag -> Order -> size -> [a] ->
   BandedHermitian offDiag size a
fromList :: forall offDiag size a.
(Natural offDiag, C size, Storable a) =>
UnaryProxy offDiag
-> Order -> size -> [a] -> BandedHermitian offDiag size a
fromList UnaryProxy offDiag
numOff Order
order size
size =
   BandedHermitian offDiag size
-> [a] -> Array (BandedHermitian offDiag size) a
forall sh a. (C sh, Storable a) => sh -> [a] -> Array sh a
CheckedArray.fromList (UnaryProxy offDiag -> Order -> size -> BandedHermitian offDiag size
forall off size.
UnaryProxy off -> Order -> size -> BandedHermitian off size
Layout.BandedHermitian UnaryProxy offDiag
numOff Order
order size
size)

identityFatOrder ::
   (Unary.Natural offDiag, Shape.C sh, Class.Floating a) =>
   Order -> sh -> BandedHermitian offDiag sh a
identityFatOrder :: forall offDiag sh a.
(Natural offDiag, C sh, Floating a) =>
Order -> sh -> BandedHermitian offDiag sh a
identityFatOrder Order
order sh
sh =
   case Order
order of
      Order
RowMajor ->
         Matrix sh (() ::+ ZeroBased offDiag) a
-> BandedHermitian offDiag sh a
forall offDiag size a.
(Natural offDiag, C size) =>
Matrix size (() ::+ ZeroBased offDiag) a
-> BandedHermitian offDiag size a
fromRowMajor (Matrix sh (() ::+ ZeroBased offDiag) a
 -> BandedHermitian offDiag sh a)
-> Matrix sh (() ::+ ZeroBased offDiag) a
-> BandedHermitian offDiag sh a
forall a b. (a -> b) -> a -> b
$
         Either Conjugation Conjugation
-> Vector sh a
-> Vector (() ::+ ZeroBased offDiag) a
-> Matrix sh (() ::+ ZeroBased offDiag) a
forall height width a.
(C height, C width, Floating a) =>
Either Conjugation Conjugation
-> Vector height a -> Vector width a -> Matrix height width a
RowMajor.tensorProduct (Conjugation -> Either Conjugation Conjugation
forall a b. a -> Either a b
Left Conjugation
NonConjugated)
            (sh -> Vector sh a
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.one sh
sh) ((() ::+ ZeroBased offDiag)
-> Index (() ::+ ZeroBased offDiag)
-> Vector (() ::+ ZeroBased offDiag) a
forall sh a.
(Indexed sh, Floating a) =>
sh -> Index sh -> Vector sh a
Vector.unit () ::+ ZeroBased offDiag
forall sh. Static sh => sh
Shape.static (Index (() ::+ ZeroBased offDiag)
 -> Vector (() ::+ ZeroBased offDiag) a)
-> Index (() ::+ ZeroBased offDiag)
-> Vector (() ::+ ZeroBased offDiag) a
forall a b. (a -> b) -> a -> b
$ () -> Either () (Index offDiag)
forall a b. a -> Either a b
Left ())
      Order
ColumnMajor ->
         Matrix sh (ZeroBased offDiag ::+ ()) a
-> BandedHermitian offDiag sh a
forall offDiag size a.
(Natural offDiag, C size) =>
Matrix size (ZeroBased offDiag ::+ ()) a
-> BandedHermitian offDiag size a
fromColumnMajor (Matrix sh (ZeroBased offDiag ::+ ()) a
 -> BandedHermitian offDiag sh a)
-> Matrix sh (ZeroBased offDiag ::+ ()) a
-> BandedHermitian offDiag sh a
forall a b. (a -> b) -> a -> b
$
         Either Conjugation Conjugation
-> Vector sh a
-> Vector (ZeroBased offDiag ::+ ()) a
-> Matrix sh (ZeroBased offDiag ::+ ()) a
forall height width a.
(C height, C width, Floating a) =>
Either Conjugation Conjugation
-> Vector height a -> Vector width a -> Matrix height width a
RowMajor.tensorProduct (Conjugation -> Either Conjugation Conjugation
forall a b. a -> Either a b
Left Conjugation
NonConjugated)
            (sh -> Vector sh a
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.one sh
sh) ((ZeroBased offDiag ::+ ())
-> Index (ZeroBased offDiag ::+ ())
-> Vector (ZeroBased offDiag ::+ ()) a
forall sh a.
(Indexed sh, Floating a) =>
sh -> Index sh -> Vector sh a
Vector.unit ZeroBased offDiag ::+ ()
forall sh. Static sh => sh
Shape.static (Index (ZeroBased offDiag ::+ ())
 -> Vector (ZeroBased offDiag ::+ ()) a)
-> Index (ZeroBased offDiag ::+ ())
-> Vector (ZeroBased offDiag ::+ ()) a
forall a b. (a -> b) -> a -> b
$ () -> Either (Index offDiag) ()
forall a b. b -> Either a b
Right ())

fromRowMajor ::
   (Unary.Natural offDiag, Shape.C size) =>
   RowMajor.Matrix size (() ::+ ShapeStatic.ZeroBased offDiag) a ->
   BandedHermitian offDiag size a
fromRowMajor :: forall offDiag size a.
(Natural offDiag, C size) =>
Matrix size (() ::+ ZeroBased offDiag) a
-> BandedHermitian offDiag size a
fromRowMajor =
   ((size, () ::+ ZeroBased offDiag) -> BandedHermitian offDiag size)
-> Array (size, () ::+ ZeroBased offDiag) a
-> Array (BandedHermitian offDiag size) a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape
      (\(size
size, () ::+ ShapeStatic.ZeroBased Proxy (Un offDiag)
k) ->
         Proxy (Un offDiag) -> Order -> size -> BandedHermitian offDiag size
forall off size.
UnaryProxy off -> Order -> size -> BandedHermitian off size
Layout.BandedHermitian Proxy (Un offDiag)
k Order
RowMajor size
size)

fromColumnMajor ::
   (Unary.Natural offDiag, Shape.C size) =>
   RowMajor.Matrix size (ShapeStatic.ZeroBased offDiag ::+ ()) a ->
   BandedHermitian offDiag size a
fromColumnMajor :: forall offDiag size a.
(Natural offDiag, C size) =>
Matrix size (ZeroBased offDiag ::+ ()) a
-> BandedHermitian offDiag size a
fromColumnMajor =
   ((size, ZeroBased offDiag ::+ ()) -> BandedHermitian offDiag size)
-> Array (size, ZeroBased offDiag ::+ ()) a
-> Array (BandedHermitian offDiag size) a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape
      (\(size
size, ShapeStatic.ZeroBased Proxy (Un offDiag)
k ::+ ()) ->
         Proxy (Un offDiag) -> Order -> size -> BandedHermitian offDiag size
forall off size.
UnaryProxy off -> Order -> size -> BandedHermitian off size
Layout.BandedHermitian Proxy (Un offDiag)
k Order
ColumnMajor size
size)

identity ::
   (Shape.C sh, Class.Floating a) => sh -> Diagonal sh a
identity :: forall sh a. (C sh, Floating a) => sh -> Diagonal sh a
identity =
   (sh -> BandedHermitian U0 sh)
-> Array sh a -> Array (BandedHermitian U0 sh) a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape (UnaryProxy U0 -> Order -> sh -> BandedHermitian U0 sh
forall off size.
UnaryProxy off -> Order -> size -> BandedHermitian off size
Layout.BandedHermitian UnaryProxy U0
forall a. Proxy a
Proxy Order
ColumnMajor) (Array sh a -> Array (BandedHermitian U0 sh) a)
-> (sh -> Array sh a) -> sh -> Array (BandedHermitian U0 sh) a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. sh -> Array sh a
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.one

diagonal ::
   (Shape.C sh, Class.Floating a) => Vector sh (RealOf a) -> Diagonal sh a
diagonal :: forall sh a.
(C sh, Floating a) =>
Vector sh (RealOf a) -> Diagonal sh a
diagonal =
   (sh -> BandedHermitian U0 sh) -> Array sh a -> Diagonal sh a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape (UnaryProxy U0 -> Order -> sh -> BandedHermitian U0 sh
forall off size.
UnaryProxy off -> Order -> size -> BandedHermitian off size
Layout.BandedHermitian UnaryProxy U0
forall a. Proxy a
Proxy Order
ColumnMajor) (Array sh a -> Diagonal sh a)
-> (Vector sh (RealOf a) -> Array sh a)
-> Vector sh (RealOf a)
-> Diagonal sh a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Vector sh (RealOf a) -> Array sh a
forall sh a.
(C sh, Floating a) =>
Vector sh (RealOf a) -> Vector sh a
Vector.fromReal

takeDiagonal ::
   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>
   BandedHermitian offDiag size a -> Vector size (RealOf a)
takeDiagonal :: forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
BandedHermitian offDiag size a -> Vector size (RealOf a)
takeDiagonal =
   TakeDiagonal (Array (BandedHermitian offDiag size)) (Array size) a
-> Array (BandedHermitian offDiag size) a -> Array size (RealOf a)
forall (f :: * -> *) (g :: * -> *) a.
TakeDiagonal f g a -> f a -> g (RealOf a)
runTakeDiagonal (TakeDiagonal (Array (BandedHermitian offDiag size)) (Array size) a
 -> Array (BandedHermitian offDiag size) a -> Array size (RealOf a))
-> TakeDiagonal
     (Array (BandedHermitian offDiag size)) (Array size) a
-> Array (BandedHermitian offDiag size) a
-> Array size (RealOf a)
forall a b. (a -> b) -> a -> b
$
   TakeDiagonal
  (Array (BandedHermitian offDiag size)) (Array size) Float
-> TakeDiagonal
     (Array (BandedHermitian offDiag size)) (Array size) Double
-> TakeDiagonal
     (Array (BandedHermitian offDiag size)) (Array size) (Complex Float)
-> TakeDiagonal
     (Array (BandedHermitian offDiag size))
     (Array size)
     (Complex Double)
-> TakeDiagonal
     (Array (BandedHermitian offDiag size)) (Array size) a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
forall (f :: * -> *).
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      ((Array (BandedHermitian offDiag size) Float
 -> Array size (RealOf Float))
-> TakeDiagonal
     (Array (BandedHermitian offDiag size)) (Array size) Float
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (RealOf a)) -> TakeDiagonal f g a
TakeDiagonal Array (BandedHermitian offDiag size) Float -> Vector size Float
Array (BandedHermitian offDiag size) Float
-> Array size (RealOf Float)
forall offDiag size a ar.
(Natural offDiag, C size, Floating a, RealOf a ~ ar, Real ar) =>
BandedHermitian offDiag size a -> Vector size ar
takeDiagonalAux) ((Array (BandedHermitian offDiag size) Double
 -> Array size (RealOf Double))
-> TakeDiagonal
     (Array (BandedHermitian offDiag size)) (Array size) Double
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (RealOf a)) -> TakeDiagonal f g a
TakeDiagonal Array (BandedHermitian offDiag size) Double -> Vector size Double
Array (BandedHermitian offDiag size) Double
-> Array size (RealOf Double)
forall offDiag size a ar.
(Natural offDiag, C size, Floating a, RealOf a ~ ar, Real ar) =>
BandedHermitian offDiag size a -> Vector size ar
takeDiagonalAux)
      ((Array (BandedHermitian offDiag size) (Complex Float)
 -> Array size (RealOf (Complex Float)))
-> TakeDiagonal
     (Array (BandedHermitian offDiag size)) (Array size) (Complex Float)
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (RealOf a)) -> TakeDiagonal f g a
TakeDiagonal Array (BandedHermitian offDiag size) (Complex Float)
-> Vector size Float
Array (BandedHermitian offDiag size) (Complex Float)
-> Array size (RealOf (Complex Float))
forall offDiag size a ar.
(Natural offDiag, C size, Floating a, RealOf a ~ ar, Real ar) =>
BandedHermitian offDiag size a -> Vector size ar
takeDiagonalAux) ((Array (BandedHermitian offDiag size) (Complex Double)
 -> Array size (RealOf (Complex Double)))
-> TakeDiagonal
     (Array (BandedHermitian offDiag size))
     (Array size)
     (Complex Double)
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (RealOf a)) -> TakeDiagonal f g a
TakeDiagonal Array (BandedHermitian offDiag size) (Complex Double)
-> Vector size Double
Array (BandedHermitian offDiag size) (Complex Double)
-> Array size (RealOf (Complex Double))
forall offDiag size a ar.
(Natural offDiag, C size, Floating a, RealOf a ~ ar, Real ar) =>
BandedHermitian offDiag size a -> Vector size ar
takeDiagonalAux)

takeDiagonalAux ::
   (Unary.Natural offDiag, Shape.C size,
    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   BandedHermitian offDiag size a -> Vector size ar
takeDiagonalAux :: forall offDiag size a ar.
(Natural offDiag, C size, Floating a, RealOf a ~ ar, Real ar) =>
BandedHermitian offDiag size a -> Vector size ar
takeDiagonalAux (Array (Layout.BandedHermitian UnaryProxy offDiag
numOff Order
order size
size) ForeignPtr a
a) =
   let k :: Int
k = UnaryProxy offDiag -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy offDiag
numOff
   in size -> (Int -> Ptr ar -> IO ()) -> Array size ar
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize size
size ((Int -> Ptr ar -> IO ()) -> Array size ar)
-> (Int -> Ptr ar -> IO ()) -> Array size ar
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr ar
yPtr -> ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
         Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
         Ptr a
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
         let xPtr :: Ptr (RealOf a)
xPtr =
               Ptr a -> Ptr (RealOf a)
forall a. Ptr a -> Ptr (RealOf a)
realPtr (Ptr a -> Ptr (RealOf a)) -> Ptr a -> Ptr (RealOf a)
forall a b. (a -> b) -> a -> b
$ Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aPtr (Int -> Ptr a) -> Int -> Ptr a
forall a b. (a -> b) -> a -> b
$
               case Order
order of
                  Order
RowMajor -> Int
0
                  Order
ColumnMajor -> Int
k
         Ptr CInt
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (Ptr a -> Int -> Int -> Int
forall a (f :: * -> *) b. Floating a => f a -> b -> b -> b
caseRealComplexFunc Ptr a
aPtr Int
1 Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1))
         Ptr CInt
incyPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
         IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ Ptr CInt -> Ptr ar -> Ptr CInt -> Ptr ar -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.copy Ptr CInt
nPtr Ptr ar
Ptr (RealOf a)
xPtr Ptr CInt
incxPtr Ptr ar
yPtr Ptr CInt
incyPtr


toHermitian ::
   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>
   BandedHermitian offDiag size a -> Hermitian size a
toHermitian :: forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
BandedHermitian offDiag size a -> Hermitian size a
toHermitian (Array (Layout.BandedHermitian UnaryProxy offDiag
numOff Order
order size
size) ForeignPtr a
a) =
   Hermitian size
-> (Int -> Ptr a -> IO ()) -> Array (Hermitian size) a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (Order -> size -> Hermitian size
forall size. Order -> size -> Hermitian size
Layout.hermitian Order
order size
size) ((Int -> Ptr a -> IO ()) -> Array (Hermitian size) a)
-> (Int -> Ptr a -> IO ()) -> Array (Hermitian size) a
forall a b. (a -> b) -> a -> b
$
   Int -> Order -> Int -> ForeignPtr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Int -> Order -> Int -> ForeignPtr a -> Int -> Ptr a -> IO ()
Mos.fromBanded (UnaryProxy offDiag -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy offDiag
numOff) Order
order (size -> Int
forall sh. C sh => sh -> Int
Shape.size size
size) ForeignPtr a
a


toBanded ::
   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>
   BandedHermitian offDiag size a ->
   Banded.Square offDiag offDiag size a
toBanded :: forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
BandedHermitian offDiag size a -> Square offDiag offDiag size a
toBanded (Array (Layout.BandedHermitian UnaryProxy offDiag
numOff Order
order size
sh) ForeignPtr a
a) =
   BandedSquare offDiag offDiag size
-> (Ptr a -> IO ()) -> Array (BandedSquare offDiag offDiag size) a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate
      ((UnaryProxy offDiag, UnaryProxy offDiag)
-> Order
-> Extent Shape Small Small size size
-> BandedSquare offDiag offDiag size
forall sub super meas vert horiz height width.
(UnaryProxy sub, UnaryProxy super)
-> Order
-> Extent meas vert horiz height width
-> Banded sub super meas vert horiz height width
Layout.Banded (UnaryProxy offDiag
numOff,UnaryProxy offDiag
numOff) Order
order (size -> Extent Shape Small Small size size
forall sh. sh -> Square sh
Extent.square size
sh)) ((Ptr a -> IO ()) -> Array (BandedSquare offDiag offDiag size) a)
-> (Ptr a -> IO ()) -> Array (BandedSquare offDiag offDiag size) 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 ->
      let n :: Int
n = size -> Int
forall sh. C sh => sh -> Int
Shape.size size
sh
          k :: Int
k = UnaryProxy offDiag -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy offDiag
numOff
          lda :: Int
lda = Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
          ldb :: Int
ldb = Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
      in case Order
order of
            Order
ColumnMajor -> do
               Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copySubMatrix Int
lda Int
n Int
lda Ptr a
aPtr Int
ldb Ptr a
bPtr
               (Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
(Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
columnToRowMajor 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
                  (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Int
k
                  Int
lda (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aPtr (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1))
                  Int
ldb (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
bPtr (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1))
               a -> Int -> Ptr a -> IO ()
forall a. Floating a => a -> Int -> Ptr a -> IO ()
fill a
forall a. Floating a => a
zero Int
k (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
bPtr (Int
ldbInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k))
            Order
RowMajor -> do
               Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copySubMatrix Int
lda Int
n Int
lda Ptr a
aPtr Int
ldb (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
bPtr Int
k)
               a -> Int -> Ptr a -> IO ()
forall a. Floating a => a -> Int -> Ptr a -> IO ()
fill a
forall a. Floating a => a
zero Int
k Ptr a
bPtr
               (Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
(Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
rowToColumnMajor 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
                  (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Int
k
                  Int
lda (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aPtr Int
1)
                  Int
ldb (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
bPtr Int
ldb)

forceOrder ::
   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>
   Order ->
   BandedHermitian offDiag size a ->
   BandedHermitian offDiag size a
forceOrder :: forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
Order
-> BandedHermitian offDiag size a -> BandedHermitian offDiag size a
forceOrder Order
newOrder BandedHermitian offDiag size a
a =
   if Order
newOrder Order -> Order -> Bool
forall a. Eq a => a -> a -> Bool
== BandedHermitian offDiag size -> Order
forall off size. BandedHermitian off size -> Order
Layout.bandedHermitianOrder (BandedHermitian offDiag size a -> BandedHermitian offDiag size
forall sh a. Array sh a -> sh
Array.shape BandedHermitian offDiag size a
a)
      then BandedHermitian offDiag size a
a
      else BandedHermitian offDiag size a -> BandedHermitian offDiag size a
forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
BandedHermitian offDiag size a -> BandedHermitian offDiag size a
flipOrder BandedHermitian offDiag size a
a

flipOrder ::
   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>
   BandedHermitian offDiag size a ->
   BandedHermitian offDiag size a
flipOrder :: forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
BandedHermitian offDiag size a -> BandedHermitian offDiag size a
flipOrder (Array (Layout.BandedHermitian UnaryProxy offDiag
numOff Order
order size
sh) ForeignPtr a
a) =
   BandedHermitian offDiag size
-> (Ptr a -> IO ()) -> Array (BandedHermitian offDiag size) a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate
      (UnaryProxy offDiag -> Order -> size -> BandedHermitian offDiag size
forall off size.
UnaryProxy off -> Order -> size -> BandedHermitian off size
Layout.BandedHermitian UnaryProxy offDiag
numOff (Order -> Order
Layout.flipOrder Order
order) size
sh) ((Ptr a -> IO ()) -> Array (BandedHermitian offDiag size) a)
-> (Ptr a -> IO ()) -> Array (BandedHermitian offDiag size) 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 ->
      let n :: Int
n = size -> Int
forall sh. C sh => sh -> Int
Shape.size size
sh
          k :: Int
k = UnaryProxy offDiag -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy offDiag
numOff Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
      in case Order
order of
            Order
ColumnMajor -> (Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
(Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
columnToRowMajor 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 Int
n Int
k Int
k Ptr a
aPtr Int
k Ptr a
bPtr
            Order
RowMajor -> (Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
(Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
rowToColumnMajor 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 Int
n Int
k Int
k Ptr a
aPtr Int
k Ptr a
bPtr

columnToRowMajor ::
   (Class.Floating a) =>
   (Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()) ->
   Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
columnToRowMajor :: forall a.
Floating a =>
(Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
columnToRowMajor Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
copy Int
n Int
k Int
lda Ptr a
aPtr Int
ldb 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
   Ptr CInt
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (Int
ldaInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
   Ptr CInt
incyPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
   Ptr CInt
inczPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
   Ptr a
zPtr <- 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
nPtr <- FortranIO () (Ptr CInt)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
   IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ [Int] -> (Int -> IO ()) -> IO ()
forall (t :: * -> *) (f :: * -> *) a b.
(Foldable t, Applicative f) =>
t a -> (a -> f b) -> f ()
for_ (Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
take Int
n [Int
0..]) ((Int -> IO ()) -> IO ()) -> (Int -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Int
i -> do
      let split :: Int
split = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
k (Int -> Int) -> Int -> Int
forall a b. (a -> b) -> a -> b
$ Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
i
      let xPtr :: Ptr a
xPtr = Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aPtr (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
lda Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
      let yPtr :: Ptr a
yPtr = Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
bPtr (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
ldb)
      Ptr CInt -> Int -> IO ()
pokeCInt Ptr CInt
nPtr Int
split
      Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
copy Ptr CInt
nPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr
      Ptr CInt -> Int -> IO ()
pokeCInt Ptr CInt
nPtr (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
split)
      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
zPtr Ptr CInt
inczPtr (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
yPtr Int
split) Ptr CInt
incyPtr

rowToColumnMajor ::
   (Class.Floating a) =>
   (Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()) ->
   Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
rowToColumnMajor :: forall a.
Floating a =>
(Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
rowToColumnMajor Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
copy Int
n Int
k Int
lda Ptr a
aPtr Int
ldb 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
   Ptr CInt
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (Int
ldaInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
   Ptr CInt
incyPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
   Ptr CInt
inczPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
   Ptr a
zPtr <- 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
nPtr <- FortranIO () (Ptr CInt)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
   IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ [Int] -> (Int -> IO ()) -> IO ()
forall (t :: * -> *) (f :: * -> *) a b.
(Foldable t, Applicative f) =>
t a -> (a -> f b) -> f ()
for_ (Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
take Int
n [Int
0..]) ((Int -> IO ()) -> IO ()) -> (Int -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Int
i -> do
      let split :: Int
split = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
      let xPtr :: Ptr a
xPtr = Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aPtr (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
lda Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (Int
splitInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k)Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
ldaInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))
      let yPtr :: Ptr a
yPtr = Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
bPtr (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
ldb)
      Ptr CInt -> Int -> IO ()
pokeCInt Ptr CInt
nPtr Int
split
      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
zPtr Ptr CInt
inczPtr Ptr a
yPtr Ptr CInt
incyPtr
      Ptr CInt -> Int -> IO ()
pokeCInt Ptr CInt
nPtr (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
split)
      Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
copy Ptr CInt
nPtr Ptr a
xPtr Ptr CInt
incxPtr (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
yPtr Int
split) Ptr CInt
incyPtr


takeTopLeft ::
   (Unary.Natural offDiag, Shape.C sh0, Shape.C sh1, Class.Floating a) =>
   BandedHermitian offDiag (sh0 ::+ sh1) a ->
   BandedHermitian offDiag sh0 a
takeTopLeft :: forall offDiag sh0 sh1 a.
(Natural offDiag, C sh0, C sh1, Floating a) =>
BandedHermitian offDiag (sh0 ::+ sh1) a
-> BandedHermitian offDiag sh0 a
takeTopLeft
   (Array (Layout.BandedHermitian UnaryProxy offDiag
numOff Order
order (sh0
sh0 ::+ sh1
_sh1)) ForeignPtr a
a) =

   BandedHermitian offDiag sh0
-> (Int -> Ptr a -> IO ()) -> Array (BandedHermitian offDiag sh0) a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (UnaryProxy offDiag -> Order -> sh0 -> BandedHermitian offDiag sh0
forall off size.
UnaryProxy off -> Order -> size -> BandedHermitian off size
Layout.BandedHermitian UnaryProxy offDiag
numOff Order
order sh0
sh0) ((Int -> Ptr a -> IO ()) -> Array (BandedHermitian offDiag sh0) a)
-> (Int -> Ptr a -> IO ()) -> Array (BandedHermitian offDiag sh0) a
forall a b. (a -> b) -> a -> b
$
      \Int
n 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 -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
copyBlock Int
n Ptr a
aPtr Ptr a
bPtr

takeBottomRight ::
   (Unary.Natural offDiag, Shape.C sh0, Shape.C sh1, Class.Floating a) =>
   BandedHermitian offDiag (sh0 ::+ sh1) a ->
   BandedHermitian offDiag sh1 a
takeBottomRight :: forall offDiag sh0 sh1 a.
(Natural offDiag, C sh0, C sh1, Floating a) =>
BandedHermitian offDiag (sh0 ::+ sh1) a
-> BandedHermitian offDiag sh1 a
takeBottomRight
   (Array (Layout.BandedHermitian UnaryProxy offDiag
numOff Order
order (sh0
sh0 ::+ sh1
sh1)) ForeignPtr a
a) =

   BandedHermitian offDiag sh1
-> (Int -> Ptr a -> IO ()) -> Array (BandedHermitian offDiag sh1) a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (UnaryProxy offDiag -> Order -> sh1 -> BandedHermitian offDiag sh1
forall off size.
UnaryProxy off -> Order -> size -> BandedHermitian off size
Layout.BandedHermitian UnaryProxy offDiag
numOff Order
order sh1
sh1) ((Int -> Ptr a -> IO ()) -> Array (BandedHermitian offDiag sh1) a)
-> (Int -> Ptr a -> IO ()) -> Array (BandedHermitian offDiag sh1) a
forall a b. (a -> b) -> a -> b
$
      \Int
n 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 ->
      Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
copyBlock Int
n
         (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aPtr (Int -> Ptr a) -> Int -> Ptr a
forall a b. (a -> b) -> a -> b
$ (UnaryProxy offDiag -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy offDiag
numOff Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
* sh0 -> Int
forall sh. C sh => sh -> Int
Shape.size sh0
sh0)
         Ptr a
bPtr


multiplyVector ::
   (Unary.Natural offDiag, Shape.C size, Eq size, Class.Floating a) =>
   Transposition -> BandedHermitian offDiag size a ->
   Vector size a -> Vector size a
multiplyVector :: forall offDiag size a.
(Natural offDiag, C size, Eq size, Floating a) =>
Transposition
-> BandedHermitian offDiag size a -> Vector size a -> Vector size a
multiplyVector Transposition
transposed
   (Array (Layout.BandedHermitian UnaryProxy offDiag
numOff Order
order size
size) ForeignPtr a
a) (Array size
sizeX ForeignPtr a
x) =
      size -> (Int -> Ptr a -> IO ()) -> Array size a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize size
size ((Int -> Ptr a -> IO ()) -> Array size a)
-> (Int -> Ptr a -> IO ()) -> Array size a
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
yPtr -> do

   String -> Bool -> IO ()
Call.assert String
"BandedHermitian.multiplyVector: shapes mismatch"
      (size
size size -> size -> Bool
forall a. Eq a => a -> a -> Bool
== size
sizeX)
   let k :: Int
k = UnaryProxy offDiag -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy offDiag
numOff
   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 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 -> FortranIO () (Ptr CInt)) -> Int -> FortranIO () (Ptr CInt)
forall a b. (a -> b) -> a -> b
$ Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
      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 a. IO a -> ContT () IO a
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 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 CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
BlasGen.hbmv Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
kPtr
            Ptr a
alphaPtr Ptr a
aPtr Ptr CInt
ldaPtr 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


gramian ::
   (Shape.C size, Eq size, Class.Floating a,
    Unary.Natural sub, Unary.Natural super) =>
   Banded.Square sub super size a ->
   BandedHermitian (sub :+: super) size a
gramian :: forall size a sub super.
(C size, Eq size, Floating a, Natural sub, Natural super) =>
Square sub super size a -> BandedHermitian (sub :+: super) size a
gramian Square sub super size a
a =
   case (UnaryProxy sub -> Nat sub, UnaryProxy super -> Nat super)
-> (UnaryProxy sub, UnaryProxy super) -> (Nat sub, Nat super)
forall a c b d. (a -> c, b -> d) -> (a, b) -> (c, d)
mapPair (UnaryProxy sub -> Nat sub
forall n. Natural n => UnaryProxy n -> Nat n
natFromProxy,UnaryProxy super -> Nat super
forall n. Natural n => UnaryProxy n -> Nat n
natFromProxy) ((UnaryProxy sub, UnaryProxy super) -> (Nat sub, Nat super))
-> (UnaryProxy sub, UnaryProxy super) -> (Nat sub, Nat super)
forall a b. (a -> b) -> a -> b
$
        Banded sub super Shape Small Small size size
-> (UnaryProxy sub, UnaryProxy super)
forall sub super meas vert horiz height width.
Banded sub super meas vert horiz height width
-> (UnaryProxy sub, UnaryProxy super)
Layout.bandedOffDiagonals (Banded sub super Shape Small Small size size
 -> (UnaryProxy sub, UnaryProxy super))
-> Banded sub super Shape Small Small size size
-> (UnaryProxy sub, UnaryProxy super)
forall a b. (a -> b) -> a -> b
$ Square sub super size a
-> Banded sub super Shape Small Small size size
forall sh a. Array sh a -> sh
Array.shape Square sub super size a
a of
      (Nat sub
sub,Nat super
super) ->
         case (Nat sub -> Nat super -> Nat (sub :+: super)
forall x y. Nat x -> Nat y -> Nat (x :+: y)
Proof.addNat Nat sub
sub Nat super
super, Nat sub -> Nat super -> AddComm sub super
forall x y. Nat x -> Nat y -> AddComm x y
Proof.addComm Nat sub
sub Nat super
super) of
            (Nat (sub :+: super)
Proof.Nat, AddComm sub super
Proof.AddComm) ->
               Square (sub :+: super) (sub :+: super) size a
-> BandedHermitian (sub :+: super) size a
forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
Square offDiag offDiag size a -> BandedHermitian offDiag size a
fromUpperPart (Square (sub :+: super) (sub :+: super) size a
 -> BandedHermitian (sub :+: super) size a)
-> Square (sub :+: super) (sub :+: super) size a
-> BandedHermitian (sub :+: super) size a
forall a b. (a -> b) -> a -> b
$ Banded super sub Shape Small Small size size a
-> Square sub super size a
-> Banded
     (super :+: sub) (super :+: sub) Shape Small Small size size a
forall subA superA subB superB subC superC meas vert horiz height
       width fuse a.
(Natural subA, Natural superA, Natural subB, Natural superB,
 (subA :+: subB) ~ subC, (superA :+: superB) ~ superC, Measure meas,
 C vert, C horiz, C height, C width, C fuse, Eq fuse, Floating a) =>
Banded subA superA meas vert horiz height fuse a
-> Banded subB superB meas vert horiz fuse width a
-> Banded subC superC meas vert horiz height width a
Banded.multiply (Square sub super size a
-> Banded super sub Shape Small Small size size a
forall sub super meas vert horiz height width a.
(Natural sub, Natural super, Measure meas, C vert, C horiz,
 C height, C width, Floating a) =>
Banded sub super meas vert horiz height width a
-> Banded super sub meas horiz vert width height a
Banded.adjoint Square sub super size a
a) Square sub super size a
a

fromUpperPart ::
   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>
   Banded.Square offDiag offDiag size a -> BandedHermitian offDiag size a
fromUpperPart :: forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
Square offDiag offDiag size a -> BandedHermitian offDiag size a
fromUpperPart (Array (Layout.Banded (UnaryProxy offDiag
sub,UnaryProxy offDiag
super) Order
order Extent Shape Small Small size size
extent) ForeignPtr a
a) =
   let sh :: size
sh = Extent Shape Small Small size size -> size
forall shape. Square shape -> shape
Extent.squareSize Extent Shape Small Small size size
extent
       n :: Int
n = size -> Int
forall sh. C sh => sh -> Int
Shape.size size
sh
       kl :: Int
kl = UnaryProxy offDiag -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy offDiag
sub
       ku :: Int
ku = UnaryProxy offDiag -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy offDiag
super
       lda :: Int
lda = Int
klInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
ku
       ldb :: Int
ldb = Int
kuInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
   in BandedHermitian offDiag size
-> (Ptr a -> IO ()) -> Array (BandedHermitian offDiag size) a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (UnaryProxy offDiag -> Order -> size -> BandedHermitian offDiag size
forall off size.
UnaryProxy off -> Order -> size -> BandedHermitian off size
Layout.BandedHermitian UnaryProxy offDiag
super Order
order size
sh) ((Ptr a -> IO ()) -> Array (BandedHermitian offDiag size) a)
-> (Ptr a -> IO ()) -> Array (BandedHermitian offDiag size) 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 ->
      case Order
order of
         Order
ColumnMajor -> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copySubMatrix Int
ldb Int
n Int
lda Ptr a
aPtr Int
ldb Ptr a
bPtr
         Order
RowMajor -> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copySubMatrix Int
ldb Int
n Int
lda (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aPtr Int
kl) Int
ldb Ptr a
bPtr


multiplyFull ::
   (Unary.Natural offDiag, Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Transposition -> BandedHermitian offDiag height a ->
   Matrix.Full meas vert horiz height width a ->
   Matrix.Full meas vert horiz height width a
multiplyFull :: forall offDiag meas vert horiz height width a.
(Natural offDiag, Measure meas, C vert, C horiz, C height,
 Eq height, C width, Floating a) =>
Transposition
-> BandedHermitian offDiag height a
-> Full meas vert horiz height width a
-> Full meas vert horiz height width a
multiplyFull Transposition
transposed BandedHermitian offDiag height a
a Full meas vert horiz height width a
b =
   case Full meas vert horiz height width -> Order
forall meas vert horiz height width.
Full meas vert horiz height width -> Order
Layout.fullOrder (Full meas vert horiz height width -> Order)
-> Full meas vert horiz height width -> Order
forall a b. (a -> b) -> a -> b
$ Full meas vert horiz height width a
-> Full meas vert horiz height width
forall sh a. Array sh a -> sh
Array.shape Full meas vert horiz height width a
b of
      Order
ColumnMajor -> Transposition
-> BandedHermitian offDiag height a
-> Full meas vert horiz height width a
-> Full meas vert horiz height width a
forall offDiag meas vert horiz height width a.
(Natural offDiag, Measure meas, C vert, C horiz, Eq height,
 C height, C width, Floating a) =>
Transposition
-> BandedHermitian offDiag height a
-> Full meas vert horiz height width a
-> Full meas vert horiz height width a
multiplyFullSpecial Transposition
transposed BandedHermitian offDiag height a
a Full meas vert horiz height width a
b
      Order
RowMajor -> Transposition
-> BandedHermitian offDiag height a
-> Full meas vert horiz height width a
-> Full meas vert horiz height width a
forall offDiag meas vert horiz height width a.
(Natural offDiag, Measure meas, C vert, C horiz, C height,
 Eq height, C width, Floating a) =>
Transposition
-> BandedHermitian offDiag height a
-> Full meas vert horiz height width a
-> Full meas vert horiz height width a
multiplyFullGeneric Transposition
transposed BandedHermitian offDiag height a
a Full meas vert horiz height width a
b

multiplyFullSpecial ::
   (Unary.Natural offDiag, Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Eq height, Shape.C height, Shape.C width, Class.Floating a) =>
   Transposition -> BandedHermitian offDiag height a ->
   Matrix.Full meas vert horiz height width a ->
   Matrix.Full meas vert horiz height width a
multiplyFullSpecial :: forall offDiag meas vert horiz height width a.
(Natural offDiag, Measure meas, C vert, C horiz, Eq height,
 C height, C width, Floating a) =>
Transposition
-> BandedHermitian offDiag height a
-> Full meas vert horiz height width a
-> Full meas vert horiz height width a
multiplyFullSpecial Transposition
transposed
      (Array (Layout.BandedHermitian UnaryProxy offDiag
numOff Order
orderA height
sizeA) ForeignPtr a
a)
      (Array (Layout.Full Order
orderB Extent meas vert horiz height width
extentB) ForeignPtr a
b) =
   Full meas vert horiz height width
-> (Ptr a -> IO ()) -> Array (Full meas vert horiz height width) a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order
-> Extent meas vert horiz height width
-> Full meas vert horiz height width
forall meas vert horiz height width.
Order
-> Extent meas vert horiz height width
-> Full meas vert horiz height width
Layout.Full Order
orderB Extent meas vert horiz height width
extentB) ((Ptr a -> IO ()) -> Array (Full meas vert horiz height width) a)
-> (Ptr a -> IO ()) -> Array (Full meas vert horiz height width) a
forall a b. (a -> b) -> a -> b
$ \Ptr a
cPtr -> do
      String -> Bool -> IO ()
Call.assert String
"BandedHermitian.multiplyFull: shapes mismatch"
         (height
sizeA height -> height -> Bool
forall a. Eq a => a -> a -> Bool
== Extent meas vert horiz height width -> height
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> height
Extent.height Extent meas vert horiz height width
extentB)
      let (height
height,width
width) = Extent meas vert horiz height width -> (height, width)
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> (height, width)
Extent.dimensions Extent meas vert horiz height width
extentB
      case Order
orderB of
         Order
ColumnMajor ->
            Transposition
-> UnaryProxy offDiag
-> (height, width)
-> Order
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> IO ()
forall offDiag height width a.
(Natural offDiag, C height, C width, Floating a) =>
Transposition
-> UnaryProxy offDiag
-> (height, width)
-> Order
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> IO ()
multiplyFullColumnMajor
               Transposition
transposed UnaryProxy offDiag
numOff (height
height,width
width) Order
orderA ForeignPtr a
a ForeignPtr a
b Ptr a
cPtr
         Order
RowMajor ->
            Transposition
-> UnaryProxy offDiag
-> (height, width)
-> Order
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> IO ()
forall offDiag height width a.
(Natural offDiag, C height, C width, Floating a) =>
Transposition
-> UnaryProxy offDiag
-> (height, width)
-> Order
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> IO ()
multiplyFullRowMajor
               Transposition
transposed UnaryProxy offDiag
numOff (height
height,width
width) Order
orderA ForeignPtr a
a ForeignPtr a
b Ptr a
cPtr

multiplyFullColumnMajor ::
   (Unary.Natural offDiag, Shape.C height, Shape.C width, Class.Floating a) =>
   Transposition -> UnaryProxy offDiag -> (height, width) ->
   Order -> ForeignPtr a -> ForeignPtr a -> Ptr a -> IO ()
multiplyFullColumnMajor :: forall offDiag height width a.
(Natural offDiag, C height, C width, Floating a) =>
Transposition
-> UnaryProxy offDiag
-> (height, width)
-> Order
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> IO ()
multiplyFullColumnMajor Transposition
transposed UnaryProxy offDiag
numOff (height
height,width
width) Order
order ForeignPtr a
a ForeignPtr a
b Ptr a
cPtr = do
   let n :: Int
n = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
   let nrhs :: Int
nrhs = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
width
   let k :: Int
k = UnaryProxy offDiag -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy offDiag
numOff
   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 -> 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 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 -> FortranIO () (Ptr CInt)) -> Int -> FortranIO () (Ptr CInt)
forall a b. (a -> b) -> a -> b
$ Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
      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
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
      let pointers :: [(Ptr a, Ptr a)]
pointers = Int -> [(Ptr a, Ptr a)] -> [(Ptr a, Ptr a)]
forall a. Int -> [a] -> [a]
take Int
nrhs ([(Ptr a, Ptr a)] -> [(Ptr a, Ptr a)])
-> [(Ptr a, Ptr a)] -> [(Ptr a, Ptr a)]
forall a b. (a -> b) -> a -> b
$ [Ptr a] -> [Ptr a] -> [(Ptr a, Ptr a)]
forall a b. [a] -> [b] -> [(a, b)]
zip (Int -> Ptr a -> [Ptr a]
forall a. Storable a => Int -> Ptr a -> [Ptr a]
pointerSeq Int
n Ptr a
bPtr) (Int -> Ptr a -> [Ptr a]
forall a. Storable a => Int -> Ptr a -> [Ptr a]
pointerSeq Int
n Ptr a
cPtr)
      case Transposition -> Order -> Order
transposeOrder Transposition
transposed Order
order of
         Order
RowMajor -> do
            Ptr a
xPtr <- Int -> FortranIO () (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
            IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ [(Ptr a, Ptr a)] -> ((Ptr a, Ptr a) -> IO ()) -> IO ()
forall (t :: * -> *) (f :: * -> *) a b.
(Foldable t, Applicative f) =>
t a -> (a -> f b) -> f ()
for_ [(Ptr a, Ptr a)]
pointers (((Ptr a, Ptr a) -> IO ()) -> IO ())
-> ((Ptr a, Ptr a) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \(Ptr a
biPtr,Ptr a
yPtr) -> 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
biPtr Ptr CInt
incxPtr Ptr a
xPtr Ptr CInt
incxPtr
               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 CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
BlasGen.hbmv Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
kPtr
                  Ptr a
alphaPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
betaPtr Ptr a
yPtr Ptr CInt
incyPtr
               Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
lacgv Ptr CInt
nPtr Ptr a
yPtr Ptr CInt
incyPtr
         Order
ColumnMajor ->
            IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ [(Ptr a, Ptr a)] -> ((Ptr a, Ptr a) -> IO ()) -> IO ()
forall (t :: * -> *) (f :: * -> *) a b.
(Foldable t, Applicative f) =>
t a -> (a -> f b) -> f ()
for_ [(Ptr a, Ptr a)]
pointers (((Ptr a, Ptr a) -> IO ()) -> IO ())
-> ((Ptr a, Ptr a) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \(Ptr a
xPtr,Ptr a
yPtr) ->
               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 CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
BlasGen.hbmv Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
kPtr
                  Ptr a
alphaPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
betaPtr Ptr a
yPtr Ptr CInt
incyPtr

multiplyFullRowMajor ::
   (Unary.Natural offDiag, Shape.C height, Shape.C width, Class.Floating a) =>
   Transposition -> UnaryProxy offDiag -> (height, width) ->
   Order -> ForeignPtr a -> ForeignPtr a -> Ptr a -> IO ()
multiplyFullRowMajor :: forall offDiag height width a.
(Natural offDiag, C height, C width, Floating a) =>
Transposition
-> UnaryProxy offDiag
-> (height, width)
-> Order
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> IO ()
multiplyFullRowMajor =
   String
-> Transposition
-> UnaryProxy offDiag
-> (height, width)
-> Order
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> IO ()
forall a. HasCallStack => String -> a
error String
"BandedHermitian.multiplyFullRowMajor: not implemented"


multiplyFullGeneric ::
   (Unary.Natural offDiag, Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Transposition -> BandedHermitian offDiag height a ->
   Matrix.Full meas vert horiz height width a ->
   Matrix.Full meas vert horiz height width a
multiplyFullGeneric :: forall offDiag meas vert horiz height width a.
(Natural offDiag, Measure meas, C vert, C horiz, C height,
 Eq height, C width, Floating a) =>
Transposition
-> BandedHermitian offDiag height a
-> Full meas vert horiz height width a
-> Full meas vert horiz height width a
multiplyFullGeneric Transposition
transposed BandedHermitian offDiag height a
a Full meas vert horiz height width a
b =
   let (Square offDiag U0 height a
lower,Square U0 offDiag height a
upper) = (BandedHermitian offDiag height a -> Square offDiag U0 height a
forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
BandedHermitian offDiag size a -> Square offDiag U0 size a
takeStrictLower BandedHermitian offDiag height a
a, BandedHermitian offDiag height a -> Square U0 offDiag height a
forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
BandedHermitian offDiag size a -> Square U0 offDiag size a
takeUpper BandedHermitian offDiag height a
a)
       (Square offDiag U0 height a
lowerT,Square U0 offDiag height a
upperT) =
         case Transposition
transposed of
            Transposition
Transposed -> (Square U0 offDiag height a -> Square offDiag U0 height a
forall meas vert horiz sub super height width a.
(Measure meas, C vert, C horiz) =>
Banded sub super meas vert horiz height width a
-> Banded super sub meas horiz vert width height a
Banded.transpose Square U0 offDiag height a
upper, Square offDiag U0 height a -> Square U0 offDiag height a
forall meas vert horiz sub super height width a.
(Measure meas, C vert, C horiz) =>
Banded sub super meas vert horiz height width a
-> Banded super sub meas horiz vert width height a
Banded.transpose Square offDiag U0 height a
lower)
            Transposition
NonTransposed -> (Square offDiag U0 height a
lower,Square U0 offDiag height a
upper)
   in Array (Unchecked (Full meas vert horiz height width)) a
-> Full meas vert horiz height width a
forall sh a. Array (Unchecked sh) a -> Array sh a
VectorPriv.recheck (Array (Unchecked (Full meas vert horiz height width)) a
 -> Full meas vert horiz height width a)
-> Array (Unchecked (Full meas vert horiz height width)) a
-> Full meas vert horiz height width a
forall a b. (a -> b) -> a -> b
$
      a
-> Array (Unchecked (Full meas vert horiz height width)) a
-> Array (Unchecked (Full meas vert horiz height width)) a
-> Array (Unchecked (Full meas vert horiz height width)) a
forall sh a.
(C sh, Eq sh, Floating a) =>
a -> Vector sh a -> Vector sh a -> Vector sh a
Vector.mac a
forall a. Floating a => a
one
         (Full meas vert horiz height width a
-> Array (Unchecked (Full meas vert horiz height width)) a
forall sh a. Array sh a -> Array (Unchecked sh) a
VectorPriv.uncheck (Full meas vert horiz height width a
 -> Array (Unchecked (Full meas vert horiz height width)) a)
-> Full meas vert horiz height width a
-> Array (Unchecked (Full meas vert horiz height width)) a
forall a b. (a -> b) -> a -> b
$
          Banded offDiag U0 meas vert horiz height height a
-> Full meas vert horiz height width a
-> Full meas vert horiz height width a
forall sub super meas vert horiz height width fuse a.
(Natural sub, Natural super, Measure meas, C vert, C horiz,
 C height, C width, C fuse, Eq fuse, Floating a) =>
Banded sub super meas vert horiz height fuse a
-> Full meas vert horiz fuse width a
-> Full meas vert horiz height width a
Banded.multiplyFull
            (Map Shape Small Small meas vert horiz height height
-> Square offDiag U0 height a
-> Banded offDiag U0 meas vert horiz height height a
forall vertA horizA vertB horizB measA measB height width sub super
       a.
(C vertA, C horizA, C vertB, C horizB) =>
Map measA vertA horizA measB vertB horizB height width
-> Banded sub super measA vertA horizA height width a
-> Banded sub super measB vertB horizB height width a
Banded.mapExtent Map Shape Small Small meas vert horiz height height
forall meas vert horiz size.
(Measure meas, C vert, C horiz) =>
Square size -> Extent meas vert horiz size size
ExtentPriv.fromSquare Square offDiag U0 height a
lowerT) Full meas vert horiz height width a
b)
         (Full meas vert horiz height width a
-> Array (Unchecked (Full meas vert horiz height width)) a
forall sh a. Array sh a -> Array (Unchecked sh) a
VectorPriv.uncheck (Full meas vert horiz height width a
 -> Array (Unchecked (Full meas vert horiz height width)) a)
-> Full meas vert horiz height width a
-> Array (Unchecked (Full meas vert horiz height width)) a
forall a b. (a -> b) -> a -> b
$
          Banded U0 offDiag meas vert horiz height height a
-> Full meas vert horiz height width a
-> Full meas vert horiz height width a
forall sub super meas vert horiz height width fuse a.
(Natural sub, Natural super, Measure meas, C vert, C horiz,
 C height, C width, C fuse, Eq fuse, Floating a) =>
Banded sub super meas vert horiz height fuse a
-> Full meas vert horiz fuse width a
-> Full meas vert horiz height width a
Banded.multiplyFull
            (Map Shape Small Small meas vert horiz height height
-> Square U0 offDiag height a
-> Banded U0 offDiag meas vert horiz height height a
forall vertA horizA vertB horizB measA measB height width sub super
       a.
(C vertA, C horizA, C vertB, C horizB) =>
Map measA vertA horizA measB vertB horizB height width
-> Banded sub super measA vertA horizA height width a
-> Banded sub super measB vertB horizB height width a
Banded.mapExtent Map Shape Small Small meas vert horiz height height
forall meas vert horiz size.
(Measure meas, C vert, C horiz) =>
Square size -> Extent meas vert horiz size size
ExtentPriv.fromSquare Square U0 offDiag height a
upperT) Full meas vert horiz height width a
b)

takeUpper ::
   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>
   BandedHermitian offDiag size a ->
   Banded.Square TypeNum.U0 offDiag size a
takeUpper :: forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
BandedHermitian offDiag size a -> Square U0 offDiag size a
takeUpper =
   (BandedHermitian offDiag size -> BandedSquare U0 offDiag size)
-> Array (BandedHermitian offDiag size) a
-> Array (BandedSquare U0 offDiag size) a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape
      (\(Layout.BandedHermitian UnaryProxy offDiag
numOff Order
order size
sh) ->
         (UnaryProxy U0, UnaryProxy offDiag)
-> Order -> size -> BandedSquare U0 offDiag size
forall sub super size.
(UnaryProxy sub, UnaryProxy super)
-> Order -> size -> BandedSquare sub super size
Layout.bandedSquare (UnaryProxy U0
forall a. Proxy a
Proxy,UnaryProxy offDiag
numOff) Order
order size
sh)

takeStrictLower ::
   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>
   BandedHermitian offDiag size a ->
   Banded.Square offDiag TypeNum.U0 size a
takeStrictLower :: forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
BandedHermitian offDiag size a -> Square offDiag U0 size a
takeStrictLower (Array (Layout.BandedHermitian UnaryProxy offDiag
numOff Order
order size
sh) ForeignPtr a
x) =
   BandedSquare offDiag U0 size
-> (Int -> Ptr a -> IO ())
-> Array (BandedSquare offDiag U0 size) a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize
      ((UnaryProxy offDiag, UnaryProxy U0)
-> Order -> size -> BandedSquare offDiag U0 size
forall sub super size.
(UnaryProxy sub, UnaryProxy super)
-> Order -> size -> BandedSquare sub super size
Layout.bandedSquare
         (UnaryProxy offDiag
numOff,UnaryProxy U0
forall a. Proxy a
Proxy) (Order -> Order
Layout.flipOrder Order
order) size
sh) ((Int -> Ptr a -> IO ()) -> Array (BandedSquare offDiag U0 size) a)
-> (Int -> Ptr a -> IO ())
-> Array (BandedSquare offDiag U0 size) a
forall a b. (a -> b) -> a -> b
$
      \Int
size Ptr a
yPtr -> 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 k :: Int
k = UnaryProxy offDiag -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy offDiag
numOff
   Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (Int -> FortranIO () (Ptr CInt)) -> Int -> FortranIO () (Ptr CInt)
forall a b. (a -> b) -> a -> b
$ size -> Int
forall sh. C sh => sh -> Int
Shape.size size
sh
   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
   Ptr CInt
sizePtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
size
   Ptr CInt
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
   Ptr CInt
incyPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
   Ptr CInt
inczPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
   Ptr CInt
ldbPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim (Int -> FortranIO () (Ptr CInt)) -> Int -> FortranIO () (Ptr CInt)
forall a b. (a -> b) -> a -> b
$ Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
   Ptr a
zPtr <- a -> ContT () IO (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
zero
   IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
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 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
sizePtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr
      let offset :: Int
offset = case Order
order of Order
ColumnMajor -> Int
k; Order
RowMajor -> Int
0
      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
zPtr Ptr CInt
inczPtr (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
yPtr Int
offset) Ptr CInt
ldbPtr


type StaticVector n = Vector (ShapeStatic.ZeroBased n)

sumRank1 ::
   (Unary.Natural k, Shape.Indexed sh, Class.Floating a) =>
   Order -> sh ->
   [(RealOf a, (Shape.Index sh, StaticVector (Unary.Succ k) a))] ->
   BandedHermitian k sh a
sumRank1 :: forall k sh a.
(Natural k, Indexed sh, Floating a) =>
Order
-> sh
-> [(RealOf a, (Index sh, StaticVector (Succ k) a))]
-> BandedHermitian k sh a
sumRank1 =
   SumRank1 k sh a
-> Order
-> sh
-> [(RealOf a, (Index sh, Array (ZeroBased (Succ k)) a))]
-> Array (BandedHermitian k sh) a
forall k sh a. SumRank1 k sh a -> SumRank1_ k sh (RealOf a) a
getSumRank1 (SumRank1 k sh a
 -> Order
 -> sh
 -> [(RealOf a, (Index sh, Array (ZeroBased (Succ k)) a))]
 -> Array (BandedHermitian k sh) a)
-> SumRank1 k sh a
-> Order
-> sh
-> [(RealOf a, (Index sh, Array (ZeroBased (Succ k)) a))]
-> Array (BandedHermitian k sh) a
forall a b. (a -> b) -> a -> b
$
   SumRank1 k sh Float
-> SumRank1 k sh Double
-> SumRank1 k sh (Complex Float)
-> SumRank1 k sh (Complex Double)
-> SumRank1 k sh a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
forall (f :: * -> *).
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (SumRank1_ k sh (RealOf Float) Float -> SumRank1 k sh Float
forall k sh a. SumRank1_ k sh (RealOf a) a -> SumRank1 k sh a
SumRank1 (SumRank1_ k sh (RealOf Float) Float -> SumRank1 k sh Float)
-> SumRank1_ k sh (RealOf Float) Float -> SumRank1 k sh Float
forall a b. (a -> b) -> a -> b
$ UnaryProxy k -> SumRank1_ k sh Float Float
forall k sh a ar.
(Natural k, Indexed sh, Floating a, RealOf a ~ ar, Real ar) =>
UnaryProxy k -> SumRank1_ k sh ar a
sumRank1Aux UnaryProxy k
forall a. Proxy a
Proxy)
      (SumRank1_ k sh (RealOf Double) Double -> SumRank1 k sh Double
forall k sh a. SumRank1_ k sh (RealOf a) a -> SumRank1 k sh a
SumRank1 (SumRank1_ k sh (RealOf Double) Double -> SumRank1 k sh Double)
-> SumRank1_ k sh (RealOf Double) Double -> SumRank1 k sh Double
forall a b. (a -> b) -> a -> b
$ UnaryProxy k -> SumRank1_ k sh Double Double
forall k sh a ar.
(Natural k, Indexed sh, Floating a, RealOf a ~ ar, Real ar) =>
UnaryProxy k -> SumRank1_ k sh ar a
sumRank1Aux UnaryProxy k
forall a. Proxy a
Proxy)
      (SumRank1_ k sh (RealOf (Complex Float)) (Complex Float)
-> SumRank1 k sh (Complex Float)
forall k sh a. SumRank1_ k sh (RealOf a) a -> SumRank1 k sh a
SumRank1 (SumRank1_ k sh (RealOf (Complex Float)) (Complex Float)
 -> SumRank1 k sh (Complex Float))
-> SumRank1_ k sh (RealOf (Complex Float)) (Complex Float)
-> SumRank1 k sh (Complex Float)
forall a b. (a -> b) -> a -> b
$ UnaryProxy k -> SumRank1_ k sh Float (Complex Float)
forall k sh a ar.
(Natural k, Indexed sh, Floating a, RealOf a ~ ar, Real ar) =>
UnaryProxy k -> SumRank1_ k sh ar a
sumRank1Aux UnaryProxy k
forall a. Proxy a
Proxy)
      (SumRank1_ k sh (RealOf (Complex Double)) (Complex Double)
-> SumRank1 k sh (Complex Double)
forall k sh a. SumRank1_ k sh (RealOf a) a -> SumRank1 k sh a
SumRank1 (SumRank1_ k sh (RealOf (Complex Double)) (Complex Double)
 -> SumRank1 k sh (Complex Double))
-> SumRank1_ k sh (RealOf (Complex Double)) (Complex Double)
-> SumRank1 k sh (Complex Double)
forall a b. (a -> b) -> a -> b
$ UnaryProxy k -> SumRank1_ k sh Double (Complex Double)
forall k sh a ar.
(Natural k, Indexed sh, Floating a, RealOf a ~ ar, Real ar) =>
UnaryProxy k -> SumRank1_ k sh ar a
sumRank1Aux UnaryProxy k
forall a. Proxy a
Proxy)

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

type SumRank1_ k sh ar a =
   Order -> sh ->
   [(ar, (Shape.Index sh, StaticVector (Unary.Succ k) a))] ->
   BandedHermitian k sh a

sumRank1Aux ::
   (Unary.Natural k, Shape.Indexed sh,
    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   UnaryProxy k -> SumRank1_ k sh ar a
sumRank1Aux :: forall k sh a ar.
(Natural k, Indexed sh, Floating a, RealOf a ~ ar, Real ar) =>
UnaryProxy k -> SumRank1_ k sh ar a
sumRank1Aux UnaryProxy k
numOff Order
order sh
size [(ar, (Index sh, StaticVector (Succ k) a))]
xs =
   BandedHermitian k sh
-> (Int -> Ptr a -> IO ()) -> Array (BandedHermitian k sh) a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize
      (UnaryProxy k -> Order -> sh -> BandedHermitian k sh
forall off size.
UnaryProxy off -> Order -> size -> BandedHermitian off size
Layout.BandedHermitian UnaryProxy k
numOff Order
order sh
size) ((Int -> Ptr a -> IO ()) -> Array (BandedHermitian k sh) a)
-> (Int -> Ptr a -> IO ()) -> Array (BandedHermitian k sh) a
forall a b. (a -> b) -> a -> b
$
         \Int
bSize 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
   let k :: Int
k = UnaryProxy k -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy k
numOff
   let n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
size
   let lda :: Int
lda = Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
   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
mPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
lda
   Ptr a
alphaPtr <- FortranIO () (Ptr a)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
   Ptr CInt
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
   Ptr CInt
kPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
k
   Ptr CInt
ldbPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
k
   Ptr CInt
bSizePtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
bSize
   IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () 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
bSize Ptr a
aPtr
      [(ar, (Index sh, StaticVector (Succ k) a))]
-> ((ar, (Index sh, StaticVector (Succ k) a)) -> IO ()) -> IO ()
forall (t :: * -> *) (f :: * -> *) a b.
(Foldable t, Applicative f) =>
t a -> (a -> f b) -> f ()
for_ [(ar, (Index sh, StaticVector (Succ k) a))]
xs (((ar, (Index sh, StaticVector (Succ k) a)) -> IO ()) -> IO ())
-> ((ar, (Index sh, StaticVector (Succ k) a)) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \(ar
alpha, (Index sh
offset, Array ZeroBased (Succ k)
_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
            let i :: Int
i = sh -> Index sh -> Int
forall sh. Indexed sh => sh -> Index sh -> Int
Shape.offset sh
size Index sh
offset
            String -> Bool -> IO ()
Call.assert String
"BandedHermitian.sumRank1: index too large" (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n)
            let bPtr :: Ptr a
bPtr = Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aPtr (Int
ldaInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
i)
            HBR_ (RealOf a) a
forall a. Floating a => HBR_ (RealOf a) a
hbr Order
order Int
k ar
RealOf a
alpha
               Ptr CChar
uploPtr Ptr CInt
mPtr Ptr CInt
kPtr Ptr a
alphaPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
bPtr Ptr CInt
incxPtr Ptr CInt
ldbPtr
      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
bSizePtr Ptr a
aPtr Ptr CInt
incxPtr


type HBR_ ar a =
   Order -> Int -> ar -> Ptr CChar -> Ptr CInt -> Ptr CInt ->
   Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ()

newtype HBR a = HBR {forall a. HBR a -> HBR_ (RealOf a) a
getHBR :: HBR_ (RealOf a) a}

hbr :: Class.Floating a => HBR_ (RealOf a) a
hbr :: forall a. Floating a => HBR_ (RealOf a) a
hbr = HBR a
-> Order
-> Int
-> RealOf a
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a. HBR a -> HBR_ (RealOf a) a
getHBR (HBR a
 -> Order
 -> Int
 -> RealOf a
 -> Ptr CChar
 -> Ptr CInt
 -> Ptr CInt
 -> Ptr a
 -> Ptr a
 -> Ptr CInt
 -> Ptr a
 -> Ptr CInt
 -> Ptr CInt
 -> IO ())
-> HBR a
-> Order
-> Int
-> RealOf a
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a b. (a -> b) -> a -> b
$ HBR Float
-> HBR Double
-> HBR (Complex Float)
-> HBR (Complex Double)
-> HBR a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
forall (f :: * -> *).
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating (HBR_ (RealOf Float) Float -> HBR Float
forall a. HBR_ (RealOf a) a -> HBR a
HBR HBR_ Float Float
HBR_ (RealOf Float) Float
forall a. Real a => HBR_ a a
syr) (HBR_ (RealOf Double) Double -> HBR Double
forall a. HBR_ (RealOf a) a -> HBR a
HBR HBR_ Double Double
HBR_ (RealOf Double) Double
forall a. Real a => HBR_ a a
syr) (HBR_ (RealOf (Complex Float)) (Complex Float)
-> HBR (Complex Float)
forall a. HBR_ (RealOf a) a -> HBR a
HBR HBR_ Float (Complex Float)
HBR_ (RealOf (Complex Float)) (Complex Float)
forall a. Real a => HBR_ a (Complex a)
her) (HBR_ (RealOf (Complex Double)) (Complex Double)
-> HBR (Complex Double)
forall a. HBR_ (RealOf a) a -> HBR a
HBR HBR_ Double (Complex Double)
HBR_ (RealOf (Complex Double)) (Complex Double)
forall a. Real a => HBR_ a (Complex a)
her)

syr :: (Class.Real a) => HBR_ a a
syr :: forall a. Real a => HBR_ a a
syr Order
order Int
k a
alpha Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
kPtr Ptr a
alphaPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
a0Ptr Ptr CInt
incaPtr Ptr CInt
ldaPtr =
   case Order
order of
      Order
ColumnMajor -> do
         let aPtr :: Ptr a
aPtr = Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
a0Ptr Int
k
         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
-> IO ()
forall a.
Real a =>
Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> IO ()
BlasReal.syr Ptr CChar
uploPtr Ptr CInt
kPtr Ptr a
alphaPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
aPtr Ptr CInt
ldaPtr
         Ptr a -> a -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr a
alphaPtr (a -> IO ()) -> (a -> a) -> a -> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a
alphaa -> a -> a
forall a. Num a => a -> a -> a
*) (a -> IO ()) -> IO a -> IO ()
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Ptr a -> Int -> IO a
forall a. Storable a => Ptr a -> Int -> IO a
peekElemOff Ptr a
xPtr Int
k
         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
xPtr Ptr CInt
incxPtr (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aPtr (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
k)) Ptr CInt
incaPtr
      Order
RowMajor -> do
         let aPtr :: Ptr a
aPtr = Ptr a
a0Ptr
         Ptr a -> a -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr a
alphaPtr (a -> IO ()) -> (a -> a) -> a -> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a
alphaa -> a -> a
forall a. Num a => a -> a -> a
*) (a -> IO ()) -> IO a -> IO ()
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Ptr a -> IO a
forall a. Storable a => Ptr a -> IO a
peek Ptr a
xPtr
         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
xPtr Ptr CInt
incxPtr Ptr a
aPtr Ptr CInt
incaPtr
         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
-> IO ()
forall a.
Real a =>
Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> IO ()
BlasReal.syr Ptr CChar
uploPtr Ptr CInt
kPtr Ptr a
alphaPtr
            (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
xPtr Int
1) Ptr CInt
incxPtr (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aPtr (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)) Ptr CInt
ldaPtr

her :: (Class.Real a) => HBR_ a (Complex a)
her :: forall a. Real a => HBR_ a (Complex a)
her Order
order Int
k a
alpha Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
kPtr Ptr (Complex a)
alphaPtr Ptr (Complex a)
xPtr Ptr CInt
incxPtr Ptr (Complex a)
a0Ptr Ptr CInt
incaPtr Ptr CInt
ldaPtr =
   case Order
order of
      Order
ColumnMajor -> do
         let aPtr :: Ptr (Complex a)
aPtr = Ptr (Complex a) -> Int -> Ptr (Complex a)
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (Complex a)
a0Ptr Int
k
         let alphaRealPtr :: Ptr (RealOf (Complex a))
alphaRealPtr = Ptr (Complex a) -> Ptr (RealOf (Complex a))
forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
alphaPtr
         Ptr a -> a -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr a
Ptr (RealOf (Complex a))
alphaRealPtr a
alpha
         Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> IO ()
BlasComplex.her Ptr CChar
uploPtr Ptr CInt
kPtr Ptr a
Ptr (RealOf (Complex a))
alphaRealPtr Ptr (Complex a)
xPtr Ptr CInt
incxPtr Ptr (Complex a)
aPtr Ptr CInt
ldaPtr
         Ptr (Complex a) -> Complex a -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr (Complex a)
alphaPtr (Complex a -> IO ())
-> (Complex a -> Complex a) -> Complex a -> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a) -> Complex a -> Complex a
forall a b. (a -> b) -> Complex a -> Complex b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a
alphaa -> a -> a
forall a. Num a => a -> a -> a
*) (Complex a -> Complex a)
-> (Complex a -> Complex a) -> Complex a -> Complex a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Complex a -> Complex a
forall a. Num a => Complex a -> Complex a
conjugate (Complex a -> IO ()) -> IO (Complex a) -> IO ()
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Ptr (Complex a) -> Int -> IO (Complex a)
forall a. Storable a => Ptr a -> Int -> IO a
peekElemOff Ptr (Complex a)
xPtr Int
k
         Ptr CInt
-> Ptr (Complex a)
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex 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 (Complex a)
alphaPtr Ptr (Complex a)
xPtr Ptr CInt
incxPtr (Ptr (Complex a) -> Int -> Ptr (Complex a)
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (Complex a)
aPtr (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
k)) Ptr CInt
incaPtr
      Order
RowMajor -> do
         let aPtr :: Ptr (Complex a)
aPtr = Ptr (Complex a)
a0Ptr
         let alphaRealPtr :: Ptr (RealOf (Complex a))
alphaRealPtr = Ptr (Complex a) -> Ptr (RealOf (Complex a))
forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
alphaPtr
         Ptr (Complex a) -> Complex a -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr (Complex a)
alphaPtr (Complex a -> IO ())
-> (Complex a -> Complex a) -> Complex a -> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a) -> Complex a -> Complex a
forall a b. (a -> b) -> Complex a -> Complex b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a
alphaa -> a -> a
forall a. Num a => a -> a -> a
*) (Complex a -> Complex a)
-> (Complex a -> Complex a) -> Complex a -> Complex a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Complex a -> Complex a
forall a. Num a => Complex a -> Complex a
conjugate (Complex a -> IO ()) -> IO (Complex a) -> IO ()
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Ptr (Complex a) -> IO (Complex a)
forall a. Storable a => Ptr a -> IO a
peek Ptr (Complex a)
xPtr
         Ptr CInt
-> Ptr (Complex a)
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex 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 (Complex a)
alphaPtr Ptr (Complex a)
xPtr Ptr CInt
incxPtr Ptr (Complex a)
aPtr Ptr CInt
incaPtr
         Ptr a -> a -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr a
Ptr (RealOf (Complex a))
alphaRealPtr a
alpha
         Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> IO ()
BlasComplex.her Ptr CChar
uploPtr Ptr CInt
kPtr Ptr a
Ptr (RealOf (Complex a))
alphaRealPtr
            (Ptr (Complex a) -> Int -> Ptr (Complex a)
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (Complex a)
xPtr Int
1) Ptr CInt
incxPtr (Ptr (Complex a) -> Int -> Ptr (Complex a)
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (Complex a)
aPtr (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)) Ptr CInt
ldaPtr