{-# LANGUAGE TypeFamilies #-}
module Numeric.BLAS.Private where

import qualified Numeric.BLAS.FFI.Real as BlasReal
import qualified Numeric.BLAS.FFI.Generic as Blas
import qualified Numeric.Netlib.Utility as Call
import qualified Numeric.Netlib.Class as Class
import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated, Conjugated))
import Numeric.BLAS.Scalar (RealOf, zero, one, minusOne, isZero)

import qualified Foreign.Marshal.Array.Guarded as ForeignArray
import qualified Foreign.Marshal.Utils as Marshal
import qualified Foreign.C.String as CStr
import Foreign.Marshal.Array (advancePtr)
import Foreign.C.Types (CChar, CInt)
import Foreign.Ptr (Ptr, castPtr)
import Foreign.Storable (Storable, peek, pokeElemOff, peekElemOff)

import Control.Monad.Trans.Cont (evalContT)
import Control.Monad.IO.Class (liftIO)
import Control.Monad (when)
import Control.Applicative (liftA2)

import qualified Data.Array.Comfort.Shape as Shape

import qualified Data.Complex as Complex
import Data.Complex (Complex((:+)))

import Prelude hiding (sum)


type ShapeInt = Shape.ZeroBased Int

shapeInt :: Int -> ShapeInt
shapeInt :: Int -> ShapeInt
shapeInt = forall n. n -> ZeroBased n
Shape.ZeroBased


realPtr :: Ptr a -> Ptr (RealOf a)
realPtr :: forall a. Ptr a -> Ptr (RealOf a)
realPtr = forall a b. Ptr a -> Ptr b
castPtr


pointerSeq :: (Storable a) => Int -> Ptr a -> [Ptr a]
pointerSeq :: forall a. Storable a => Int -> Ptr a -> [Ptr a]
pointerSeq Int
k Ptr a
ptr = forall a. (a -> a) -> a -> [a]
iterate (forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Int
k) Ptr a
ptr


fill :: (Class.Floating a) => a -> Int -> Ptr a -> IO ()
fill :: forall a. Floating a => a -> Int -> Ptr a -> IO ()
fill a
a Int
n Ptr a
dstPtr = forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
   Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
   Ptr a
srcPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
a
   Ptr CInt
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
   Ptr CInt
incyPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
   forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy Ptr CInt
nPtr Ptr a
srcPtr Ptr CInt
incxPtr Ptr a
dstPtr Ptr CInt
incyPtr


copyConjugate ::
   (Class.Floating a) =>
   Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
copyConjugate :: 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
yPtr Ptr CInt
incyPtr = do
   forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy Ptr CInt
nPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr
   forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
lacgv Ptr CInt
nPtr Ptr a
yPtr Ptr CInt
incyPtr



newtype Sum a = Sum {forall a. Sum a -> Int -> Ptr a -> Int -> IO a
runSum :: Int -> Ptr a -> Int -> IO a}

sum :: Class.Floating a => Int -> Ptr a -> Int -> IO a
sum :: forall a. Floating a => Int -> Ptr a -> Int -> IO a
sum =
   forall a. Sum a -> Int -> Ptr a -> Int -> IO a
runSum forall a b. (a -> b) -> a -> b
$
   forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (forall a. (Int -> Ptr a -> Int -> IO a) -> Sum a
Sum forall a. Real a => Int -> Ptr a -> Int -> IO a
sumReal)
      (forall a. (Int -> Ptr a -> Int -> IO a) -> Sum a
Sum forall a. Real a => Int -> Ptr a -> Int -> IO a
sumReal)
      (forall a. (Int -> Ptr a -> Int -> IO a) -> Sum a
Sum forall a. Real a => Int -> Ptr (Complex a) -> Int -> IO (Complex a)
sumComplex)
      (forall a. (Int -> Ptr a -> Int -> IO a) -> Sum a
Sum forall a. Real a => Int -> Ptr (Complex a) -> Int -> IO (Complex a)
sumComplex)

sumReal :: Class.Real a => Int -> Ptr a -> Int -> IO a
sumReal :: forall a. Real a => Int -> Ptr a -> Int -> IO a
sumReal Int
n Ptr a
xPtr Int
incx =
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
incx
      Ptr a
yPtr <- forall a r. Real a => a -> FortranIO r (Ptr a)
Call.real forall a. Floating a => a
one
      Ptr CInt
incyPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.dot Ptr CInt
nPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr

sumComplex, sumComplexAlt ::
   Class.Real a => Int -> Ptr (Complex a) -> Int -> IO (Complex a)
sumComplex :: forall a. Real a => Int -> Ptr (Complex a) -> Int -> IO (Complex a)
sumComplex Int
n Ptr (Complex a)
xPtr Int
incx =
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      let sxPtr :: Ptr (RealOf (Complex a))
sxPtr = forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
xPtr
      Ptr CInt
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (Int
2forall a. Num a => a -> a -> a
*Int
incx)
      Ptr a
yPtr <- forall a r. Real a => a -> FortranIO r (Ptr a)
Call.real forall a. Floating a => a
one
      Ptr CInt
incyPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$
         forall (f :: * -> *) a b c.
Applicative f =>
(a -> b -> c) -> f a -> f b -> f c
liftA2 forall a. a -> a -> Complex a
(Complex.:+)
            (forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.dot Ptr CInt
nPtr Ptr (RealOf (Complex a))
sxPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr)
            (forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.dot Ptr CInt
nPtr (forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (RealOf (Complex a))
sxPtr Int
1) Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr)

sumComplexAlt :: forall a. Real a => Int -> Ptr (Complex a) -> Int -> IO (Complex a)
sumComplexAlt Int
n Ptr (Complex a)
aPtr Int
inca =
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr CChar
transPtr <- forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'N'
      Ptr CInt
mPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
2
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
onePtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number forall a. Floating a => a
one
      Ptr CInt
inc0Ptr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
      let saPtr :: Ptr (RealOf (Complex a))
saPtr = forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
aPtr
      Ptr CInt
ldaPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim (Int
2forall a. Num a => a -> a -> a
*Int
inca)
      Ptr a
sxPtr <- forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      Ptr CInt
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      Ptr a
betaPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number forall a. Floating a => a
zero
      Ptr (Complex a)
yPtr <- forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      let syPtr :: Ptr (RealOf (Complex a))
syPtr = forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
yPtr
      Ptr CInt
incyPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ do
         forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy Ptr CInt
nPtr Ptr a
onePtr Ptr CInt
inc0Ptr Ptr a
sxPtr Ptr CInt
incxPtr
         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 ()
gemv
            Ptr CChar
transPtr Ptr CInt
mPtr Ptr CInt
nPtr Ptr a
onePtr Ptr (RealOf (Complex a))
saPtr Ptr CInt
ldaPtr
            Ptr a
sxPtr Ptr CInt
incxPtr Ptr a
betaPtr Ptr (RealOf (Complex a))
syPtr Ptr CInt
incyPtr
         forall a. Storable a => Ptr a -> IO a
peek Ptr (Complex a)
yPtr


mul ::
   (Class.Floating a) =>
   Conjugation -> Int -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
mul :: forall a.
Floating a =>
Conjugation
-> Int -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
mul Conjugation
conj Int
n Ptr a
aPtr Int
inca Ptr a
xPtr Int
incx Ptr a
yPtr Int
incy =
   forall a.
Floating a =>
Conjugation
-> Int
-> Ptr a
-> Int
-> Ptr a
-> Int
-> a
-> Ptr a
-> Int
-> IO ()
mulAdd Conjugation
conj Int
n Ptr a
aPtr Int
inca Ptr a
xPtr Int
incx forall a. Floating a => a
zero Ptr a
yPtr Int
incy

mulAdd ::
   (Class.Floating a) =>
   Conjugation ->
   Int -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()
mulAdd :: forall a.
Floating a =>
Conjugation
-> Int
-> Ptr a
-> Int
-> Ptr a
-> Int
-> a
-> Ptr a
-> Int
-> IO ()
mulAdd Conjugation
conj Int
n Ptr a
aPtr Int
inca Ptr a
xPtr Int
incx a
beta Ptr a
yPtr Int
incy = forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
   Ptr CChar
transPtr <- forall r. Char -> FortranIO r (Ptr CChar)
Call.char forall a b. (a -> b) -> a -> b
$ case Conjugation
conj of Conjugation
NonConjugated -> Char
'N'; Conjugation
Conjugated -> Char
'C'
   Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
   Ptr CInt
klPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
   Ptr CInt
kuPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
   Ptr a
alphaPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number forall a. Floating a => a
one
   Ptr CInt
ldaPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
inca
   Ptr CInt
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
incx
   Ptr a
betaPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
beta
   Ptr CInt
incyPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
incy
   forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$
      forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
Blas.gbmv Ptr CChar
transPtr
         Ptr CInt
nPtr Ptr CInt
nPtr Ptr CInt
klPtr Ptr CInt
kuPtr Ptr a
alphaPtr Ptr a
aPtr Ptr CInt
ldaPtr
         Ptr a
xPtr Ptr CInt
incxPtr Ptr a
betaPtr Ptr a
yPtr Ptr CInt
incyPtr

{- |
Use the foldBalanced trick.
-}
product :: (Class.Floating a) => Int -> Ptr a -> Int -> IO a
product :: forall a. Floating a => Int -> Ptr a -> Int -> IO a
product Int
n Ptr a
aPtr Int
inca =
   case forall a. Ord a => a -> a -> Ordering
compare Int
n Int
1 of
      Ordering
LT -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a. Floating a => a
one
      Ordering
EQ -> forall a. Storable a => Ptr a -> IO a
peek Ptr a
aPtr
      Ordering
GT -> let n2 :: Int
n2 = forall a. Integral a => a -> a -> a
div Int
n Int
2; new :: Int
new = Int
nforall a. Num a => a -> a -> a
-Int
n2
            in forall a b. Storable a => Int -> (Ptr a -> IO b) -> IO b
ForeignArray.alloca (Int
2forall a. Num a => a -> a -> a
*Int
newforall a. Num a => a -> a -> a
-Int
1) forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> do
         forall a.
Floating a =>
Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
mulPairs Int
n2 Ptr a
aPtr Int
inca Ptr a
xPtr Int
1
         forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (forall a. Integral a => a -> Bool
odd Int
n) forall a b. (a -> b) -> a -> b
$ forall a. Storable a => Ptr a -> Int -> a -> IO ()
pokeElemOff Ptr a
xPtr Int
n2 forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< forall a. Storable a => Ptr a -> Int -> IO a
peekElemOff Ptr a
aPtr ((Int
nforall a. Num a => a -> a -> a
-Int
1)forall a. Num a => a -> a -> a
*Int
inca)
         forall a. Floating a => Int -> Ptr a -> IO a
productLoop Int
new Ptr a
xPtr

{- |
If 'mul' would be based on a scalar loop
we would not need to cut the vector into chunks.

The invariance is:
When calling @productLoop n xPtr@,
starting from xPtr there is storage allocated for 2*n-1 elements.
-}
productLoop :: (Class.Floating a) => Int -> Ptr a -> IO a
productLoop :: forall a. Floating a => Int -> Ptr a -> IO a
productLoop Int
n Ptr a
xPtr =
   if Int
nforall a. Eq a => a -> a -> Bool
==Int
1
      then forall a. Storable a => Ptr a -> IO a
peek Ptr a
xPtr
      else do
         let n2 :: Int
n2 = forall a. Integral a => a -> a -> a
div Int
n Int
2
         forall a.
Floating a =>
Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
mulPairs Int
n2 Ptr a
xPtr Int
1 (forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
xPtr Int
n) Int
1
         forall a. Floating a => Int -> Ptr a -> IO a
productLoop (Int
nforall a. Num a => a -> a -> a
-Int
n2) (forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
xPtr (Int
2forall a. Num a => a -> a -> a
*Int
n2))

mulPairs ::
   (Class.Floating a) =>
   Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
mulPairs :: forall a.
Floating a =>
Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
mulPairs Int
n Ptr a
aPtr Int
inca Ptr a
xPtr Int
incx =
   let inca2 :: Int
inca2 = Int
2forall a. Num a => a -> a -> a
*Int
inca
   in forall a.
Floating a =>
Conjugation
-> Int -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
mul Conjugation
NonConjugated Int
n Ptr a
aPtr Int
inca2 (forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aPtr Int
inca) Int
inca2 Ptr a
xPtr Int
incx


newtype LACGV a = LACGV {forall a. LACGV a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
getLACGV :: Ptr CInt -> Ptr a -> Ptr CInt -> IO ()}

lacgv :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
lacgv :: forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
lacgv =
   forall a. LACGV a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
getLACGV forall a b. (a -> b) -> a -> b
$
   forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (forall a. (Ptr CInt -> Ptr a -> Ptr CInt -> IO ()) -> LACGV a
LACGV forall a b. (a -> b) -> a -> b
$ forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *) a. Monad m => a -> m a
return ())
      (forall a. (Ptr CInt -> Ptr a -> Ptr CInt -> IO ()) -> LACGV a
LACGV forall a b. (a -> b) -> a -> b
$ forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *) a. Monad m => a -> m a
return ())
      (forall a. (Ptr CInt -> Ptr a -> Ptr CInt -> IO ()) -> LACGV a
LACGV forall a.
Real a =>
Ptr CInt -> Ptr (Complex a) -> Ptr CInt -> IO ()
clacgv)
      (forall a. (Ptr CInt -> Ptr a -> Ptr CInt -> IO ()) -> LACGV a
LACGV forall a.
Real a =>
Ptr CInt -> Ptr (Complex a) -> Ptr CInt -> IO ()
clacgv)

clacgv :: Class.Real a => Ptr CInt -> Ptr (Complex a) -> Ptr CInt -> IO ()
clacgv :: forall a.
Real a =>
Ptr CInt -> Ptr (Complex a) -> Ptr CInt -> IO ()
clacgv Ptr CInt
nPtr Ptr (Complex a)
xPtr Ptr CInt
incxPtr =
   forall a b. Storable a => a -> (Ptr a -> IO b) -> IO b
Marshal.with forall a. Floating a => a
minusOne forall a b. (a -> b) -> a -> b
$ \Ptr a
saPtr -> do
      CInt
incx <- forall a. Storable a => Ptr a -> IO a
peek Ptr CInt
incxPtr
      forall a b. Storable a => a -> (Ptr a -> IO b) -> IO b
Marshal.with (CInt
2forall a. Num a => a -> a -> a
*CInt
incx) forall a b. (a -> b) -> a -> b
$ \Ptr CInt
incyPtr ->
         forall a. Real a => Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> IO ()
BlasReal.scal Ptr CInt
nPtr Ptr a
saPtr (forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr (forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
xPtr) Int
1) Ptr CInt
incyPtr


{-
Work around an inconsistency of BLAS.
In case of a zero-column matrix
BLAS's gemv and gbmv do not initialize the target vector.
In contrast, these work-arounds do.
-}
{-# INLINE gemv #-}
gemv ::
   (Class.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 ()
gemv :: 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 ()
gemv Ptr CChar
transPtr Ptr CInt
mPtr Ptr CInt
nPtr
      Ptr a
alphaPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
betaPtr Ptr a
yPtr Ptr CInt
incyPtr = do
   forall a.
Floating a =>
Ptr CChar
-> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> IO ()
initializeMV Ptr CChar
transPtr Ptr CInt
mPtr Ptr CInt
nPtr Ptr a
betaPtr Ptr a
yPtr Ptr CInt
incyPtr
   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 ()
Blas.gemv Ptr CChar
transPtr Ptr CInt
mPtr Ptr CInt
nPtr
      Ptr a
alphaPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
betaPtr Ptr a
yPtr Ptr CInt
incyPtr

{-# INLINE gbmv #-}
gbmv ::
   (Class.Floating a) =>
   Ptr CChar -> Ptr CInt -> Ptr CInt -> Ptr CInt -> Ptr CInt ->
   Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->
   Ptr a -> Ptr a -> Ptr CInt -> IO ()
gbmv :: forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
gbmv Ptr CChar
transPtr Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
klPtr Ptr CInt
kuPtr
      Ptr a
alphaPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
betaPtr Ptr a
yPtr Ptr CInt
incyPtr = do
   forall a.
Floating a =>
Ptr CChar
-> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> IO ()
initializeMV Ptr CChar
transPtr Ptr CInt
mPtr Ptr CInt
nPtr Ptr a
betaPtr Ptr a
yPtr Ptr CInt
incyPtr
   forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
Blas.gbmv Ptr CChar
transPtr Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
klPtr Ptr CInt
kuPtr
      Ptr a
alphaPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
betaPtr Ptr a
yPtr Ptr CInt
incyPtr

initializeMV ::
   Class.Floating a =>
   Ptr CChar -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> IO ()
initializeMV :: forall a.
Floating a =>
Ptr CChar
-> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> IO ()
initializeMV Ptr CChar
transPtr Ptr CInt
mPtr Ptr CInt
nPtr Ptr a
betaPtr Ptr a
yPtr Ptr CInt
incyPtr = do
   CChar
trans <- forall a. Storable a => Ptr a -> IO a
peek Ptr CChar
transPtr
   let (Ptr CInt
mtPtr,Ptr CInt
ntPtr) =
         if CChar
trans forall a. Eq a => a -> a -> Bool
== Char -> CChar
CStr.castCharToCChar Char
'N'
            then (Ptr CInt
mPtr,Ptr CInt
nPtr) else (Ptr CInt
nPtr,Ptr CInt
mPtr)
   CInt
n <- forall a. Storable a => Ptr a -> IO a
peek Ptr CInt
ntPtr
   a
beta <- forall a. Storable a => Ptr a -> IO a
peek Ptr a
betaPtr
   forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (CInt
n forall a. Eq a => a -> a -> Bool
== CInt
0 Bool -> Bool -> Bool
&& forall a. Floating a => a -> Bool
isZero a
beta) forall a b. (a -> b) -> a -> b
$
      forall a b. Storable a => a -> (Ptr a -> IO b) -> IO b
Marshal.with CInt
0 forall a b. (a -> b) -> a -> b
$ \Ptr CInt
incbPtr ->
      forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy Ptr CInt
mtPtr Ptr a
betaPtr Ptr CInt
incbPtr Ptr a
yPtr Ptr CInt
incyPtr


{-
ToDo:

type ComplexShape =
         Shape.NestedTuple Shape.TupleAccessor (Complex Shape.Element)

This would allow the use of Complex.realPart as accessor,
but it requires GHC>7.6.3 or so, where realPart has no RealFloat constraint.
-}
type ComplexShape = Shape.NestedTuple Shape.TupleIndex (Complex Shape.Element)

ixReal, ixImaginary :: Shape.ElementIndex (Complex Shape.Element)
ElementIndex (Complex Element)
ixReal :+ ElementIndex (Complex Element)
ixImaginary =
   forall tuple.
ElementTuple tuple =>
NestedTuple TupleIndex tuple
-> DataTuple tuple (ElementIndex tuple)
Shape.indexTupleFromShape (forall sh. Static sh => sh
Shape.static :: ComplexShape)