    {-# LANGUAGE CPP #-}
    {-# LANGUAGE FlexibleContexts #-}
    {-# LANGUAGE FlexibleInstances #-}
    {-# LANGUAGE MultiParamTypeClasses #-}
    {-# LANGUAGE FunctionalDependencies #-}
    {-# LANGUAGE ScopedTypeVariables #-}
    {-# LANGUAGE ViewPatterns #-}

    {-# OPTIONS_GHC -fno-warn-missing-signatures #-}
    {-# OPTIONS_GHC -fno-warn-orphans #-}

    -----------------------------------------------------------------------------
    {- |
    Module      :  Internal.Util
    Copyright   :  (c) Alberto Ruiz 2013
    License     :  BSD3
    Maintainer  :  Alberto Ruiz
    Stability   :  provisional

    -}
    -----------------------------------------------------------------------------

    module Internal.Util(

        -- * Convenience functions
        vector, matrix,
        disp,
        formatSparse,
        approxInt,
        dispDots,
        dispBlanks,
        formatShort,
        dispShort,
        zeros, ones,
        diagl,
        row,
        col,
        (&), (¦), (|||), (——), (===),
        (?), (¿),
        Indexable(..), size,
        Numeric,
        rand, randn,
        cross,
        norm,
        ℕ,ℤ,ℝ,ℂ,iC,
        Normed(..), norm_Frob, norm_nuclear,
        magnit,
        normalize,
        mt,
        (~!~),
        pairwiseD2,
        rowOuters,
        null1,
        null1sym,
        -- * Convolution
        -- ** 1D
        corr, conv, corrMin,
        -- ** 2D
        corr2, conv2, separable,
        block2x2,block3x3,view1,unView1,foldMatrix,
        gaussElim_1, gaussElim_2, gaussElim,
        luST, luSolve', luSolve'', luPacked', luPacked'',
        invershur
    ) where

    import Internal.Vector
    import Internal.Matrix hiding (size)
    import Internal.Numeric
    import Internal.Element
    import Internal.Container
    import Internal.Vectorized
    import Internal.IO
    import Internal.Algorithms hiding (Normed,linearSolve',luSolve', luPacked')
    import Numeric.Matrix()
    import Numeric.Vector()
    import Internal.Random
    import Internal.Convolution
    import Control.Monad(when,forM_)
    import Text.Printf
    import Data.List.Split(splitOn)
    import Data.List(intercalate,sortBy,foldl')
    import Control.Arrow((&&&),(***))
    import Data.Complex
    import Data.Function(on)
    import Internal.ST
    #if MIN_VERSION_base(4,11,0)
    import Prelude hiding ((<>))
    #endif

    type ℝ = Double
    type ℕ = Int
    type ℤ = Int
    type ℂ = Complex Double

    -- | imaginary unit
    iC :: C
    iC :: Complex Double
    iC = Double
    0Double -> Double -> Complex Double
    forall a. a -> a -> Complex a
    :+Double
    1

    {- | Create a real vector.

    >>> vector [1..5]
    [1.0,2.0,3.0,4.0,5.0]
    it :: Vector R

    -}
    vector :: [R] -> Vector R
    vector :: [Double] -> Vector Double
    vector = [Double] -> Vector Double
    forall a. Storable a => [a] -> Vector a
    fromList

    {- | Create a real matrix.

    >>> matrix 5 [1..15]
    (3><5)
     [  1.0,  2.0,  3.0,  4.0,  5.0
     ,  6.0,  7.0,  8.0,  9.0, 10.0
     , 11.0, 12.0, 13.0, 14.0, 15.0 ]

    -}
    matrix
      :: Int -- ^ number of columns
      -> [R] -- ^ elements in row order
      -> Matrix R
    matrix :: Int -> [Double] -> Matrix Double
    matrix Int
    c = Int -> Vector Double -> Matrix Double
    forall t. Storable t => Int -> Vector t -> Matrix t
    reshape Int
    c (Vector Double -> Matrix Double)
    -> ([Double] -> Vector Double) -> [Double] -> Matrix Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [Double] -> Vector Double
    forall a. Storable a => [a] -> Vector a
    fromList

    {- | print a real matrix with given number of digits after the decimal point

    >>> disp 5 $ ident 2 / 3
    2x2
    0.33333  0.00000
    0.00000  0.33333

    -}
    disp :: Int -> Matrix Double -> IO ()

    disp :: Int -> Matrix Double -> IO ()
    disp Int
    n = String -> IO ()
    putStr (String -> IO ())
    -> (Matrix Double -> String) -> Matrix Double -> IO ()
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Int -> Matrix Double -> String
    dispf Int
    n

    {- | create a real diagonal matrix from a list

    >>> diagl [1,2,3]
    (3><3)
     [ 1.0, 0.0, 0.0
     , 0.0, 2.0, 0.0
     , 0.0, 0.0, 3.0 ]

    -}
    diagl :: [Double] -> Matrix Double
    diagl :: [Double] -> Matrix Double
    diagl = Vector Double -> Matrix Double
    forall a. (Num a, Element a) => Vector a -> Matrix a
    diag (Vector Double -> Matrix Double)
    -> ([Double] -> Vector Double) -> [Double] -> Matrix Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [Double] -> Vector Double
    forall a. Storable a => [a] -> Vector a
    fromList

    -- | a real matrix of zeros
    zeros :: Int -- ^ rows
          -> Int -- ^ columns
          -> Matrix Double
    zeros :: Int -> Int -> Matrix Double
    zeros Int
    r Int
    c = Double -> (Int, Int) -> Matrix Double
    forall e d (c :: * -> *). Konst e d c => e -> d -> c e
    konst Double
    0 (Int
    r,Int
    c)

    -- | a real matrix of ones
    ones :: Int -- ^ rows
         -> Int -- ^ columns
         -> Matrix Double
    ones :: Int -> Int -> Matrix Double
    ones Int
    r Int
    c = Double -> (Int, Int) -> Matrix Double
    forall e d (c :: * -> *). Konst e d c => e -> d -> c e
    konst Double
    1 (Int
    r,Int
    c)

    -- | concatenation of real vectors
    infixl 3 &
    (&) :: Vector Double -> Vector Double -> Vector Double
    Vector Double
    a & :: Vector Double -> Vector Double -> Vector Double
    & Vector Double
    b = [Vector Double] -> Vector Double
    forall t. Storable t => [Vector t] -> Vector t
    vjoin [Vector Double
    a,Vector Double
    b]

    {- | horizontal concatenation

    >>> ident 3 ||| konst 7 (3,4)
    (3><7)
     [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0
     , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0
     , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ]

    -}
    infixl 3 |||
    (|||) :: Element t => Matrix t -> Matrix t -> Matrix t
    Matrix t
    a ||| :: forall t. Element t => Matrix t -> Matrix t -> Matrix t
    ||| Matrix t
    b = [[Matrix t]] -> Matrix t
    forall t. Element t => [[Matrix t]] -> Matrix t
    fromBlocks [[Matrix t
    a,Matrix t
    b]]

    -- | a synonym for ('|||') (unicode 0x00a6, broken bar)
    infixl 3 ¦
    (¦) :: Matrix Double -> Matrix Double -> Matrix Double
    ¦ :: Matrix Double -> Matrix Double -> Matrix Double
    (¦) = Matrix Double -> Matrix Double -> Matrix Double
    forall t. Element t => Matrix t -> Matrix t -> Matrix t
    (|||)

    -- | vertical concatenation
    --
    (===) :: Element t => Matrix t -> Matrix t -> Matrix t
    infixl 2 ===
    Matrix t
    a === :: forall t. Element t => Matrix t -> Matrix t -> Matrix t
    === Matrix t
    b = [[Matrix t]] -> Matrix t
    forall t. Element t => [[Matrix t]] -> Matrix t
    fromBlocks [[Matrix t
    a],[Matrix t
    b]]

    -- | a synonym for ('===') (unicode 0x2014, em dash)
    (——) :: Matrix Double -> Matrix Double -> Matrix Double
    infixl 2 ——
    —— :: Matrix Double -> Matrix Double -> Matrix Double
    (——) = Matrix Double -> Matrix Double -> Matrix Double
    forall t. Element t => Matrix t -> Matrix t -> Matrix t
    (===)

    -- | create a single row real matrix from a list
    --
    -- >>> row [2,3,1,8]
    -- (1><4)
    --  [ 2.0, 3.0, 1.0, 8.0 ]
    --
    row :: [Double] -> Matrix Double
    row :: [Double] -> Matrix Double
    row = Vector Double -> Matrix Double
    forall a. Storable a => Vector a -> Matrix a
    asRow (Vector Double -> Matrix Double)
    -> ([Double] -> Vector Double) -> [Double] -> Matrix Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [Double] -> Vector Double
    forall a. Storable a => [a] -> Vector a
    fromList

    -- | create a single column real matrix from a list
    --
    -- >>> col [7,-2,4]
    -- (3><1)
    --  [  7.0
    --  , -2.0
    --  ,  4.0 ]
    --
    col :: [Double] -> Matrix Double
    col :: [Double] -> Matrix Double
    col = Vector Double -> Matrix Double
    forall a. Storable a => Vector a -> Matrix a
    asColumn (Vector Double -> Matrix Double)
    -> ([Double] -> Vector Double) -> [Double] -> Matrix Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [Double] -> Vector Double
    forall a. Storable a => [a] -> Vector a
    fromList

    {- | extract rows

    >>> (20><4) [1..] ? [2,1,1]
    (3><4)
     [ 9.0, 10.0, 11.0, 12.0
     , 5.0,  6.0,  7.0,  8.0
     , 5.0,  6.0,  7.0,  8.0 ]

    -}
    infixl 9 ?
    (?) :: Element t => Matrix t -> [Int] -> Matrix t
    ? :: forall t. Element t => Matrix t -> [Int] -> Matrix t
    (?) = ([Int] -> Matrix t -> Matrix t) -> Matrix t -> [Int] -> Matrix t
    forall a b c. (a -> b -> c) -> b -> a -> c
    flip [Int] -> Matrix t -> Matrix t
    forall t. Element t => [Int] -> Matrix t -> Matrix t
    extractRows

    {- | extract columns

    (unicode 0x00bf, inverted question mark, Alt-Gr ?)

    >>> (3><4) [1..] ¿ [3,0]
    (3><2)
     [  4.0, 1.0
     ,  8.0, 5.0
     , 12.0, 9.0 ]

    -}
    infixl 9 ¿
    (¿) :: Element t => Matrix t -> [Int] -> Matrix t
    ¿ :: forall t. Element t => Matrix t -> [Int] -> Matrix t
    (¿)= ([Int] -> Matrix t -> Matrix t) -> Matrix t -> [Int] -> Matrix t
    forall a b c. (a -> b -> c) -> b -> a -> c
    flip [Int] -> Matrix t -> Matrix t
    forall t. Element t => [Int] -> Matrix t -> Matrix t
    extractColumns

    cross :: Product t => Vector t -> Vector t -> Vector t
    -- ^ cross product (for three-element vectors)
    cross :: forall t. Product t => Vector t -> Vector t -> Vector t
    cross Vector t
    x Vector t
    y | Vector t -> Int
    forall t. Storable t => Vector t -> Int
    dim Vector t
    x Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    == Int
    3 Bool -> Bool -> Bool
    && Vector t -> Int
    forall t. Storable t => Vector t -> Int
    dim Vector t
    y Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    == Int
    3 = [t] -> Vector t
    forall a. Storable a => [a] -> Vector a
    fromList [t
    z1,t
    z2,t
    z3]
              | Bool
    otherwise = String -> Vector t
    forall a. HasCallStack => String -> a
    error (String -> Vector t) -> String -> Vector t
    forall a b. (a -> b) -> a -> b
    $ String
    "the cross product requires 3-element vectors (sizes given: "
                                    String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++Int -> String
    forall a. Show a => a -> String
    show (Vector t -> Int
    forall t. Storable t => Vector t -> Int
    dim Vector t
    x)String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++String
    " and "String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++Int -> String
    forall a. Show a => a -> String
    show (Vector t -> Int
    forall t. Storable t => Vector t -> Int
    dim Vector t
    y)String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++String
    ")"
      where
        [t
    x1,t
    x2,t
    x3] = Vector t -> [t]
    forall a. Storable a => Vector a -> [a]
    toList Vector t
    x
        [t
    y1,t
    y2,t
    y3] = Vector t -> [t]
    forall a. Storable a => Vector a -> [a]
    toList Vector t
    y
        z1 :: t
    z1 = t
    x2t -> t -> t
    forall a. Num a => a -> a -> a
    *t
    y3t -> t -> t
    forall a. Num a => a -> a -> a
    -t
    x3t -> t -> t
    forall a. Num a => a -> a -> a
    *t
    y2
        z2 :: t
    z2 = t
    x3t -> t -> t
    forall a. Num a => a -> a -> a
    *t
    y1t -> t -> t
    forall a. Num a => a -> a -> a
    -t
    x1t -> t -> t
    forall a. Num a => a -> a -> a
    *t
    y3
        z3 :: t
    z3 = t
    x1t -> t -> t
    forall a. Num a => a -> a -> a
    *t
    y2t -> t -> t
    forall a. Num a => a -> a -> a
    -t
    x2t -> t -> t
    forall a. Num a => a -> a -> a
    *t
    y1

    {-# SPECIALIZE cross :: Vector Double -> Vector Double -> Vector Double #-}
    {-# SPECIALIZE cross :: Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double) #-}

    norm :: Vector Double -> Double
    -- ^ 2-norm of real vector
    norm :: Vector Double -> Double
    norm = NormType -> Vector Double -> RealOf Double
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm2

    -- | p-norm for vectors, operator norm for matrices
    class Normed a
      where
        norm_0   :: a -> R
        norm_1   :: a -> R
        norm_2   :: a -> R
        norm_Inf :: a -> R

    instance Normed (Vector R)
      where
        norm_0 :: Vector Double -> Double
    norm_0 Vector Double
    v = Vector Double -> Double
    forall (c :: * -> *) e. Container c e => c e -> e
    sumElements (Vector Double -> Vector Double
    forall e (c :: * -> *). (Ord e, Container c e) => c e -> c e
    step (Vector Double -> Vector Double
    forall a. Num a => a -> a
    abs Vector Double
    v Vector Double -> Vector Double -> Vector Double
    forall a. Num a => a -> a -> a
    - Double -> Vector Double
    forall (c :: * -> *) e. Container c e => e -> c e
    scalar (Double
    epsDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    *Vector Double -> RealOf Double
    forall e. Product e => Vector e -> RealOf e
    normInf Vector Double
    v)))
        norm_1 :: Vector Double -> Double
    norm_1 = NormType -> Vector Double -> RealOf Double
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1
        norm_2 :: Vector Double -> Double
    norm_2 = NormType -> Vector Double -> RealOf Double
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm2
        norm_Inf :: Vector Double -> Double
    norm_Inf = NormType -> Vector Double -> RealOf Double
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    Infinity

    instance Normed (Vector C)
      where
        norm_0 :: Vector (Complex Double) -> Double
    norm_0 Vector (Complex Double)
    v = Vector Double -> Double
    forall (c :: * -> *) e. Container c e => c e -> e
    sumElements (Vector Double -> Vector Double
    forall e (c :: * -> *). (Ord e, Container c e) => c e -> c e
    step ((Vector Double, Vector Double) -> Vector Double
    forall a b. (a, b) -> a
    fst (Vector (Complex Double) -> (Vector Double, Vector Double)
    forall t (c :: * -> *).
    (Convert t, Complexable c, RealElement t) =>
    c (Complex t) -> (c t, c t)
    forall (c :: * -> *).
    (Complexable c, RealElement Double) =>
    c (Complex Double) -> (c Double, c Double)
    fromComplex (Vector (Complex Double) -> Vector (Complex Double)
    forall a. Num a => a -> a
    abs Vector (Complex Double)
    v)) Vector Double -> Vector Double -> Vector Double
    forall a. Num a => a -> a -> a
    - Double -> Vector Double
    forall (c :: * -> *) e. Container c e => e -> c e
    scalar (Double
    epsDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    *Vector (Complex Double) -> RealOf (Complex Double)
    forall e. Product e => Vector e -> RealOf e
    normInf Vector (Complex Double)
    v)))
        norm_1 :: Vector (Complex Double) -> Double
    norm_1 = NormType -> Vector (Complex Double) -> RealOf (Complex Double)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1
        norm_2 :: Vector (Complex Double) -> Double
    norm_2 = NormType -> Vector (Complex Double) -> RealOf (Complex Double)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm2
        norm_Inf :: Vector (Complex Double) -> Double
    norm_Inf = NormType -> Vector (Complex Double) -> RealOf (Complex Double)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    Infinity

    instance Normed (Matrix R)
      where
        norm_0 :: Matrix Double -> Double
    norm_0 = Vector Double -> Double
    forall a. Normed a => a -> Double
    norm_0 (Vector Double -> Double)
    -> (Matrix Double -> Vector Double) -> Matrix Double -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Double -> Vector Double
    forall t. Element t => Matrix t -> Vector t
    flatten
        norm_1 :: Matrix Double -> Double
    norm_1 = NormType -> Matrix Double -> RealOf Double
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1
        norm_2 :: Matrix Double -> Double
    norm_2 = NormType -> Matrix Double -> RealOf Double
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm2
        norm_Inf :: Matrix Double -> Double
    norm_Inf = NormType -> Matrix Double -> RealOf Double
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    Infinity

    instance Normed (Matrix C)
      where
        norm_0 :: Matrix (Complex Double) -> Double
    norm_0 = Vector (Complex Double) -> Double
    forall a. Normed a => a -> Double
    norm_0 (Vector (Complex Double) -> Double)
    -> (Matrix (Complex Double) -> Vector (Complex Double))
    -> Matrix (Complex Double)
    -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix (Complex Double) -> Vector (Complex Double)
    forall t. Element t => Matrix t -> Vector t
    flatten
        norm_1 :: Matrix (Complex Double) -> Double
    norm_1 = NormType -> Matrix (Complex Double) -> RealOf (Complex Double)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm1
        norm_2 :: Matrix (Complex Double) -> Double
    norm_2 = NormType -> Matrix (Complex Double) -> RealOf (Complex Double)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    PNorm2
        norm_Inf :: Matrix (Complex Double) -> Double
    norm_Inf = NormType -> Matrix (Complex Double) -> RealOf (Complex Double)
    forall (c :: * -> *) t. Normed c t => NormType -> c t -> RealOf t
    pnorm NormType
    Infinity

    instance Normed (Vector I)
      where
        norm_0 :: Vector I -> Double
    norm_0 = I -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (I -> Double) -> (Vector I -> I) -> Vector I -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector I -> I
    forall (c :: * -> *) e. Container c e => c e -> e
    sumElements (Vector I -> I) -> (Vector I -> Vector I) -> Vector I -> I
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector I -> Vector I
    forall e (c :: * -> *). (Ord e, Container c e) => c e -> c e
    step (Vector I -> Vector I)
    -> (Vector I -> Vector I) -> Vector I -> Vector I
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector I -> Vector I
    forall a. Num a => a -> a
    abs
        norm_1 :: Vector I -> Double
    norm_1 = I -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (I -> Double) -> (Vector I -> I) -> Vector I -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector I -> I
    Vector I -> RealOf I
    forall e. Product e => Vector e -> RealOf e
    norm1
        norm_2 :: Vector I -> Double
    norm_2 Vector I
    v = Double -> Double
    forall a. Floating a => a -> a
    sqrt (Double -> Double) -> (I -> Double) -> I -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . I -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (I -> Double) -> I -> Double
    forall a b. (a -> b) -> a -> b
    $ Vector I -> Vector I -> I
    forall t. Numeric t => Vector t -> Vector t -> t
    dot Vector I
    v Vector I
    v
        norm_Inf :: Vector I -> Double
    norm_Inf = I -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (I -> Double) -> (Vector I -> I) -> Vector I -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector I -> I
    Vector I -> RealOf I
    forall e. Product e => Vector e -> RealOf e
    normInf

    instance Normed (Vector Z)
      where
        norm_0 :: Vector Z -> Double
    norm_0 = Z -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Z -> Double) -> (Vector Z -> Z) -> Vector Z -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector Z -> Z
    forall (c :: * -> *) e. Container c e => c e -> e
    sumElements (Vector Z -> Z) -> (Vector Z -> Vector Z) -> Vector Z -> Z
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector Z -> Vector Z
    forall e (c :: * -> *). (Ord e, Container c e) => c e -> c e
    step (Vector Z -> Vector Z)
    -> (Vector Z -> Vector Z) -> Vector Z -> Vector Z
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector Z -> Vector Z
    forall a. Num a => a -> a
    abs
        norm_1 :: Vector Z -> Double
    norm_1 = Z -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Z -> Double) -> (Vector Z -> Z) -> Vector Z -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector Z -> Z
    Vector Z -> RealOf Z
    forall e. Product e => Vector e -> RealOf e
    norm1
        norm_2 :: Vector Z -> Double
    norm_2 Vector Z
    v = Double -> Double
    forall a. Floating a => a -> a
    sqrt (Double -> Double) -> (Z -> Double) -> Z -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Z -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Z -> Double) -> Z -> Double
    forall a b. (a -> b) -> a -> b
    $ Vector Z -> Vector Z -> Z
    forall t. Numeric t => Vector t -> Vector t -> t
    dot Vector Z
    v Vector Z
    v
        norm_Inf :: Vector Z -> Double
    norm_Inf = Z -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Z -> Double) -> (Vector Z -> Z) -> Vector Z -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector Z -> Z
    Vector Z -> RealOf Z
    forall e. Product e => Vector e -> RealOf e
    normInf

    instance Normed (Vector Float)
      where
        norm_0 :: Vector Float -> Double
    norm_0 = Vector Double -> Double
    forall a. Normed a => a -> Double
    norm_0 (Vector Double -> Double)
    -> (Vector Float -> Vector Double) -> Vector Float -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector Float -> Vector Double
    Vector Float -> Vector (DoubleOf Float)
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c t -> c (DoubleOf t)
    forall (c :: * -> *).
    Complexable c =>
    c Float -> c (DoubleOf Float)
    double
        norm_1 :: Vector Float -> Double
    norm_1 = Vector Double -> Double
    forall a. Normed a => a -> Double
    norm_1 (Vector Double -> Double)
    -> (Vector Float -> Vector Double) -> Vector Float -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector Float -> Vector Double
    Vector Float -> Vector (DoubleOf Float)
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c t -> c (DoubleOf t)
    forall (c :: * -> *).
    Complexable c =>
    c Float -> c (DoubleOf Float)
    double
        norm_2 :: Vector Float -> Double
    norm_2 = Vector Double -> Double
    forall a. Normed a => a -> Double
    norm_2 (Vector Double -> Double)
    -> (Vector Float -> Vector Double) -> Vector Float -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector Float -> Vector Double
    Vector Float -> Vector (DoubleOf Float)
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c t -> c (DoubleOf t)
    forall (c :: * -> *).
    Complexable c =>
    c Float -> c (DoubleOf Float)
    double
        norm_Inf :: Vector Float -> Double
    norm_Inf = Vector Double -> Double
    forall a. Normed a => a -> Double
    norm_Inf (Vector Double -> Double)
    -> (Vector Float -> Vector Double) -> Vector Float -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector Float -> Vector Double
    Vector Float -> Vector (DoubleOf Float)
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c t -> c (DoubleOf t)
    forall (c :: * -> *).
    Complexable c =>
    c Float -> c (DoubleOf Float)
    double

    instance Normed (Vector (Complex Float))
      where
        norm_0 :: Vector (Complex Float) -> Double
    norm_0 = Vector (Complex Double) -> Double
    forall a. Normed a => a -> Double
    norm_0 (Vector (Complex Double) -> Double)
    -> (Vector (Complex Float) -> Vector (Complex Double))
    -> Vector (Complex Float)
    -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector (Complex Float) -> Vector (Complex Double)
    Vector (Complex Float) -> Vector (DoubleOf (Complex Float))
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c t -> c (DoubleOf t)
    forall (c :: * -> *).
    Complexable c =>
    c (Complex Float) -> c (DoubleOf (Complex Float))
    double
        norm_1 :: Vector (Complex Float) -> Double
    norm_1 = Vector (Complex Double) -> Double
    forall a. Normed a => a -> Double
    norm_1 (Vector (Complex Double) -> Double)
    -> (Vector (Complex Float) -> Vector (Complex Double))
    -> Vector (Complex Float)
    -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector (Complex Float) -> Vector (Complex Double)
    Vector (Complex Float) -> Vector (DoubleOf (Complex Float))
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c t -> c (DoubleOf t)
    forall (c :: * -> *).
    Complexable c =>
    c (Complex Float) -> c (DoubleOf (Complex Float))
    double
        norm_2 :: Vector (Complex Float) -> Double
    norm_2 = Vector (Complex Double) -> Double
    forall a. Normed a => a -> Double
    norm_2 (Vector (Complex Double) -> Double)
    -> (Vector (Complex Float) -> Vector (Complex Double))
    -> Vector (Complex Float)
    -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector (Complex Float) -> Vector (Complex Double)
    Vector (Complex Float) -> Vector (DoubleOf (Complex Float))
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c t -> c (DoubleOf t)
    forall (c :: * -> *).
    Complexable c =>
    c (Complex Float) -> c (DoubleOf (Complex Float))
    double
        norm_Inf :: Vector (Complex Float) -> Double
    norm_Inf = Vector (Complex Double) -> Double
    forall a. Normed a => a -> Double
    norm_Inf (Vector (Complex Double) -> Double)
    -> (Vector (Complex Float) -> Vector (Complex Double))
    -> Vector (Complex Float)
    -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector (Complex Float) -> Vector (Complex Double)
    Vector (Complex Float) -> Vector (DoubleOf (Complex Float))
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c t -> c (DoubleOf t)
    forall (c :: * -> *).
    Complexable c =>
    c (Complex Float) -> c (DoubleOf (Complex Float))
    double

    -- | Frobenius norm (Schatten p-norm with p=2)
    norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R
    norm_Frob :: forall t. (Normed (Vector t), Element t) => Matrix t -> Double
    norm_Frob = Vector t -> Double
    forall a. Normed a => a -> Double
    norm_2 (Vector t -> Double)
    -> (Matrix t -> Vector t) -> Matrix t -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> Vector t
    forall t. Element t => Matrix t -> Vector t
    flatten

    -- | Sum of singular values (Schatten p-norm with p=1)
    norm_nuclear :: Field t => Matrix t -> R
    norm_nuclear :: forall t. Field t => Matrix t -> Double
    norm_nuclear = Vector Double -> Double
    forall (c :: * -> *) e. Container c e => c e -> e
    sumElements (Vector Double -> Double)
    -> (Matrix t -> Vector Double) -> Matrix t -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> Vector Double
    forall t. Field t => Matrix t -> Vector Double
    singularValues

    {- | Check if the absolute value or complex magnitude is greater than a given threshold

    >>> magnit 1E-6 (1E-12 :: R)
    False
    >>> magnit 1E-6 (3+iC :: C)
    True
    >>> magnit 0 (3 :: I ./. 5)
    True

    -}
    magnit :: (Element t, Normed (Vector t)) => R -> t -> Bool
    magnit :: forall t. (Element t, Normed (Vector t)) => Double -> t -> Bool
    magnit Double
    e t
    x = Vector t -> Double
    forall a. Normed a => a -> Double
    norm_1 ([t] -> Vector t
    forall a. Storable a => [a] -> Vector a
    fromList [t
    x]) Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    > Double
    e

    -- | Obtains a vector in the same direction with 2-norm=1
    normalize :: (Normed (Vector t), Num (Vector t), Field t) => Vector t -> Vector t
    normalize :: forall t.
    (Normed (Vector t), Num (Vector t), Field t) =>
    Vector t -> Vector t
    normalize Vector t
    v = Vector t
    v Vector t -> Vector t -> Vector t
    forall a. Fractional a => a -> a -> a
    / Vector (RealOf t) -> Vector t
    forall t (c :: * -> *).
    (Convert t, Complexable c) =>
    c (RealOf t) -> c t
    forall (c :: * -> *). Complexable c => c (RealOf t) -> c t
    real (Double -> Vector Double
    forall (c :: * -> *) e. Container c e => e -> c e
    scalar (Vector t -> Double
    forall a. Normed a => a -> Double
    norm_2 Vector t
    v))

    -- | trans . inv
    mt :: Matrix Double -> Matrix Double
    mt :: Matrix Double -> Matrix Double
    mt = Matrix Double -> Matrix Double
    forall t. Matrix t -> Matrix t
    trans (Matrix Double -> Matrix Double)
    -> (Matrix Double -> Matrix Double)
    -> Matrix Double
    -> Matrix Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Double -> Matrix Double
    forall t. Field t => Matrix t -> Matrix t
    inv

    --------------------------------------------------------------------------------
    {- |

    >>> size $ vector [1..10]
    10
    >>> size $ (2><5)[1..10::Double]
    (2,5)

    -}
    size :: Container c t => c t -> IndexOf c
    size :: forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    size = c t -> IndexOf c
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    size'

    {- | Alternative indexing function.

    >>> vector [1..10] ! 3
    4.0

    On a matrix it gets the k-th row as a vector:

    >>> matrix 5 [1..15] ! 1
    [6.0,7.0,8.0,9.0,10.0]
    it :: Vector Double

    >>> matrix 5 [1..15] ! 1 ! 3
    9.0

    -}
    class Indexable c t | c -> t , t -> c
      where
        infixl 9 !
        (!) :: c -> Int -> t

    instance Indexable (Vector Double) Double
      where
        ! :: Vector Double -> Int -> Double
    (!) = Vector Double -> Int -> Double
    forall t. Storable t => Vector t -> Int -> t
    (@>)

    instance Indexable (Vector Float) Float
      where
        ! :: Vector Float -> Int -> Float
    (!) = Vector Float -> Int -> Float
    forall t. Storable t => Vector t -> Int -> t
    (@>)

    instance Indexable (Vector I) I
      where
        ! :: Vector I -> Int -> I
    (!) = Vector I -> Int -> I
    forall t. Storable t => Vector t -> Int -> t
    (@>)

    instance Indexable (Vector Z) Z
      where
        ! :: Vector Z -> Int -> Z
    (!) = Vector Z -> Int -> Z
    forall t. Storable t => Vector t -> Int -> t
    (@>)

    instance Indexable (Vector (Complex Double)) (Complex Double)
      where
        ! :: Vector (Complex Double) -> Int -> Complex Double
    (!) = Vector (Complex Double) -> Int -> Complex Double
    forall t. Storable t => Vector t -> Int -> t
    (@>)

    instance Indexable (Vector (Complex Float)) (Complex Float)
      where
        ! :: Vector (Complex Float) -> Int -> Complex Float
    (!) = Vector (Complex Float) -> Int -> Complex Float
    forall t. Storable t => Vector t -> Int -> t
    (@>)

    instance Element t => Indexable (Matrix t) (Vector t)
      where
        Matrix t
    m ! :: Matrix t -> Int -> Vector t
    ! Int
    j = Int -> Int -> Vector t -> Vector t
    forall t. Storable t => Int -> Int -> Vector t -> Vector t
    subVector (Int
    jInt -> Int -> Int
    forall a. Num a => a -> a -> a
    *Int
    c) Int
    c (Matrix t -> Vector t
    forall t. Element t => Matrix t -> Vector t
    flatten Matrix t
    m)
          where
            c :: Int
    c = Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    m

    --------------------------------------------------------------------------------

    -- | Matrix of pairwise squared distances of row vectors
    -- (using the matrix product trick in blog.smola.org)
    pairwiseD2 :: Matrix Double -> Matrix Double -> Matrix Double
    pairwiseD2 :: Matrix Double -> Matrix Double -> Matrix Double
    pairwiseD2 Matrix Double
    x Matrix Double
    y | Bool
    ok = Vector Double
    x2 Vector Double -> Vector Double -> Matrix Double
    forall t. Product t => Vector t -> Vector t -> Matrix t
    `outer` Vector Double
    oy Matrix Double -> Matrix Double -> Matrix Double
    forall a. Num a => a -> a -> a
    + Vector Double
    ox Vector Double -> Vector Double -> Matrix Double
    forall t. Product t => Vector t -> Vector t -> Matrix t
    `outer` Vector Double
    y2 Matrix Double -> Matrix Double -> Matrix Double
    forall a. Num a => a -> a -> a
    - Matrix Double
    2Matrix Double -> Matrix Double -> Matrix Double
    forall a. Num a => a -> a -> a
    * Matrix Double
    x Matrix Double -> Matrix Double -> Matrix Double
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Matrix Double -> Matrix Double
    forall t. Matrix t -> Matrix t
    trans Matrix Double
    y
                   | Bool
    otherwise = String -> Matrix Double
    forall a. HasCallStack => String -> a
    error (String -> Matrix Double) -> String -> Matrix Double
    forall a b. (a -> b) -> a -> b
    $ String
    "pairwiseD2 with different number of columns: "
                                       String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++ (Int, Int) -> String
    forall a. Show a => a -> String
    show (Matrix Double -> IndexOf Matrix
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    size Matrix Double
    x) String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++ String
    ", " String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++ (Int, Int) -> String
    forall a. Show a => a -> String
    show (Matrix Double -> IndexOf Matrix
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    size Matrix Double
    y)
      where
        ox :: Vector Double
    ox = Int -> Vector Double
    forall {e} {d} {c :: * -> *}. (Konst e d c, Num e) => d -> c e
    one (Matrix Double -> Int
    forall t. Matrix t -> Int
    rows Matrix Double
    x)
        oy :: Vector Double
    oy = Int -> Vector Double
    forall {e} {d} {c :: * -> *}. (Konst e d c, Num e) => d -> c e
    one (Matrix Double -> Int
    forall t. Matrix t -> Int
    rows Matrix Double
    y)
        oc :: Vector Double
    oc = Int -> Vector Double
    forall {e} {d} {c :: * -> *}. (Konst e d c, Num e) => d -> c e
    one (Matrix Double -> Int
    forall t. Matrix t -> Int
    cols Matrix Double
    x)
        one :: d -> c e
    one d
    k = e -> d -> c e
    forall e d (c :: * -> *). Konst e d c => e -> d -> c e
    konst e
    1 d
    k
        x2 :: Vector Double
    x2 = Matrix Double
    x Matrix Double -> Matrix Double -> Matrix Double
    forall a. Num a => a -> a -> a
    * Matrix Double
    x Matrix Double -> Vector Double -> Vector Double
    forall t. Product t => Matrix t -> Vector t -> Vector t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Vector Double
    oc
        y2 :: Vector Double
    y2 = Matrix Double
    y Matrix Double -> Matrix Double -> Matrix Double
    forall a. Num a => a -> a -> a
    * Matrix Double
    y Matrix Double -> Vector Double -> Vector Double
    forall t. Product t => Matrix t -> Vector t -> Vector t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Vector Double
    oc
        ok :: Bool
    ok = Matrix Double -> Int
    forall t. Matrix t -> Int
    cols Matrix Double
    x Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    == Matrix Double -> Int
    forall t. Matrix t -> Int
    cols Matrix Double
    y

    --------------------------------------------------------------------------------

    {- | outer products of rows

    >>> a
    (3><2)
     [   1.0,   2.0
     ,  10.0,  20.0
     , 100.0, 200.0 ]
    >>> b
    (3><3)
     [ 1.0, 2.0, 3.0
     , 4.0, 5.0, 6.0
     , 7.0, 8.0, 9.0 ]

    >>> rowOuters a (b ||| 1)
    (3><8)
     [   1.0,   2.0,   3.0,   1.0,    2.0,    4.0,    6.0,   2.0
     ,  40.0,  50.0,  60.0,  10.0,   80.0,  100.0,  120.0,  20.0
     , 700.0, 800.0, 900.0, 100.0, 1400.0, 1600.0, 1800.0, 200.0 ]

    -}
    rowOuters :: Matrix Double -> Matrix Double -> Matrix Double
    rowOuters :: Matrix Double -> Matrix Double -> Matrix Double
    rowOuters Matrix Double
    a Matrix Double
    b = Matrix Double
    a' Matrix Double -> Matrix Double -> Matrix Double
    forall a. Num a => a -> a -> a
    * Matrix Double
    b'
      where
        a' :: Matrix Double
    a' = Matrix Double -> Matrix Double -> Matrix Double
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    kronecker Matrix Double
    a (Int -> Int -> Matrix Double
    ones Int
    1 (Matrix Double -> Int
    forall t. Matrix t -> Int
    cols Matrix Double
    b))
        b' :: Matrix Double
    b' = Matrix Double -> Matrix Double -> Matrix Double
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    kronecker (Int -> Int -> Matrix Double
    ones Int
    1 (Matrix Double -> Int
    forall t. Matrix t -> Int
    cols Matrix Double
    a)) Matrix Double
    b

    --------------------------------------------------------------------------------

    -- | solution of overconstrained homogeneous linear system
    null1 :: Matrix R -> Vector R
    null1 :: Matrix Double -> Vector Double
    null1 = [Vector Double] -> Vector Double
    forall a. HasCallStack => [a] -> a
    last ([Vector Double] -> Vector Double)
    -> (Matrix Double -> [Vector Double])
    -> Matrix Double
    -> Vector Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Double -> [Vector Double]
    forall t. Element t => Matrix t -> [Vector t]
    toColumns (Matrix Double -> [Vector Double])
    -> (Matrix Double -> Matrix Double)
    -> Matrix Double
    -> [Vector Double]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Vector Double, Matrix Double) -> Matrix Double
    forall a b. (a, b) -> b
    snd ((Vector Double, Matrix Double) -> Matrix Double)
    -> (Matrix Double -> (Vector Double, Matrix Double))
    -> Matrix Double
    -> Matrix Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Double -> (Vector Double, Matrix Double)
    forall t. Field t => Matrix t -> (Vector Double, Matrix t)
    rightSV

    -- | solution of overconstrained homogeneous symmetric linear system
    null1sym :: Herm R -> Vector R
    null1sym :: Herm Double -> Vector Double
    null1sym = [Vector Double] -> Vector Double
    forall a. HasCallStack => [a] -> a
    last ([Vector Double] -> Vector Double)
    -> (Herm Double -> [Vector Double]) -> Herm Double -> Vector Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix Double -> [Vector Double]
    forall t. Element t => Matrix t -> [Vector t]
    toColumns (Matrix Double -> [Vector Double])
    -> (Herm Double -> Matrix Double) -> Herm Double -> [Vector Double]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Vector Double, Matrix Double) -> Matrix Double
    forall a b. (a, b) -> b
    snd ((Vector Double, Matrix Double) -> Matrix Double)
    -> (Herm Double -> (Vector Double, Matrix Double))
    -> Herm Double
    -> Matrix Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Herm Double -> (Vector Double, Matrix Double)
    forall t. Field t => Herm t -> (Vector Double, Matrix t)
    eigSH

    --------------------------------------------------------------------------------

    infixl 0 ~!~
    Bool
    c ~!~ :: Bool -> String -> f ()
    ~!~ String
    msg = Bool -> f () -> f ()
    forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
    when Bool
    c (String -> f ()
    forall a. HasCallStack => String -> a
    error String
    msg)

    --------------------------------------------------------------------------------

    formatSparse :: String -> String -> String -> Int -> Matrix Double -> String

    formatSparse :: String -> String -> String -> Int -> Matrix Double -> String
    formatSparse String
    zeroI String
    _zeroF String
    sep Int
    _ (Matrix Double -> Maybe (Matrix Double)
    approxInt -> Just Matrix Double
    m) = String -> (Double -> String) -> Matrix Double -> String
    forall t.
    Element t =>
    String -> (t -> String) -> Matrix t -> String
    format String
    sep Double -> String
    forall {t}. (Eq t, Num t, PrintfArg t) => t -> String
    f Matrix Double
    m
      where
        f :: t -> String
    f t
    0 = String
    zeroI
        f t
    x = String -> t -> String
    forall r. PrintfType r => String -> r
    printf String
    "%.0f" t
    x

    formatSparse String
    zeroI String
    zeroF String
    sep Int
    n Matrix Double
    m = String -> (Double -> String) -> Matrix Double -> String
    forall t.
    Element t =>
    String -> (t -> String) -> Matrix t -> String
    format String
    sep Double -> String
    f Matrix Double
    m
      where
        f :: Double -> String
    f Double
    x | Double -> Double
    forall a. Num a => a -> a
    abs (Double
    x::Double) Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    < Double
    2Double -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    forall x. RealFloat x => x
    peps = String
    zeroIString -> String -> String
    forall a. [a] -> [a] -> [a]
    ++String
    zeroF
            | Double -> Double
    forall a. Num a => a -> a
    abs (Int -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Double -> Int
    forall b. Integral b => Double -> b
    forall a b. (RealFrac a, Integral b) => a -> b
    round Double
    x::Int) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    x) Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double -> Double
    forall a. Num a => a -> a
    abs Double
    x Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    < Double
    2Double -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    forall x. RealFloat x => x
    peps
                = String -> Double -> String
    forall r. PrintfType r => String -> r
    printf (String
    "%.0f."String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++Int -> Char -> String
    forall a. Int -> a -> [a]
    replicate Int
    n Char
    ' ') Double
    x
            | Bool
    otherwise = String -> Double -> String
    forall r. PrintfType r => String -> r
    printf (String
    "%."String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++Int -> String
    forall a. Show a => a -> String
    show Int
    nString -> String -> String
    forall a. [a] -> [a] -> [a]
    ++String
    "f") Double
    x

    approxInt :: Matrix Double -> Maybe (Matrix Double)
    approxInt Matrix Double
    m
        | Vector Double -> Double
    forall a. Normed a => a -> Double
    norm_Inf (Vector Double
    v Vector Double -> Vector Double -> Vector Double
    forall a. Num a => a -> a -> a
    - Vector Double
    vi) Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    < Double
    2Double -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    forall x. RealFloat x => x
    peps Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Vector Double -> Double
    forall a. Normed a => a -> Double
    norm_Inf Vector Double
    v = Matrix Double -> Maybe (Matrix Double)
    forall a. a -> Maybe a
    Just (Int -> Vector Double -> Matrix Double
    forall t. Storable t => Int -> Vector t -> Matrix t
    reshape (Matrix Double -> Int
    forall t. Matrix t -> Int
    cols Matrix Double
    m) Vector Double
    vi)
        | Bool
    otherwise = Maybe (Matrix Double)
    forall a. Maybe a
    Nothing
      where
        v :: Vector Double
    v = Matrix Double -> Vector Double
    forall t. Element t => Matrix t -> Vector t
    flatten Matrix Double
    m
        vi :: Vector Double
    vi = Vector Double -> Vector Double
    roundVector Vector Double
    v

    dispDots :: Int -> Matrix Double -> IO ()
    dispDots Int
    n = String -> IO ()
    putStr (String -> IO ())
    -> (Matrix Double -> String) -> Matrix Double -> IO ()
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . String -> String -> String -> Int -> Matrix Double -> String
    formatSparse String
    "." (Int -> Char -> String
    forall a. Int -> a -> [a]
    replicate Int
    n Char
    ' ') String
    "  " Int
    n

    dispBlanks :: Int -> Matrix Double -> IO ()
    dispBlanks Int
    n = String -> IO ()
    putStr (String -> IO ())
    -> (Matrix Double -> String) -> Matrix Double -> IO ()
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . String -> String -> String -> Int -> Matrix Double -> String
    formatSparse String
    "" String
    "" String
    "  " Int
    n

    formatShort :: String -> (t -> String) -> Int -> Int -> Matrix t -> String
    formatShort String
    sep t -> String
    fmt Int
    maxr Int
    maxc Matrix t
    m = String
    auxm4
      where
        (Int
    rm,Int
    cm) = Matrix t -> IndexOf Matrix
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    size Matrix t
    m
        (Int
    r1,Int
    r2,Int
    r3)
            | Int
    rm Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Int
    maxr = (Int
    rm,Int
    0,Int
    0)
            | Bool
    otherwise  = (Int
    maxrInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    3,Int
    rmInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    maxrInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1,Int
    2)
        (Int
    c1,Int
    c2,Int
    c3)
            | Int
    cm Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Int
    maxc = (Int
    cm,Int
    0,Int
    0)
            | Bool
    otherwise  = (Int
    maxcInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    3,Int
    cmInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    maxcInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1,Int
    2)
        [ [Matrix t
    a,Matrix t
    _,Matrix t
    b]
         ,[Matrix t
    _,Matrix t
    _,Matrix t
    _]
         ,[Matrix t
    c,Matrix t
    _,Matrix t
    d]] = [Int] -> [Int] -> Matrix t -> [[Matrix t]]
    forall t. Element t => [Int] -> [Int] -> Matrix t -> [[Matrix t]]
    toBlocks [Int
    r1,Int
    r2,Int
    r3]
                              [Int
    c1,Int
    c2,Int
    c3] Matrix t
    m
        auxm :: Matrix t
    auxm = [[Matrix t]] -> Matrix t
    forall t. Element t => [[Matrix t]] -> Matrix t
    fromBlocks [[Matrix t
    a,Matrix t
    b],[Matrix t
    c,Matrix t
    d]]
        auxm2 :: String
    auxm2
            | Int
    cm Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    > Int
    maxc = String -> (t -> String) -> Matrix t -> String
    forall t.
    Element t =>
    String -> (t -> String) -> Matrix t -> String
    format String
    "|" t -> String
    fmt Matrix t
    auxm
            | Bool
    otherwise = String -> (t -> String) -> Matrix t -> String
    forall t.
    Element t =>
    String -> (t -> String) -> Matrix t -> String
    format String
    sep t -> String
    fmt Matrix t
    auxm
        auxm3 :: [String]
    auxm3
            | Int
    cm Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    > Int
    maxc = (String -> String) -> [String] -> [String]
    forall a b. (a -> b) -> [a] -> [b]
    map ([String] -> String
    f ([String] -> String) -> (String -> [String]) -> String -> String
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . String -> String -> [String]
    forall a. Eq a => [a] -> [a] -> [[a]]
    splitOn String
    "|") (String -> [String]
    lines String
    auxm2)
            | Bool
    otherwise = (String -> [String]
    lines String
    auxm2)
        f :: [String] -> String
    f [String]
    items = String -> [String] -> String
    forall a. [a] -> [[a]] -> [a]
    intercalate String
    sep (Int -> [String] -> [String]
    forall a. Int -> [a] -> [a]
    take (Int
    maxcInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    3) [String]
    items) String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++ String
    "  .. " String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++
                  String -> [String] -> String
    forall a. [a] -> [[a]] -> [a]
    intercalate String
    sep (Int -> [String] -> [String]
    forall a. Int -> [a] -> [a]
    drop (Int
    maxcInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    3) [String]
    items)
        auxm4 :: String
    auxm4
            | Int
    rm Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    > Int
    maxr = [String] -> String
    unlines (Int -> [String] -> [String]
    forall a. Int -> [a] -> [a]
    take (Int
    maxrInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    3) [String]
    auxm3 [String] -> [String] -> [String]
    forall a. [a] -> [a] -> [a]
    ++ String
    vsep String -> [String] -> [String]
    forall a. a -> [a] -> [a]
    : Int -> [String] -> [String]
    forall a. Int -> [a] -> [a]
    drop (Int
    maxrInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    3) [String]
    auxm3)
            | Bool
    otherwise = [String] -> String
    unlines [String]
    auxm3
        vsep :: String
    vsep = (Char -> Char) -> String -> String
    forall a b. (a -> b) -> [a] -> [b]
    map Char -> Char
    g ([String] -> String
    forall a. HasCallStack => [a] -> a
    head [String]
    auxm3)
        g :: Char -> Char
    g Char
    '.' = Char
    ':'
        g Char
    _ = Char
    ' '

    dispShort :: Int -> Int -> Int -> Matrix Double -> IO ()
    dispShort :: Int -> Int -> Int -> Matrix Double -> IO ()
    dispShort Int
    maxr Int
    maxc Int
    dec Matrix Double
    m =
        String -> Int -> Int -> String -> IO ()
    forall r. PrintfType r => String -> r
    printf String
    "%dx%d\n%s" (Matrix Double -> Int
    forall t. Matrix t -> Int
    rows Matrix Double
    m) (Matrix Double -> Int
    forall t. Matrix t -> Int
    cols Matrix Double
    m) (String
    -> (Double -> String) -> Int -> Int -> Matrix Double -> String
    forall {t}.
    (Num t, Container Vector t) =>
    String -> (t -> String) -> Int -> Int -> Matrix t -> String
    formatShort String
    "  " Double -> String
    fmt Int
    maxr Int
    maxc Matrix Double
    m)
      where
        fmt :: Double -> String
    fmt = String -> Double -> String
    forall r. PrintfType r => String -> r
    printf (String
    "%."String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++Int -> String
    forall a. Show a => a -> String
    show Int
    dec String -> String -> String
    forall a. [a] -> [a] -> [a]
    ++String
    "f")

    --------------------------------------------------------------------------------

    -- matrix views

    block2x2 :: Int -> Int -> Matrix t -> [[Matrix t]]
    block2x2 Int
    r Int
    c Matrix t
    m = [[Matrix t
    m11,Matrix t
    m12],[Matrix t
    m21,Matrix t
    m22]]
      where
        m11 :: Matrix t
    m11 = Matrix t
    m Matrix t -> (Extractor, Extractor) -> Matrix t
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Int -> Extractor
    Take Int
    r, Int -> Extractor
    Take Int
    c)
        m12 :: Matrix t
    m12 = Matrix t
    m Matrix t -> (Extractor, Extractor) -> Matrix t
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Int -> Extractor
    Take Int
    r, Int -> Extractor
    Drop Int
    c)
        m21 :: Matrix t
    m21 = Matrix t
    m Matrix t -> (Extractor, Extractor) -> Matrix t
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Int -> Extractor
    Drop Int
    r, Int -> Extractor
    Take Int
    c)
        m22 :: Matrix t
    m22 = Matrix t
    m Matrix t -> (Extractor, Extractor) -> Matrix t
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Int -> Extractor
    Drop Int
    r, Int -> Extractor
    Drop Int
    c)

    block3x3 :: Int -> Int -> Int -> Int -> Matrix t -> [[Matrix t]]
    block3x3 Int
    r Int
    nr Int
    c Int
    nc Matrix t
    m = [[Matrix t
    m Matrix t -> (Extractor, Extractor) -> Matrix t
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? ([Extractor]
    er [Extractor] -> Int -> Extractor
    forall a. HasCallStack => [a] -> Int -> a
    !! Int
    i, [Extractor]
    ec [Extractor] -> Int -> Extractor
    forall a. HasCallStack => [a] -> Int -> a
    !! Int
    j) | Int
    j <- [Int
    0..Int
    2] ] | Int
    i <- [Int
    0..Int
    2] ]
      where
        er :: [Extractor]
    er = [ Int -> Int -> Int -> Extractor
    Range Int
    0 Int
    1 (Int
    rInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1), Int -> Int -> Int -> Extractor
    Range Int
    r Int
    1 (Int
    rInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    nrInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1), Int -> Extractor
    Drop (Int
    nrInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    r) ]
        ec :: [Extractor]
    ec = [ Int -> Int -> Int -> Extractor
    Range Int
    0 Int
    1 (Int
    cInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1), Int -> Int -> Int -> Extractor
    Range Int
    c Int
    1 (Int
    cInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    ncInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1), Int -> Extractor
    Drop (Int
    ncInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    c) ]

    view1 :: Numeric t => Matrix t -> Maybe (View1 t)
    view1 :: forall t. Numeric t => Matrix t -> Maybe (View1 t)
    view1 Matrix t
    m
        | Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    > Int
    0 Bool -> Bool -> Bool
    && Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    m Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    > Int
    0 = View1 t -> Maybe (View1 t)
    forall a. a -> Maybe a
    Just (t
    e, Matrix t -> Vector t
    forall t. Element t => Matrix t -> Vector t
    flatten Matrix t
    m12, Matrix t -> Vector t
    forall t. Element t => Matrix t -> Vector t
    flatten Matrix t
    m21 , Matrix t
    m22)
        | Bool
    otherwise = Maybe (View1 t)
    forall a. Maybe a
    Nothing
      where
        [[Matrix t
    m11,Matrix t
    m12],[Matrix t
    m21,Matrix t
    m22]] = Int -> Int -> Matrix t -> [[Matrix t]]
    forall {t}. Element t => Int -> Int -> Matrix t -> [[Matrix t]]
    block2x2 Int
    1 Int
    1 Matrix t
    m
        e :: t
    e = Matrix t
    m11 Matrix t -> IndexOf Matrix -> t
    forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
    `atIndex` (Int
    0, Int
    0)

    unView1 :: Numeric t => View1 t -> Matrix t
    unView1 :: forall t. Numeric t => View1 t -> Matrix t
    unView1 (t
    e,Vector t
    r,Vector t
    c,Matrix t
    m) = [[Matrix t]] -> Matrix t
    forall t. Element t => [[Matrix t]] -> Matrix t
    fromBlocks [[t -> Matrix t
    forall (c :: * -> *) e. Container c e => e -> c e
    scalar t
    e, Vector t -> Matrix t
    forall a. Storable a => Vector a -> Matrix a
    asRow Vector t
    r],[Vector t -> Matrix t
    forall a. Storable a => Vector a -> Matrix a
    asColumn Vector t
    c, Matrix t
    m]]

    type View1 t = (t, Vector t, Vector t, Matrix t)

    foldMatrix :: Numeric t => (Matrix t -> Matrix t) -> (View1 t -> View1 t) -> (Matrix t -> Matrix t)
    foldMatrix :: forall t.
    Numeric t =>
    (Matrix t -> Matrix t)
    -> (View1 t -> View1 t) -> Matrix t -> Matrix t
    foldMatrix Matrix t -> Matrix t
    g View1 t -> View1 t
    f ( (View1 t -> View1 t
    f (View1 t -> View1 t) -> Maybe (View1 t) -> Maybe (View1 t)
    forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
    <$>) (Maybe (View1 t) -> Maybe (View1 t))
    -> (Matrix t -> Maybe (View1 t)) -> Matrix t -> Maybe (View1 t)
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> Maybe (View1 t)
    forall t. Numeric t => Matrix t -> Maybe (View1 t)
    view1 (Matrix t -> Maybe (View1 t))
    -> (Matrix t -> Matrix t) -> Matrix t -> Maybe (View1 t)
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> Matrix t
    g -> Just (t
    e,Vector t
    r,Vector t
    c,Matrix t
    m)) = View1 t -> Matrix t
    forall t. Numeric t => View1 t -> Matrix t
    unView1 (t
    e, Vector t
    r, Vector t
    c, (Matrix t -> Matrix t)
    -> (View1 t -> View1 t) -> Matrix t -> Matrix t
    forall t.
    Numeric t =>
    (Matrix t -> Matrix t)
    -> (View1 t -> View1 t) -> Matrix t -> Matrix t
    foldMatrix Matrix t -> Matrix t
    g View1 t -> View1 t
    f Matrix t
    m)
    foldMatrix Matrix t -> Matrix t
    _ View1 t -> View1 t
    _ Matrix t
    m = Matrix t
    m

    swapMax :: Int -> Matrix e -> (IndexOf Vector, Matrix e)
    swapMax Int
    k Matrix e
    m
        | Matrix e -> Int
    forall t. Matrix t -> Int
    rows Matrix e
    m Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    > Int
    0 Bool -> Bool -> Bool
    && IndexOf Vector
    jIndexOf Vector -> IndexOf Vector -> Bool
    forall a. Ord a => a -> a -> Bool
    >IndexOf Vector
    0 = (IndexOf Vector
    j, Matrix e
    m Matrix e -> (Extractor, Extractor) -> Matrix e
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Vector I -> Extractor
    Pos ([Int] -> Vector I
    idxs [Int]
    swapped), Extractor
    All))
        | Bool
    otherwise  = (IndexOf Vector
    0,Matrix e
    m)
      where
        j :: IndexOf Vector
    j = Vector e -> IndexOf Vector
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    maxIndex (Vector e -> IndexOf Vector) -> Vector e -> IndexOf Vector
    forall a b. (a -> b) -> a -> b
    $ Vector e -> Vector e
    forall a. Num a => a -> a
    abs (Matrix e -> Matrix e
    forall m mt. Transposable m mt => m -> mt
    tr Matrix e
    m Matrix e -> Int -> Vector e
    forall c t. Indexable c t => c -> Int -> t
    ! Int
    k)
        swapped :: [Int]
    swapped = Int
    IndexOf Vector
    jInt -> [Int] -> [Int]
    forall a. a -> [a] -> [a]
    :[Int
    1..Int
    IndexOf Vector
    jInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1] [Int] -> [Int] -> [Int]
    forall a. [a] -> [a] -> [a]
    ++ Int
    0Int -> [Int] -> [Int]
    forall a. a -> [a] -> [a]
    :[Int
    IndexOf Vector
    jInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1..Matrix e -> Int
    forall t. Matrix t -> Int
    rows Matrix e
    mInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1]

    down :: (Matrix t -> Matrix t) -> Matrix t -> Matrix t
    down Matrix t -> Matrix t
    g Matrix t
    a = (Matrix t -> Matrix t)
    -> (View1 t -> View1 t) -> Matrix t -> Matrix t
    forall t.
    Numeric t =>
    (Matrix t -> Matrix t)
    -> (View1 t -> View1 t) -> Matrix t -> Matrix t
    foldMatrix Matrix t -> Matrix t
    g View1 t -> View1 t
    forall {t} {a} {c}.
    (Eq t, Product t, Fractional t, Container Vector t, Num a, Num c,
     Num (Vector t)) =>
    (t, Vector t, Vector t, Matrix t) -> (a, Vector t, c, Matrix t)
    f Matrix t
    a
      where
        f :: (t, Vector t, Vector t, Matrix t) -> (a, Vector t, c, Matrix t)
    f (t
    e,Vector t
    r,Vector t
    c,Matrix t
    m)
            | t
    e t -> t -> Bool
    forall a. Eq a => a -> a -> Bool
    /= t
    0    = (a
    1, Vector t
    r', c
    0, Matrix t
    m Matrix t -> Matrix t -> Matrix t
    forall a. Num a => a -> a -> a
    - Vector t -> Vector t -> Matrix t
    forall t. Product t => Vector t -> Vector t -> Matrix t
    outer Vector t
    c Vector t
    r')
            | Bool
    otherwise = String -> (a, Vector t, c, Matrix t)
    forall a. HasCallStack => String -> a
    error String
    "singular!"
          where
            r' :: Vector t
    r' = Vector t
    r Vector t -> Vector t -> Vector t
    forall a. Fractional a => a -> a -> a
    / t -> Vector t
    forall (c :: * -> *) e. Container c e => e -> c e
    scalar t
    e

    -- | generic reference implementation of gaussian elimination
    --
    -- @a <> gaussElim a b = b@
    --
    gaussElim_2
      :: (Eq t, Fractional t, Num (Vector t), Numeric t)
      => Matrix t -> Matrix t -> Matrix t

    gaussElim_2 :: forall t.
    (Eq t, Fractional t, Num (Vector t), Numeric t) =>
    Matrix t -> Matrix t -> Matrix t
    gaussElim_2 Matrix t
    a Matrix t
    b = Matrix t -> Matrix t
    flipudrl Matrix t
    r
      where
        flipudrl :: Matrix t -> Matrix t
    flipudrl = Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t
    flipud (Matrix t -> Matrix t)
    -> (Matrix t -> Matrix t) -> Matrix t -> Matrix t
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t
    fliprl
        splitColsAt :: Int -> Matrix t -> (Matrix t, Matrix t)
    splitColsAt Int
    n = (Int -> Matrix t -> Matrix t
    forall t. Element t => Int -> Matrix t -> Matrix t
    takeColumns Int
    n (Matrix t -> Matrix t)
    -> (Matrix t -> Matrix t) -> Matrix t -> (Matrix t, Matrix t)
    forall b c c'. (b -> c) -> (b -> c') -> b -> (c, c')
    forall (a :: * -> * -> *) b c c'.
    Arrow a =>
    a b c -> a b c' -> a b (c, c')
    &&& Int -> Matrix t -> Matrix t
    forall t. Element t => Int -> Matrix t -> Matrix t
    dropColumns Int
    n)
        go :: (Matrix t -> Matrix t)
    -> Matrix t -> Matrix t -> (Matrix t, Matrix t)
    go Matrix t -> Matrix t
    f Matrix t
    x Matrix t
    y = Int -> Matrix t -> (Matrix t, Matrix t)
    forall {t}. Element t => Int -> Matrix t -> (Matrix t, Matrix t)
    splitColsAt (Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    a) ((Matrix t -> Matrix t) -> Matrix t -> Matrix t
    forall {t}.
    (Numeric t, Eq t, Fractional t, Num (Vector t)) =>
    (Matrix t -> Matrix t) -> Matrix t -> Matrix t
    down Matrix t -> Matrix t
    f (Matrix t -> Matrix t) -> Matrix t -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ Matrix t
    x Matrix t -> Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t -> Matrix t
    ||| Matrix t
    y)
        (Matrix t
    a1,Matrix t
    b1) = (Matrix t -> Matrix t)
    -> Matrix t -> Matrix t -> (Matrix t, Matrix t)
    forall {t}.
    (Numeric t, Eq t, Fractional t, Num (Vector t)) =>
    (Matrix t -> Matrix t)
    -> Matrix t -> Matrix t -> (Matrix t, Matrix t)
    go ((Int, Matrix t) -> Matrix t
    forall a b. (a, b) -> b
    snd ((Int, Matrix t) -> Matrix t)
    -> (Matrix t -> (Int, Matrix t)) -> Matrix t -> Matrix t
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Int -> Matrix t -> (Int, Matrix t)
    forall {e}.
    (Container Vector e, Num (Vector e), CTrans e) =>
    Int -> Matrix e -> (Int, Matrix e)
    swapMax Int
    0) Matrix t
    a Matrix t
    b
        ( Matrix t
    _, Matrix t
    r) = (Matrix t -> Matrix t)
    -> Matrix t -> Matrix t -> (Matrix t, Matrix t)
    forall {t}.
    (Numeric t, Eq t, Fractional t, Num (Vector t)) =>
    (Matrix t -> Matrix t)
    -> Matrix t -> Matrix t -> (Matrix t, Matrix t)
    go Matrix t -> Matrix t
    forall a. a -> a
    id (Matrix t -> Matrix t
    flipudrl (Matrix t -> Matrix t) -> Matrix t -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ Matrix t
    a1) (Matrix t -> Matrix t
    flipudrl (Matrix t -> Matrix t) -> Matrix t -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ Matrix t
    b1)

    --------------------------------------------------------------------------------

    gaussElim_1
      :: (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)
      => Matrix t -> Matrix t -> Matrix t

    gaussElim_1 :: forall t.
    (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
     Numeric t) =>
    Matrix t -> Matrix t -> Matrix t
    gaussElim_1 Matrix t
    x Matrix t
    y = Int -> Matrix t -> Matrix t
    forall t. Element t => Int -> Matrix t -> Matrix t
    dropColumns (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    x) (Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t
    flipud (Matrix t -> Matrix t) -> Matrix t -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ [Vector t] -> Matrix t
    forall t. Element t => [Vector t] -> Matrix t
    fromRows [Vector t]
    s2)
      where
        rs :: [Vector t]
    rs = Matrix t -> [Vector t]
    forall t. Element t => Matrix t -> [Vector t]
    toRows (Matrix t -> [Vector t]) -> Matrix t -> [Vector t]
    forall a b. (a -> b) -> a -> b
    $ Matrix t
    x Matrix t -> Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t -> Matrix t
    ||| Matrix t
    y
        s1 :: Matrix t
    s1 = [Vector t] -> Matrix t
    forall t. Element t => [Vector t] -> Matrix t
    fromRows ([Vector t] -> Matrix t) -> [Vector t] -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ Int -> Int -> [Vector t] -> [Vector t]
    forall t.
    (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
     Numeric t) =>
    Int -> Int -> [Vector t] -> [Vector t]
    pivotDown (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    x) Int
    0 [Vector t]
    rs      -- interesting
        s2 :: [Vector t]
    s2 = Int -> [Vector t] -> [Vector t]
    forall t.
    (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
     Numeric t) =>
    Int -> [Vector t] -> [Vector t]
    pivotUp (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    xInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1) (Matrix t -> [Vector t]
    forall t. Element t => Matrix t -> [Vector t]
    toRows (Matrix t -> [Vector t]) -> Matrix t -> [Vector t]
    forall a b. (a -> b) -> a -> b
    $ Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t
    flipud Matrix t
    s1)

    pivotDown
      :: forall t . (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)
      => Int -> Int -> [Vector t] -> [Vector t]
    pivotDown :: forall t.
    (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
     Numeric t) =>
    Int -> Int -> [Vector t] -> [Vector t]
    pivotDown Int
    t Int
    n [Vector t]
    xs
        | Int
    t Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    == Int
    n    = []
        | Bool
    otherwise = Vector t
    y Vector t -> [Vector t] -> [Vector t]
    forall a. a -> [a] -> [a]
    : Int -> Int -> [Vector t] -> [Vector t]
    forall t.
    (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
     Numeric t) =>
    Int -> Int -> [Vector t] -> [Vector t]
    pivotDown Int
    t (Int
    nInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1) [Vector t]
    ys
      where
        Vector t
    y:[Vector t]
    ys = (Int, [Vector t]) -> [Vector t]
    redu (Int -> [Vector t] -> (Int, [Vector t])
    forall {b} {a}.
    (Ord b, Num b, Indexable a b) =>
    Int -> [a] -> (Int, [a])
    pivot Int
    n [Vector t]
    xs)

        pivot :: Int -> [a] -> (Int, [a])
    pivot Int
    k = (Int -> [a] -> Int
    forall a b. a -> b -> a
    const Int
    k ([a] -> Int) -> ([a] -> [a]) -> [a] -> (Int, [a])
    forall b c c'. (b -> c) -> (b -> c') -> b -> (c, c')
    forall (a :: * -> * -> *) b c c'.
    Arrow a =>
    a b c -> a b c' -> a b (c, c')
    &&& [a] -> [a]
    forall a. a -> a
    id)
                ([a] -> (Int, [a])) -> ([a] -> [a]) -> [a] -> (Int, [a])
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (a -> a -> Ordering) -> [a] -> [a]
    forall a. (a -> a -> Ordering) -> [a] -> [a]
    sortBy ((b -> b -> Ordering) -> b -> b -> Ordering
    forall a b c. (a -> b -> c) -> b -> a -> c
    flip b -> b -> Ordering
    forall a. Ord a => a -> a -> Ordering
    compare (b -> b -> Ordering) -> (a -> b) -> a -> a -> Ordering
    forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
    `on` (b -> b
    forall a. Num a => a -> a
    abs(b -> b) -> (a -> b) -> a -> b
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (a -> Int -> b
    forall c t. Indexable c t => c -> Int -> t
    ! Int
    k)))

        redu :: (Int, [Vector t]) -> [Vector t]
        redu :: (Int, [Vector t]) -> [Vector t]
    redu (Int
    k,Vector t
    x:[Vector t]
    zs)
            | t
    p t -> t -> Bool
    forall a. Eq a => a -> a -> Bool
    == t
    0 = String -> [Vector t]
    forall a. HasCallStack => String -> a
    error String
    "gauss: singular!"  -- FIXME
            | Bool
    otherwise = Vector t
    u Vector t -> [Vector t] -> [Vector t]
    forall a. a -> [a] -> [a]
    : (Vector t -> Vector t) -> [Vector t] -> [Vector t]
    forall a b. (a -> b) -> [a] -> [b]
    map Vector t -> Vector t
    f [Vector t]
    zs
          where
            p :: t
    p = Vector t
    xVector t -> Int -> t
    forall c t. Indexable c t => c -> Int -> t
    !Int
    k
            u :: Vector t
    u = t -> Vector t -> Vector t
    forall t (c :: * -> *). Linear t c => t -> c t -> c t
    scale (t -> t
    forall a. Fractional a => a -> a
    recip (Vector t
    xVector t -> Int -> t
    forall c t. Indexable c t => c -> Int -> t
    !Int
    k)) Vector t
    x
            f :: Vector t -> Vector t
    f Vector t
    z = Vector t
    z Vector t -> Vector t -> Vector t
    forall a. Num a => a -> a -> a
    - t -> Vector t -> Vector t
    forall t (c :: * -> *). Linear t c => t -> c t -> c t
    scale (Vector t
    zVector t -> Int -> t
    forall c t. Indexable c t => c -> Int -> t
    !Int
    k) Vector t
    u
        redu (Int
    _,[]) = []

    pivotUp
      :: forall t . (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)
      => Int -> [Vector t] -> [Vector t]
    pivotUp :: forall t.
    (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
     Numeric t) =>
    Int -> [Vector t] -> [Vector t]
    pivotUp Int
    n [Vector t]
    xs
        | Int
    n Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    == -Int
    1 = []
        | Bool
    otherwise = Vector t
    y Vector t -> [Vector t] -> [Vector t]
    forall a. a -> [a] -> [a]
    : Int -> [Vector t] -> [Vector t]
    forall t.
    (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t,
     Numeric t) =>
    Int -> [Vector t] -> [Vector t]
    pivotUp (Int
    nInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1) [Vector t]
    ys
      where
        Vector t
    y:[Vector t]
    ys = (Int, [Vector t]) -> [Vector t]
    redu' (Int
    n,[Vector t]
    xs)

        redu' :: (Int, [Vector t]) -> [Vector t]
        redu' :: (Int, [Vector t]) -> [Vector t]
    redu' (Int
    k,Vector t
    x:[Vector t]
    zs) = Vector t
    u Vector t -> [Vector t] -> [Vector t]
    forall a. a -> [a] -> [a]
    : (Vector t -> Vector t) -> [Vector t] -> [Vector t]
    forall a b. (a -> b) -> [a] -> [b]
    map Vector t -> Vector t
    f [Vector t]
    zs
          where
            u :: Vector t
    u = Vector t
    x
            f :: Vector t -> Vector t
    f Vector t
    z = Vector t
    z Vector t -> Vector t -> Vector t
    forall a. Num a => a -> a -> a
    - t -> Vector t -> Vector t
    forall t (c :: * -> *). Linear t c => t -> c t -> c t
    scale (Vector t
    zVector t -> Int -> t
    forall c t. Indexable c t => c -> Int -> t
    !Int
    k) Vector t
    u
        redu' (Int
    _,[]) = []

    --------------------------------------------------------------------------------

    gaussElim :: Matrix t -> Matrix t -> Matrix t
    gaussElim Matrix t
    a Matrix t
    b = Int -> Matrix t -> Matrix t
    forall t. Element t => Int -> Matrix t -> Matrix t
    dropColumns (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    a) (Matrix t -> Matrix t) -> Matrix t -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ (Matrix t, ()) -> Matrix t
    forall a b. (a, b) -> a
    fst ((Matrix t, ()) -> Matrix t) -> (Matrix t, ()) -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ (forall s. (Int, Int) -> STMatrix s t -> ST s ())
    -> Matrix t -> (Matrix t, ())
    forall t u.
    Element t =>
    (forall s. (Int, Int) -> STMatrix s t -> ST s u)
    -> Matrix t -> (Matrix t, u)
    mutable (Int, Int) -> STMatrix s t -> ST s ()
    forall s. (Int, Int) -> STMatrix s t -> ST s ()
    forall {t} {b} {s}.
    (Container Vector t, Num (Vector t), Eq t, Fractional t) =>
    (Int, b) -> STMatrix s t -> ST s ()
    gaussST (Matrix t
    a Matrix t -> Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t -> Matrix t
    ||| Matrix t
    b)

    gaussST :: (Int, b) -> STMatrix s t -> ST s ()
    gaussST (Int
    r,b
    _) STMatrix s t
    x = do
        let n :: Int
    n = Int
    rInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1
            axpy :: STMatrix s t -> t -> Int -> Int -> ST s ()
    axpy STMatrix s t
    m t
    a Int
    i Int
    j = RowOper t -> STMatrix s t -> ST s ()
    forall t s.
    (Num t, Element t) =>
    RowOper t -> STMatrix s t -> ST s ()
    rowOper (t -> Int -> Int -> ColRange -> RowOper t
    forall t. t -> Int -> Int -> ColRange -> RowOper t
    AXPY t
    a Int
    i Int
    j ColRange
    AllCols)     STMatrix s t
    m
            swap :: STMatrix s t -> Int -> Int -> ST s ()
    swap STMatrix s t
    m Int
    i Int
    j   = RowOper t -> STMatrix s t -> ST s ()
    forall t s.
    (Num t, Element t) =>
    RowOper t -> STMatrix s t -> ST s ()
    rowOper (Int -> Int -> ColRange -> RowOper t
    forall t. Int -> Int -> ColRange -> RowOper t
    SWAP Int
    i Int
    j ColRange
    AllCols)       STMatrix s t
    m
            scal :: STMatrix s t -> t -> Int -> ST s ()
    scal STMatrix s t
    m t
    a Int
    i   = RowOper t -> STMatrix s t -> ST s ()
    forall t s.
    (Num t, Element t) =>
    RowOper t -> STMatrix s t -> ST s ()
    rowOper (t -> RowRange -> ColRange -> RowOper t
    forall t. t -> RowRange -> ColRange -> RowOper t
    SCAL t
    a (Int -> RowRange
    Row Int
    i) ColRange
    AllCols) STMatrix s t
    m
        [Int] -> (Int -> ST s ()) -> ST s ()
    forall (t :: * -> *) (m :: * -> *) a b.
    (Foldable t, Monad m) =>
    t a -> (a -> m b) -> m ()
    forM_ [Int
    0..Int
    n] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
    forall a b. (a -> b) -> a -> b
    $ \Int
    i -> do
            c <- Vector t -> Int
    Vector t -> IndexOf Vector
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    maxIndex (Vector t -> Int) -> (Matrix t -> Vector t) -> Matrix t -> Int
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector t -> Vector t
    forall a. Num a => a -> a
    abs (Vector t -> Vector t)
    -> (Matrix t -> Vector t) -> Matrix t -> Vector t
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> Vector t
    forall t. Element t => Matrix t -> Vector t
    flatten (Matrix t -> Int) -> ST s (Matrix t) -> ST s Int
    forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
    <$> STMatrix s t -> RowRange -> ColRange -> ST s (Matrix t)
    forall a t s.
    Element a =>
    STMatrix t a -> RowRange -> ColRange -> ST s (Matrix a)
    extractMatrix STMatrix s t
    x (Int -> RowRange
    FromRow Int
    i) (Int -> ColRange
    Col Int
    i)
            swap x i (i+c)
            a <- readMatrix x i i
            when (a == 0) $ error "singular!"
            scal x (recip a) i
            forM_ [i+1..n] $ \Int
    j -> do
                b <- STMatrix s t -> Int -> Int -> ST s t
    forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
    readMatrix STMatrix s t
    x Int
    j Int
    i
                axpy x (-b) i j
        [Int] -> (Int -> ST s ()) -> ST s ()
    forall (t :: * -> *) (m :: * -> *) a b.
    (Foldable t, Monad m) =>
    t a -> (a -> m b) -> m ()
    forM_ [Int
    n,Int
    nInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1..Int
    1] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
    forall a b. (a -> b) -> a -> b
    $ \Int
    i -> do
            [Int] -> (Int -> ST s ()) -> ST s ()
    forall (t :: * -> *) (m :: * -> *) a b.
    (Foldable t, Monad m) =>
    t a -> (a -> m b) -> m ()
    forM_ [Int
    iInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1,Int
    iInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    2..Int
    0] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
    forall a b. (a -> b) -> a -> b
    $ \Int
    j -> do
                b <- STMatrix s t -> Int -> Int -> ST s t
    forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
    readMatrix STMatrix s t
    x Int
    j Int
    i
                axpy x (-b) i j

    luST :: (t -> Bool) -> (Int, b) -> STMatrix s t -> ST s [Int]
    luST t -> Bool
    ok (Int
    r,b
    _) STMatrix s t
    x = do
        let axpy :: STMatrix s t -> t -> Int -> Int -> ST s ()
    axpy STMatrix s t
    m t
    a Int
    i Int
    j = RowOper t -> STMatrix s t -> ST s ()
    forall t s.
    (Num t, Element t) =>
    RowOper t -> STMatrix s t -> ST s ()
    rowOper (t -> Int -> Int -> ColRange -> RowOper t
    forall t. t -> Int -> Int -> ColRange -> RowOper t
    AXPY t
    a Int
    i Int
    j (Int -> ColRange
    FromCol (Int
    iInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1))) STMatrix s t
    m
            swap :: STMatrix s t -> Int -> Int -> ST s ()
    swap STMatrix s t
    m Int
    i Int
    j   = RowOper t -> STMatrix s t -> ST s ()
    forall t s.
    (Num t, Element t) =>
    RowOper t -> STMatrix s t -> ST s ()
    rowOper (Int -> Int -> ColRange -> RowOper t
    forall t. Int -> Int -> ColRange -> RowOper t
    SWAP Int
    i Int
    j ColRange
    AllCols)           STMatrix s t
    m
        p <- Int -> ST s (STVector s Int)
    forall t s. Storable t => Int -> ST s (STVector s t)
    newUndefinedVector Int
    r
        forM_ [0..r-1] $ \Int
    i -> do
            k <- Vector t -> Int
    Vector t -> IndexOf Vector
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    maxIndex (Vector t -> Int) -> (Matrix t -> Vector t) -> Matrix t -> Int
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector t -> Vector t
    forall a. Num a => a -> a
    abs (Vector t -> Vector t)
    -> (Matrix t -> Vector t) -> Matrix t -> Vector t
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix t -> Vector t
    forall t. Element t => Matrix t -> Vector t
    flatten (Matrix t -> Int) -> ST s (Matrix t) -> ST s Int
    forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
    <$> STMatrix s t -> RowRange -> ColRange -> ST s (Matrix t)
    forall a t s.
    Element a =>
    STMatrix t a -> RowRange -> ColRange -> ST s (Matrix a)
    extractMatrix STMatrix s t
    x (Int -> RowRange
    FromRow Int
    i) (Int -> ColRange
    Col Int
    i)
            writeVector p i (k+i)
            swap x i (i+k)
            a <- readMatrix x i i
            when (ok a) $ do
                forM_ [i+1..r-1] $ \Int
    j -> do
                    b <- (t -> t -> t
    forall a. Fractional a => a -> a -> a
    /t
    a) (t -> t) -> ST s t -> ST s t
    forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
    <$> STMatrix s t -> Int -> Int -> ST s t
    forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
    readMatrix STMatrix s t
    x Int
    j Int
    i
                    axpy x (-b) i j
                    writeMatrix x j i b
        v <- unsafeFreezeVector p
        return (toList v)

    {- | Experimental implementation of 'luPacked'
         for any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.

    >>> let m = ident 5 + (5><5) [0..] :: Matrix (Z ./. 17)
    (5><5)
     [  1,  1,  2,  3,  4
     ,  5,  7,  7,  8,  9
     , 10, 11, 13, 13, 14
     , 15, 16,  0,  2,  2
     ,  3,  4,  5,  6,  8 ]

    >>> let (l,u,p,s) = luFact $ luPacked' m
    >>> l
    (5><5)
     [  1,  0, 0,  0, 0
     ,  6,  1, 0,  0, 0
     , 12,  7, 1,  0, 0
     ,  7, 10, 7,  1, 0
     ,  8,  2, 6, 11, 1 ]
    >>> u
    (5><5)
     [ 15, 16,  0,  2,  2
     ,  0, 13,  7, 13, 14
     ,  0,  0, 15,  0, 11
     ,  0,  0,  0, 15, 15
     ,  0,  0,  0,  0,  1 ]

    -}
    luPacked' :: Matrix t -> LU t
    luPacked' Matrix t
    x = Matrix t -> [Int] -> LU t
    forall t. Matrix t -> [Int] -> LU t
    LU Matrix t
    m [Int]
    p
      where
        (Matrix t
    m,[Int]
    p) = (forall s. (Int, Int) -> STMatrix s t -> ST s [Int])
    -> Matrix t -> (Matrix t, [Int])
    forall t u.
    Element t =>
    (forall s. (Int, Int) -> STMatrix s t -> ST s u)
    -> Matrix t -> (Matrix t, u)
    mutable ((t -> Bool) -> (Int, Int) -> STMatrix s t -> ST s [Int]
    forall {t} {b} {s}.
    (Container Vector t, Num (Vector t), Fractional t) =>
    (t -> Bool) -> (Int, b) -> STMatrix s t -> ST s [Int]
    luST (Double -> t -> Bool
    forall t. (Element t, Normed (Vector t)) => Double -> t -> Bool
    magnit Double
    0)) Matrix t
    x

    --------------------------------------------------------------------------------

    scalS :: t -> Slice s t -> ST s ()
    scalS t
    a (Slice STMatrix s t
    x Int
    r0 Int
    c0 Int
    nr Int
    nc) = RowOper t -> STMatrix s t -> ST s ()
    forall t s.
    (Num t, Element t) =>
    RowOper t -> STMatrix s t -> ST s ()
    rowOper (t -> RowRange -> ColRange -> RowOper t
    forall t. t -> RowRange -> ColRange -> RowOper t
    SCAL t
    a (Int -> Int -> RowRange
    RowRange Int
    r0 (Int
    r0Int -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    nrInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1)) (Int -> Int -> ColRange
    ColRange Int
    c0 (Int
    c0Int -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    ncInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1))) STMatrix s t
    x

    view :: STMatrix s a
    -> Int -> Int -> ST s (a, Slice s a, Slice s a, Slice s a)
    view STMatrix s a
    x Int
    k Int
    r = do
        d <- STMatrix s a -> Int -> Int -> ST s a
    forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
    readMatrix STMatrix s a
    x Int
    k Int
    k
        let rr = Int
    rInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1Int -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    k
            o  = if Int
    k Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    < Int
    rInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1 then Int
    1 else Int
    0
            s = STMatrix s a -> Int -> Int -> Int -> Int -> Slice s a
    forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    Slice STMatrix s a
    x (Int
    kInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1) (Int
    kInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1) Int
    rr Int
    rr
            u = STMatrix s a -> Int -> Int -> Int -> Int -> Slice s a
    forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    Slice STMatrix s a
    x Int
    k     (Int
    kInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1) Int
    o  Int
    rr
            l = STMatrix s a -> Int -> Int -> Int -> Int -> Slice s a
    forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    Slice STMatrix s a
    x (Int
    kInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1) Int
    k     Int
    rr Int
    o
        return (d,u,l,s)

    withVec :: Int
    -> (t -> t -> STVector s t -> ST s a) -> t -> t -> ST s (Vector t)
    withVec Int
    r t -> t -> STVector s t -> ST s a
    f = \t
    s t
    x -> do
        p <- Int -> ST s (STVector s t)
    forall t s. Storable t => Int -> ST s (STVector s t)
    newUndefinedVector Int
    r
        _ <- f s x p
        v <- unsafeFreezeVector p
        return v

    luPacked'' :: Matrix t -> (Matrix t, [Int])
    luPacked'' Matrix t
    m = (Matrix t -> Matrix t
    forall a. a -> a
    id (Matrix t -> Matrix t)
    -> (Vector Int -> [Int])
    -> (Matrix t, Vector Int)
    -> (Matrix t, [Int])
    forall b c b' c'. (b -> c) -> (b' -> c') -> (b, b') -> (c, c')
    forall (a :: * -> * -> *) b c b' c'.
    Arrow a =>
    a b c -> a b' c' -> a (b, b') (c, c')
    *** Vector Int -> [Int]
    forall a. Storable a => Vector a -> [a]
    toList) ((forall s. (Int, Int) -> STMatrix s t -> ST s (Vector Int))
    -> Matrix t -> (Matrix t, Vector Int)
    forall t u.
    Element t =>
    (forall s. (Int, Int) -> STMatrix s t -> ST s u)
    -> Matrix t -> (Matrix t, u)
    mutable (Int
    -> ((Int, Int) -> STMatrix s t -> STVector s Int -> ST s ())
    -> (Int, Int)
    -> STMatrix s t
    -> ST s (Vector Int)
    forall {t} {t} {t} {s} {a}.
    Storable t =>
    Int
    -> (t -> t -> STVector s t -> ST s a) -> t -> t -> ST s (Vector t)
    withVec (Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m) (Int, Int) -> STMatrix s t -> STVector s Int -> ST s ()
    forall {t} {b} {s}.
    (Container Vector t, Num (Vector t), Normed (Vector t),
     Fractional t) =>
    (Int, b) -> STMatrix s t -> STVector s Int -> ST s ()
    lu2) Matrix t
    m)
      where
        lu2 :: (Int, b) -> STMatrix s t -> STVector s Int -> ST s ()
    lu2 (Int
    r,b
    _) STMatrix s t
    x STVector s Int
    p = do
            [Int] -> (Int -> ST s ()) -> ST s ()
    forall (t :: * -> *) (m :: * -> *) a b.
    (Foldable t, Monad m) =>
    t a -> (a -> m b) -> m ()
    forM_ [Int
    0..Int
    rInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1] ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
    forall a b. (a -> b) -> a -> b
    $ \Int
    k -> do
                STMatrix s t -> STVector s Int -> Int -> ST s ()
    forall {e} {s}.
    (Container Vector e, Num e, Num (Vector e)) =>
    STMatrix s e -> STVector s Int -> Int -> ST s ()
    pivot STMatrix s t
    x STVector s Int
    p Int
    k
                (d,u,l,s) <- STMatrix s t
    -> Int -> Int -> ST s (t, Slice s t, Slice s t, Slice s t)
    forall {a} {s}.
    Storable a =>
    STMatrix s a
    -> Int -> Int -> ST s (a, Slice s a, Slice s a, Slice s a)
    view STMatrix s t
    x Int
    k Int
    r
                when (magnit 0 d) $ do
                    scalS (recip d) l
                    gemmm 1 s (-1) l u

        pivot :: STMatrix s e -> STVector s Int -> Int -> ST s ()
    pivot STMatrix s e
    x STVector s Int
    p Int
    k = do
            j <- Vector e -> Int
    Vector e -> IndexOf Vector
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    maxIndex (Vector e -> Int) -> (Matrix e -> Vector e) -> Matrix e -> Int
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector e -> Vector e
    forall a. Num a => a -> a
    abs (Vector e -> Vector e)
    -> (Matrix e -> Vector e) -> Matrix e -> Vector e
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Matrix e -> Vector e
    forall t. Element t => Matrix t -> Vector t
    flatten (Matrix e -> Int) -> ST s (Matrix e) -> ST s Int
    forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
    <$> STMatrix s e -> RowRange -> ColRange -> ST s (Matrix e)
    forall a t s.
    Element a =>
    STMatrix t a -> RowRange -> ColRange -> ST s (Matrix a)
    extractMatrix STMatrix s e
    x (Int -> RowRange
    FromRow Int
    k) (Int -> ColRange
    Col Int
    k)
            writeVector p k (j+k)
            swap k (k+j)
          where
            swap :: Int -> Int -> ST s ()
    swap Int
    i Int
    j = RowOper e -> STMatrix s e -> ST s ()
    forall t s.
    (Num t, Element t) =>
    RowOper t -> STMatrix s t -> ST s ()
    rowOper (Int -> Int -> ColRange -> RowOper e
    forall t. Int -> Int -> ColRange -> RowOper t
    SWAP Int
    i Int
    j ColRange
    AllCols) STMatrix s e
    x

    --------------------------------------------------------------------------------

    rowRange :: Matrix t -> [Int]
    rowRange Matrix t
    m = [Int
    0..Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m Int -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1]

    at :: Int -> Extractor
    at Int
    k = Vector I -> Extractor
    Pos ([Int] -> Vector I
    idxs[Int
    k])

    backSust' :: Matrix t -> Matrix t -> Matrix t
    backSust' Matrix t
    lup Matrix t
    rhs = (Matrix t -> (Matrix t, Matrix t, Matrix t) -> Matrix t)
    -> Matrix t -> [(Matrix t, Matrix t, Matrix t)] -> Matrix t
    forall b a. (b -> a -> b) -> b -> [a] -> b
    forall (t :: * -> *) b a.
    Foldable t =>
    (b -> a -> b) -> b -> t a -> b
    foldl' Matrix t -> (Matrix t, Matrix t, Matrix t) -> Matrix t
    forall {t} {a :: * -> *}.
    (Container Vector t, Fractional t, Num (Vector t),
     Mul a Matrix Matrix, Product t) =>
    Matrix t -> (Matrix t, a t, Matrix t) -> Matrix t
    f (Matrix t
    rhsMatrix t -> [Int] -> Matrix t
    forall t. Element t => Matrix t -> [Int] -> Matrix t
    ?[]) ([(Matrix t, Matrix t, Matrix t)]
    -> [(Matrix t, Matrix t, Matrix t)]
    forall a. [a] -> [a]
    reverse [(Matrix t, Matrix t, Matrix t)]
    ls)
      where
        ls :: [(Matrix t, Matrix t, Matrix t)]
    ls  = [ (Int -> Matrix t
    d Int
    k , Int -> Matrix t
    u Int
    k , Int -> Matrix t
    b Int
    k) | Int
    k <- Matrix t -> [Int]
    forall {t}. Matrix t -> [Int]
    rowRange Matrix t
    lup ]
          where
            d :: Int -> Matrix t
    d Int
    k = Matrix t
    lup Matrix t -> (Extractor, Extractor) -> Matrix t
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Int -> Extractor
    at Int
    k, Int -> Extractor
    at Int
    k)
            u :: Int -> Matrix t
    u Int
    k = Matrix t
    lup Matrix t -> (Extractor, Extractor) -> Matrix t
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Int -> Extractor
    at Int
    k, Int -> Extractor
    Drop (Int
    kInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1))
            b :: Int -> Matrix t
    b Int
    k = Matrix t
    rhs Matrix t -> (Extractor, Extractor) -> Matrix t
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Int -> Extractor
    at Int
    k, Extractor
    All)

        f :: Matrix t -> (Matrix t, a t, Matrix t) -> Matrix t
    f Matrix t
    x (Matrix t
    d,a t
    u,Matrix t
    b) = (Matrix t
    b Matrix t -> Matrix t -> Matrix t
    forall a. Num a => a -> a -> a
    - a t
    ua t -> Matrix t -> Matrix t
    forall t. Product t => a t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <>Matrix t
    x) Matrix t -> Matrix t -> Matrix t
    forall a. Fractional a => a -> a -> a
    / Matrix t
    d
                           Matrix t -> Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t -> Matrix t
    ===
                            Matrix t
    x

    forwSust' :: Matrix t -> Matrix t -> Matrix t
    forwSust' Matrix t
    lup Matrix t
    rhs = (Matrix t -> (Matrix t, Matrix t) -> Matrix t)
    -> Matrix t -> [(Matrix t, Matrix t)] -> Matrix t
    forall b a. (b -> a -> b) -> b -> [a] -> b
    forall (t :: * -> *) b a.
    Foldable t =>
    (b -> a -> b) -> b -> t a -> b
    foldl' Matrix t -> (Matrix t, Matrix t) -> Matrix t
    forall {t} {a :: * -> *}.
    (Container Vector t, Num (Vector t), Mul a Matrix Matrix,
     Product t) =>
    Matrix t -> (a t, Matrix t) -> Matrix t
    f (Matrix t
    rhsMatrix t -> [Int] -> Matrix t
    forall t. Element t => Matrix t -> [Int] -> Matrix t
    ?[]) [(Matrix t, Matrix t)]
    ls
      where
        ls :: [(Matrix t, Matrix t)]
    ls  = [ (Int -> Matrix t
    l Int
    k , Int -> Matrix t
    b Int
    k) | Int
    k <- Matrix t -> [Int]
    forall {t}. Matrix t -> [Int]
    rowRange Matrix t
    lup ]
          where
            l :: Int -> Matrix t
    l Int
    k = Matrix t
    lup Matrix t -> (Extractor, Extractor) -> Matrix t
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Int -> Extractor
    at Int
    k, Int -> Extractor
    Take Int
    k)
            b :: Int -> Matrix t
    b Int
    k = Matrix t
    rhs Matrix t -> (Extractor, Extractor) -> Matrix t
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Int -> Extractor
    at Int
    k, Extractor
    All)

        f :: Matrix t -> (a t, Matrix t) -> Matrix t
    f Matrix t
    x (a t
    l,Matrix t
    b) =     Matrix t
    x
                       Matrix t -> Matrix t -> Matrix t
    forall t. Element t => Matrix t -> Matrix t -> Matrix t
    ===
                    (Matrix t
    b Matrix t -> Matrix t -> Matrix t
    forall a. Num a => a -> a -> a
    - a t
    la t -> Matrix t -> Matrix t
    forall t. Product t => a t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <>Matrix t
    x)

    luSolve'' :: LU t -> Matrix t -> Matrix t
    luSolve'' (LU Matrix t
    lup [Int]
    p) Matrix t
    b = Matrix t -> Matrix t -> Matrix t
    forall {t}.
    (Container Vector t, Fractional t, Product t, Num (Vector t)) =>
    Matrix t -> Matrix t -> Matrix t
    backSust' Matrix t
    lup (Matrix t -> Matrix t -> Matrix t
    forall {t}.
    (Container Vector t, Product t, Num (Vector t)) =>
    Matrix t -> Matrix t -> Matrix t
    forwSust' Matrix t
    lup Matrix t
    pb)
      where
        pb :: Matrix t
    pb = Matrix t
    b Matrix t -> (Extractor, Extractor) -> Matrix t
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Vector I -> Extractor
    Pos ([Int] -> Vector I
    fixPerm' [Int]
    p), Extractor
    All)

    --------------------------------------------------------------------------------

    forwSust :: Matrix t -> Matrix t -> Matrix t
    forwSust Matrix t
    lup Matrix t
    rhs = (Matrix t, ()) -> Matrix t
    forall a b. (a, b) -> a
    fst ((Matrix t, ()) -> Matrix t) -> (Matrix t, ()) -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ (forall s. (Int, Int) -> STMatrix s t -> ST s ())
    -> Matrix t -> (Matrix t, ())
    forall t u.
    Element t =>
    (forall s. (Int, Int) -> STMatrix s t -> ST s u)
    -> Matrix t -> (Matrix t, u)
    mutable (Int, Int) -> STMatrix s t -> ST s ()
    forall s. (Int, Int) -> STMatrix s t -> ST s ()
    f Matrix t
    rhs
      where
        f :: (Int, Int) -> STMatrix s t -> ST s ()
    f (Int
    r,Int
    c) STMatrix s t
    x = do
            l <- Matrix t -> ST s (STMatrix s t)
    forall t s. Storable t => Matrix t -> ST s (STMatrix s t)
    unsafeThawMatrix Matrix t
    lup
            let go Int
    k = t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
    forall t s.
    Element t =>
    t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
    gemmm t
    1 (STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    Slice STMatrix s t
    x Int
    k Int
    0 Int
    1 Int
    c) (-t
    1) (STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    Slice STMatrix s t
    l Int
    k Int
    0 Int
    1 Int
    k) (STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    Slice STMatrix s t
    x Int
    0 Int
    0 Int
    k Int
    c)
            mapM_ go [0..r-1]

    backSust :: Matrix t -> Matrix t -> Matrix t
    backSust Matrix t
    lup Matrix t
    rhs = (Matrix t, ()) -> Matrix t
    forall a b. (a, b) -> a
    fst ((Matrix t, ()) -> Matrix t) -> (Matrix t, ()) -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ (forall s. (Int, Int) -> STMatrix s t -> ST s ())
    -> Matrix t -> (Matrix t, ())
    forall t u.
    Element t =>
    (forall s. (Int, Int) -> STMatrix s t -> ST s u)
    -> Matrix t -> (Matrix t, u)
    mutable (Int, Int) -> STMatrix s t -> ST s ()
    forall s. (Int, Int) -> STMatrix s t -> ST s ()
    f Matrix t
    rhs
      where
        f :: (Int, Int) -> STMatrix s t -> ST s ()
    f (Int
    r,Int
    c) STMatrix s t
    m = do
            l <- Matrix t -> ST s (STMatrix s t)
    forall t s. Storable t => Matrix t -> ST s (STMatrix s t)
    unsafeThawMatrix Matrix t
    lup
            let d Int
    k = t -> t
    forall a. Fractional a => a -> a
    recip (Matrix t
    lup Matrix t -> IndexOf Matrix -> t
    forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
    `atIndex` (Int
    k,Int
    k))
                u Int
    k = STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    Slice STMatrix s t
    l Int
    k (Int
    kInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1) Int
    1 (Int
    rInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1Int -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    k)
                b Int
    k = STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    Slice STMatrix s t
    m Int
    k Int
    0 Int
    1 Int
    c
                x Int
    k = STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    forall s t. STMatrix s t -> Int -> Int -> Int -> Int -> Slice s t
    Slice STMatrix s t
    m (Int
    kInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1) Int
    0 (Int
    rInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1Int -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    k) Int
    c
                scal Int
    k = RowOper t -> STMatrix s t -> ST s ()
    forall t s.
    (Num t, Element t) =>
    RowOper t -> STMatrix s t -> ST s ()
    rowOper (t -> RowRange -> ColRange -> RowOper t
    forall t. t -> RowRange -> ColRange -> RowOper t
    SCAL (Int -> t
    d Int
    k) (Int -> RowRange
    Row Int
    k) ColRange
    AllCols) STMatrix s t
    m

                go Int
    k = t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
    forall t s.
    Element t =>
    t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
    gemmm t
    1 (Int -> Slice s t
    b Int
    k) (-t
    1) (Int -> Slice s t
    u Int
    k) (Int -> Slice s t
    x Int
    k) ST s () -> ST s () -> ST s ()
    forall a b. ST s a -> ST s b -> ST s b
    forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
    >> Int -> ST s ()
    scal Int
    k
            mapM_ go [r-1,r-2..0]

    {- | Experimental implementation of 'luSolve' for any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.

    >>> let a = (2><2) [1,2,3,5] :: Matrix (Z ./. 13)
    (2><2)
     [ 1, 2
     , 3, 5 ]
    >>> b
    (2><3)
     [ 5, 1, 3
     , 8, 6, 3 ]

    >>> luSolve' (luPacked' a) b
    (2><3)
     [ 4,  7, 4
     , 7, 10, 6 ]

    -}
    luSolve' :: LU t -> Matrix t -> Matrix t
    luSolve' (LU Matrix t
    lup [Int]
    p) Matrix t
    b = Matrix t -> Matrix t -> Matrix t
    forall {t}.
    (Fractional t, Container Vector t) =>
    Matrix t -> Matrix t -> Matrix t
    backSust Matrix t
    lup (Matrix t -> Matrix t -> Matrix t
    forall {t}. (Element t, Num t) => Matrix t -> Matrix t -> Matrix t
    forwSust Matrix t
    lup Matrix t
    pb)
      where
        pb :: Matrix t
    pb = Matrix t
    b Matrix t -> (Extractor, Extractor) -> Matrix t
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Vector I -> Extractor
    Pos ([Int] -> Vector I
    fixPerm' [Int]
    p), Extractor
    All)

    --------------------------------------------------------------------------------

    data MatrixView t b
        = Elem t
        | Block b b b b
      deriving Int -> MatrixView t b -> String -> String
    [MatrixView t b] -> String -> String
    MatrixView t b -> String
    (Int -> MatrixView t b -> String -> String)
    -> (MatrixView t b -> String)
    -> ([MatrixView t b] -> String -> String)
    -> Show (MatrixView t b)
    forall a.
    (Int -> a -> String -> String)
    -> (a -> String) -> ([a] -> String -> String) -> Show a
    forall t b.
    (Show t, Show b) =>
    Int -> MatrixView t b -> String -> String
    forall t b.
    (Show t, Show b) =>
    [MatrixView t b] -> String -> String
    forall t b. (Show t, Show b) => MatrixView t b -> String
    $cshowsPrec :: forall t b.
    (Show t, Show b) =>
    Int -> MatrixView t b -> String -> String
    showsPrec :: Int -> MatrixView t b -> String -> String
    $cshow :: forall t b. (Show t, Show b) => MatrixView t b -> String
    show :: MatrixView t b -> String
    $cshowList :: forall t b.
    (Show t, Show b) =>
    [MatrixView t b] -> String -> String
    showList :: [MatrixView t b] -> String -> String
    Show

    viewBlock' :: Int -> Int -> Matrix t -> MatrixView t (Matrix t)
    viewBlock' Int
    r Int
    c Matrix t
    m
        | (Int
    rt,Int
    ct) (Int, Int) -> (Int, Int) -> Bool
    forall a. Eq a => a -> a -> Bool
    == (Int
    1,Int
    1) = t -> MatrixView t (Matrix t)
    forall t b. t -> MatrixView t b
    Elem (Matrix t -> Int -> Int -> t
    forall t. Storable t => Matrix t -> Int -> Int -> t
    atM' Matrix t
    m Int
    0 Int
    0)
        | Bool
    otherwise        = Matrix t
    -> Matrix t -> Matrix t -> Matrix t -> MatrixView t (Matrix t)
    forall t b. b -> b -> b -> b -> MatrixView t b
    Block Matrix t
    m11 Matrix t
    m12 Matrix t
    m21 Matrix t
    m22
      where
        (Int
    rt,Int
    ct) = Matrix t -> IndexOf Matrix
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    size Matrix t
    m
        m11 :: Matrix t
    m11 = (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
    subm (Int
    0,Int
    0) (Int
    r,Int
    c)       Matrix t
    m
        m12 :: Matrix t
    m12 = (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
    subm (Int
    0,Int
    c) (Int
    r,Int
    ctInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    c)    Matrix t
    m
        m21 :: Matrix t
    m21 = (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
    subm (Int
    r,Int
    0) (Int
    rtInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    r,Int
    c)    Matrix t
    m
        m22 :: Matrix t
    m22 = (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
    subm (Int
    r,Int
    c) (Int
    rtInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    r,Int
    ctInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    c) Matrix t
    m
        subm :: (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
    subm = (Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
    forall a.
    Element a =>
    (Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
    subMatrix

    viewBlock :: Matrix t -> MatrixView t (Matrix t)
    viewBlock Matrix t
    m = Int -> Int -> Matrix t -> MatrixView t (Matrix t)
    forall {t}.
    (Num t, Container Vector t) =>
    Int -> Int -> Matrix t -> MatrixView t (Matrix t)
    viewBlock' Int
    n Int
    n Matrix t
    m
      where
        n :: Int
    n = Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m Int -> Int -> Int
    forall a. Integral a => a -> a -> a
    `div` Int
    2

    invershur :: Matrix t -> Matrix t
    invershur (Matrix t -> MatrixView t (Matrix t)
    forall {t}.
    (Num t, Container Vector t) =>
    Matrix t -> MatrixView t (Matrix t)
    viewBlock -> Block Matrix t
    a Matrix t
    b Matrix t
    c Matrix t
    d) = [[Matrix t]] -> Matrix t
    forall t. Element t => [[Matrix t]] -> Matrix t
    fromBlocks [[Matrix t
    a',Matrix t
    b'],[Matrix t
    c',Matrix t
    d']]
      where
        r1 :: Matrix t
    r1 = Matrix t -> Matrix t
    invershur Matrix t
    a
        r2 :: Matrix t
    r2 = Matrix t
    c Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Matrix t
    r1
        r3 :: Matrix t
    r3 = Matrix t
    r1 Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Matrix t
    b
        r4 :: Matrix t
    r4 = Matrix t
    c Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Matrix t
    r3
        r5 :: Matrix t
    r5 = Matrix t
    r4Matrix t -> Matrix t -> Matrix t
    forall a. Num a => a -> a -> a
    -Matrix t
    d
        r6 :: Matrix t
    r6 = Matrix t -> Matrix t
    invershur Matrix t
    r5
        b' :: Matrix t
    b' = Matrix t
    r3 Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Matrix t
    r6
        c' :: Matrix t
    c' = Matrix t
    r6 Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Matrix t
    r2
        r7 :: Matrix t
    r7 = Matrix t
    r3 Matrix t -> Matrix t -> Matrix t
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Matrix t
    c'
        a' :: Matrix t
    a' = Matrix t
    r1Matrix t -> Matrix t -> Matrix t
    forall a. Num a => a -> a -> a
    -Matrix t
    r7
        d' :: Matrix t
    d' = -Matrix t
    r6

    invershur Matrix t
    x = Matrix t -> Matrix t
    forall a. Fractional a => a -> a
    recip Matrix t
    x

    --------------------------------------------------------------------------------

    instance Testable (Matrix I) where
       checkT :: Matrix I -> (Bool, IO ())
    checkT Matrix I
    _ = (Bool, IO ())
    test

    test :: (Bool, IO())
    test :: (Bool, IO ())
    test = ([Bool] -> Bool
    forall (t :: * -> *). Foldable t => t Bool -> Bool
    and [Bool]
    ok, () -> IO ()
    forall a. a -> IO a
    forall (m :: * -> *) a. Monad m => a -> m a
    return ())
      where
        m :: Matrix I
    m  = (Int
    3Int -> Int -> [I] -> Matrix I
    forall a. Storable a => Int -> Int -> [a] -> Matrix a
    ><Int
    4) [I
    1..I
    12] :: Matrix I
        r :: Matrix I
    r  = (Int
    2Int -> Int -> [I] -> Matrix I
    forall a. Storable a => Int -> Int -> [a] -> Matrix a
    ><Int
    3) [I
    1,I
    2,I
    3,I
    4,I
    3,I
    2]
        c :: Matrix I
    c  = (Int
    3Int -> Int -> [I] -> Matrix I
    forall a. Storable a => Int -> Int -> [a] -> Matrix a
    ><Int
    2) [I
    0,I
    4,I
    4,I
    1,I
    2,I
    3]
        p :: Matrix I
    p  = (Int
    9Int -> Int -> [I] -> Matrix I
    forall a. Storable a => Int -> Int -> [a] -> Matrix a
    ><Int
    10) [I
    0..I
    89] :: Matrix I
        ep :: Matrix I
    ep = (Int
    2Int -> Int -> [I] -> Matrix I
    forall a. Storable a => Int -> Int -> [a] -> Matrix a
    ><Int
    3) [I
    10,I
    24,I
    32,I
    44,I
    31,I
    23]
        md :: Matrix Double
    md = Matrix I -> Matrix Double
    forall (c :: * -> *) e. Container c e => c I -> c e
    fromInt Matrix I
    m      :: Matrix Double
        ok :: [Bool]
    ok = [ Matrix I -> Matrix I
    forall m mt. Transposable m mt => m -> mt
    tr Matrix I
    m Matrix I -> Matrix I -> Matrix I
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Matrix I
    m Matrix I -> Matrix I -> Bool
    forall a. Eq a => a -> a -> Bool
    == Matrix Double -> Matrix I
    forall (c :: * -> *) e. Container c e => c e -> c I
    toInt (Matrix Double -> Matrix Double
    forall m mt. Transposable m mt => m -> mt
    tr Matrix Double
    md Matrix Double -> Matrix Double -> Matrix Double
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Matrix Double
    md)
             , Matrix I
    m Matrix I -> Matrix I -> Matrix I
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Matrix I -> Matrix I
    forall m mt. Transposable m mt => m -> mt
    tr Matrix I
    m Matrix I -> Matrix I -> Bool
    forall a. Eq a => a -> a -> Bool
    == Matrix Double -> Matrix I
    forall (c :: * -> *) e. Container c e => c e -> c I
    toInt (Matrix Double
    md Matrix Double -> Matrix Double -> Matrix Double
    forall t. Product t => Matrix t -> Matrix t -> Matrix t
    forall (a :: * -> *) (b :: * -> *) (c :: * -> *) t.
    (Mul a b c, Product t) =>
    a t -> b t -> c t
    <> Matrix Double -> Matrix Double
    forall m mt. Transposable m mt => m -> mt
    tr Matrix Double
    md)
             , Matrix I
    m Matrix I -> (Extractor, Extractor) -> Matrix I
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Int -> Extractor
    Take Int
    2, Int -> Extractor
    Take Int
    3) Matrix I -> Matrix I -> Bool
    forall a. Eq a => a -> a -> Bool
    == Matrix I -> Matrix I -> Matrix I -> Matrix I
    forall t. Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t
    remap (Vector I -> Matrix I
    forall a. Storable a => Vector a -> Matrix a
    asColumn (Int -> Vector I
    range Int
    2)) (Vector I -> Matrix I
    forall a. Storable a => Vector a -> Matrix a
    asRow (Int -> Vector I
    range Int
    3)) Matrix I
    m
             , Matrix I -> Matrix I -> Matrix I -> Matrix I
    forall t. Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t
    remap Matrix I
    r (Matrix I -> Matrix I
    forall m mt. Transposable m mt => m -> mt
    tr Matrix I
    c) Matrix I
    p Matrix I -> Matrix I -> Bool
    forall a. Eq a => a -> a -> Bool
    == Matrix I
    ep
             , Matrix I -> Matrix I
    forall m mt. Transposable m mt => m -> mt
    tr Matrix I
    p Matrix I -> (Extractor, Extractor) -> Matrix I
    forall t.
    Element t =>
    Matrix t -> (Extractor, Extractor) -> Matrix t
    ?? (Vector I -> Extractor
    PosCyc ([Int] -> Vector I
    idxs[-Int
    5,Int
    13]), Vector I -> Extractor
    Pos ([Int] -> Vector I
    idxs[Int
    3,Int
    7,Int
    1])) Matrix I -> Matrix I -> Bool
    forall a. Eq a => a -> a -> Bool
    == (Int
    2Int -> Int -> [I] -> Matrix I
    forall a. Storable a => Int -> Int -> [a] -> Matrix a
    ><Int
    3) [I
    35,I
    75,I
    15,I
    33,I
    73,I
    13]
             ]