{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE ConstraintKinds #-}
module Numeric.LAPACK.Matrix.Triangular.Private where

import qualified Numeric.LAPACK.Matrix.Private as Matrix
import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
import qualified Numeric.LAPACK.Matrix.Shape.Box as Box
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import Numeric.LAPACK.Matrix.Shape.Private
         (Order(RowMajor,ColumnMajor), flipOrder, uploFromOrder,
          Empty, Filled, NonUnit)
import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated))
import Numeric.LAPACK.Matrix.Private (Full)
import Numeric.LAPACK.Scalar (zero)
import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))
import Numeric.LAPACK.Private
         (pointerSeq, copyBlock, copyCondConjugateToTemp,
          pokeCInt, fill, withInfo, errorCodeMsg)

import qualified Numeric.LAPACK.FFI.Generic as LapackGen
import qualified Numeric.Netlib.Utility as Call
import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Storable.Unchecked as Array
import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Storable.Unchecked (Array(Array))
import Data.Array.Comfort.Shape ((:+:)((:+:)))

import Foreign.Marshal.Alloc (alloca)
import Foreign.Marshal.Array (advancePtr)
import Foreign.C.Types (CInt)
import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
import Foreign.Ptr (Ptr)
import Foreign.Storable (Storable)

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

import Data.Foldable (forM_)


diagonalPointers :: (Storable a) => Order -> Int -> Ptr a -> [Ptr a]
diagonalPointers :: Order -> Int -> Ptr a -> [Ptr a]
diagonalPointers Order
order Int
n Ptr a
aPtr =
   Int -> [Ptr a] -> [Ptr a]
forall a. Int -> [a] -> [a]
take Int
n ([Ptr a] -> [Ptr a]) -> [Ptr a] -> [Ptr a]
forall a b. (a -> b) -> a -> b
$ (Ptr a -> Int -> Ptr a) -> Ptr a -> [Int] -> [Ptr a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl 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 -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate Int -> Int
forall a. Enum a => a -> a
pred Int
n
      Order
ColumnMajor -> (Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate Int -> Int
forall a. Enum a => a -> a
succ Int
2

diagonalPointerPairs ::
   (Storable a, Storable b) =>
   Order -> Int -> Ptr a -> Ptr b -> [(Ptr a, Ptr b)]
diagonalPointerPairs :: Order -> Int -> Ptr a -> Ptr b -> [(Ptr a, Ptr b)]
diagonalPointerPairs Order
order Int
n Ptr a
aPtr Ptr b
bPtr =
   [Ptr a] -> [Ptr b] -> [(Ptr a, Ptr b)]
forall a b. [a] -> [b] -> [(a, b)]
zip (Int -> Ptr a -> [Ptr a]
forall a. Storable a => Int -> Ptr a -> [Ptr a]
pointerSeq Int
1 Ptr a
aPtr) ([Ptr b] -> [(Ptr a, Ptr b)]) -> [Ptr b] -> [(Ptr a, Ptr b)]
forall a b. (a -> b) -> a -> b
$ Order -> Int -> Ptr b -> [Ptr b]
forall a. Storable a => Order -> Int -> Ptr a -> [Ptr a]
diagonalPointers Order
order Int
n Ptr b
bPtr


columnMajorPointers ::
   (Storable a) => Int -> Ptr a -> Ptr a -> [(Int, ((Ptr a, Ptr a), Ptr a))]
columnMajorPointers :: Int -> Ptr a -> Ptr a -> [(Int, ((Ptr a, Ptr a), Ptr a))]
columnMajorPointers Int
n Ptr a
fullPtr Ptr a
packedPtr =
   let ds :: [Int]
ds = (Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate Int -> Int
forall a. Enum a => a -> a
succ Int
1
   in  Int
-> [(Int, ((Ptr a, Ptr a), Ptr a))]
-> [(Int, ((Ptr a, Ptr a), Ptr a))]
forall a. Int -> [a] -> [a]
take Int
n ([(Int, ((Ptr a, Ptr a), Ptr a))]
 -> [(Int, ((Ptr a, Ptr a), Ptr a))])
-> [(Int, ((Ptr a, Ptr a), Ptr a))]
-> [(Int, ((Ptr a, Ptr a), Ptr a))]
forall a b. (a -> b) -> a -> b
$ [Int]
-> [((Ptr a, Ptr a), Ptr a)] -> [(Int, ((Ptr a, Ptr a), Ptr a))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
ds ([((Ptr a, Ptr a), Ptr a)] -> [(Int, ((Ptr a, Ptr a), Ptr a))])
-> [((Ptr a, Ptr a), Ptr a)] -> [(Int, ((Ptr a, Ptr a), Ptr a))]
forall a b. (a -> b) -> a -> b
$
       [(Ptr a, Ptr a)] -> [Ptr a] -> [((Ptr a, Ptr a), Ptr a)]
forall a b. [a] -> [b] -> [(a, b)]
zip
         ([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
1 Ptr a
fullPtr) (Int -> Ptr a -> [Ptr a]
forall a. Storable a => Int -> Ptr a -> [Ptr a]
pointerSeq Int
n Ptr a
fullPtr))
         ((Ptr a -> Int -> Ptr a) -> Ptr a -> [Int] -> [Ptr a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
packedPtr [Int]
ds)

rowMajorPointers ::
   (Storable a) => Int -> Ptr a -> Ptr a -> [(Int, (Ptr a, Ptr a))]
rowMajorPointers :: Int -> Ptr a -> Ptr a -> [(Int, (Ptr a, Ptr a))]
rowMajorPointers Int
n Ptr a
fullPtr Ptr a
packedPtr =
   let ds :: [Int]
ds = (Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate Int -> Int
forall a. Enum a => a -> a
pred Int
n
   in  Int -> [(Int, (Ptr a, Ptr a))] -> [(Int, (Ptr a, Ptr a))]
forall a. Int -> [a] -> [a]
take Int
n ([(Int, (Ptr a, Ptr a))] -> [(Int, (Ptr a, Ptr a))])
-> [(Int, (Ptr a, Ptr a))] -> [(Int, (Ptr a, Ptr a))]
forall a b. (a -> b) -> a -> b
$ [Int] -> [(Ptr a, Ptr a)] -> [(Int, (Ptr a, Ptr a))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
ds ([(Ptr a, Ptr a)] -> [(Int, (Ptr a, Ptr a))])
-> [(Ptr a, Ptr a)] -> [(Int, (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
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Ptr a
fullPtr) ((Ptr a -> Int -> Ptr a) -> Ptr a -> [Int] -> [Ptr a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
packedPtr [Int]
ds)


forPointers :: [(Int, a)] -> (Ptr CInt -> a -> IO ()) -> IO ()
forPointers :: [(Int, a)] -> (Ptr CInt -> a -> IO ()) -> IO ()
forPointers [(Int, a)]
xs Ptr CInt -> a -> IO ()
act =
   (Ptr CInt -> IO ()) -> IO ()
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr CInt
nPtr ->
   [(Int, a)] -> ((Int, a) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [(Int, a)]
xs (((Int, a) -> IO ()) -> IO ()) -> ((Int, a) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \(Int
d,a
ptrs) -> do
      Ptr CInt -> Int -> IO ()
pokeCInt Ptr CInt
nPtr Int
d
      Ptr CInt -> a -> IO ()
act Ptr CInt
nPtr a
ptrs


copyTriangleToTemp ::
   Class.Floating a =>
   Conjugation -> Order -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyTriangleToTemp :: Conjugation -> Order -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyTriangleToTemp Conjugation
conj Order
order =
   Conjugation -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
forall a r.
Floating a =>
Conjugation -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyCondConjugateToTemp (Conjugation -> Int -> ForeignPtr a -> ContT r IO (Ptr a))
-> Conjugation -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
forall a b. (a -> b) -> a -> b
$
   case Order
order of
      Order
RowMajor -> Conjugation
conj
      Order
ColumnMajor -> Conjugation
NonConjugated


unpackToTemp ::
   Storable a =>
   (Int -> Ptr a -> Ptr a -> IO ()) ->
   Int -> ForeignPtr a -> ContT r IO (Ptr a)
unpackToTemp :: (Int -> Ptr a -> Ptr a -> IO ())
-> Int -> ForeignPtr a -> ContT r IO (Ptr a)
unpackToTemp Int -> Ptr a -> Ptr a -> IO ()
f Int
n ForeignPtr a
a = do
   Ptr a
apPtr <- ((Ptr a -> IO r) -> IO r) -> ContT r IO (Ptr a)
forall k (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO r) -> IO r) -> ContT r IO (Ptr a))
-> ((Ptr a -> IO r) -> IO r) -> ContT r IO (Ptr a)
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> (Ptr a -> IO r) -> IO r
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
a
   Ptr a
aPtr <- Int -> ContT r IO (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)
   IO () -> ContT r IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT r IO ()) -> IO () -> ContT r IO ()
forall a b. (a -> b) -> a -> b
$ Int -> Ptr a -> Ptr a -> IO ()
f Int
n Ptr a
apPtr Ptr a
aPtr
   Ptr a -> ContT r IO (Ptr a)
forall (m :: * -> *) a. Monad m => a -> m a
return Ptr a
aPtr


unpack :: Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
unpack :: Order -> Int -> Ptr a -> Ptr a -> IO ()
unpack Order
order Int
n Ptr a
packedPtr Ptr a
fullPtr =
   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
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim 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
$ String -> String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
errorCodeMsg String
"tpttr" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CChar
-> Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ()
LapackGen.tpttr Ptr CChar
uploPtr Ptr CInt
nPtr Ptr a
packedPtr Ptr a
fullPtr Ptr CInt
ldaPtr

pack :: Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
pack :: Order -> Int -> Ptr a -> Ptr a -> IO ()
pack Order
order Int
n = Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
packRect Order
order Int
n Int
n

packRect :: Class.Floating a => Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
packRect :: Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
packRect Order
order Int
n Int
ld Ptr a
fullPtr Ptr a
packedPtr =
   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
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
ld
      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
$ String -> String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
errorCodeMsg String
"trttp" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CChar
-> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
LapackGen.trttp Ptr CChar
uploPtr Ptr CInt
nPtr Ptr a
fullPtr Ptr CInt
ldaPtr Ptr a
packedPtr


unpackZero, _unpackZero ::
   Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
_unpackZero :: Order -> Int -> Ptr a -> Ptr a -> IO ()
_unpackZero Order
order Int
n Ptr a
packedPtr Ptr a
fullPtr = do
   a -> Int -> Ptr a -> IO ()
forall a. Floating a => a -> Int -> Ptr a -> IO ()
fill a
forall a. Floating a => a
zero (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n) Ptr a
fullPtr
   Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
unpack Order
order Int
n Ptr a
packedPtr Ptr a
fullPtr

unpackZero :: Order -> Int -> Ptr a -> Ptr a -> IO ()
unpackZero Order
order Int
n Ptr a
packedPtr Ptr a
fullPtr = do
   a -> Order -> Int -> Ptr a -> IO ()
forall a. Floating a => a -> Order -> Int -> Ptr a -> IO ()
fillTriangle a
forall a. Floating a => a
zero (Order -> Order
flipOrder Order
order) Int
n Ptr a
fullPtr
   Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
unpack Order
order Int
n Ptr a
packedPtr Ptr a
fullPtr

fillTriangle :: Class.Floating a => a -> Order -> Int -> Ptr a -> IO ()
fillTriangle :: a -> Order -> Int -> Ptr a -> IO ()
fillTriangle a
z Order
order Int
n 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
   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 a
zPtr <- a -> FortranIO () (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
z
   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 CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
LapackGen.laset Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
nPtr Ptr a
zPtr Ptr a
zPtr Ptr a
aPtr Ptr CInt
nPtr



uncheck :: Triangular lo diag up sh a -> Triangular lo diag up (Unchecked sh) a
uncheck :: Triangular lo diag up sh a
-> Triangular lo diag up (Unchecked sh) a
uncheck =
   (Triangular lo diag up sh -> Triangular lo diag up (Unchecked sh))
-> Triangular lo diag up sh a
-> Triangular lo diag up (Unchecked sh) a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape ((Triangular lo diag up sh -> Triangular lo diag up (Unchecked sh))
 -> Triangular lo diag up sh a
 -> Triangular lo diag up (Unchecked sh) a)
-> (Triangular lo diag up sh
    -> Triangular lo diag up (Unchecked sh))
-> Triangular lo diag up sh a
-> Triangular lo diag up (Unchecked sh) a
forall a b. (a -> b) -> a -> b
$
      \(MatrixShape.Triangular diag
diag (lo, up)
uplo Order
order sh
sh) ->
         diag
-> (lo, up)
-> Order
-> Unchecked sh
-> Triangular lo diag up (Unchecked sh)
forall lo diag up size.
diag -> (lo, up) -> Order -> size -> Triangular lo diag up size
MatrixShape.Triangular diag
diag (lo, up)
uplo Order
order (sh -> Unchecked sh
forall sh. sh -> Unchecked sh
Unchecked sh
sh)

recheck :: Triangular lo diag up (Unchecked sh) a -> Triangular lo diag up sh a
recheck :: Triangular lo diag up (Unchecked sh) a
-> Triangular lo diag up sh a
recheck =
   (Triangular lo diag up (Unchecked sh) -> Triangular lo diag up sh)
-> Triangular lo diag up (Unchecked sh) a
-> Triangular lo diag up sh a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape ((Triangular lo diag up (Unchecked sh) -> Triangular lo diag up sh)
 -> Triangular lo diag up (Unchecked sh) a
 -> Triangular lo diag up sh a)
-> (Triangular lo diag up (Unchecked sh)
    -> Triangular lo diag up sh)
-> Triangular lo diag up (Unchecked sh) a
-> Triangular lo diag up sh a
forall a b. (a -> b) -> a -> b
$
      \(MatrixShape.Triangular diag
diag (lo, up)
uplo Order
order (Unchecked sh
sh)) ->
         diag -> (lo, up) -> Order -> sh -> Triangular lo diag up sh
forall lo diag up size.
diag -> (lo, up) -> Order -> size -> Triangular lo diag up size
MatrixShape.Triangular diag
diag (lo, up)
uplo Order
order sh
sh


stack ::
   (Box.Box sh0, Box.HeightOf sh0 ~ height, Shape.C height, Eq height,
    Box.Box sh1, Box.WidthOf sh1 ~ width, Shape.C width, Eq width,
    Shape.C sh2, Class.Floating a) =>
   String -> (height:+:width -> sh2) ->
   Array sh0 a -> Matrix.General height width a -> Array sh1 a -> Array sh2 a
stack :: String
-> ((height :+: width) -> sh2)
-> Array sh0 a
-> General height width a
-> Array sh1 a
-> Array sh2 a
stack String
name (height :+: width) -> sh2
consShape
      (Array sh0
sha ForeignPtr a
a) (Array (MatrixShape.Full Order
order Extent Big Big height width
extent) ForeignPtr a
b) (Array sh1
shc ForeignPtr a
c) =
   let (height
height,width
width) = Extent Big Big height width -> (height, width)
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> (height, width)
Extent.dimensions Extent Big Big height width
extent
   in sh2 -> (Ptr a -> IO ()) -> Array sh2 a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate ((height :+: width) -> sh2
consShape (height
height height -> width -> height :+: width
forall sh0 sh1. sh0 -> sh1 -> sh0 :+: sh1
:+: width
width)) ((Ptr a -> IO ()) -> Array sh2 a)
-> (Ptr a -> IO ()) -> Array sh2 a
forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> do
      String -> Bool -> IO ()
Call.assert (String
nameString -> String -> String
forall a. [a] -> [a] -> [a]
++String
".stack: height shapes mismatch") (Bool -> IO ()) -> Bool -> IO ()
forall a b. (a -> b) -> a -> b
$
         height
height height -> height -> Bool
forall a. Eq a => a -> a -> Bool
== sh0 -> HeightOf sh0
forall shape. Box shape => shape -> HeightOf shape
Box.height sh0
sha
      String -> Bool -> IO ()
Call.assert (String
nameString -> String -> String
forall a. [a] -> [a] -> [a]
++String
".stack: width shapes mismatch") (Bool -> IO ()) -> Bool -> IO ()
forall a b. (a -> b) -> a -> b
$
         width
width width -> width -> Bool
forall a. Eq a => a -> a -> Bool
== sh1 -> WidthOf sh1
forall shape. Box shape => shape -> WidthOf shape
Box.width sh1
shc
      let m :: Int
m = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
      let n :: Int
n = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
width
      ForeignPtr a -> (Ptr a -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
a ((Ptr a -> IO ()) -> IO ()) -> (Ptr a -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr a
aPtr -> (Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
(Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyTriangleA Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
copyBlock Order
order Int
m Int
n Ptr a
aPtr Ptr a
xPtr
      ForeignPtr a -> (Ptr a -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
b ((Ptr a -> IO ()) -> IO ()) -> (Ptr a -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr a
bPtr -> (Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
(Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyRectangle Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
copyBlock Order
order Int
m Int
n Ptr a
bPtr Ptr a
xPtr
      ForeignPtr a -> (Ptr a -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
c ((Ptr a -> IO ()) -> IO ()) -> (Ptr a -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr a
cPtr -> (Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
(Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyTriangleC Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
copyBlock Order
order Int
m Int
n Ptr a
cPtr Ptr a
xPtr

takeTopRight ::
   (Shape.C sh, Shape.C height, Shape.C width, Class.Floating a) =>
   (sh -> (MatrixShape.Order, height:+:width)) ->
   Array sh a -> Matrix.General height width a
takeTopRight :: (sh -> (Order, height :+: width))
-> Array sh a -> General height width a
takeTopRight sh -> (Order, height :+: width)
getShapes (Array sh
sh ForeignPtr a
x) =
   let (Order
order, height
height:+:width
width) = sh -> (Order, height :+: width)
getShapes sh
sh
   in General height width -> (Ptr a -> IO ()) -> General height width a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> height -> width -> General height width
forall height width.
Order -> height -> width -> General height width
MatrixShape.general Order
order height
height width
width) ((Ptr a -> IO ()) -> General height width a)
-> (Ptr a -> IO ()) -> General height width a
forall a b. (a -> b) -> a -> b
$ \Ptr a
bPtr -> do
      let m :: Int
m = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
      let n :: Int
n = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
width
      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
$ (Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
(Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyRectangle ((Ptr a -> Ptr a -> IO ()) -> Ptr a -> Ptr a -> IO ()
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Ptr a -> Ptr a -> IO ()) -> Ptr a -> Ptr a -> IO ())
-> (Int -> Ptr a -> Ptr a -> IO ())
-> Int
-> Ptr a
-> Ptr a
-> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
copyBlock) Order
order Int
m Int
n Ptr a
bPtr

takeTopLeft ::
   (Shape.C sh, Shape.C sha, Shape.C height, Shape.C width, Class.Floating a) =>
   (sh -> (sha, (MatrixShape.Order, height:+:width))) ->
   Array sh a -> Array sha a
takeTopLeft :: (sh -> (sha, (Order, height :+: width)))
-> Array sh a -> Array sha a
takeTopLeft sh -> (sha, (Order, height :+: width))
getShapes (Array sh
sh ForeignPtr a
x) =
   let (sha
sha, (Order
order, height
height:+:width
width)) = sh -> (sha, (Order, height :+: width))
getShapes sh
sh
   in sha -> (Ptr a -> IO ()) -> Array sha a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate sha
sha ((Ptr a -> IO ()) -> Array sha a)
-> (Ptr a -> IO ()) -> Array sha a
forall a b. (a -> b) -> a -> b
$ \Ptr a
aPtr -> do
      let m :: Int
m = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
      let n :: Int
n = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
width
      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
$ (Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
(Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyTriangleA ((Ptr a -> Ptr a -> IO ()) -> Ptr a -> Ptr a -> IO ()
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Ptr a -> Ptr a -> IO ()) -> Ptr a -> Ptr a -> IO ())
-> (Int -> Ptr a -> Ptr a -> IO ())
-> Int
-> Ptr a
-> Ptr a
-> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
copyBlock) Order
order Int
m Int
n Ptr a
aPtr

takeBottomRight ::
   (Shape.C sh, Shape.C shc, Shape.C height, Shape.C width, Class.Floating a) =>
   (sh -> (shc, (MatrixShape.Order, height:+:width))) ->
   Array sh a -> Array shc a
takeBottomRight :: (sh -> (shc, (Order, height :+: width)))
-> Array sh a -> Array shc a
takeBottomRight sh -> (shc, (Order, height :+: width))
getShapes (Array sh
sh ForeignPtr a
x) =
   let (shc
shc, (Order
order, height
height:+:width
width)) = sh -> (shc, (Order, height :+: width))
getShapes sh
sh
   in shc -> (Ptr a -> IO ()) -> Array shc a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate shc
shc ((Ptr a -> IO ()) -> Array shc a)
-> (Ptr a -> IO ()) -> Array shc a
forall a b. (a -> b) -> a -> b
$ \Ptr a
cPtr -> do
      let m :: Int
m = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
      let n :: Int
n = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
width
      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
$ (Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
(Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyTriangleC ((Ptr a -> Ptr a -> IO ()) -> Ptr a -> Ptr a -> IO ()
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Ptr a -> Ptr a -> IO ()) -> Ptr a -> Ptr a -> IO ())
-> (Int -> Ptr a -> Ptr a -> IO ())
-> Int
-> Ptr a
-> Ptr a
-> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
copyBlock) Order
order Int
m Int
n Ptr a
cPtr

{-# INLINE copyTriangleA #-}
copyTriangleA ::
   (Class.Floating a) =>
   (Int -> Ptr a -> Ptr a -> IO ()) ->
   Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyTriangleA :: (Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyTriangleA Int -> Ptr a -> Ptr a -> IO ()
copy Order
order Int
m Int
n Ptr a
aPtr Ptr a
xPtr =
   case Order
order of
      Order
ColumnMajor -> Int -> Ptr a -> Ptr a -> IO ()
copy (Int -> Int
Shape.triangleSize Int
m) Ptr a
aPtr Ptr a
xPtr
      Order
RowMajor ->
         [(Int, (Ptr a, Ptr a))]
-> ((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ ([Int] -> [(Ptr a, Ptr a)] -> [(Int, (Ptr a, Ptr a))]
forall a b. [a] -> [b] -> [(a, b)]
zip ((Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate Int -> Int
forall a. Enum a => a -> a
pred Int
m) ([(Ptr a, Ptr a)] -> [(Int, (Ptr a, Ptr a))])
-> [(Ptr a, Ptr a)] -> [(Int, (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 (Order -> Int -> Ptr a -> [Ptr a]
forall a. Storable a => Order -> Int -> Ptr a -> [Ptr a]
diagonalPointers Order
order Int
m Ptr a
aPtr)
                    (Order -> Int -> Ptr a -> [Ptr a]
forall a. Storable a => Order -> Int -> Ptr a -> [Ptr a]
diagonalPointers Order
order (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
n) Ptr a
xPtr)) (((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ())
-> ((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
            \(Int
k,(Ptr a
aiPtr,Ptr a
xiPtr)) -> Int -> Ptr a -> Ptr a -> IO ()
copy Int
k Ptr a
aiPtr Ptr a
xiPtr

{-# INLINE copyTriangleC #-}
copyTriangleC ::
   (Class.Floating a) =>
   (Int -> Ptr a -> Ptr a -> IO ()) ->
   Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyTriangleC :: (Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyTriangleC Int -> Ptr a -> Ptr a -> IO ()
copy Order
order Int
m Int
n Ptr a
cPtr Ptr a
xPtr =
   case Order
order of
      Order
RowMajor ->
         let triSize :: Int
triSize = Int -> Int
Shape.triangleSize Int
n
         in Int -> Ptr a -> Ptr a -> IO ()
copy Int
triSize Ptr a
cPtr
               (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
xPtr (Int -> Ptr a) -> Int -> Ptr a
forall a b. (a -> b) -> a -> b
$ Int -> Int
Shape.triangleSize (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
n) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
triSize)
      Order
ColumnMajor ->
         [(Int, (Ptr a, Ptr a))]
-> ((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ ([Int] -> [(Ptr a, Ptr a)] -> [(Int, (Ptr a, Ptr a))]
forall a b. [a] -> [b] -> [(a, b)]
zip ((Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate Int -> Int
forall a. Enum a => a -> a
succ Int
0) ([(Ptr a, Ptr a)] -> [(Int, (Ptr a, Ptr a))])
-> [(Ptr a, Ptr a)] -> [(Int, (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 (Order -> Int -> Ptr a -> [Ptr a]
forall a. Storable a => Order -> Int -> Ptr a -> [Ptr a]
diagonalPointers Order
order Int
n Ptr a
cPtr)
                    (Int -> [Ptr a] -> [Ptr a]
forall a. Int -> [a] -> [a]
drop Int
m ([Ptr a] -> [Ptr a]) -> [Ptr a] -> [Ptr a]
forall a b. (a -> b) -> a -> b
$ Order -> Int -> Ptr a -> [Ptr a]
forall a. Storable a => Order -> Int -> Ptr a -> [Ptr a]
diagonalPointers Order
order (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
n) Ptr a
xPtr)) (((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ())
-> ((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
            \(Int
k,(Ptr a
aiPtr,Ptr a
xiPtr)) ->
               Int -> Ptr a -> Ptr a -> IO ()
copy (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
aiPtr (-Int
k)) (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
xiPtr (-Int
k))

{-# INLINE copyRectangle #-}
copyRectangle ::
   (Class.Floating a) =>
   (Int -> Ptr a -> Ptr a -> IO ()) ->
   Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyRectangle :: (Int -> Ptr a -> Ptr a -> IO ())
-> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyRectangle Int -> Ptr a -> Ptr a -> IO ()
copy Order
order Int
m Int
n Ptr a
bPtr Ptr a
xPtr =
   case Order
order of
      Order
RowMajor ->
         [(Int, (Ptr a, Ptr a))]
-> ((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ (Int -> [(Int, (Ptr a, Ptr a))] -> [(Int, (Ptr a, Ptr a))]
forall a. Int -> [a] -> [a]
take Int
m ([(Int, (Ptr a, Ptr a))] -> [(Int, (Ptr a, Ptr a))])
-> [(Int, (Ptr a, Ptr a))] -> [(Int, (Ptr a, Ptr a))]
forall a b. (a -> b) -> a -> b
$ [Int] -> [(Ptr a, Ptr a)] -> [(Int, (Ptr a, Ptr a))]
forall a b. [a] -> [b] -> [(a, b)]
zip ((Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate Int -> Int
forall a. Enum a => a -> a
pred Int
m) ([(Ptr a, Ptr a)] -> [(Int, (Ptr a, Ptr a))])
-> [(Ptr a, Ptr a)] -> [(Int, (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) (Order -> Int -> Ptr a -> [Ptr a]
forall a. Storable a => Order -> Int -> Ptr a -> [Ptr a]
diagonalPointers Order
order (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
n) Ptr a
xPtr)) (((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ())
-> ((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
            \(Int
k,(Ptr a
biPtr,Ptr a
xiPtr)) -> Int -> Ptr a -> Ptr a -> IO ()
copy Int
n Ptr a
biPtr (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
xiPtr Int
k)
      Order
ColumnMajor ->
         [(Int, (Ptr a, Ptr a))]
-> ((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ (Int -> [(Int, (Ptr a, Ptr a))] -> [(Int, (Ptr a, Ptr a))]
forall a. Int -> [a] -> [a]
take Int
n ([(Int, (Ptr a, Ptr a))] -> [(Int, (Ptr a, Ptr a))])
-> [(Int, (Ptr a, Ptr a))] -> [(Int, (Ptr a, Ptr a))]
forall a b. (a -> b) -> a -> b
$ [Int] -> [(Ptr a, Ptr a)] -> [(Int, (Ptr a, Ptr a))]
forall a b. [a] -> [b] -> [(a, b)]
zip ((Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate Int -> Int
forall a. Enum a => a -> a
succ Int
m) ([(Ptr a, Ptr a)] -> [(Int, (Ptr a, Ptr a))])
-> [(Ptr a, Ptr a)] -> [(Int, (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
m Ptr a
bPtr)
                    (Int -> [Ptr a] -> [Ptr a]
forall a. Int -> [a] -> [a]
drop Int
m ([Ptr a] -> [Ptr a]) -> [Ptr a] -> [Ptr a]
forall a b. (a -> b) -> a -> b
$ Order -> Int -> Ptr a -> [Ptr a]
forall a. Storable a => Order -> Int -> Ptr a -> [Ptr a]
diagonalPointers Order
order (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
n) Ptr a
xPtr)) (((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ())
-> ((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
            \(Int
k,(Ptr a
biPtr,Ptr a
xiPtr)) -> Int -> Ptr a -> Ptr a -> IO ()
copy Int
m Ptr a
biPtr (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
xiPtr (-Int
k))



type Triangular lo diag up sh = Array (MatrixShape.Triangular lo diag up sh)

type FlexDiagonal diag sh =
         Triangular MatrixShape.Empty diag MatrixShape.Empty sh

newtype MultiplyRight diag sh a b lo up =
   MultiplyRight {MultiplyRight diag sh a b lo up -> Triangular lo diag up sh a -> b
getMultiplyRight :: Triangular lo diag up sh a -> b}

newtype Map diag sh0 sh1 a lo up =
   Map {Map diag sh0 sh1 a lo up
-> Triangular lo diag up sh0 a -> Triangular lo diag up sh1 a
getMap :: Triangular lo diag up sh0 a -> Triangular lo diag up sh1 a}

newtype Power diag sh a lo up =
   Power {
      Power diag sh a lo up
-> Triangular lo diag up sh a
-> Triangular lo (PowerDiag lo up diag) up sh a
getPower ::
         Triangular lo diag up sh a ->
         Triangular lo (PowerDiag lo up diag) up sh a
   }

type family PowerDiag lo up diag
type instance PowerDiag Empty up diag = diag
type instance PowerDiag Filled Empty diag = diag
type instance PowerDiag Filled Filled diag = NonUnit

type PowerContentDiag lo diag up =
      (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag,
       PowerDiag lo up diag ~ diag, PowerDiag up lo diag ~ diag)


fromBanded ::
   (Class.Floating a) =>
   Int -> Order -> Int -> ForeignPtr a -> Int -> Ptr a -> IO ()
fromBanded :: Int -> Order -> Int -> ForeignPtr a -> Int -> Ptr a -> IO ()
fromBanded Int
k Order
order Int
n ForeignPtr a
a Int
bSize 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 -> 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
bPtr
      let lda :: Int
lda = Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
      let pointers :: [(Int, (Ptr a, Ptr a))]
pointers =
            [Int] -> [(Ptr a, Ptr a)] -> [(Int, (Ptr a, Ptr a))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0..] ([(Ptr a, Ptr a)] -> [(Int, (Ptr a, Ptr a))])
-> [(Ptr a, Ptr a)] -> [(Int, (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
lda Ptr a
aPtr) ([Ptr a] -> [(Ptr a, Ptr a)]) -> [Ptr a] -> [(Ptr a, Ptr a)]
forall a b. (a -> b) -> a -> b
$
            Order -> Int -> Ptr a -> [Ptr a]
forall a. Storable a => Order -> Int -> Ptr a -> [Ptr a]
diagonalPointers Order
order Int
n Ptr a
bPtr
      case Order
order of
         Order
ColumnMajor ->
            [(Int, (Ptr a, Ptr a))]
-> ((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [(Int, (Ptr a, Ptr a))]
pointers (((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ())
-> ((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \(Int
i,(Ptr a
xPtr,Ptr a
yPtr)) ->
               let j :: Int
j = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
i Int
k
               in Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
copyBlock (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
xPtr (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
j)) (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
yPtr (-Int
j))
         Order
RowMajor ->
            [(Int, (Ptr a, Ptr a))]
-> ((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [(Int, (Ptr a, Ptr a))]
pointers (((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ())
-> ((Int, (Ptr a, Ptr a)) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \(Int
i,(Ptr a
xPtr,Ptr a
yPtr)) ->
               Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
copyBlock (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
lda (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
i)) Ptr a
xPtr Ptr a
yPtr


type FlexLower diag sh = Array (MatrixShape.LowerTriangular diag sh)

takeLower ::
   (MatrixShape.TriDiag diag,
    Extent.C horiz, Shape.C height, Shape.C width, Class.Floating a) =>
   (diag, Order -> Int -> Ptr a -> IO ()) ->
   Full Extent.Small horiz height width a -> FlexLower diag height a
takeLower :: (diag, Order -> Int -> Ptr a -> IO ())
-> Full Small horiz height width a -> FlexLower diag height a
takeLower (diag
diag, Order -> Int -> Ptr a -> IO ()
fillDiag) (Array (MatrixShape.Full Order
order Extent Small horiz height width
extent) ForeignPtr a
a) =
   let (height
height,width
width) = Extent Small horiz height width -> (height, width)
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> (height, width)
Extent.dimensions Extent Small horiz height width
extent
       m :: Int
m = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
       n :: Int
n = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
width
       k :: Int
k = case Order
order of Order
RowMajor -> Int
n; Order
ColumnMajor -> Int
m
   in Triangular Filled diag Empty height
-> (Ptr a -> IO ()) -> FlexLower diag height a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate
         (diag
-> (Filled, Empty)
-> Order
-> height
-> Triangular Filled diag Empty height
forall lo diag up size.
diag -> (lo, up) -> Order -> size -> Triangular lo diag up size
MatrixShape.Triangular diag
diag (Filled, Empty)
MatrixShape.lower Order
order height
height) ((Ptr a -> IO ()) -> FlexLower diag height a)
-> (Ptr a -> IO ()) -> FlexLower diag height a
forall a b. (a -> b) -> a -> b
$ \Ptr a
lPtr ->
      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 -> do
         let dstOrder :: Order
dstOrder = Order -> Order
flipOrder Order
order
         Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
packRect Order
dstOrder Int
m Int
k Ptr a
aPtr Ptr a
lPtr
         Order -> Int -> Ptr a -> IO ()
fillDiag Order
dstOrder Int
m Ptr a
lPtr


fromUpperPart ::
   (Extent.C vert, Shape.C height, Shape.C width, Shape.C shape,
    Class.Floating a) =>
   (Order -> width -> shape) ->
   Full vert Extent.Small height width a -> Array shape a
fromUpperPart :: (Order -> width -> shape)
-> Full vert Small height width a -> Array shape a
fromUpperPart Order -> width -> shape
shape (Array (MatrixShape.Full Order
order Extent vert Small height width
extent) ForeignPtr a
a) =
   let (height
height,width
width) = Extent vert Small height width -> (height, width)
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> (height, width)
Extent.dimensions Extent vert Small height width
extent
       m :: Int
m = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
       n :: Int
n = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
width
       k :: Int
k = case Order
order of Order
RowMajor -> Int
n; Order
ColumnMajor -> Int
m
   in shape -> (Ptr a -> IO ()) -> Array shape a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> width -> shape
shape Order
order width
width) ((Ptr a -> IO ()) -> Array shape a)
-> (Ptr a -> IO ()) -> Array shape 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 -> Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
packRect Order
order Int
n Int
k Ptr a
aPtr Ptr a
bPtr