{-# LANGUAGE TypeOperators #-}
module Numeric.LAPACK.Matrix.Square.Basic (
   Square,
   mapSize,
   toFull,
   fromGeneral,
   fromScalar,
   toScalar,
   fromList,
   autoFromList,

   transpose,
   adjoint,

   identity,
   identityOrder,
   identityFrom,
   identityFromWidth,
   identityFromHeight,
   diagonal,
   takeDiagonal,
   trace,

   stack,

   square,
   power,
   congruence,
   congruenceAdjoint,
   ) where


import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
import qualified Numeric.LAPACK.Matrix.Extent.Private as ExtentPriv
import qualified Numeric.LAPACK.Matrix.Extent as Extent
import qualified Numeric.LAPACK.Matrix.Basic as Basic
import qualified Numeric.LAPACK.Matrix.Private as Matrix
import qualified Numeric.LAPACK.Vector as Vector
import qualified Numeric.LAPACK.Private as Private
import Numeric.LAPACK.Matrix.Shape.Private
         (Order(RowMajor, ColumnMajor), swapOnRowMajor)
import Numeric.LAPACK.Matrix.Private
         (Full, mapExtent,
          General, argGeneral, Square, argSquare, ShapeInt, shapeInt)
import Numeric.LAPACK.Vector (Vector)
import Numeric.LAPACK.Scalar (zero, one)
import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))
import Numeric.LAPACK.Private (pokeCInt)

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

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

import Foreign.ForeignPtr (withForeignPtr)
import Foreign.Storable (Storable, peek, poke)

import System.IO.Unsafe (unsafePerformIO)

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

import Data.Function.HT (powerAssociative)


mapSize :: (sh0 -> sh1) -> Square sh0 a -> Square sh1 a
mapSize :: (sh0 -> sh1) -> Square sh0 a -> Square sh1 a
mapSize sh0 -> sh1
f =
   (Full Small Small sh0 sh0 -> Full Small Small sh1 sh1)
-> Square sh0 a -> Square sh1 a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape
      (\(MatrixShape.Full Order
order Extent Small Small sh0 sh0
extent) ->
         Order -> Extent Small Small sh1 sh1 -> Full Small Small sh1 sh1
forall vert horiz height width.
Order
-> Extent vert horiz height width -> Full vert horiz height width
MatrixShape.Full Order
order (Extent Small Small sh1 sh1 -> Full Small Small sh1 sh1)
-> Extent Small Small sh1 sh1 -> Full Small Small sh1 sh1
forall a b. (a -> b) -> a -> b
$ (sh0 -> sh1)
-> Extent Small Small sh0 sh0 -> Extent Small Small sh1 sh1
forall shA shB. (shA -> shB) -> Square shA -> Square shB
ExtentPriv.mapSquareSize sh0 -> sh1
f Extent Small Small sh0 sh0
extent)

toFull ::
   (Extent.C vert, Extent.C horiz) => Square sh a -> Full vert horiz sh sh a
toFull :: Square sh a -> Full vert horiz sh sh a
toFull = Map Small Small vert horiz sh sh
-> Square sh a -> Full vert horiz sh sh a
forall vertA horizA vertB horizB height width a.
(C vertA, C horizA, C vertB, C horizB) =>
Map vertA horizA vertB horizB height width
-> Full vertA horizA height width a
-> Full vertB horizB height width a
mapExtent Map Small Small vert horiz sh sh
forall vert horiz size.
(C vert, C horiz) =>
Map Small Small vert horiz size size
Extent.fromSquare

fromGeneral :: (Eq sh) => General sh sh a -> Square sh a
fromGeneral :: General sh sh a -> Square sh a
fromGeneral = Map Big Big Small Small sh sh -> General sh sh a -> Square sh a
forall vertA horizA vertB horizB height width a.
(C vertA, C horizA, C vertB, C horizB) =>
Map vertA horizA vertB horizB height width
-> Full vertA horizA height width a
-> Full vertB horizB height width a
mapExtent ((Extent Big Big sh sh -> Extent Small Small sh sh)
-> Map Big Big Small Small sh sh
forall vertA horizA vertB horizB height width.
(Extent vertA horizA height width
 -> Extent vertB horizB height width)
-> Map vertA horizA vertB horizB height width
ExtentPriv.Map Extent Big Big sh sh -> Extent Small Small sh sh
forall vert horiz size.
(C vert, C horiz, Eq size) =>
Extent vert horiz size size -> Square size
ExtentPriv.squareFromGeneral)


fromScalar :: (Storable a) => a -> Square () a
fromScalar :: a -> Square () a
fromScalar a
a =
   Square () -> (Ptr a -> IO ()) -> Square () a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> () -> Square ()
forall sh. Order -> sh -> Square sh
MatrixShape.square Order
RowMajor ()) ((Ptr a -> IO ()) -> Square () a)
-> (Ptr a -> IO ()) -> Square () a
forall a b. (a -> b) -> a -> b
$ (Ptr a -> a -> IO ()) -> a -> Ptr a -> IO ()
forall a b c. (a -> b -> c) -> b -> a -> c
flip Ptr a -> a -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke a
a

toScalar :: (Storable a) => Square () a -> a
toScalar :: Square () a -> a
toScalar = (Order -> () -> ForeignPtr a -> a) -> Square () a -> a
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> () -> ForeignPtr a -> a) -> Square () a -> a)
-> (Order -> () -> ForeignPtr a -> a) -> Square () a -> a
forall a b. (a -> b) -> a -> b
$ \Order
_ () ForeignPtr a
a ->
   IO a -> a
forall a. IO a -> a
unsafePerformIO (IO a -> a) -> IO a -> a
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> (Ptr a -> IO a) -> IO a
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
a Ptr a -> IO a
forall a. Storable a => Ptr a -> IO a
peek

fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Square sh a
fromList :: sh -> [a] -> Square sh a
fromList sh
sh =
   Square sh -> [a] -> Square sh a
forall sh a. (C sh, Storable a) => sh -> [a] -> Array sh a
CheckedArray.fromList (Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
MatrixShape.square Order
RowMajor sh
sh)

autoFromList :: (Storable a) => [a] -> Square ShapeInt a
autoFromList :: [a] -> Square ShapeInt a
autoFromList [a]
xs =
   let n :: Int
n = [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs
       m :: Int
m = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n :: Double)
   in if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
m
        then ShapeInt -> [a] -> Square ShapeInt a
forall sh a. (C sh, Storable a) => sh -> [a] -> Square sh a
fromList (Int -> ShapeInt
shapeInt Int
m) [a]
xs
        else [Char] -> Square ShapeInt a
forall a. HasCallStack => [Char] -> a
error [Char]
"Square.autoFromList: no quadratic number of elements"


transpose :: Square sh a -> Square sh a
transpose :: Square sh a -> Square sh a
transpose = (Full Small Small sh sh -> Full Small Small sh sh)
-> Square sh a -> Square sh a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape Full Small Small sh sh -> Full Small Small sh sh
forall vert horiz height width.
(C vert, C horiz) =>
Full vert horiz height width -> Full horiz vert width height
MatrixShape.transpose

adjoint :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a
adjoint :: Square sh a -> Square sh a
adjoint = Square sh a -> Square sh a
forall sh a. Square sh a -> Square sh a
transpose (Square sh a -> Square sh a)
-> (Square sh a -> Square sh a) -> Square sh a -> Square sh a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Square sh a -> Square sh a
forall sh a. (C sh, Floating a) => Vector sh a -> Vector sh a
Vector.conjugate


identity :: (Shape.C sh, Class.Floating a) => sh -> Square sh a
identity :: sh -> Square sh a
identity = Order -> sh -> Square sh a
forall sh a. (C sh, Floating a) => Order -> sh -> Square sh a
identityOrder Order
RowMajor

identityFrom :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a
identityFrom :: Square sh a -> Square sh a
identityFrom = (Order -> sh -> ForeignPtr a -> Square sh a)
-> Square sh a -> Square sh a
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr a -> Square sh a)
 -> Square sh a -> Square sh a)
-> (Order -> sh -> ForeignPtr a -> Square sh a)
-> Square sh a
-> Square sh a
forall a b. (a -> b) -> a -> b
$ \Order
order sh
sh ForeignPtr a
_ -> Order -> sh -> Square sh a
forall sh a. (C sh, Floating a) => Order -> sh -> Square sh a
identityOrder Order
order sh
sh

identityFromWidth ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   General height width a -> Square width a
identityFromWidth :: General height width a -> Square width a
identityFromWidth =
   (Order -> height -> width -> ForeignPtr a -> Square width a)
-> General height width a -> Square width a
forall height width a b.
(Order -> height -> width -> ForeignPtr a -> b)
-> General height width a -> b
argGeneral ((Order -> height -> width -> ForeignPtr a -> Square width a)
 -> General height width a -> Square width a)
-> (Order -> height -> width -> ForeignPtr a -> Square width a)
-> General height width a
-> Square width a
forall a b. (a -> b) -> a -> b
$ \Order
order height
_ width
width ForeignPtr a
_ -> Order -> width -> Square width a
forall sh a. (C sh, Floating a) => Order -> sh -> Square sh a
identityOrder Order
order width
width

identityFromHeight ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   General height width a -> Square height a
identityFromHeight :: General height width a -> Square height a
identityFromHeight =
   (Order -> height -> width -> ForeignPtr a -> Square height a)
-> General height width a -> Square height a
forall height width a b.
(Order -> height -> width -> ForeignPtr a -> b)
-> General height width a -> b
argGeneral ((Order -> height -> width -> ForeignPtr a -> Square height a)
 -> General height width a -> Square height a)
-> (Order -> height -> width -> ForeignPtr a -> Square height a)
-> General height width a
-> Square height a
forall a b. (a -> b) -> a -> b
$ \Order
order height
height width
_ ForeignPtr a
_ -> Order -> height -> Square height a
forall sh a. (C sh, Floating a) => Order -> sh -> Square sh a
identityOrder Order
order height
height

identityOrder, _identityOrder ::
   (Shape.C sh, Class.Floating a) => Order -> sh -> Square sh a
identityOrder :: Order -> sh -> Square sh a
identityOrder Order
order sh
sh =
   Square sh -> (Ptr a -> IO ()) -> Square sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
MatrixShape.square Order
order sh
sh) ((Ptr a -> IO ()) -> Square sh a)
-> (Ptr a -> IO ()) -> Square sh a
forall a b. (a -> b) -> a -> b
$ \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
'A'
      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
$ sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
      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
zero
      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
one
      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
alphaPtr Ptr a
betaPtr Ptr a
aPtr Ptr CInt
nPtr

_identityOrder :: Order -> sh -> Square sh a
_identityOrder Order
order sh
sh =
   Square sh -> (Int -> Ptr a -> IO ()) -> Square sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
MatrixShape.square Order
order sh
sh) ((Int -> Ptr a -> IO ()) -> Square sh a)
-> (Int -> Ptr a -> IO ()) -> Square sh a
forall a b. (a -> b) -> a -> b
$ \Int
blockSize 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
      Ptr CInt
nPtr <- FortranIO () (Ptr CInt)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      Ptr a
xPtr <- 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
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
      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 CInt -> Int -> IO ()
pokeCInt Ptr CInt
nPtr Int
blockSize
         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
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr
         let n :: CInt
n = Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CInt) -> Int -> CInt
forall a b. (a -> b) -> a -> b
$ sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
         Ptr CInt -> CInt -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr CInt
nPtr CInt
n
         Ptr a -> a -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr a
xPtr a
forall a. Floating a => a
one
         Ptr CInt -> CInt -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr CInt
incyPtr (CInt
nCInt -> CInt -> CInt
forall a. Num a => a -> a -> a
+CInt
1)
         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
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr

diagonal :: (Shape.C sh, Class.Floating a) => Vector sh a -> Square sh a
diagonal :: Vector sh a -> Square sh a
diagonal (Array sh
sh ForeignPtr a
x) =
   Square sh -> (Int -> Ptr a -> IO ()) -> Square sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
MatrixShape.square Order
ColumnMajor sh
sh) ((Int -> Ptr a -> IO ()) -> Square sh a)
-> (Int -> Ptr a -> IO ()) -> Square sh a
forall a b. (a -> b) -> a -> b
$
      \Int
blockSize 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
      Ptr CInt
nPtr <- FortranIO () (Ptr CInt)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      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 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
      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
      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 -> Int -> IO ()
pokeCInt Ptr CInt
nPtr Int
blockSize
         Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.copy Ptr CInt
nPtr Ptr a
zPtr Ptr CInt
inczPtr Ptr a
yPtr Ptr CInt
incyPtr
         let n :: CInt
n = Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CInt) -> Int -> CInt
forall a b. (a -> b) -> a -> b
$ sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
         Ptr CInt -> CInt -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr CInt
nPtr CInt
n
         Ptr CInt -> CInt -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr CInt
incyPtr (CInt
nCInt -> CInt -> CInt
forall a. Num a => a -> a -> a
+CInt
1)
         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
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr

takeDiagonal :: (Shape.C sh, Class.Floating a) => Square sh a -> Vector sh a
takeDiagonal :: Square sh a -> Vector sh a
takeDiagonal = (Order -> sh -> ForeignPtr a -> Vector sh a)
-> Square sh a -> Vector sh a
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr a -> Vector sh a)
 -> Square sh a -> Vector sh a)
-> (Order -> sh -> ForeignPtr a -> Vector sh a)
-> Square sh a
-> Vector sh a
forall a b. (a -> b) -> a -> b
$ \Order
_ sh
sh ForeignPtr a
x ->
   sh -> (Int -> Ptr a -> IO ()) -> Vector sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
sh ((Int -> Ptr a -> IO ()) -> Vector sh a)
-> (Int -> Ptr a -> IO ()) -> Vector sh a
forall a b. (a -> b) -> a -> b
$ \Int
n 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
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      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
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (Int
nInt -> 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 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
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr

trace :: (Shape.C sh, Class.Floating a) => Square sh a -> a
trace :: Square sh a -> a
trace = (Order -> sh -> ForeignPtr a -> a) -> Square sh a -> a
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr a -> a) -> Square sh a -> a)
-> (Order -> sh -> ForeignPtr a -> a) -> Square sh a -> a
forall a b. (a -> b) -> a -> b
$ \Order
_ sh
sh ForeignPtr a
x -> IO a -> a
forall a. IO a -> a
unsafePerformIO (IO a -> a) -> IO a -> a
forall a b. (a -> b) -> a -> b
$ do
   let n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
   ForeignPtr a -> (Ptr a -> IO a) -> IO a
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
x ((Ptr a -> IO a) -> IO a) -> (Ptr a -> IO a) -> IO a
forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> Int -> Ptr a -> Int -> IO a
forall a. Floating a => Int -> Ptr a -> Int -> IO a
Private.sum Int
n Ptr a
xPtr (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)


stack ::
   (Extent.C vert, Extent.C horiz,
    Shape.C sizeA, Eq sizeA, Shape.C sizeB, Eq sizeB, Class.Floating a) =>
   Square sizeA a -> Full vert horiz sizeA sizeB a ->
   Full horiz vert sizeB sizeA a -> Square sizeB a ->
   Square (sizeA:+:sizeB) a
stack :: Square sizeA a
-> Full vert horiz sizeA sizeB a
-> Full horiz vert sizeB sizeA a
-> Square sizeB a
-> Square (sizeA :+: sizeB) a
stack Square sizeA a
a Full vert horiz sizeA sizeB a
b Full horiz vert sizeB sizeA a
c Square sizeB a
d = Square sizeA a
-> General sizeA sizeB a
-> General sizeB sizeA a
-> Square sizeB a
-> Square (sizeA :+: sizeB) a
forall vert horiz heightA heightB widthA widthB a.
(C vert, C horiz, C heightA, Eq heightA, C heightB, Eq heightB,
 C widthA, Eq widthA, C widthB, Eq widthB, Floating a) =>
Full vert horiz heightA widthA a
-> General heightA widthB a
-> General heightB widthA a
-> Full vert horiz heightB widthB a
-> Full vert horiz (heightA :+: heightB) (widthA :+: widthB) a
Basic.stack Square sizeA a
a (Full vert horiz sizeA sizeB a -> General sizeA sizeB a
forall vert horiz height width a.
(C vert, C horiz) =>
Full vert horiz height width a -> General height width a
Matrix.fromFull Full vert horiz sizeA sizeB a
b) (Full horiz vert sizeB sizeA a -> General sizeB sizeA a
forall vert horiz height width a.
(C vert, C horiz) =>
Full vert horiz height width a -> General height width a
Matrix.fromFull Full horiz vert sizeB sizeA a
c) Square sizeB a
d


square :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a
square :: Square sh a -> Square sh a
square Square sh a
a = Square sh a -> Square sh a -> Square sh a
forall sh a.
(C sh, Floating a) =>
Square sh a -> Square sh a -> Square sh a
multiplyCommutativeUnchecked Square sh a
a Square sh a
a

power ::
   (Shape.C sh, Class.Floating a) =>
   Integer -> Square sh a -> Square sh a
power :: Integer -> Square sh a -> Square sh a
power Integer
n Square sh a
a =
   (Square sh a -> Square sh a -> Square sh a)
-> Square sh a -> Square sh a -> Integer -> Square sh a
forall a. (a -> a -> a) -> a -> a -> Integer -> a
powerAssociative Square sh a -> Square sh a -> Square sh a
forall sh a.
(C sh, Floating a) =>
Square sh a -> Square sh a -> Square sh a
multiplyCommutativeUnchecked (Square sh a -> Square sh a
forall sh a. (C sh, Floating a) => Square sh a -> Square sh a
identityFrom Square sh a
a) Square sh a
a Integer
n

congruence ::
   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Square height a -> General height width a -> Square width a
congruence :: Square height a -> General height width a -> Square width a
congruence Square height a
b General height width a
a0 =
   let a :: Full Big Big height (Unchecked width) a
a = (width -> Unchecked width)
-> General height width a
-> Full Big Big height (Unchecked width) a
forall vert horiz widthA widthB height a.
(GeneralTallWide vert horiz, GeneralTallWide horiz vert) =>
(widthA -> widthB)
-> Full vert horiz height widthA a
-> Full vert horiz height widthB a
Basic.mapWidth width -> Unchecked width
forall sh. sh -> Unchecked sh
Unchecked General height width a
a0
   in (Unchecked width -> width)
-> Square (Unchecked width) a -> Square width a
forall sh0 sh1 a. (sh0 -> sh1) -> Square sh0 a -> Square sh1 a
mapSize (\(Unchecked width
sh) -> width
sh) (Square (Unchecked width) a -> Square width a)
-> Square (Unchecked width) a -> Square width a
forall a b. (a -> b) -> a -> b
$ General (Unchecked width) (Unchecked width) a
-> Square (Unchecked width) a
forall sh a. Eq sh => General sh sh a -> Square sh a
fromGeneral (General (Unchecked width) (Unchecked width) a
 -> Square (Unchecked width) a)
-> General (Unchecked width) (Unchecked width) a
-> Square (Unchecked width) a
forall a b. (a -> b) -> a -> b
$
      Full Big Big (Unchecked width) height a
-> Full Big Big height (Unchecked width) a
-> General (Unchecked width) (Unchecked width) a
forall vert horiz height fuse width a.
(C vert, C horiz, C height, C fuse, Eq fuse, C width,
 Floating a) =>
Full vert horiz height fuse a
-> Full vert horiz fuse width a -> Full vert horiz height width a
Basic.multiply (Full Big Big height (Unchecked width) a
-> Full Big Big (Unchecked width) height a
forall vert horiz height width a.
(C vert, C horiz, C height, C width, Floating a) =>
Full vert horiz height width a -> Full horiz vert width height a
Basic.adjoint Full Big Big height (Unchecked width) a
a) (Full Big Big height (Unchecked width) a
 -> General (Unchecked width) (Unchecked width) a)
-> Full Big Big height (Unchecked width) a
-> General (Unchecked width) (Unchecked width) a
forall a b. (a -> b) -> a -> b
$ Full Big Big height height a
-> Full Big Big height (Unchecked width) a
-> Full Big Big height (Unchecked width) a
forall vert horiz height fuse width a.
(C vert, C horiz, C height, C fuse, Eq fuse, C width,
 Floating a) =>
Full vert horiz height fuse a
-> Full vert horiz fuse width a -> Full vert horiz height width a
Basic.multiply (Square height a -> Full Big Big height height a
forall vert horiz sh a.
(C vert, C horiz) =>
Square sh a -> Full vert horiz sh sh a
toFull Square height a
b) Full Big Big height (Unchecked width) a
a

congruenceAdjoint ::
   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
   General height width a -> Square width a -> Square height a
congruenceAdjoint :: General height width a -> Square width a -> Square height a
congruenceAdjoint General height width a
a0 Square width a
b =
   let a :: Full Big Big (Unchecked height) width a
a = (height -> Unchecked height)
-> General height width a
-> Full Big Big (Unchecked height) width a
forall vert horiz heightA heightB width a.
(GeneralTallWide vert horiz, GeneralTallWide horiz vert) =>
(heightA -> heightB)
-> Full vert horiz heightA width a
-> Full vert horiz heightB width a
Basic.mapHeight height -> Unchecked height
forall sh. sh -> Unchecked sh
Unchecked General height width a
a0
   in (Unchecked height -> height)
-> Square (Unchecked height) a -> Square height a
forall sh0 sh1 a. (sh0 -> sh1) -> Square sh0 a -> Square sh1 a
mapSize (\(Unchecked height
sh) -> height
sh) (Square (Unchecked height) a -> Square height a)
-> Square (Unchecked height) a -> Square height a
forall a b. (a -> b) -> a -> b
$ General (Unchecked height) (Unchecked height) a
-> Square (Unchecked height) a
forall sh a. Eq sh => General sh sh a -> Square sh a
fromGeneral (General (Unchecked height) (Unchecked height) a
 -> Square (Unchecked height) a)
-> General (Unchecked height) (Unchecked height) a
-> Square (Unchecked height) a
forall a b. (a -> b) -> a -> b
$
      Full Big Big (Unchecked height) width a
-> Full Big Big width (Unchecked height) a
-> General (Unchecked height) (Unchecked height) a
forall vert horiz height fuse width a.
(C vert, C horiz, C height, C fuse, Eq fuse, C width,
 Floating a) =>
Full vert horiz height fuse a
-> Full vert horiz fuse width a -> Full vert horiz height width a
Basic.multiply Full Big Big (Unchecked height) width a
a (Full Big Big width (Unchecked height) a
 -> General (Unchecked height) (Unchecked height) a)
-> Full Big Big width (Unchecked height) a
-> General (Unchecked height) (Unchecked height) a
forall a b. (a -> b) -> a -> b
$ Full Big Big width width a
-> Full Big Big width (Unchecked height) a
-> Full Big Big width (Unchecked height) a
forall vert horiz height fuse width a.
(C vert, C horiz, C height, C fuse, Eq fuse, C width,
 Floating a) =>
Full vert horiz height fuse a
-> Full vert horiz fuse width a -> Full vert horiz height width a
Basic.multiply (Square width a -> Full Big Big width width a
forall vert horiz sh a.
(C vert, C horiz) =>
Square sh a -> Full vert horiz sh sh a
toFull Square width a
b) (Full Big Big width (Unchecked height) a
 -> Full Big Big width (Unchecked height) a)
-> Full Big Big width (Unchecked height) a
-> Full Big Big width (Unchecked height) a
forall a b. (a -> b) -> a -> b
$ Full Big Big (Unchecked height) width a
-> Full Big Big width (Unchecked height) a
forall vert horiz height width a.
(C vert, C horiz, C height, C width, Floating a) =>
Full vert horiz height width a -> Full horiz vert width height a
Basic.adjoint Full Big Big (Unchecked height) width a
a

{-
orderA and orderB must be equal but this is not checked.
-}
multiplyCommutativeUnchecked ::
   (Shape.C sh, Class.Floating a) =>
   Square sh a -> Square sh a -> Square sh a
multiplyCommutativeUnchecked :: Square sh a -> Square sh a -> Square sh a
multiplyCommutativeUnchecked
   (Array shape :: Square sh
shape@(MatrixShape.Full Order
order Extent Small Small sh sh
extent) ForeignPtr a
a)
   (Array Square sh
_ ForeignPtr a
b) =
      Square sh -> (Ptr a -> IO ()) -> Square sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate Square sh
shape ((Ptr a -> IO ()) -> Square sh a)
-> (Ptr a -> IO ()) -> Square sh a
forall a b. (a -> b) -> a -> b
$ \Ptr a
cPtr ->
   let n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size (sh -> Int) -> sh -> Int
forall a b. (a -> b) -> a -> b
$ Extent Small Small sh sh -> sh
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> height
Extent.height Extent Small Small sh sh
extent
       (ForeignPtr a
at,ForeignPtr a
bt) = Order
-> (ForeignPtr a, ForeignPtr a) -> (ForeignPtr a, ForeignPtr a)
forall a. Order -> (a, a) -> (a, a)
swapOnRowMajor Order
order (ForeignPtr a
a,ForeignPtr a
b)
   in  Order
-> Order
-> Int
-> Int
-> Int
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> IO ()
forall a.
Floating a =>
Order
-> Order
-> Int
-> Int
-> Int
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> IO ()
Private.multiplyMatrix Order
ColumnMajor Order
ColumnMajor Int
n Int
n Int
n ForeignPtr a
at ForeignPtr a
bt Ptr a
cPtr