module Numeric.LAPACK.Matrix.Lazy.UpperTriangular where

import qualified Numeric.LAPACK.Matrix.Array.Mosaic as Mosaic
import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
import qualified Numeric.LAPACK.Matrix.Quadratic as Quad
import Numeric.LAPACK.Matrix.Array.Indexed ((#!))
import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))

import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Storable as StArray
import qualified Data.Array.Comfort.Boxed.Unchecked as Array
import qualified Data.Array.Comfort.Boxed as CheckedArray
import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Boxed ((!))

import Prelude hiding (sqrt)



type Upper sh = Array.Array (Shape.Triangular Shape.Upper (Shape.Deferred sh))
type Vector sh = Array.Array (Shape.Deferred sh)


sample :: Shape.Indexed sh => sh -> (Shape.Index sh -> b) -> Array.Array sh b
sample :: sh -> (Index sh -> b) -> Array sh b
sample sh
shape Index sh -> b
f = (Index sh -> b) -> Array sh (Index sh) -> Array sh b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Index sh -> b
f (Array sh (Index sh) -> Array sh b)
-> Array sh (Index sh) -> Array sh b
forall a b. (a -> b) -> a -> b
$ sh -> Array sh (Index sh)
forall sh. Indexed sh => sh -> Array sh (Index sh)
CheckedArray.indices sh
shape


fromStorable ::
   (Shape.C sh, Class.Floating a) =>
   Mosaic.FlexUpperP pack diag sh a -> Upper sh a
fromStorable :: FlexUpperP pack diag sh a -> Upper sh a
fromStorable FlexUpperP pack diag sh a
a0 =
   let a :: Quadratic pack diag Empty Filled (Deferred sh) a
a = (sh -> Deferred sh)
-> FlexUpperP pack diag sh a
-> Quadratic pack diag Empty Filled (Deferred sh) a
forall shA shB pack property lower upper a.
(C shA, C shB) =>
(shA -> shB)
-> Quadratic pack property lower upper shA a
-> Quadratic pack property lower upper shB a
Quad.mapSize sh -> Deferred sh
forall sh. sh -> Deferred sh
Shape.Deferred FlexUpperP pack diag sh a
a0
   in Triangular Upper (Deferred sh)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall sh b. Indexed sh => sh -> (Index sh -> b) -> Array sh b
sample (Upper -> Deferred sh -> Triangular Upper (Deferred sh)
forall part size. part -> size -> Triangular part size
Shape.Triangular Upper
Shape.Upper (Quadratic pack diag Empty Filled (Deferred sh) a -> Deferred sh
forall pack property lower upper sh a.
Quadratic pack property lower upper sh a -> sh
Quad.size Quadratic pack diag Empty Filled (Deferred sh) a
a)) (Quadratic pack diag Empty Filled (Deferred sh) a
aQuadratic pack diag Empty Filled (Deferred sh) a
-> (Index (Deferred sh), Index (Deferred sh)) -> a
forall meas vert horiz height width a pack property lower upper.
(Measure meas, C vert, C horiz, Indexed height, Indexed width,
 Floating a) =>
ArrayMatrix
  pack property lower upper meas vert horiz height width a
-> (Index height, Index width) -> a
#!)

toStorable ::
   (Shape.C sh, Class.Floating a) => Upper sh a -> Mosaic.Upper sh a
toStorable :: Upper sh a -> Upper sh a
toStorable =
   (Deferred sh -> sh)
-> Quadratic Packed Arbitrary Empty Filled (Deferred sh) a
-> Upper sh a
forall shA shB pack property lower upper a.
(C shA, C shB) =>
(shA -> shB)
-> Quadratic pack property lower upper shA a
-> Quadratic pack property lower upper shB a
Quad.mapSize (\(Shape.Deferred sh
sh) -> sh
sh) (Quadratic Packed Arbitrary Empty Filled (Deferred sh) a
 -> Upper sh a)
-> (Upper sh a
    -> Quadratic Packed Arbitrary Empty Filled (Deferred sh) a)
-> Upper sh a
-> Upper sh a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Array (Triangular Upper (Deferred sh)) a
-> Quadratic Packed Arbitrary Empty Filled (Deferred sh) a
forall sh a.
(C sh, Floating a) =>
Array (Triangular Upper sh) a -> Upper sh a
Triangular.fromUpperRowMajor (Array (Triangular Upper (Deferred sh)) a
 -> Quadratic Packed Arbitrary Empty Filled (Deferred sh) a)
-> (Upper sh a -> Array (Triangular Upper (Deferred sh)) a)
-> Upper sh a
-> Quadratic Packed Arbitrary Empty Filled (Deferred sh) a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Upper sh a -> Array (Triangular Upper (Deferred sh)) a
forall sh a. (C sh, Storable a) => Array sh a -> Array sh a
StArray.fromBoxed


scaleColumns :: (Shape.C sh, Num a) => Vector sh a -> Upper sh a -> Upper sh a
scaleColumns :: Vector sh a -> Upper sh a -> Upper sh a
scaleColumns Vector sh a
d Upper sh a
a = Triangular Upper (Deferred sh)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall sh b. Indexed sh => sh -> (Index sh -> b) -> Array sh b
sample (Upper sh a -> Triangular Upper (Deferred sh)
forall sh a. Array sh a -> sh
Array.shape Upper sh a
a) ((Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall a b. (a -> b) -> a -> b
$ \(i,j) -> Upper sh a
aUpper sh a -> Index (Triangular Upper (Deferred sh)) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!(DeferredIndex sh
i,DeferredIndex sh
j) a -> a -> a
forall a. Num a => a -> a -> a
* Vector sh a
dVector sh a -> Index (Deferred sh) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!Index (Deferred sh)
DeferredIndex sh
j

scaleRows :: (Shape.C sh, Num a) => Vector sh a -> Upper sh a -> Upper sh a
scaleRows :: Vector sh a -> Upper sh a -> Upper sh a
scaleRows Vector sh a
d Upper sh a
a = Triangular Upper (Deferred sh)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall sh b. Indexed sh => sh -> (Index sh -> b) -> Array sh b
sample (Upper sh a -> Triangular Upper (Deferred sh)
forall sh a. Array sh a -> sh
Array.shape Upper sh a
a) ((Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall a b. (a -> b) -> a -> b
$ \(i,j) -> Upper sh a
aUpper sh a -> Index (Triangular Upper (Deferred sh)) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!(DeferredIndex sh
i,DeferredIndex sh
j) a -> a -> a
forall a. Num a => a -> a -> a
* Vector sh a
dVector sh a -> Index (Deferred sh) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!Index (Deferred sh)
DeferredIndex sh
i


{- |
It is an unchecked error if the shapes mismatch.
-}
multiply :: (Shape.C sh, Num a) => Upper sh a -> Upper sh a -> Upper sh a
multiply :: Upper sh a -> Upper sh a -> Upper sh a
multiply Upper sh a
a Upper sh a
b =
   Triangular Upper (Deferred sh)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall sh b. Indexed sh => sh -> (Index sh -> b) -> Array sh b
sample (Upper sh a -> Triangular Upper (Deferred sh)
forall sh a. Array sh a -> sh
Array.shape Upper sh a
a) ((Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall a b. (a -> b) -> a -> b
$
      \(i@(Shape.DeferredIndex di), j@(Shape.DeferredIndex dj)) ->
         [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([a] -> a) -> [a] -> a
forall a b. (a -> b) -> a -> b
$ (DeferredIndex sh -> a) -> [DeferredIndex sh] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (\DeferredIndex sh
k -> Upper sh a
aUpper sh a -> Index (Triangular Upper (Deferred sh)) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!(DeferredIndex sh
i,DeferredIndex sh
k)a -> a -> a
forall a. Num a => a -> a -> a
*Upper sh a
bUpper sh a -> Index (Triangular Upper (Deferred sh)) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!(DeferredIndex sh
k,DeferredIndex sh
j)) ([DeferredIndex sh] -> [a]) -> [DeferredIndex sh] -> [a]
forall a b. (a -> b) -> a -> b
$ (Int -> DeferredIndex sh) -> [Int] -> [DeferredIndex sh]
forall a b. (a -> b) -> [a] -> [b]
map Int -> DeferredIndex sh
forall sh. Int -> DeferredIndex sh
Shape.DeferredIndex [Int
di..Int
dj]

{- |
@multiplyStrictPart a b@
is almost @multiply (clearDiagonal a) (clearDiagonal b)@,
but it is more lazy
since it does not access the elements that are multiplied with zero.
-}
multiplyStrictPart ::
   (Shape.C sh, Num a) => Upper sh a -> Upper sh a -> Upper sh a
multiplyStrictPart :: Upper sh a -> Upper sh a -> Upper sh a
multiplyStrictPart Upper sh a
a Upper sh a
b =
   Triangular Upper (Deferred sh)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall sh b. Indexed sh => sh -> (Index sh -> b) -> Array sh b
sample (Upper sh a -> Triangular Upper (Deferred sh)
forall sh a. Array sh a -> sh
Array.shape Upper sh a
a) ((Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall a b. (a -> b) -> a -> b
$
      \(i@(Shape.DeferredIndex di), j@(Shape.DeferredIndex dj)) ->
         [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([a] -> a) -> [a] -> a
forall a b. (a -> b) -> a -> b
$ (DeferredIndex sh -> a) -> [DeferredIndex sh] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (\DeferredIndex sh
k -> Upper sh a
aUpper sh a -> Index (Triangular Upper (Deferred sh)) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!(DeferredIndex sh
i,DeferredIndex sh
k)a -> a -> a
forall a. Num a => a -> a -> a
*Upper sh a
bUpper sh a -> Index (Triangular Upper (Deferred sh)) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!(DeferredIndex sh
k,DeferredIndex sh
j)) ([DeferredIndex sh] -> [a]) -> [DeferredIndex sh] -> [a]
forall a b. (a -> b) -> a -> b
$
            (Int -> DeferredIndex sh) -> [Int] -> [DeferredIndex sh]
forall a b. (a -> b) -> [a] -> [b]
map Int -> DeferredIndex sh
forall sh. Int -> DeferredIndex sh
Shape.DeferredIndex [Int
diInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1..Int
djInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]

takeDiagonal :: (Shape.C sh, Num a) => Upper sh a -> Vector sh a
takeDiagonal :: Upper sh a -> Vector sh a
takeDiagonal Upper sh a
a =
   Deferred sh -> (Index (Deferred sh) -> a) -> Vector sh a
forall sh b. Indexed sh => sh -> (Index sh -> b) -> Array sh b
sample (Triangular Upper (Deferred sh) -> Deferred sh
forall part size. Triangular part size -> size
Shape.triangularSize (Triangular Upper (Deferred sh) -> Deferred sh)
-> Triangular Upper (Deferred sh) -> Deferred sh
forall a b. (a -> b) -> a -> b
$ Upper sh a -> Triangular Upper (Deferred sh)
forall sh a. Array sh a -> sh
Array.shape Upper sh a
a) ((Index (Deferred sh) -> a) -> Vector sh a)
-> (Index (Deferred sh) -> a) -> Vector sh a
forall a b. (a -> b) -> a -> b
$ \Index (Deferred sh)
i -> Upper sh a
aUpper sh a -> Index (Triangular Upper (Deferred sh)) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!(Index (Deferred sh)
DeferredIndex sh
i,Index (Deferred sh)
DeferredIndex sh
i)

replaceDiagonal ::
   (Shape.C sh, Num a) => Vector sh a -> Upper sh a -> Upper sh a
replaceDiagonal :: Vector sh a -> Upper sh a -> Upper sh a
replaceDiagonal Vector sh a
d Upper sh a
a =
   Triangular Upper (Deferred sh)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall sh b. Indexed sh => sh -> (Index sh -> b) -> Array sh b
sample (Upper sh a -> Triangular Upper (Deferred sh)
forall sh a. Array sh a -> sh
Array.shape Upper sh a
a) ((Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall a b. (a -> b) -> a -> b
$ \(i,j) -> if DeferredIndex sh
iDeferredIndex sh -> DeferredIndex sh -> Bool
forall a. Eq a => a -> a -> Bool
==DeferredIndex sh
j then Vector sh a
dVector sh a -> Index (Deferred sh) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!Index (Deferred sh)
DeferredIndex sh
i else Upper sh a
aUpper sh a -> Index (Triangular Upper (Deferred sh)) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!(DeferredIndex sh
i,DeferredIndex sh
j)

rank2Part :: (Shape.C sh, Num a) => Vector sh a -> Upper sh a
rank2Part :: Vector sh a -> Upper sh a
rank2Part Vector sh a
d =
   Triangular Upper (Deferred sh)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall sh b. Indexed sh => sh -> (Index sh -> b) -> Array sh b
sample (Upper -> Deferred sh -> Triangular Upper (Deferred sh)
forall part size. part -> size -> Triangular part size
Shape.Triangular Upper
Shape.Upper (Deferred sh -> Triangular Upper (Deferred sh))
-> Deferred sh -> Triangular Upper (Deferred sh)
forall a b. (a -> b) -> a -> b
$ Vector sh a -> Deferred sh
forall sh a. Array sh a -> sh
Array.shape Vector sh a
d) ((Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall a b. (a -> b) -> a -> b
$ \(i,j) -> Vector sh a
dVector sh a -> Index (Deferred sh) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!Index (Deferred sh)
DeferredIndex sh
i a -> a -> a
forall a. Num a => a -> a -> a
+ Vector sh a
dVector sh a -> Index (Deferred sh) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!Index (Deferred sh)
DeferredIndex sh
j

rank2DiffPart :: (Shape.C sh, Num a) => Vector sh a -> Upper sh a
rank2DiffPart :: Vector sh a -> Upper sh a
rank2DiffPart Vector sh a
d =
   Triangular Upper (Deferred sh)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall sh b. Indexed sh => sh -> (Index sh -> b) -> Array sh b
sample (Upper -> Deferred sh -> Triangular Upper (Deferred sh)
forall part size. part -> size -> Triangular part size
Shape.Triangular Upper
Shape.Upper (Deferred sh -> Triangular Upper (Deferred sh))
-> Deferred sh -> Triangular Upper (Deferred sh)
forall a b. (a -> b) -> a -> b
$ Vector sh a -> Deferred sh
forall sh a. Array sh a -> sh
Array.shape Vector sh a
d) ((Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a)
-> (Index (Triangular Upper (Deferred sh)) -> a) -> Upper sh a
forall a b. (a -> b) -> a -> b
$ \(i,j) -> Vector sh a
dVector sh a -> Index (Deferred sh) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!Index (Deferred sh)
DeferredIndex sh
i a -> a -> a
forall a. Num a => a -> a -> a
- Vector sh a
dVector sh a -> Index (Deferred sh) -> a
forall sh a. Indexed sh => Array sh a -> Index sh -> a
!Index (Deferred sh)
DeferredIndex sh
j


{- |
Lazy implicit solver
-}
{-
A = (D+U)*(D+U) = D*D + D*U + U*D + U*U

let U = divide (A-D*D - U*U) (toColumn(D)*1^T+1*toRow(D))
-}
sqrt :: (Shape.C sh, Fractional a) => (a -> a) -> Upper sh a -> Upper sh a
sqrt :: (a -> a) -> Upper sh a -> Upper sh a
sqrt a -> a
sqrtF = (Upper (Unchecked sh) a -> Upper (Unchecked sh) a)
-> Upper sh a -> Upper sh a
forall sh a.
C sh =>
(Upper (Unchecked sh) a -> Upper (Unchecked sh) a)
-> Upper sh a -> Upper sh a
applyUnchecked ((Upper (Unchecked sh) a -> Upper (Unchecked sh) a)
 -> Upper sh a -> Upper sh a)
-> (Upper (Unchecked sh) a -> Upper (Unchecked sh) a)
-> Upper sh a
-> Upper sh a
forall a b. (a -> b) -> a -> b
$ \Upper (Unchecked sh) a
a ->
   let d :: Array (Deferred (Unchecked sh)) a
d = (a -> a)
-> Array (Deferred (Unchecked sh)) a
-> Array (Deferred (Unchecked sh)) a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> a
sqrtF (Array (Deferred (Unchecked sh)) a
 -> Array (Deferred (Unchecked sh)) a)
-> Array (Deferred (Unchecked sh)) a
-> Array (Deferred (Unchecked sh)) a
forall a b. (a -> b) -> a -> b
$ Upper (Unchecked sh) a -> Array (Deferred (Unchecked sh)) a
forall sh a. (C sh, Num a) => Upper sh a -> Vector sh a
takeDiagonal Upper (Unchecked sh) a
a
       u :: Upper (Unchecked sh) a
u =
         Triangular Upper (Deferred (Unchecked sh))
-> Upper (Unchecked sh) a -> Upper (Unchecked sh) a
forall sh1 sh0 a. sh1 -> Array sh0 a -> Array sh1 a
Array.reshape (Upper (Unchecked sh) a
-> Triangular Upper (Deferred (Unchecked sh))
forall sh a. Array sh a -> sh
Array.shape Upper (Unchecked sh) a
a) (Upper (Unchecked sh) a -> Upper (Unchecked sh) a)
-> Upper (Unchecked sh) a -> Upper (Unchecked sh) a
forall a b. (a -> b) -> a -> b
$
         (a -> a -> a)
-> Upper (Unchecked sh) a
-> Upper (Unchecked sh) a
-> Upper (Unchecked sh) a
forall sh a b c.
C sh =>
(a -> b -> c) -> Array sh a -> Array sh b -> Array sh c
Array.zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/)
            ((a -> a -> a)
-> Upper (Unchecked sh) a
-> Upper (Unchecked sh) a
-> Upper (Unchecked sh) a
forall sh a b c.
C sh =>
(a -> b -> c) -> Array sh a -> Array sh b -> Array sh c
Array.zipWith (-) Upper (Unchecked sh) a
a (Upper (Unchecked sh) a
-> Upper (Unchecked sh) a -> Upper (Unchecked sh) a
forall sh a.
(C sh, Num a) =>
Upper sh a -> Upper sh a -> Upper sh a
multiplyStrictPart Upper (Unchecked sh) a
u Upper (Unchecked sh) a
u))
            (Array (Deferred (Unchecked sh)) a -> Upper (Unchecked sh) a
forall sh a. (C sh, Num a) => Vector sh a -> Upper sh a
rank2Part Array (Deferred (Unchecked sh)) a
d)
   in Array (Deferred (Unchecked sh)) a
-> Upper (Unchecked sh) a -> Upper (Unchecked sh) a
forall sh a.
(C sh, Num a) =>
Vector sh a -> Upper sh a -> Upper sh a
replaceDiagonal Array (Deferred (Unchecked sh)) a
d Upper (Unchecked sh) a
u


{- |
Parlett recursion for lifting a scalar function to an upper triangular matrix.

Given A and the diagonal of f(A) it solves A*f(A) = f(A)*A.

Requires distinct values on the diagonal,
where even almost close values can produce dramatic errors.
But it admits for a nice lazy implicit implementation.
-}
{-
/g h i\   /a b c\   /  x y\   /  x y\   /a b c\
|  j k| = |  d e| * |    z| - |    z| * |  d e|
\    l/   \    f/   \     /   \     /   \    f/

h = a*x - x*d = (a-d)*x
k = d*z - z*f = (d-f)*z
i = a*y+b*z - (x*e+y*f) = b*z-x*e + (a-f)*y


0 = A*f(A) - f(A)*A

A = E+V
f(A) = D+U

0 = A*(D+U) - (D+U)*A
D*A-A*D = (E+V)*U - U*(E+V)
D*A-A*D + U*V-V*U = E*U-U*E

U = divide (D*A-A*D + U*V-V*U) (toColumn(E)*1^T - 1*toRow(E))
-}
parlett :: (Shape.C sh, Fractional a) => (a -> a) -> Upper sh a -> Upper sh a
parlett :: (a -> a) -> Upper sh a -> Upper sh a
parlett a -> a
f = (Upper (Unchecked sh) a -> Upper (Unchecked sh) a)
-> Upper sh a -> Upper sh a
forall sh a.
C sh =>
(Upper (Unchecked sh) a -> Upper (Unchecked sh) a)
-> Upper sh a -> Upper sh a
applyUnchecked ((Upper (Unchecked sh) a -> Upper (Unchecked sh) a)
 -> Upper sh a -> Upper sh a)
-> (Upper (Unchecked sh) a -> Upper (Unchecked sh) a)
-> Upper sh a
-> Upper sh a
forall a b. (a -> b) -> a -> b
$ \Upper (Unchecked sh) a
a ->
   let e :: Vector (Unchecked sh) a
e = Upper (Unchecked sh) a -> Vector (Unchecked sh) a
forall sh a. (C sh, Num a) => Upper sh a -> Vector sh a
takeDiagonal Upper (Unchecked sh) a
a
       d :: Vector (Unchecked sh) a
d = (a -> a) -> Vector (Unchecked sh) a -> Vector (Unchecked sh) a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> a
f Vector (Unchecked sh) a
e
       u :: Upper (Unchecked sh) a
u =
         Triangular Upper (Deferred (Unchecked sh))
-> Upper (Unchecked sh) a -> Upper (Unchecked sh) a
forall sh1 sh0 a. sh1 -> Array sh0 a -> Array sh1 a
Array.reshape (Upper (Unchecked sh) a
-> Triangular Upper (Deferred (Unchecked sh))
forall sh a. Array sh a -> sh
Array.shape Upper (Unchecked sh) a
a) (Upper (Unchecked sh) a -> Upper (Unchecked sh) a)
-> Upper (Unchecked sh) a -> Upper (Unchecked sh) a
forall a b. (a -> b) -> a -> b
$
         (a -> a -> a)
-> Upper (Unchecked sh) a
-> Upper (Unchecked sh) a
-> Upper (Unchecked sh) a
forall sh a b c.
C sh =>
(a -> b -> c) -> Array sh a -> Array sh b -> Array sh c
Array.zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/)
            ((a -> a -> a)
-> Upper (Unchecked sh) a
-> Upper (Unchecked sh) a
-> Upper (Unchecked sh) a
forall sh a b c.
C sh =>
(a -> b -> c) -> Array sh a -> Array sh b -> Array sh c
Array.zipWith (-)
               ((a -> a -> a)
-> Upper (Unchecked sh) a
-> Upper (Unchecked sh) a
-> Upper (Unchecked sh) a
forall sh a b c.
C sh =>
(a -> b -> c) -> Array sh a -> Array sh b -> Array sh c
Array.zipWith a -> a -> a
forall a. Num a => a -> a -> a
(+) (Vector (Unchecked sh) a
-> Upper (Unchecked sh) a -> Upper (Unchecked sh) a
forall sh a.
(C sh, Num a) =>
Vector sh a -> Upper sh a -> Upper sh a
scaleRows    Vector (Unchecked sh) a
d Upper (Unchecked sh) a
a) (Upper (Unchecked sh) a
-> Upper (Unchecked sh) a -> Upper (Unchecked sh) a
forall sh a.
(C sh, Num a) =>
Upper sh a -> Upper sh a -> Upper sh a
multiplyStrictPart Upper (Unchecked sh) a
u Upper (Unchecked sh) a
a))
               ((a -> a -> a)
-> Upper (Unchecked sh) a
-> Upper (Unchecked sh) a
-> Upper (Unchecked sh) a
forall sh a b c.
C sh =>
(a -> b -> c) -> Array sh a -> Array sh b -> Array sh c
Array.zipWith a -> a -> a
forall a. Num a => a -> a -> a
(+) (Vector (Unchecked sh) a
-> Upper (Unchecked sh) a -> Upper (Unchecked sh) a
forall sh a.
(C sh, Num a) =>
Vector sh a -> Upper sh a -> Upper sh a
scaleColumns Vector (Unchecked sh) a
d Upper (Unchecked sh) a
a) (Upper (Unchecked sh) a
-> Upper (Unchecked sh) a -> Upper (Unchecked sh) a
forall sh a.
(C sh, Num a) =>
Upper sh a -> Upper sh a -> Upper sh a
multiplyStrictPart Upper (Unchecked sh) a
a Upper (Unchecked sh) a
u)))
            (Vector (Unchecked sh) a -> Upper (Unchecked sh) a
forall sh a. (C sh, Num a) => Vector sh a -> Upper sh a
rank2DiffPart Vector (Unchecked sh) a
e)
   in Vector (Unchecked sh) a
-> Upper (Unchecked sh) a -> Upper (Unchecked sh) a
forall sh a.
(C sh, Num a) =>
Vector sh a -> Upper sh a -> Upper sh a
replaceDiagonal Vector (Unchecked sh) a
d Upper (Unchecked sh) a
u


applyUnchecked ::
   (Shape.C sh) =>
   (Upper (Unchecked sh) a -> Upper (Unchecked sh) a) ->
   Upper sh a -> Upper sh a
applyUnchecked :: (Upper (Unchecked sh) a -> Upper (Unchecked sh) a)
-> Upper sh a -> Upper sh a
applyUnchecked Upper (Unchecked sh) a -> Upper (Unchecked sh) a
f =
   (Triangular Upper (Deferred (Unchecked sh))
 -> Triangular Upper (Deferred sh))
-> Upper (Unchecked sh) a -> Upper sh a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape
      (\(Shape.Triangular Upper
uplo (Shape.Deferred (Unchecked sh
sh))) ->
         Upper -> Deferred sh -> Triangular Upper (Deferred sh)
forall part size. part -> size -> Triangular part size
Shape.Triangular Upper
uplo (sh -> Deferred sh
forall sh. sh -> Deferred sh
Shape.Deferred sh
sh)) (Upper (Unchecked sh) a -> Upper sh a)
-> (Upper sh a -> Upper (Unchecked sh) a)
-> Upper sh a
-> Upper sh a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Upper (Unchecked sh) a -> Upper (Unchecked sh) a
f (Upper (Unchecked sh) a -> Upper (Unchecked sh) a)
-> (Upper sh a -> Upper (Unchecked sh) a)
-> Upper sh a
-> Upper (Unchecked sh) a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (Triangular Upper (Deferred sh)
 -> Triangular Upper (Deferred (Unchecked sh)))
-> Upper sh a -> Upper (Unchecked sh) a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape
      (\(Shape.Triangular Upper
uplo (Shape.Deferred sh
sh)) ->
         Upper
-> Deferred (Unchecked sh)
-> Triangular Upper (Deferred (Unchecked sh))
forall part size. part -> size -> Triangular part size
Shape.Triangular Upper
uplo (Unchecked sh -> Deferred (Unchecked sh)
forall sh. sh -> Deferred sh
Shape.Deferred (sh -> Unchecked sh
forall sh. sh -> Unchecked sh
Unchecked sh
sh)))