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

import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
import qualified Numeric.LAPACK.Matrix.Extent as Extent
import qualified Numeric.LAPACK.Matrix.Banded.Basic as Banded
import qualified Numeric.LAPACK.Matrix.Triangular.Private as TriangularPriv
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 Numeric.LAPACK.ShapeStatic as ShapeStatic
import Numeric.LAPACK.Matrix.Hermitian.Private (TakeDiagonal(..))
import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian)
import Numeric.LAPACK.Matrix.Shape.Private
         (Order(RowMajor,ColumnMajor), flipOrder, uploFromOrder,
          UnaryProxy, natFromProxy)
import Numeric.LAPACK.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,
          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 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 (MatrixShape.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 :: UnaryProxy offDiag
-> Order -> size -> [a] -> BandedHermitian offDiag size a
fromList UnaryProxy offDiag
numOff Order
order size
size =
   BandedHermitian offDiag size
-> [a] -> 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
MatrixShape.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 :: 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 (() Shape.:+: ShapeStatic.ZeroBased offDiag) a ->
   BandedHermitian offDiag size a
fromRowMajor :: Matrix size (() :+: ZeroBased offDiag) a
-> BandedHermitian offDiag size a
fromRowMajor =
   ((size, () :+: ZeroBased offDiag) -> BandedHermitian offDiag size)
-> Matrix size (() :+: ZeroBased offDiag) a
-> BandedHermitian offDiag size a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape
      (\(size
size, () Shape.:+: ShapeStatic.ZeroBased UnaryProxy offDiag
k) ->
         UnaryProxy offDiag -> Order -> size -> BandedHermitian offDiag size
forall off size.
UnaryProxy off -> Order -> size -> BandedHermitian off size
MatrixShape.BandedHermitian UnaryProxy offDiag
k Order
RowMajor size
size)

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

identity ::
   (Shape.C sh, Class.Floating a) => sh -> Diagonal sh a
identity :: sh -> Diagonal sh a
identity =
   (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
MatrixShape.BandedHermitian UnaryProxy U0
forall a. Proxy a
Proxy Order
ColumnMajor) (Array sh a -> Diagonal sh a)
-> (sh -> Array sh a) -> sh -> Diagonal 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 :: 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
MatrixShape.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 :: BandedHermitian offDiag size a -> Vector size (RealOf a)
takeDiagonal =
   TakeDiagonal (Array (BandedHermitian offDiag size)) (Array size) a
-> BandedHermitian offDiag size a -> Vector 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
 -> BandedHermitian offDiag size a -> Vector size (RealOf a))
-> TakeDiagonal
     (Array (BandedHermitian offDiag size)) (Array size) a
-> BandedHermitian offDiag size a
-> Vector 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
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
-> 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
-> 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)
-> 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)
-> 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 :: BandedHermitian offDiag size a -> Vector size ar
takeDiagonalAux (Array (MatrixShape.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 ()) -> Vector 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 ()) -> Vector size ar)
-> (Int -> Ptr ar -> IO ()) -> Vector 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 (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 :: BandedHermitian offDiag size a -> Hermitian size a
toHermitian (Array (MatrixShape.BandedHermitian UnaryProxy offDiag
numOff Order
order size
size) ForeignPtr a
a) =
   Hermitian size -> (Int -> Ptr a -> IO ()) -> 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
MatrixShape.Hermitian Order
order size
size) ((Int -> Ptr a -> IO ()) -> Hermitian size a)
-> (Int -> Ptr a -> IO ()) -> 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 ()
TriangularPriv.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 :: BandedHermitian offDiag size a -> Square offDiag offDiag size a
toBanded (Array (MatrixShape.BandedHermitian UnaryProxy offDiag
numOff Order
order size
sh) ForeignPtr a
a) =
   Banded offDiag offDiag Small Small size size
-> (Ptr a -> IO ()) -> Square 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 Small Small size size
-> Banded offDiag offDiag Small Small size size
forall sub super vert horiz height width.
(UnaryProxy sub, UnaryProxy super)
-> Order
-> Extent vert horiz height width
-> Banded sub super vert horiz height width
MatrixShape.Banded (UnaryProxy offDiag
numOff,UnaryProxy offDiag
numOff) Order
order (size -> Extent Small Small size size
forall sh. sh -> Square sh
Extent.square size
sh)) ((Ptr a -> IO ()) -> Square offDiag offDiag size a)
-> (Ptr a -> IO ()) -> Square 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 ->
      case Order
order of
         Order
ColumnMajor -> UnaryProxy offDiag -> size -> Ptr a -> Ptr a -> IO ()
forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
UnaryProxy offDiag -> size -> Ptr a -> Ptr a -> IO ()
toBandedColumnMajor UnaryProxy offDiag
numOff size
sh Ptr a
aPtr Ptr a
bPtr
         Order
RowMajor -> UnaryProxy offDiag -> size -> Ptr a -> Ptr a -> IO ()
forall offDiag size a.
(Natural offDiag, C size, Floating a) =>
UnaryProxy offDiag -> size -> Ptr a -> Ptr a -> IO ()
toBandedRowMajor UnaryProxy offDiag
numOff size
sh Ptr a
aPtr Ptr a
bPtr

toBandedColumnMajor ::
   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>
   UnaryProxy offDiag -> size -> Ptr a -> Ptr a -> IO ()
toBandedColumnMajor :: UnaryProxy offDiag -> size -> Ptr a -> Ptr a -> IO ()
toBandedColumnMajor UnaryProxy offDiag
numOff size
size Ptr a
aPtr Ptr a
bPtr = do
   let n :: Int
n = size -> Int
forall sh. C sh => sh -> Int
Shape.size size
size
   let k :: Int
k = UnaryProxy offDiag -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy offDiag
numOff
   let lda0 :: Int
lda0 = Int
k
   let lda :: Int
lda = Int
lda0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
   let ldb0 :: Int
ldb0 = Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
k
   let ldb :: Int
ldb = Int
ldb0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
   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
   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
lda0
      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 (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 top :: Int
top = Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
         let bottom :: Int
bottom = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
n (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
kInt -> 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
1)Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
lda0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
topInt -> 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
ldb0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k)
         Ptr CInt -> Int -> IO ()
pokeCInt Ptr CInt
nPtr (Int
bottomInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
top)
         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
xPtr Ptr CInt
incxPtr (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
yPtr Int
top) Ptr CInt
incyPtr
         Ptr CInt -> Int -> IO ()
pokeCInt Ptr CInt
nPtr (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
bottom)
         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
bottom) Ptr CInt
incyPtr

toBandedRowMajor ::
   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>
   UnaryProxy offDiag -> size -> Ptr a -> Ptr a -> IO ()
toBandedRowMajor :: UnaryProxy offDiag -> size -> Ptr a -> Ptr a -> IO ()
toBandedRowMajor UnaryProxy offDiag
numOff size
size Ptr a
aPtr Ptr a
bPtr = do
   let n :: Int
n = size -> Int
forall sh. C sh => sh -> Int
Shape.size size
size
   let k :: Int
k = UnaryProxy offDiag -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy offDiag
numOff
   let lda0 :: Int
lda0 = Int
k
   let lda :: Int
lda = Int
lda0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
   let ldb0 :: Int
ldb0 = Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
k
   let ldb :: Int
ldb = Int
ldb0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
   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)
   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
lda0
      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 (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 left :: Int
left = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k)
         let xPtr :: Ptr a
xPtr = Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aPtr (Int
leftInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
lda0Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
i)
         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
ldb0)
         Ptr CInt -> Int -> IO ()
pokeCInt Ptr CInt
nPtr (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
left)
         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
i) Ptr CInt
incyPtr
         Ptr CInt -> Int -> IO ()
pokeCInt Ptr CInt
nPtr (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
left)
         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
xPtr Ptr CInt
incxPtr (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
yPtr (Int
leftInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k)) Ptr CInt
incyPtr


multiplyVector ::
   (Unary.Natural offDiag, Shape.C size, Eq size, Class.Floating a) =>
   Transposition -> BandedHermitian offDiag size a ->
   Vector size a -> Vector size a
multiplyVector :: Transposition
-> BandedHermitian offDiag size a -> Vector size a -> Vector size a
multiplyVector Transposition
transposed
   (Array (MatrixShape.BandedHermitian UnaryProxy offDiag
numOff Order
order size
size) ForeignPtr a
a) (Array size
sizeX ForeignPtr a
x) =
      size -> (Int -> Ptr a -> IO ()) -> Vector 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 ()) -> Vector size a)
-> (Int -> Ptr a -> IO ()) -> Vector 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 (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 :: 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 Small Small size size
-> (UnaryProxy sub, UnaryProxy super)
forall sub super vert horiz height width.
Banded sub super vert horiz height width
-> (UnaryProxy sub, UnaryProxy super)
MatrixShape.bandedOffDiagonals (Banded sub super Small Small size size
 -> (UnaryProxy sub, UnaryProxy super))
-> Banded sub super Small Small size size
-> (UnaryProxy sub, UnaryProxy super)
forall a b. (a -> b) -> a -> b
$ Square sub super size a -> Banded sub super 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 Small Small size size a
-> Square sub super size a
-> Square (sub :+: super) (sub :+: super) size a
forall subA superA subB superB subC superC vert horiz height width
       fuse a.
(Natural subA, Natural superA, Natural subB, Natural superB,
 (subA :+: subB) ~ subC, (superA :+: superB) ~ superC, C vert,
 C horiz, C height, C width, C fuse, Eq fuse, Floating a) =>
Banded subA superA vert horiz height fuse a
-> Banded subB superB vert horiz fuse width a
-> Banded subC superC vert horiz height width a
Banded.multiply (Square sub super size a -> Banded super sub Small Small size size a
forall sub super vert horiz height width a.
(Natural sub, Natural super, C vert, C horiz, C height, C width,
 Floating a) =>
Banded sub super vert horiz height width a
-> Banded super sub 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 :: Square offDiag offDiag size a -> BandedHermitian offDiag size a
fromUpperPart (Array (MatrixShape.Banded (UnaryProxy offDiag
sub,UnaryProxy offDiag
super) Order
order Extent Small Small size size
extent) ForeignPtr a
a) =
   let sh :: size
sh = Extent Small Small size size -> size
forall height width. Extent Small Small height width -> height
Extent.squareSize Extent 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 ()) -> 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
MatrixShape.BandedHermitian UnaryProxy offDiag
super Order
order size
sh) ((Ptr a -> IO ()) -> BandedHermitian offDiag size a)
-> (Ptr a -> IO ()) -> 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.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Transposition -> BandedHermitian offDiag height a ->
   Matrix.Full vert horiz height width a ->
   Matrix.Full vert horiz height width a
multiplyFull :: Transposition
-> BandedHermitian offDiag height a
-> Full vert horiz height width a
-> Full vert horiz height width a
multiplyFull Transposition
transposed BandedHermitian offDiag height a
a Full vert horiz height width a
b =
   case Full vert horiz height width -> Order
forall vert horiz height width.
Full vert horiz height width -> Order
MatrixShape.fullOrder (Full vert horiz height width -> Order)
-> Full vert horiz height width -> Order
forall a b. (a -> b) -> a -> b
$ Full vert horiz height width a -> Full vert horiz height width
forall sh a. Array sh a -> sh
Array.shape Full vert horiz height width a
b of
      Order
ColumnMajor -> Transposition
-> BandedHermitian offDiag height a
-> Full vert horiz height width a
-> Full vert horiz height width a
forall offDiag vert horiz height width a.
(Natural offDiag, C vert, C horiz, Eq height, C height, C width,
 Floating a) =>
Transposition
-> BandedHermitian offDiag height a
-> Full vert horiz height width a
-> Full vert horiz height width a
multiplyFullSpecial Transposition
transposed BandedHermitian offDiag height a
a Full vert horiz height width a
b
      Order
RowMajor -> Transposition
-> BandedHermitian offDiag height a
-> Full vert horiz height width a
-> Full vert horiz height width a
forall offDiag vert horiz height width a.
(Natural offDiag, C vert, C horiz, C height, Eq height, C width,
 Floating a) =>
Transposition
-> BandedHermitian offDiag height a
-> Full vert horiz height width a
-> Full vert horiz height width a
multiplyFullGeneric Transposition
transposed BandedHermitian offDiag height a
a Full vert horiz height width a
b

multiplyFullSpecial ::
   (Unary.Natural offDiag, Extent.C vert, Extent.C horiz,
    Eq height, Shape.C height, Shape.C width, Class.Floating a) =>
   Transposition -> BandedHermitian offDiag height a ->
   Matrix.Full vert horiz height width a ->
   Matrix.Full vert horiz height width a
multiplyFullSpecial :: Transposition
-> BandedHermitian offDiag height a
-> Full vert horiz height width a
-> Full vert horiz height width a
multiplyFullSpecial Transposition
transposed
      (Array (MatrixShape.BandedHermitian UnaryProxy offDiag
numOff Order
orderA height
sizeA) ForeignPtr a
a)
      (Array (MatrixShape.Full Order
orderB Extent vert horiz height width
extentB) ForeignPtr a
b) =
   Full vert horiz height width
-> (Ptr a -> IO ()) -> Full vert horiz height width a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order
-> Extent vert horiz height width -> Full vert horiz height width
forall vert horiz height width.
Order
-> Extent vert horiz height width -> Full vert horiz height width
MatrixShape.Full Order
orderB Extent vert horiz height width
extentB) ((Ptr a -> IO ()) -> Full vert horiz height width a)
-> (Ptr a -> IO ()) -> Full vert horiz height width a
forall a b. (a -> b) -> a -> b
$ \Ptr a
cPtr -> do
      String -> Bool -> IO ()
Call.assert String
"BandedHermitian.multiplyFull: shapes mismatch"
         (height
sizeA height -> height -> Bool
forall a. Eq a => a -> a -> Bool
== Extent vert horiz height width -> height
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> height
Extent.height Extent vert horiz height width
extentB)
      let (height
height,width
width) = Extent vert horiz height width -> (height, width)
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> (height, width)
Extent.dimensions Extent vert horiz height width
extentB
      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 :: 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 (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 (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 :: 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.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Transposition -> BandedHermitian offDiag height a ->
   Matrix.Full vert horiz height width a ->
   Matrix.Full vert horiz height width a
multiplyFullGeneric :: Transposition
-> BandedHermitian offDiag height a
-> Full vert horiz height width a
-> Full vert horiz height width a
multiplyFullGeneric Transposition
transposed BandedHermitian offDiag height a
a Full 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 vert horiz sub super height width a.
(C vert, C horiz) =>
Banded sub super vert horiz height width a
-> Banded super sub horiz vert width height a
Banded.transpose Square U0 offDiag height a
upper, Square offDiag U0 height a -> Square U0 offDiag height a
forall vert horiz sub super height width a.
(C vert, C horiz) =>
Banded sub super vert horiz height width a
-> Banded super sub 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 a
-> Full vert horiz height width a
-> Full vert horiz height width a
-> Full vert horiz height width a
forall sh a.
(C sh, Floating a) =>
a -> Array sh a -> Array sh a -> Array sh a
VectorPriv.mac a
forall a. Floating a => a
one
         (Banded offDiag U0 vert horiz height height a
-> Full vert horiz height width a -> Full vert horiz height width a
forall sub super vert horiz height width fuse a.
(Natural sub, Natural super, C vert, C horiz, C height, C width,
 C fuse, Eq fuse, Floating a) =>
Banded sub super vert horiz height fuse a
-> Full vert horiz fuse width a -> Full vert horiz height width a
Banded.multiplyFull (Map Small Small vert horiz height height
-> Square offDiag U0 height a
-> Banded offDiag U0 vert horiz height height a
forall vertA horizA vertB horizB height width super sub a.
(C vertA, C horizA, C vertB, C horizB) =>
Map vertA horizA vertB horizB height width
-> Banded super sub vertA horizA height width a
-> Banded super sub vertB horizB height width a
Banded.mapExtent Map Small Small vert horiz height height
forall vert horiz size.
(C vert, C horiz) =>
Map Small Small vert horiz size size
Extent.fromSquare Square offDiag U0 height a
lowerT) Full vert horiz height width a
b)
         (Banded U0 offDiag vert horiz height height a
-> Full vert horiz height width a -> Full vert horiz height width a
forall sub super vert horiz height width fuse a.
(Natural sub, Natural super, C vert, C horiz, C height, C width,
 C fuse, Eq fuse, Floating a) =>
Banded sub super vert horiz height fuse a
-> Full vert horiz fuse width a -> Full vert horiz height width a
Banded.multiplyFull (Map Small Small vert horiz height height
-> Square U0 offDiag height a
-> Banded U0 offDiag vert horiz height height a
forall vertA horizA vertB horizB height width super sub a.
(C vertA, C horizA, C vertB, C horizB) =>
Map vertA horizA vertB horizB height width
-> Banded super sub vertA horizA height width a
-> Banded super sub vertB horizB height width a
Banded.mapExtent Map Small Small vert horiz height height
forall vert horiz size.
(C vert, C horiz) =>
Map Small Small vert horiz size size
Extent.fromSquare Square U0 offDiag height a
upperT) Full 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 :: BandedHermitian offDiag size a -> Square U0 offDiag size a
takeUpper =
   (BandedHermitian offDiag size
 -> Banded U0 offDiag Small Small size size)
-> BandedHermitian offDiag size a -> Square U0 offDiag size a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape
      (\(MatrixShape.BandedHermitian UnaryProxy offDiag
numOff Order
order size
sh) ->
         (UnaryProxy U0, UnaryProxy offDiag)
-> Order -> size -> Banded U0 offDiag Small Small size size
forall sub super size.
(UnaryProxy sub, UnaryProxy super)
-> Order -> size -> Banded sub super Small Small size size
MatrixShape.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 :: BandedHermitian offDiag size a -> Square offDiag U0 size a
takeStrictLower (Array (MatrixShape.BandedHermitian UnaryProxy offDiag
numOff Order
order size
sh) ForeignPtr a
x) =
   Banded offDiag U0 Small Small size size
-> (Int -> Ptr a -> IO ()) -> Square 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 -> Banded offDiag U0 Small Small size size
forall sub super size.
(UnaryProxy sub, UnaryProxy super)
-> Order -> size -> Banded sub super Small Small size size
MatrixShape.bandedSquare (UnaryProxy offDiag
numOff,UnaryProxy U0
forall a. Proxy a
Proxy) (Order -> Order
flipOrder Order
order) size
sh) ((Int -> Ptr a -> IO ()) -> Square offDiag U0 size a)
-> (Int -> Ptr a -> IO ()) -> Square 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 (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 :: 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, StaticVector (Succ k) a))]
-> 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, StaticVector (Succ k) a))]
 -> BandedHermitian k sh a)
-> SumRank1 k sh a
-> Order
-> sh
-> [(RealOf a, (Index sh, StaticVector (Succ k) a))]
-> 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
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 {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 :: 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
MatrixShape.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 (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 {HBR a -> HBR_ (RealOf a) a
getHBR :: HBR_ (RealOf a) a}

hbr :: Class.Floating a => HBR_ (RealOf a) a
hbr :: HBR_ (RealOf a) a
hbr = HBR a -> HBR_ (RealOf a) a
forall a. HBR a -> HBR_ (RealOf a) a
getHBR (HBR a -> HBR_ (RealOf a) a) -> HBR a -> HBR_ (RealOf a) a
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
Class.switchFloating (HBR_ (RealOf Float) Float -> HBR Float
forall a. HBR_ (RealOf a) a -> HBR a
HBR 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_ (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_ (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_ (RealOf (Complex Double)) (Complex Double)
forall a. Real a => HBR_ a (Complex a)
her)

syr :: (Class.Real a) => HBR_ a a
syr :: 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 :: 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 (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 (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