{-# LANGUAGE TypeFamilies #-}
module Numeric.LAPACK.Linear.LowerUpper (
   LowerUpper,
   Plain.Square,
   Plain.Tall,
   Plain.Wide,
   Plain.Transposition(..),
   Plain.Conjugation(..),
   Plain.Inversion(..),
   Plain.mapExtent,
   fromMatrix,
   toMatrix,
   solve,
   multiplyFull,

   Plain.determinant,

   extractP,
   multiplyP,

   extractL,
   wideExtractL,
   wideMultiplyL,
   wideSolveL,

   extractU,
   tallExtractU,
   tallMultiplyU,
   tallSolveU,

   Plain.caseTallWide,
   ) where

import qualified Numeric.LAPACK.Linear.Plain as Plain
import Numeric.LAPACK.Linear.Plain (LowerUpper)

import qualified Numeric.LAPACK.Matrix.Array.Triangular as Tri
import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Matrix.Permutation as PermMatrix
import qualified Numeric.LAPACK.Matrix as Matrix
import Numeric.LAPACK.Matrix.Array (Full)
import Numeric.LAPACK.Matrix.Modifier (Transposition, Conjugation, Inversion)

import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Shape as Shape


{- |
@LowerUpper.fromMatrix a@
computes the LU decomposition of matrix @a@ with row pivotisation.

You can reconstruct @a@ from @lu@ depending on whether @a@ is tall or wide.

> LU.multiplyP NonInverted lu $ LU.extractL lu ##*# LU.tallExtractU lu
> LU.multiplyP NonInverted lu $ LU.wideExtractL lu #*## LU.extractU lu
-}
fromMatrix ::
   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,
    Class.Floating a) =>
   Full vert horiz height width a ->
   LowerUpper vert horiz height width a
fromMatrix :: Full vert horiz height width a
-> LowerUpper vert horiz height width a
fromMatrix = Full vert horiz height width a
-> LowerUpper vert horiz 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
-> LowerUpper vert horiz height width a
Plain.fromMatrix (Full vert horiz height width a
 -> LowerUpper vert horiz height width a)
-> (Full vert horiz height width a
    -> Full vert horiz height width a)
-> Full vert horiz height width a
-> LowerUpper vert horiz height width a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Full vert horiz height width a -> Full vert horiz height width a
forall sh a. ArrayMatrix sh a -> Array sh a
ArrMatrix.toVector

solve ::
   (Extent.C vert, Extent.C horiz, Eq height, Shape.C height, Shape.C width,
    Class.Floating a) =>
   Plain.Square height a ->
   Full vert horiz height width a ->
   Full vert horiz height width a
solve :: Square height a
-> Full vert horiz height width a -> Full vert horiz height width a
solve = (Array (Full vert horiz height width) a
 -> Array (Full vert horiz height width) a)
-> Full vert horiz height width a -> Full vert horiz height width a
forall shA a shB b.
(Array shA a -> Array shB b)
-> ArrayMatrix shA a -> ArrayMatrix shB b
ArrMatrix.lift1 ((Array (Full vert horiz height width) a
  -> Array (Full vert horiz height width) a)
 -> Full vert horiz height width a
 -> Full vert horiz height width a)
-> (Square height a
    -> Array (Full vert horiz height width) a
    -> Array (Full vert horiz height width) a)
-> Square height a
-> Full vert horiz height width a
-> Full vert horiz height width a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Square height a
-> Array (Full vert horiz height width) a
-> Array (Full vert horiz height width) a
forall vert horiz height width a.
(C vert, C horiz, Eq height, C height, C width, Floating a) =>
Square height a
-> Full vert horiz height width a -> Full vert horiz height width a
Plain.solve


extractP ::
   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>
   Inversion -> LowerUpper vert horiz height width a ->
   Matrix.Permutation height a
extractP :: Inversion
-> LowerUpper vert horiz height width a -> Permutation height a
extractP Inversion
inverted = Permutation height -> Permutation height a
forall sh a. C sh => Permutation sh -> Matrix (Permutation sh) a
PermMatrix.fromPermutation (Permutation height -> Permutation height a)
-> (LowerUpper vert horiz height width a -> Permutation height)
-> LowerUpper vert horiz height width a
-> Permutation height a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Inversion
-> LowerUpper vert horiz height width a -> Permutation height
forall vert horiz height width a.
(C vert, C horiz, C height, C width) =>
Inversion
-> LowerUpper vert horiz height width a -> Permutation height
Plain.extractP Inversion
inverted

multiplyP ::
   (Extent.C vertA, Extent.C horizA, Extent.C vertB, Extent.C horizB,
    Eq height, Shape.C height, Shape.C widthA, Shape.C widthB,
    Class.Floating a) =>
   Inversion ->
   LowerUpper vertA horizA height widthA a ->
   Full vertB horizB height widthB a ->
   Full vertB horizB height widthB a
multiplyP :: Inversion
-> LowerUpper vertA horizA height widthA a
-> Full vertB horizB height widthB a
-> Full vertB horizB height widthB a
multiplyP Inversion
inverted = (Array (Full vertB horizB height widthB) a
 -> Array (Full vertB horizB height widthB) a)
-> Full vertB horizB height widthB a
-> Full vertB horizB height widthB a
forall shA a shB b.
(Array shA a -> Array shB b)
-> ArrayMatrix shA a -> ArrayMatrix shB b
ArrMatrix.lift1 ((Array (Full vertB horizB height widthB) a
  -> Array (Full vertB horizB height widthB) a)
 -> Full vertB horizB height widthB a
 -> Full vertB horizB height widthB a)
-> (LowerUpper vertA horizA height widthA a
    -> Array (Full vertB horizB height widthB) a
    -> Array (Full vertB horizB height widthB) a)
-> LowerUpper vertA horizA height widthA a
-> Full vertB horizB height widthB a
-> Full vertB horizB height widthB a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Inversion
-> LowerUpper vertA horizA height widthA a
-> Array (Full vertB horizB height widthB) a
-> Array (Full vertB horizB height widthB) a
forall vertA horizA vertB horizB height widthA widthB a.
(C vertA, C horizA, C vertB, C horizB, Eq height, C height,
 C widthA, C widthB, Floating a) =>
Inversion
-> LowerUpper vertA horizA height widthA a
-> Full vertB horizB height widthB a
-> Full vertB horizB height widthB a
Plain.multiplyP Inversion
inverted



extractL ::
   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,
    Class.Floating a) =>
   LowerUpper vert horiz height width a ->
   Full vert horiz height width a
extractL :: LowerUpper vert horiz height width a
-> Full vert horiz height width a
extractL = Array (Full vert horiz height width) a
-> Full vert horiz height width a
forall shA a. Array shA a -> ArrayMatrix shA a
ArrMatrix.lift0 (Array (Full vert horiz height width) a
 -> Full vert horiz height width a)
-> (LowerUpper vert horiz height width a
    -> Array (Full vert horiz height width) a)
-> LowerUpper vert horiz height width a
-> Full vert horiz height width a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. LowerUpper vert horiz height width a
-> Array (Full vert horiz height width) a
forall vert horiz height width a.
(C vert, C horiz, C height, C width, Floating a) =>
LowerUpper vert horiz height width a
-> Full vert horiz height width a
Plain.extractL

wideExtractL ::
   (Extent.C horiz, Shape.C height, Shape.C width, Class.Floating a) =>
   LowerUpper Extent.Small horiz height width a -> Tri.UnitLower height a
wideExtractL :: LowerUpper Small horiz height width a -> UnitLower height a
wideExtractL = Array (LowerTriangular Unit height) a -> UnitLower height a
forall shA a. Array shA a -> ArrayMatrix shA a
ArrMatrix.lift0 (Array (LowerTriangular Unit height) a -> UnitLower height a)
-> (LowerUpper Small horiz height width a
    -> Array (LowerTriangular Unit height) a)
-> LowerUpper Small horiz height width a
-> UnitLower height a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. LowerUpper Small horiz height width a
-> Array (LowerTriangular Unit height) a
forall horiz height width a.
(C horiz, C height, C width, Floating a) =>
LowerUpper Small horiz height width a -> UnitLower height a
Plain.wideExtractL

{- |
@wideMultiplyL transposed lu a@ multiplies the square part of @lu@
containing the lower triangular matrix with @a@.

> wideMultiplyL NonTransposed lu a == wideExtractL lu #*## a
> wideMultiplyL Transposed lu a == Tri.transpose (wideExtractL lu) #*## a
-}
wideMultiplyL ::
   (Extent.C horizA, Extent.C vert, Extent.C horiz, Shape.C height, Eq height,
    Shape.C widthA, Shape.C widthB, Class.Floating a) =>
   Transposition ->
   LowerUpper Extent.Small horizA height widthA a ->
   Full vert horiz height widthB a ->
   Full vert horiz height widthB a
wideMultiplyL :: Transposition
-> LowerUpper Small horizA height widthA a
-> Full vert horiz height widthB a
-> Full vert horiz height widthB a
wideMultiplyL Transposition
transposed = (Array (Full vert horiz height widthB) a
 -> Array (Full vert horiz height widthB) a)
-> Full vert horiz height widthB a
-> Full vert horiz height widthB a
forall shA a shB b.
(Array shA a -> Array shB b)
-> ArrayMatrix shA a -> ArrayMatrix shB b
ArrMatrix.lift1 ((Array (Full vert horiz height widthB) a
  -> Array (Full vert horiz height widthB) a)
 -> Full vert horiz height widthB a
 -> Full vert horiz height widthB a)
-> (LowerUpper Small horizA height widthA a
    -> Array (Full vert horiz height widthB) a
    -> Array (Full vert horiz height widthB) a)
-> LowerUpper Small horizA height widthA a
-> Full vert horiz height widthB a
-> Full vert horiz height widthB a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Transposition
-> LowerUpper Small horizA height widthA a
-> Array (Full vert horiz height widthB) a
-> Array (Full vert horiz height widthB) a
forall horizA vert horiz height widthA widthB a.
(C horizA, C vert, C horiz, C height, Eq height, C widthA,
 C widthB, Floating a) =>
Transposition
-> LowerUpper Small horizA height widthA a
-> Full vert horiz height widthB a
-> Full vert horiz height widthB a
Plain.wideMultiplyL Transposition
transposed

wideSolveL ::
   (Extent.C horizA, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   Transposition -> Conjugation ->
   LowerUpper Extent.Small horizA height width a ->
   Full vert horiz height nrhs a -> Full vert horiz height nrhs a
wideSolveL :: Transposition
-> Conjugation
-> LowerUpper Small horizA height width a
-> Full vert horiz height nrhs a
-> Full vert horiz height nrhs a
wideSolveL Transposition
transposed Conjugation
conjugated =
   (Array (Full vert horiz height nrhs) a
 -> Array (Full vert horiz height nrhs) a)
-> Full vert horiz height nrhs a -> Full vert horiz height nrhs a
forall shA a shB b.
(Array shA a -> Array shB b)
-> ArrayMatrix shA a -> ArrayMatrix shB b
ArrMatrix.lift1 ((Array (Full vert horiz height nrhs) a
  -> Array (Full vert horiz height nrhs) a)
 -> Full vert horiz height nrhs a -> Full vert horiz height nrhs a)
-> (LowerUpper Small horizA height width a
    -> Array (Full vert horiz height nrhs) a
    -> Array (Full vert horiz height nrhs) a)
-> LowerUpper Small horizA height width a
-> Full vert horiz height nrhs a
-> Full vert horiz height nrhs a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Transposition
-> Conjugation
-> LowerUpper Small horizA height width a
-> Array (Full vert horiz height nrhs) a
-> Array (Full vert horiz height nrhs) a
forall horizA vert horiz height width nrhs a.
(C horizA, C vert, C horiz, C height, Eq height, C width, C nrhs,
 Floating a) =>
Transposition
-> Conjugation
-> LowerUpper Small horizA height width a
-> Full vert horiz height nrhs a
-> Full vert horiz height nrhs a
Plain.wideSolveL Transposition
transposed Conjugation
conjugated


extractU ::
   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,
    Class.Floating a) =>
   LowerUpper vert horiz height width a ->
   Full vert horiz height width a
extractU :: LowerUpper vert horiz height width a
-> Full vert horiz height width a
extractU = Array (Full vert horiz height width) a
-> Full vert horiz height width a
forall shA a. Array shA a -> ArrayMatrix shA a
ArrMatrix.lift0 (Array (Full vert horiz height width) a
 -> Full vert horiz height width a)
-> (LowerUpper vert horiz height width a
    -> Array (Full vert horiz height width) a)
-> LowerUpper vert horiz height width a
-> Full vert horiz height width a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. LowerUpper vert horiz height width a
-> Array (Full vert horiz height width) a
forall vert horiz height width a.
(C vert, C horiz, C height, C width, Floating a) =>
LowerUpper vert horiz height width a
-> Full vert horiz height width a
Plain.extractU

tallExtractU ::
   (Extent.C vert, Shape.C height, Shape.C width, Class.Floating a) =>
   LowerUpper vert Extent.Small height width a -> Tri.Upper width a
tallExtractU :: LowerUpper vert Small height width a -> Upper width a
tallExtractU = Array (UpperTriangular NonUnit width) a -> Upper width a
forall shA a. Array shA a -> ArrayMatrix shA a
ArrMatrix.lift0 (Array (UpperTriangular NonUnit width) a -> Upper width a)
-> (LowerUpper vert Small height width a
    -> Array (UpperTriangular NonUnit width) a)
-> LowerUpper vert Small height width a
-> Upper width a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. LowerUpper vert Small height width a
-> Array (UpperTriangular NonUnit width) a
forall vert height width a.
(C vert, C height, C width, Floating a) =>
LowerUpper vert Small height width a -> Upper width a
Plain.tallExtractU

{- |
@tallMultiplyU transposed lu a@ multiplies the square part of @lu@
containing the upper triangular matrix with @a@.

> tallMultiplyU NonTransposed lu a == tallExtractU lu #*## a
> tallMultiplyU Transposed lu a == Tri.transpose (tallExtractU lu) #*## a
-}
tallMultiplyU ::
   (Extent.C vertA, Extent.C vert, Extent.C horiz, Shape.C height, Eq height,
    Shape.C heightA, Shape.C widthB, Class.Floating a) =>
   Transposition ->
   LowerUpper vertA Extent.Small heightA height a ->
   Full vert horiz height widthB a ->
   Full vert horiz height widthB a
tallMultiplyU :: Transposition
-> LowerUpper vertA Small heightA height a
-> Full vert horiz height widthB a
-> Full vert horiz height widthB a
tallMultiplyU Transposition
transposed = (Array (Full vert horiz height widthB) a
 -> Array (Full vert horiz height widthB) a)
-> Full vert horiz height widthB a
-> Full vert horiz height widthB a
forall shA a shB b.
(Array shA a -> Array shB b)
-> ArrayMatrix shA a -> ArrayMatrix shB b
ArrMatrix.lift1 ((Array (Full vert horiz height widthB) a
  -> Array (Full vert horiz height widthB) a)
 -> Full vert horiz height widthB a
 -> Full vert horiz height widthB a)
-> (LowerUpper vertA Small heightA height a
    -> Array (Full vert horiz height widthB) a
    -> Array (Full vert horiz height widthB) a)
-> LowerUpper vertA Small heightA height a
-> Full vert horiz height widthB a
-> Full vert horiz height widthB a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Transposition
-> LowerUpper vertA Small heightA height a
-> Array (Full vert horiz height widthB) a
-> Array (Full vert horiz height widthB) a
forall vertA vert horiz height heightA widthB a.
(C vertA, C vert, C horiz, C height, Eq height, C heightA,
 C widthB, Floating a) =>
Transposition
-> LowerUpper vertA Small heightA height a
-> Full vert horiz height widthB a
-> Full vert horiz height widthB a
Plain.tallMultiplyU Transposition
transposed

tallSolveU ::
   (Extent.C vertA, Extent.C vert, Extent.C horiz,
    Shape.C height, Shape.C width, Eq width, Shape.C nrhs, Class.Floating a) =>
   Transposition -> Conjugation ->
   LowerUpper vertA Extent.Small height width a ->
   Full vert horiz width nrhs a -> Full vert horiz width nrhs a
tallSolveU :: Transposition
-> Conjugation
-> LowerUpper vertA Small height width a
-> Full vert horiz width nrhs a
-> Full vert horiz width nrhs a
tallSolveU Transposition
transposed Conjugation
conjugated =
   (Array (Full vert horiz width nrhs) a
 -> Array (Full vert horiz width nrhs) a)
-> Full vert horiz width nrhs a -> Full vert horiz width nrhs a
forall shA a shB b.
(Array shA a -> Array shB b)
-> ArrayMatrix shA a -> ArrayMatrix shB b
ArrMatrix.lift1 ((Array (Full vert horiz width nrhs) a
  -> Array (Full vert horiz width nrhs) a)
 -> Full vert horiz width nrhs a -> Full vert horiz width nrhs a)
-> (LowerUpper vertA Small height width a
    -> Array (Full vert horiz width nrhs) a
    -> Array (Full vert horiz width nrhs) a)
-> LowerUpper vertA Small height width a
-> Full vert horiz width nrhs a
-> Full vert horiz width nrhs a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Transposition
-> Conjugation
-> LowerUpper vertA Small height width a
-> Array (Full vert horiz width nrhs) a
-> Array (Full vert horiz width nrhs) a
forall vertA vert horiz height width nrhs a.
(C vertA, C vert, C horiz, C height, C width, Eq width, C nrhs,
 Floating a) =>
Transposition
-> Conjugation
-> LowerUpper vertA Small height width a
-> Full vert horiz width nrhs a
-> Full vert horiz width nrhs a
Plain.tallSolveU Transposition
transposed Conjugation
conjugated



toMatrix ::
   (Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
   LowerUpper vert horiz height width a ->
   Full vert horiz height width a
toMatrix :: LowerUpper vert horiz height width a
-> Full vert horiz height width a
toMatrix = Array (Full vert horiz height width) a
-> Full vert horiz height width a
forall shA a. Array shA a -> ArrayMatrix shA a
ArrMatrix.lift0 (Array (Full vert horiz height width) a
 -> Full vert horiz height width a)
-> (LowerUpper vert horiz height width a
    -> Array (Full vert horiz height width) a)
-> LowerUpper vert horiz height width a
-> Full vert horiz height width a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. LowerUpper vert horiz height width a
-> Array (Full vert horiz height width) a
forall vert horiz height width a.
(C vert, C horiz, C height, Eq height, C width, Eq width,
 Floating a) =>
LowerUpper vert horiz height width a
-> Full vert horiz height width a
Plain.toMatrix


multiplyFull ::
   (Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C fuse, Eq fuse,
    Class.Floating a) =>
   LowerUpper vert horiz height fuse a ->
   Full vert horiz fuse width a ->
   Full vert horiz height width a
multiplyFull :: LowerUpper vert horiz height fuse a
-> Full vert horiz fuse width a -> Full vert horiz height width a
multiplyFull = (Array (Full vert horiz fuse width) a
 -> Array (Full vert horiz height width) a)
-> Full vert horiz fuse width a -> Full vert horiz height width a
forall shA a shB b.
(Array shA a -> Array shB b)
-> ArrayMatrix shA a -> ArrayMatrix shB b
ArrMatrix.lift1 ((Array (Full vert horiz fuse width) a
  -> Array (Full vert horiz height width) a)
 -> Full vert horiz fuse width a -> Full vert horiz height width a)
-> (LowerUpper vert horiz height fuse a
    -> Array (Full vert horiz fuse width) a
    -> Array (Full vert horiz height width) a)
-> LowerUpper vert horiz height fuse a
-> Full vert horiz fuse width a
-> Full vert horiz height width a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. LowerUpper vert horiz height fuse a
-> Array (Full vert horiz fuse width) a
-> Array (Full vert horiz height width) a
forall vert horiz height width fuse a.
(C vert, C horiz, C height, Eq height, C width, C fuse, Eq fuse,
 Floating a) =>
LowerUpper vert horiz height fuse a
-> Full vert horiz fuse width a -> Full vert horiz height width a
Plain.multiplyFull