{-# LANGUAGE TypeFamilies #-}
module Numeric.LAPACK.Orthogonal.Plain (
   leastSquares,
   minimumNorm,
   leastSquaresMinimumNormRCond,
   pseudoInverseRCond,

   leastSquaresConstraint,
   gaussMarkovLinearModel,

   determinantAbsolute,
   complement,
   extractComplement,
   ) where

import qualified Numeric.LAPACK.Orthogonal.Private as HH

import qualified Numeric.LAPACK.Matrix.Square.Basic as Square
import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Matrix.Basic as Basic
import qualified Numeric.LAPACK.Vector as Vector
import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor,ColumnMajor))
import Numeric.LAPACK.Matrix.Private
         (Full, General, Tall, Wide, ShapeInt, shapeInt)
import Numeric.LAPACK.Vector (Vector)
import Numeric.LAPACK.Scalar (RealOf, zero, absolute)
import Numeric.LAPACK.Private
         (lacgv, peekCInt,
          copySubMatrix, copyToTemp,
          copyToColumnMajorTemp, copyToSubColumnMajor,
          withAutoWorkspaceInfo, rankMsg, errorCodeMsg, createHigherArray)

import qualified Numeric.LAPACK.FFI.Generic as LapackGen
import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
import qualified Numeric.LAPACK.FFI.Real as LapackReal
import qualified Numeric.Netlib.Utility as Call
import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
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 System.IO.Unsafe (unsafePerformIO)

import Foreign.Marshal.Array (pokeArray)
import Foreign.C.Types (CInt, CChar)
import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
import Foreign.Ptr (Ptr)

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

import Data.Complex (Complex)
import Data.Tuple.HT (mapSnd)


leastSquares ::
   (Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   Full horiz Extent.Small height width a ->
   Full vert horiz height nrhs a ->
   Full vert horiz width nrhs a
leastSquares :: Full horiz Small height width a
-> Full vert horiz height nrhs a -> Full vert horiz width nrhs a
leastSquares
   (Array shapeA :: Full horiz Small height width
shapeA@(MatrixShape.Full Order
orderA Extent horiz Small height width
extentA) ForeignPtr a
a)
   (Array        (MatrixShape.Full Order
orderB Extent vert horiz height nrhs
extentB) ForeignPtr a
b) =

 case Extent vert horiz width height
-> Extent vert horiz height nrhs
-> Maybe (Extent vert horiz width nrhs)
forall vert horiz fuse height width.
(C vert, C horiz, Eq fuse) =>
Extent vert horiz height fuse
-> Extent vert horiz fuse width
-> Maybe (Extent vert horiz height width)
Extent.fuse (Extent Small horiz width height -> Extent vert horiz width height
forall vert horiz height width.
(C vert, C horiz) =>
Extent Small horiz height width -> Extent vert horiz height width
Extent.generalizeWide (Extent Small horiz width height -> Extent vert horiz width height)
-> Extent Small horiz width height
-> Extent vert horiz width height
forall a b. (a -> b) -> a -> b
$ Extent horiz Small height width -> Extent Small horiz width height
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> Extent horiz vert width height
Extent.transpose Extent horiz Small height width
extentA) Extent vert horiz height nrhs
extentB of
  Maybe (Extent vert horiz width nrhs)
Nothing -> [Char] -> Full vert horiz width nrhs a
forall a. HasCallStack => [Char] -> a
error [Char]
"leastSquares: height shapes mismatch"
  Just Extent vert horiz width nrhs
extent ->
      Full vert horiz width nrhs
-> (Ptr a -> IO ()) -> Full vert horiz width nrhs a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> Extent vert horiz width nrhs -> Full vert horiz width nrhs
forall vert horiz height width.
Order
-> Extent vert horiz height width -> Full vert horiz height width
MatrixShape.Full Order
ColumnMajor Extent vert horiz width nrhs
extent) ((Ptr a -> IO ()) -> Full vert horiz width nrhs a)
-> (Ptr a -> IO ()) -> Full vert horiz width nrhs a
forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> do

   let widthA :: width
widthA = Extent horiz Small height width -> width
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> width
Extent.width Extent horiz Small height width
extentA
   let (height
height,nrhs
widthB) = Extent vert horiz height nrhs -> (height, nrhs)
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> (height, width)
Extent.dimensions Extent vert horiz height nrhs
extentB
   let (Int
m,Int
n) = Full horiz Small height width -> (Int, Int)
forall vert horiz height width.
(C vert, C horiz, C height, C width) =>
Full vert horiz height width -> (Int, Int)
MatrixShape.dimensions Full horiz Small height width
shapeA
   let lda :: Int
lda = Int
m
   let nrhs :: Int
nrhs = nrhs -> Int
forall sh. C sh => sh -> Int
Shape.size nrhs
widthB
   let ldb :: Int
ldb = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
   let ldx :: Int
ldx = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
widthA
   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
mPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
nrhsPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
nrhs
      (Ptr CChar
transPtr,Ptr a
aPtr) <- Order -> Int -> ForeignPtr a -> ContT () IO (Ptr CChar, Ptr a)
forall a r.
Floating a =>
Order -> Int -> ForeignPtr a -> ContT r IO (Ptr CChar, Ptr a)
adjointA Order
orderA (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n) ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
lda
      Ptr CInt
ldbPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
ldb
      Ptr a
bPtr <- Order -> Int -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
orderB Int
ldb Int
nrhs ForeignPtr a
b
      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
$ [Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
[Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo [Char]
rankMsg [Char]
"gels" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.gels Ptr CChar
transPtr
            Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
nrhsPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
bPtr Ptr CInt
ldbPtr
      IO () -> ContT () IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copySubMatrix Int
ldx Int
nrhs Int
ldb Ptr a
bPtr Int
ldx Ptr a
xPtr

minimumNorm ::
   (Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   Full Extent.Small vert height width a ->
   Full vert horiz height nrhs a ->
   Full vert horiz width nrhs a
minimumNorm :: Full Small vert height width a
-> Full vert horiz height nrhs a -> Full vert horiz width nrhs a
minimumNorm
   (Array shapeA :: Full Small vert height width
shapeA@(MatrixShape.Full Order
orderA Extent Small vert height width
extentA) ForeignPtr a
a)
   (Array        (MatrixShape.Full Order
orderB Extent vert horiz height nrhs
extentB) ForeignPtr a
b) =

 case Extent vert horiz width height
-> Extent vert horiz height nrhs
-> Maybe (Extent vert horiz width nrhs)
forall vert horiz fuse height width.
(C vert, C horiz, Eq fuse) =>
Extent vert horiz height fuse
-> Extent vert horiz fuse width
-> Maybe (Extent vert horiz height width)
Extent.fuse (Extent vert Small width height -> Extent vert horiz width height
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert Small height width -> Extent vert horiz height width
Extent.generalizeTall (Extent vert Small width height -> Extent vert horiz width height)
-> Extent vert Small width height -> Extent vert horiz width height
forall a b. (a -> b) -> a -> b
$ Extent Small vert height width -> Extent vert Small width height
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> Extent horiz vert width height
Extent.transpose Extent Small vert height width
extentA) Extent vert horiz height nrhs
extentB of
  Maybe (Extent vert horiz width nrhs)
Nothing -> [Char] -> Full vert horiz width nrhs a
forall a. HasCallStack => [Char] -> a
error [Char]
"minimumNorm: height shapes mismatch"
  Just Extent vert horiz width nrhs
extent ->
      Full vert horiz width nrhs
-> (Ptr a -> IO ()) -> Full vert horiz width nrhs a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> Extent vert horiz width nrhs -> Full vert horiz width nrhs
forall vert horiz height width.
Order
-> Extent vert horiz height width -> Full vert horiz height width
MatrixShape.Full Order
ColumnMajor Extent vert horiz width nrhs
extent) ((Ptr a -> IO ()) -> Full vert horiz width nrhs a)
-> (Ptr a -> IO ()) -> Full vert horiz width nrhs a
forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> do

   let widthA :: width
widthA = Extent Small vert height width -> width
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> width
Extent.width Extent Small vert height width
extentA
   let (height
height,nrhs
widthB) = Extent vert horiz height nrhs -> (height, nrhs)
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> (height, width)
Extent.dimensions Extent vert horiz height nrhs
extentB
   let (Int
m,Int
n) = Full Small vert height width -> (Int, Int)
forall vert horiz height width.
(C vert, C horiz, C height, C width) =>
Full vert horiz height width -> (Int, Int)
MatrixShape.dimensions Full Small vert height width
shapeA
   let lda :: Int
lda = Int
m
   let nrhs :: Int
nrhs = nrhs -> Int
forall sh. C sh => sh -> Int
Shape.size nrhs
widthB
   let ldb :: Int
ldb = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
   let ldx :: Int
ldx = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
widthA
   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
mPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
nrhsPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
nrhs
      (Ptr CChar
transPtr,Ptr a
aPtr) <- Order -> Int -> ForeignPtr a -> ContT () IO (Ptr CChar, Ptr a)
forall a r.
Floating a =>
Order -> Int -> ForeignPtr a -> ContT r IO (Ptr CChar, Ptr a)
adjointA Order
orderA (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n) ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
lda
      Ptr a
bPtr <- ((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
b
      Ptr CInt
ldxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
ldx
      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
$ Order -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Order -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copyToSubColumnMajor Order
orderB Int
ldb Int
nrhs Ptr a
bPtr Int
ldx Ptr a
xPtr
      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
$ [Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
[Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo [Char]
rankMsg [Char]
"gels" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.gels Ptr CChar
transPtr
            Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
nrhsPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
xPtr Ptr CInt
ldxPtr


adjointA ::
   Class.Floating a =>
   Order -> Int -> ForeignPtr a -> ContT r IO (Ptr CChar, Ptr a)
adjointA :: Order -> Int -> ForeignPtr a -> ContT r IO (Ptr CChar, Ptr a)
adjointA Order
order Int
size ForeignPtr a
a = do
   Ptr a
aPtr <- Int -> ForeignPtr a -> ContT r IO (Ptr a)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToTemp Int
size ForeignPtr a
a
   Char
trans <-
      case Order
order of
         Order
RowMajor -> do
            Ptr CInt
sizePtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
size
            Ptr CInt
incPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
            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
$ Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
lacgv Ptr CInt
sizePtr Ptr a
aPtr Ptr CInt
incPtr
            Char -> ContT r IO Char
forall (m :: * -> *) a. Monad m => a -> m a
return (Char -> ContT r IO Char) -> Char -> ContT r IO Char
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> Char
forall a (f :: * -> *). Floating a => f a -> Char
HH.invChar ForeignPtr a
a
         Order
ColumnMajor -> Char -> ContT r IO Char
forall (m :: * -> *) a. Monad m => a -> m a
return Char
'N'
   Ptr CChar
transPtr <- Char -> FortranIO r (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
trans
   (Ptr CChar, Ptr a) -> ContT r IO (Ptr CChar, Ptr a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Ptr CChar
transPtr, Ptr a
aPtr)


leastSquaresMinimumNormRCond ::
   (Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   RealOf a ->
   Full horiz vert height width a ->
   Full vert horiz height nrhs a ->
   (Int, Full vert horiz width nrhs a)
leastSquaresMinimumNormRCond :: RealOf a
-> Full horiz vert height width a
-> Full vert horiz height nrhs a
-> (Int, Full vert horiz width nrhs a)
leastSquaresMinimumNormRCond RealOf a
rcond
      (Array (MatrixShape.Full Order
orderA Extent horiz vert height width
extentA) ForeignPtr a
a)
      (Array (MatrixShape.Full Order
orderB Extent vert horiz height nrhs
extentB) ForeignPtr a
b) =
   case Extent vert horiz width height
-> Extent vert horiz height nrhs
-> Maybe (Extent vert horiz width nrhs)
forall vert horiz fuse height width.
(C vert, C horiz, Eq fuse) =>
Extent vert horiz height fuse
-> Extent vert horiz fuse width
-> Maybe (Extent vert horiz height width)
Extent.fuse (Extent horiz vert height width -> Extent vert horiz width height
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> Extent horiz vert width height
Extent.transpose Extent horiz vert height width
extentA) Extent vert horiz height nrhs
extentB of
      Maybe (Extent vert horiz width nrhs)
Nothing -> [Char] -> (Int, Full vert horiz width nrhs a)
forall a. HasCallStack => [Char] -> a
error [Char]
"leastSquaresMinimumNormRCond: height shapes mismatch"
      Just Extent vert horiz width nrhs
extent ->
         let widthA :: width
widthA = Extent horiz vert height width -> width
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> width
Extent.width Extent horiz vert height width
extentA
             (height
height,nrhs
widthB) = Extent vert horiz height nrhs -> (height, nrhs)
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> (height, width)
Extent.dimensions Extent vert horiz height nrhs
extentB
             shapeX :: Full vert horiz width nrhs
shapeX = Order -> Extent vert horiz width nrhs -> Full vert horiz width nrhs
forall vert horiz height width.
Order
-> Extent vert horiz height width -> Full vert horiz height width
MatrixShape.Full Order
ColumnMajor Extent vert horiz width nrhs
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
widthA
             nrhs :: Int
nrhs = nrhs -> Int
forall sh. C sh => sh -> Int
Shape.size nrhs
widthB
         in  if Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                then (Int
0, Full vert horiz width nrhs -> Full vert horiz width nrhs a
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.zero Full vert horiz width nrhs
shapeX)
                else
                  if Int
nrhs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                     then
                        ((Int, Array (General width ()) a) -> Int
forall a b. (a, b) -> a
fst ((Int, Array (General width ()) a) -> Int)
-> (Int, Array (General width ()) a) -> Int
forall a b. (a -> b) -> a -> b
$ IO (Int, Array (General width ()) a)
-> (Int, Array (General width ()) a)
forall a. IO a -> a
unsafePerformIO (IO (Int, Array (General width ()) a)
 -> (Int, Array (General width ()) a))
-> IO (Int, Array (General width ()) a)
-> (Int, Array (General width ()) a)
forall a b. (a -> b) -> a -> b
$
                         case height -> Vector height a
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.zero height
height of
                           Array height
_ ForeignPtr a
b1 ->
                              RealOf a
-> General width ()
-> Order
-> ForeignPtr a
-> Order
-> ForeignPtr a
-> Int
-> Int
-> Int
-> IO (Int, Array (General width ()) a)
forall sh a.
(C sh, Floating a) =>
RealOf a
-> sh
-> Order
-> ForeignPtr a
-> Order
-> ForeignPtr a
-> Int
-> Int
-> Int
-> IO (Int, Array sh a)
leastSquaresMinimumNormIO RealOf a
rcond
                                 (Order -> width -> () -> General width ()
forall height width.
Order -> height -> width -> General height width
MatrixShape.general Order
ColumnMajor width
widthA ())
                                 Order
orderA ForeignPtr a
a Order
orderB ForeignPtr a
b1 Int
m Int
n Int
1,
                         Full vert horiz width nrhs -> Full vert horiz width nrhs a
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.zero Full vert horiz width nrhs
shapeX)
                     else
                        IO (Int, Full vert horiz width nrhs a)
-> (Int, Full vert horiz width nrhs a)
forall a. IO a -> a
unsafePerformIO (IO (Int, Full vert horiz width nrhs a)
 -> (Int, Full vert horiz width nrhs a))
-> IO (Int, Full vert horiz width nrhs a)
-> (Int, Full vert horiz width nrhs a)
forall a b. (a -> b) -> a -> b
$
                        RealOf a
-> Full vert horiz width nrhs
-> Order
-> ForeignPtr a
-> Order
-> ForeignPtr a
-> Int
-> Int
-> Int
-> IO (Int, Full vert horiz width nrhs a)
forall sh a.
(C sh, Floating a) =>
RealOf a
-> sh
-> Order
-> ForeignPtr a
-> Order
-> ForeignPtr a
-> Int
-> Int
-> Int
-> IO (Int, Array sh a)
leastSquaresMinimumNormIO RealOf a
rcond Full vert horiz width nrhs
shapeX
                           Order
orderA ForeignPtr a
a Order
orderB ForeignPtr a
b Int
m Int
n Int
nrhs

leastSquaresMinimumNormIO ::
   (Shape.C sh, Class.Floating a) =>
   RealOf a -> sh ->
   Order -> ForeignPtr a ->
   Order -> ForeignPtr a ->
   Int -> Int -> Int -> IO (Int, Array sh a)
leastSquaresMinimumNormIO :: RealOf a
-> sh
-> Order
-> ForeignPtr a
-> Order
-> ForeignPtr a
-> Int
-> Int
-> Int
-> IO (Int, Array sh a)
leastSquaresMinimumNormIO RealOf a
rcond sh
shapeX Order
orderA ForeignPtr a
a Order
orderB ForeignPtr a
b Int
m Int
n Int
nrhs =
   sh
-> Int
-> Int
-> Int
-> ((Ptr a, Int) -> IO Int)
-> IO (Int, Array sh a)
forall sh a rank.
(C sh, Floating a) =>
sh
-> Int
-> Int
-> Int
-> ((Ptr a, Int) -> IO rank)
-> IO (rank, Array sh a)
createHigherArray sh
shapeX Int
m Int
n Int
nrhs (((Ptr a, Int) -> IO Int) -> IO (Int, Array sh a))
-> ((Ptr a, Int) -> IO Int) -> IO (Int, Array sh a)
forall a b. (a -> b) -> a -> b
$ \(Ptr a
tmpPtr,Int
ldtmp) -> do

   let lda :: Int
lda = Int
m
   ContT Int IO Int -> IO Int
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT Int IO Int -> IO Int) -> ContT Int IO Int -> IO Int
forall a b. (a -> b) -> a -> b
$ do
      Ptr a
aPtr <- Order -> Int -> Int -> ForeignPtr a -> ContT Int IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
orderA Int
m Int
n ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO Int (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
lda
      Ptr CInt
ldtmpPtr <- Int -> FortranIO Int (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
ldtmp
      Ptr a
bPtr <- ((Ptr a -> IO Int) -> IO Int) -> ContT Int IO (Ptr a)
forall k (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO Int) -> IO Int) -> ContT Int IO (Ptr a))
-> ((Ptr a -> IO Int) -> IO Int) -> ContT Int IO (Ptr a)
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> (Ptr a -> IO Int) -> IO Int
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
b
      IO () -> ContT Int IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT Int IO ()) -> IO () -> ContT Int IO ()
forall a b. (a -> b) -> a -> b
$ Order -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Order -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copyToSubColumnMajor Order
orderB Int
m Int
nrhs Ptr a
bPtr Int
ldtmp Ptr a
tmpPtr
      Ptr CInt
jpvtPtr <- Int -> FortranIO Int (Ptr CInt)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      IO () -> ContT Int IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT Int IO ()) -> IO () -> ContT Int IO ()
forall a b. (a -> b) -> a -> b
$ Ptr CInt -> [CInt] -> IO ()
forall a. Storable a => Ptr a -> [a] -> IO ()
pokeArray Ptr CInt
jpvtPtr (Int -> CInt -> [CInt]
forall a. Int -> a -> [a]
replicate Int
n CInt
0)
      Ptr CInt
rankPtr <- FortranIO Int (Ptr CInt)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      GELSY_ Int (RealOf a) a
forall a r. Floating a => GELSY_ r (RealOf a) a
gelsy Int
m Int
n Int
nrhs Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
tmpPtr Ptr CInt
ldtmpPtr Ptr CInt
jpvtPtr RealOf a
rcond Ptr CInt
rankPtr
      IO Int -> ContT Int IO Int
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Int -> ContT Int IO Int) -> IO Int -> ContT Int IO Int
forall a b. (a -> b) -> a -> b
$ Ptr CInt -> IO Int
peekCInt Ptr CInt
rankPtr


type GELSY_ r ar a =
   Int -> Int -> Int -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->
   Ptr CInt -> ar -> Ptr CInt -> ContT r IO ()

newtype GELSY r a = GELSY {GELSY r a -> GELSY_ r (RealOf a) a
getGELSY :: GELSY_ r (RealOf a) a}

gelsy :: (Class.Floating a) => GELSY_ r (RealOf a) a
gelsy :: GELSY_ r (RealOf a) a
gelsy =
   GELSY r a -> GELSY_ r (RealOf a) a
forall r a. GELSY r a -> GELSY_ r (RealOf a) a
getGELSY (GELSY r a -> GELSY_ r (RealOf a) a)
-> GELSY r a -> GELSY_ r (RealOf a) a
forall a b. (a -> b) -> a -> b
$
   GELSY r Float
-> GELSY r Double
-> GELSY r (Complex Float)
-> GELSY r (Complex Double)
-> GELSY r a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (GELSY_ r (RealOf Float) Float -> GELSY r Float
forall r a. GELSY_ r (RealOf a) a -> GELSY r a
GELSY GELSY_ r (RealOf Float) Float
forall a r. Real a => GELSY_ r a a
gelsyReal)
      (GELSY_ r (RealOf Double) Double -> GELSY r Double
forall r a. GELSY_ r (RealOf a) a -> GELSY r a
GELSY GELSY_ r (RealOf Double) Double
forall a r. Real a => GELSY_ r a a
gelsyReal)
      (GELSY_ r (RealOf (Complex Float)) (Complex Float)
-> GELSY r (Complex Float)
forall r a. GELSY_ r (RealOf a) a -> GELSY r a
GELSY GELSY_ r (RealOf (Complex Float)) (Complex Float)
forall a r. Real a => GELSY_ r a (Complex a)
gelsyComplex)
      (GELSY_ r (RealOf (Complex Double)) (Complex Double)
-> GELSY r (Complex Double)
forall r a. GELSY_ r (RealOf a) a -> GELSY r a
GELSY GELSY_ r (RealOf (Complex Double)) (Complex Double)
forall a r. Real a => GELSY_ r a (Complex a)
gelsyComplex)

gelsyReal :: (Class.Real a) => GELSY_ r a a
gelsyReal :: GELSY_ r a a
gelsyReal Int
m Int
n Int
nrhs Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
bPtr Ptr CInt
ldbPtr Ptr CInt
jpvtPtr a
rcond Ptr CInt
rankPtr = do
   Ptr CInt
mPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
   Ptr CInt
nPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
   Ptr CInt
nrhsPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
nrhs
   Ptr a
rcondPtr <- a -> FortranIO r (Ptr a)
forall a r. Real a => a -> FortranIO r (Ptr a)
Call.real a
rcond
   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
$ [Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
[Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo [Char]
errorCodeMsg [Char]
"gelsy" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
      Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackReal.gelsy Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
nrhsPtr
         Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
bPtr Ptr CInt
ldbPtr Ptr CInt
jpvtPtr Ptr a
rcondPtr Ptr CInt
rankPtr

gelsyComplex :: (Class.Real a) => GELSY_ r a (Complex a)
gelsyComplex :: GELSY_ r a (Complex a)
gelsyComplex Int
m Int
n Int
nrhs Ptr (Complex a)
aPtr Ptr CInt
ldaPtr Ptr (Complex a)
bPtr Ptr CInt
ldbPtr Ptr CInt
jpvtPtr a
rcond Ptr CInt
rankPtr = do
   Ptr CInt
mPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
   Ptr CInt
nPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
   Ptr CInt
nrhsPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
nrhs
   Ptr a
rcondPtr <- a -> FortranIO r (Ptr a)
forall a r. Real a => a -> FortranIO r (Ptr a)
Call.real a
rcond
   Ptr a
rworkPtr <- Int -> FortranIO r (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray (Int
2Int -> 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
$
      [Char]
-> [Char]
-> (Ptr (Complex a) -> Ptr CInt -> Ptr CInt -> IO ())
-> IO ()
forall a.
Floating a =>
[Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo [Char]
errorCodeMsg [Char]
"gelsy" ((Ptr (Complex a) -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr (Complex a) -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr (Complex a)
workPtr Ptr CInt
lworkPtr Ptr CInt
infoPtr ->
      Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> IO ()
LapackComplex.gelsy Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
nrhsPtr
         Ptr (Complex a)
aPtr Ptr CInt
ldaPtr Ptr (Complex a)
bPtr Ptr CInt
ldbPtr Ptr CInt
jpvtPtr Ptr a
rcondPtr Ptr CInt
rankPtr
         Ptr (Complex a)
workPtr Ptr CInt
lworkPtr Ptr a
rworkPtr Ptr CInt
infoPtr


pseudoInverseRCond ::
   (Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
   RealOf a ->
   Full vert horiz height width a ->
   (Int, Full horiz vert width height a)
pseudoInverseRCond :: RealOf a
-> Full vert horiz height width a
-> (Int, Full horiz vert width height a)
pseudoInverseRCond RealOf a
rcond Full vert horiz height width a
a =
   case Full vert horiz height width a
-> Either (Tall height width a) (Wide height width a)
forall vert horiz height width a.
(C vert, C horiz, C height, C width) =>
Full vert horiz height width a
-> Either (Tall height width a) (Wide height width a)
Basic.caseTallWide Full vert horiz height width a
a of
      Left Tall height width a
_ ->
         (Full vert horiz height width a -> Full horiz vert width height a)
-> (Int, Full vert horiz height width a)
-> (Int, Full horiz vert width height a)
forall b c a. (b -> c) -> (a, b) -> (a, c)
mapSnd Full vert horiz height width a -> Full horiz vert width height a
forall vert horiz height width a.
(C vert, C horiz) =>
Full vert horiz height width a -> Full horiz vert width height a
Basic.transpose ((Int, Full vert horiz height width a)
 -> (Int, Full horiz vert width height a))
-> (Int, Full vert horiz height width a)
-> (Int, Full horiz vert width height a)
forall a b. (a -> b) -> a -> b
$
         RealOf a
-> Full horiz vert width height a
-> Full vert horiz width width a
-> (Int, Full vert horiz height width a)
forall vert horiz height width nrhs a.
(C vert, C horiz, C height, Eq height, C width, C nrhs,
 Floating a) =>
RealOf a
-> Full horiz vert height width a
-> Full vert horiz height nrhs a
-> (Int, Full vert horiz width nrhs a)
leastSquaresMinimumNormRCond RealOf a
rcond (Full vert horiz height width a -> Full horiz vert width height a
forall vert horiz height width a.
(C vert, C horiz) =>
Full vert horiz height width a -> Full horiz vert width height a
Basic.transpose Full vert horiz height width a
a) (Full vert horiz width width a
 -> (Int, Full vert horiz height width a))
-> Full vert horiz width width a
-> (Int, Full vert horiz height width a)
forall a b. (a -> b) -> a -> b
$
         Square width a -> Full vert horiz width width a
forall vert horiz sh a.
(C vert, C horiz) =>
Square sh a -> Full vert horiz sh sh a
Square.toFull (Square width a -> Full vert horiz width width a)
-> Square width a -> Full vert horiz width width a
forall a b. (a -> b) -> a -> b
$ width -> Square width a
forall sh a. (C sh, Floating a) => sh -> Square sh a
Square.identity (width -> Square width a) -> width -> Square width a
forall a b. (a -> b) -> a -> b
$
         Full vert horiz height width -> width
forall vert horiz height width.
(C vert, C horiz) =>
Full vert horiz height width -> width
MatrixShape.fullWidth (Full vert horiz height width -> width)
-> Full vert horiz height width -> width
forall a b. (a -> b) -> a -> b
$ Full vert horiz height width a -> Full vert horiz height width
forall sh a. Array sh a -> sh
Array.shape Full vert horiz height width a
a
      Right Wide height width a
_ ->
         RealOf a
-> Full vert horiz height width a
-> Full horiz vert height height a
-> (Int, Full horiz vert width height a)
forall vert horiz height width nrhs a.
(C vert, C horiz, C height, Eq height, C width, C nrhs,
 Floating a) =>
RealOf a
-> Full horiz vert height width a
-> Full vert horiz height nrhs a
-> (Int, Full vert horiz width nrhs a)
leastSquaresMinimumNormRCond RealOf a
rcond Full vert horiz height width a
a (Full horiz vert height height a
 -> (Int, Full horiz vert width height a))
-> Full horiz vert height height a
-> (Int, Full horiz vert width height a)
forall a b. (a -> b) -> a -> b
$
         Square height a -> Full horiz vert height height a
forall vert horiz sh a.
(C vert, C horiz) =>
Square sh a -> Full vert horiz sh sh a
Square.toFull (Square height a -> Full horiz vert height height a)
-> Square height a -> Full horiz vert height height a
forall a b. (a -> b) -> a -> b
$ height -> Square height a
forall sh a. (C sh, Floating a) => sh -> Square sh a
Square.identity (height -> Square height a) -> height -> Square height a
forall a b. (a -> b) -> a -> b
$
         Full vert horiz height width -> height
forall vert horiz height width.
(C vert, C horiz) =>
Full vert horiz height width -> height
MatrixShape.fullHeight (Full vert horiz height width -> height)
-> Full vert horiz height width -> height
forall a b. (a -> b) -> a -> b
$ Full vert horiz height width a -> Full vert horiz height width
forall sh a. Array sh a -> sh
Array.shape Full vert horiz height width a
a


leastSquaresConstraint ::
   (Shape.C height, Eq height,
    Shape.C width, Eq width,
    Shape.C constraints, Eq constraints, Class.Floating a) =>
   General height width a -> Vector height a ->
   Wide constraints width a -> Vector constraints a ->
   Vector width a
leastSquaresConstraint :: General height width a
-> Vector height a
-> Wide constraints width a
-> Vector constraints a
-> Vector width a
leastSquaresConstraint
   (Array (MatrixShape.Full Order
orderA Extent Big Big height width
extentA) ForeignPtr a
a) Vector height a
c
   (Array (MatrixShape.Full Order
orderB Extent Small Big constraints width
extentB) ForeignPtr a
b) Vector constraints a
d =

 let sameShape :: [Char] -> p -> p -> p
sameShape [Char]
name p
shape0 p
shape1 =
      if p
shape0 p -> p -> Bool
forall a. Eq a => a -> a -> Bool
== p
shape1
         then p
shape0
         else [Char] -> p
forall a. HasCallStack => [Char] -> a
error ([Char] -> p) -> [Char] -> p
forall a b. (a -> b) -> a -> b
$ [Char]
"leastSquaresConstraint: " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
name [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
" shapes mismatch"
     width :: width
width = [Char] -> width -> width -> width
forall p. Eq p => [Char] -> p -> p -> p
sameShape [Char]
"width" (Extent Big Big height width -> width
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> width
Extent.width Extent Big Big height width
extentA) (Extent Small Big constraints width -> width
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> width
Extent.width Extent Small Big constraints width
extentB)
 in
   width -> (Ptr a -> IO ()) -> Vector width a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate width
width ((Ptr a -> IO ()) -> Vector width a)
-> (Ptr a -> IO ()) -> Vector width a
forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> do

   let height :: height
height = [Char] -> height -> height -> height
forall p. Eq p => [Char] -> p -> p -> p
sameShape [Char]
"height" (Extent Big Big height width -> height
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> height
Extent.height Extent Big Big height width
extentA) (Vector height a -> height
forall sh a. Array sh a -> sh
Array.shape Vector height a
c)
   let constraints :: constraints
constraints =
         [Char] -> constraints -> constraints -> constraints
forall p. Eq p => [Char] -> p -> p -> p
sameShape [Char]
"constraints" (Extent Small Big constraints width -> constraints
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> height
Extent.height Extent Small Big constraints width
extentB) (Vector constraints a -> constraints
forall sh a. Array sh a -> sh
Array.shape Vector constraints a
d)
   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
   let p :: Int
p = constraints -> Int
forall sh. C sh => sh -> Int
Shape.size constraints
constraints
   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
mPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
pPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
p
      Ptr a
aPtr <- Order -> Int -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
orderA Int
m Int
n ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
m
      Ptr a
bPtr <- Order -> Int -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
orderB Int
p Int
n ForeignPtr a
b
      Ptr CInt
ldbPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
p
      Ptr a
cPtr <- Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToTemp Int
m (Vector height a -> ForeignPtr a
forall sh a. Array sh a -> ForeignPtr a
Array.buffer Vector height a
c)
      Ptr a
dPtr <- Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToTemp Int
p (Vector constraints a -> ForeignPtr a
forall sh a. Array sh a -> ForeignPtr a
Array.buffer Vector constraints a
d)
      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
$ [Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
[Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo [Char]
rankMsg [Char]
"gglse" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.gglse
            Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
pPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
bPtr Ptr CInt
ldbPtr Ptr a
cPtr Ptr a
dPtr Ptr a
xPtr

gaussMarkovLinearModel ::
   (Shape.C height, Eq height,
    Shape.C width, Eq width,
    Shape.C opt, Eq opt, Class.Floating a) =>
   Tall height width a -> General height opt a -> Vector height a ->
   (Vector width a, Vector opt a)
gaussMarkovLinearModel :: Tall height width a
-> General height opt a
-> Vector height a
-> (Vector width a, Vector opt a)
gaussMarkovLinearModel
   (Array (MatrixShape.Full Order
orderA Extent Big Small height width
extentA) ForeignPtr a
a)
   (Array (MatrixShape.Full Order
orderB Extent Big Big height opt
extentB) ForeignPtr a
b) Vector height a
d =

   let width :: width
width = Extent Big Small height width -> width
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> width
Extent.width Extent Big Small height width
extentA in
   let opt :: opt
opt = Extent Big Big height opt -> opt
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> width
Extent.width Extent Big Big height opt
extentB in
   width
-> (Int -> Ptr a -> IO (Vector opt a))
-> (Vector width a, Vector opt a)
forall sh a b.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO b) -> (Array sh a, b)
Array.unsafeCreateWithSizeAndResult width
width ((Int -> Ptr a -> IO (Vector opt a))
 -> (Vector width a, Vector opt a))
-> (Int -> Ptr a -> IO (Vector opt a))
-> (Vector width a, Vector opt a)
forall a b. (a -> b) -> a -> b
$ \Int
m Ptr a
xPtr -> do
   opt -> (Int -> Ptr a -> IO ()) -> IO (Vector opt a)
forall (m :: * -> *) sh a.
(PrimMonad m, C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> m (Array sh a)
ArrayIO.unsafeCreateWithSize opt
opt ((Int -> Ptr a -> IO ()) -> IO (Vector opt a))
-> (Int -> Ptr a -> IO ()) -> IO (Vector opt a)
forall a b. (a -> b) -> a -> b
$ \Int
p Ptr a
yPtr -> do

   let sameHeight :: p -> p -> p
sameHeight p
shape0 p
shape1 =
         if p
shape0 p -> p -> Bool
forall a. Eq a => a -> a -> Bool
== p
shape1
            then p
shape0
            else [Char] -> p
forall a. HasCallStack => [Char] -> a
error ([Char] -> p) -> [Char] -> p
forall a b. (a -> b) -> a -> b
$ [Char]
"gaussMarkovLinearModel: height shapes mismatch"
       height :: height
height =
         height -> height -> height
forall p. Eq p => p -> p -> p
sameHeight (Vector height a -> height
forall sh a. Array sh a -> sh
Array.shape Vector height a
d) (height -> height) -> height -> height
forall a b. (a -> b) -> a -> b
$
         height -> height -> height
forall p. Eq p => p -> p -> p
sameHeight (Extent Big Small height width -> height
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> height
Extent.height Extent Big Small height width
extentA) (Extent Big Big height opt -> height
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> height
Extent.height Extent Big Big height opt
extentB)
   let n :: Int
n = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
   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
mPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
pPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
p
      Ptr a
aPtr <- Order -> Int -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
orderA Int
n Int
m ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr a
bPtr <- Order -> Int -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
orderB Int
n Int
p ForeignPtr a
b
      Ptr CInt
ldbPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr a
dPtr <- Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToTemp Int
n (Vector height a -> ForeignPtr a
forall sh a. Array sh a -> ForeignPtr a
Array.buffer Vector height a
d)
      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
$ [Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
[Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo [Char]
rankMsg [Char]
"ggglm" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.ggglm
            Ptr CInt
nPtr Ptr CInt
mPtr Ptr CInt
pPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
bPtr Ptr CInt
ldbPtr Ptr a
dPtr Ptr a
xPtr Ptr a
yPtr


determinantAbsolute ::
   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,
    Class.Floating a) =>
   Full vert horiz height width a -> RealOf a
determinantAbsolute :: Full vert horiz height width a -> RealOf a
determinantAbsolute =
   a -> RealOf a
forall a. Floating a => a -> RealOf a
absolute (a -> RealOf a)
-> (Full vert horiz height width a -> a)
-> Full vert horiz height width a
-> RealOf a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (Full Big Small height width a -> a)
-> (Wide height width a -> a)
-> Either (Full Big Small height width a) (Wide height width a)
-> a
forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either (Householder Big Small height width a -> a
forall vert height width a.
(C vert, C height, C width, Floating a) =>
Householder vert Small height width a -> a
HH.determinantR (Householder Big Small height width a -> a)
-> (Full Big Small height width a
    -> Householder Big Small height width a)
-> Full Big Small height width a
-> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Full Big Small height width a
-> Householder Big Small height width a
forall vert horiz height width a.
(C vert, C horiz, C height, C width, Floating a) =>
Full vert horiz height width a
-> Householder vert horiz height width a
HH.fromMatrix) (a -> Wide height width a -> a
forall a b. a -> b -> a
const a
forall a. Floating a => a
zero) (Either (Full Big Small height width a) (Wide height width a) -> a)
-> (Full vert horiz height width a
    -> Either (Full Big Small height width a) (Wide height width a))
-> Full vert horiz height width a
-> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Full vert horiz height width a
-> Either (Full Big Small height width a) (Wide height width a)
forall vert horiz height width a.
(C vert, C horiz, C height, C width) =>
Full vert horiz height width a
-> Either (Tall height width a) (Wide height width a)
Basic.caseTallWide


complement ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   Tall height width a -> Tall height ShapeInt a
complement :: Tall height width a -> Tall height ShapeInt a
complement = Tall height width a -> Tall height ShapeInt a
forall height width a.
(C height, C width, Floating a) =>
Tall height width a -> Tall height ShapeInt a
extractComplement (Tall height width a -> Tall height ShapeInt a)
-> (Tall height width a -> Tall height width a)
-> Tall height width a
-> Tall height ShapeInt a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Tall height width a -> Tall height width a
forall vert horiz height width a.
(C vert, C horiz, C height, C width, Floating a) =>
Full vert horiz height width a
-> Householder vert horiz height width a
HH.fromMatrix

extractComplement ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   HH.Tall height width a -> Tall height ShapeInt a
extractComplement :: Tall height width a -> Tall height ShapeInt a
extractComplement Tall height width a
qr =
   Int -> Tall height ShapeInt a -> Tall height ShapeInt a
forall horiz height a.
(C horiz, C height, Floating a) =>
Int
-> Full Big horiz height ShapeInt a
-> Full Big horiz height ShapeInt a
Basic.dropColumns
      (width -> Int
forall sh. C sh => sh -> Int
Shape.size (width -> Int) -> width -> Int
forall a b. (a -> b) -> a -> b
$ Split Reflector Big Small height width -> width
forall vert horiz lower height width.
(C vert, C horiz) =>
Split lower vert horiz height width -> width
MatrixShape.splitWidth (Split Reflector Big Small height width -> width)
-> Split Reflector Big Small height width -> width
forall a b. (a -> b) -> a -> b
$ Array (Split Reflector Big Small height width) a
-> Split Reflector Big Small height width
forall sh a. Array sh a -> sh
Array.shape (Array (Split Reflector Big Small height width) a
 -> Split Reflector Big Small height width)
-> Array (Split Reflector Big Small height width) a
-> Split Reflector Big Small height width
forall a b. (a -> b) -> a -> b
$ Tall height width a
-> Array (Split Reflector Big Small height width) a
forall vert horiz height width a.
Matrix (Hh vert horiz height width) a
-> Array (Split Reflector vert horiz height width) a
HH.split_ Tall height width a
qr) (Tall height ShapeInt a -> Tall height ShapeInt a)
-> Tall height ShapeInt a -> Tall height ShapeInt a
forall a b. (a -> b) -> a -> b
$
   (height -> ShapeInt)
-> Full Big Small height height a -> Tall height ShapeInt 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 (Int -> ShapeInt
shapeInt (Int -> ShapeInt) -> (height -> Int) -> height -> ShapeInt
forall b c a. (b -> c) -> (a -> b) -> a -> c
. height -> Int
forall sh. C sh => sh -> Int
Shape.size) (Full Big Small height height a -> Tall height ShapeInt a)
-> Full Big Small height height a -> Tall height ShapeInt a
forall a b. (a -> b) -> a -> b
$ Square height a -> Full Big Small height height a
forall vert horiz sh a.
(C vert, C horiz) =>
Square sh a -> Full vert horiz sh sh a
Square.toFull (Square height a -> Full Big Small height height a)
-> Square height a -> Full Big Small height height a
forall a b. (a -> b) -> a -> b
$
   Tall height width a -> Square height a
forall vert horiz height width a.
(C vert, C horiz, C height, C width, Floating a) =>
Householder vert horiz height width a -> Square height a
HH.extractQ Tall height width a
qr