lapack 0.2.4 → 0.3

raw patch · 100 files changed

+12191/−5311 lines, 100 filesdep +blaze-htmldep +hyperdep +textdep ~basedep ~comfort-arraydep ~monoid-transformerPVP ok

version bump matches the API change (PVP)

Dependencies added: blaze-html, hyper, text

Dependency ranges changed: base, comfort-array, monoid-transformer, semigroups, transformers, utility-ht

API changes (from Hackage documentation)

Files

• Changes.md +15/−0
• Makefile +5/−1
• ReadMe.md +225/−0
• example/EconomicAllocation.hs +0/−71
• lapack.cabal +48/−26
• src/Numeric/LAPACK/Example/EconomicAllocation.hs +81/−0
• src/Numeric/LAPACK/Format.hs +28/−322
• src/Numeric/LAPACK/Linear/LowerUpper.hs +48/−300
• src/Numeric/LAPACK/Linear/Plain.hs +522/−0
• src/Numeric/LAPACK/Matrix.hs +765/−923
• src/Numeric/LAPACK/Matrix/Array.hs +237/−0
• src/Numeric/LAPACK/Matrix/Array/Banded.hs +24/−0
• src/Numeric/LAPACK/Matrix/Array/Basic.hs +22/−0
• src/Numeric/LAPACK/Matrix/Array/Triangular.hs +26/−0
• src/Numeric/LAPACK/Matrix/Banded.hs +162/−6
• src/Numeric/LAPACK/Matrix/Banded/Basic.hs +85/−5
• src/Numeric/LAPACK/Matrix/Banded/Linear.hs +2/−2
• src/Numeric/LAPACK/Matrix/BandedHermitian.hs +103/−6
• src/Numeric/LAPACK/Matrix/BandedHermitian/Basic.hs +65/−33
• src/Numeric/LAPACK/Matrix/BandedHermitian/Eigen.hs +3/−3
• src/Numeric/LAPACK/Matrix/BandedHermitianPositiveDefinite.hs +49/−2
• src/Numeric/LAPACK/Matrix/BandedHermitianPositiveDefinite/Linear.hs +2/−8
• src/Numeric/LAPACK/Matrix/Basic.hs +391/−52
• src/Numeric/LAPACK/Matrix/Class.hs +107/−0
• src/Numeric/LAPACK/Matrix/Divide.hs +101/−67
• src/Numeric/LAPACK/Matrix/Extent.hs +9/−2
• src/Numeric/LAPACK/Matrix/Extent/Private.hs +150/−0
• src/Numeric/LAPACK/Matrix/Hermitian.hs +290/−9
• src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs +324/−110
• src/Numeric/LAPACK/Matrix/Hermitian/Eigen.hs +6/−6
• src/Numeric/LAPACK/Matrix/Hermitian/Linear.hs +9/−12
• src/Numeric/LAPACK/Matrix/Hermitian/Private.hs +3/−5
• src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite.hs +45/−2
• src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs +6/−13
• src/Numeric/LAPACK/Matrix/Indexed.hs +21/−69
• src/Numeric/LAPACK/Matrix/Inverse.hs +56/−0
• src/Numeric/LAPACK/Matrix/Modifier.hs +49/−0
• src/Numeric/LAPACK/Matrix/Multiply.hs +147/−564
• src/Numeric/LAPACK/Matrix/Permutation.hs +71/−0
• src/Numeric/LAPACK/Matrix/Plain.hs +629/−0
• src/Numeric/LAPACK/Matrix/Plain/Class.hs +298/−0
• src/Numeric/LAPACK/Matrix/Plain/Divide.hs +162/−0
• src/Numeric/LAPACK/Matrix/Plain/Format.hs +348/−0
• src/Numeric/LAPACK/Matrix/Plain/Indexed.hs +86/−0
• src/Numeric/LAPACK/Matrix/Plain/Multiply.hs +608/−0
• src/Numeric/LAPACK/Matrix/Private.hs +12/−20
• src/Numeric/LAPACK/Matrix/RowMajor.hs +212/−0
• src/Numeric/LAPACK/Matrix/Shape.hs +1/−0
• src/Numeric/LAPACK/Matrix/Shape/Private.hs +11/−1
• src/Numeric/LAPACK/Matrix/Special.hs +11/−0
• src/Numeric/LAPACK/Matrix/Square.hs +248/−20
• src/Numeric/LAPACK/Matrix/Square/Basic.hs +50/−20
• src/Numeric/LAPACK/Matrix/Square/Eigen.hs +17/−30
• src/Numeric/LAPACK/Matrix/Square/Linear.hs +3/−4
• src/Numeric/LAPACK/Matrix/Symmetric/Private.hs +4/−6
• src/Numeric/LAPACK/Matrix/Triangular.hs +425/−14
• src/Numeric/LAPACK/Matrix/Triangular/Basic.hs +163/−126
• src/Numeric/LAPACK/Matrix/Triangular/Linear.hs +11/−24
• src/Numeric/LAPACK/Matrix/Triangular/Private.hs +96/−62
• src/Numeric/LAPACK/Matrix/Type.hs +117/−0
• src/Numeric/LAPACK/Orthogonal.hs +34/−225
• src/Numeric/LAPACK/Orthogonal/Householder.hs +129/−13
• src/Numeric/LAPACK/Orthogonal/Plain.hs +278/−0
• src/Numeric/LAPACK/Orthogonal/Private.hs +120/−65
• src/Numeric/LAPACK/Output.hs +107/−0
• src/Numeric/LAPACK/Permutation.hs +35/−10
• src/Numeric/LAPACK/Permutation/Private.hs +201/−103
• src/Numeric/LAPACK/Private.hs +132/−24
• src/Numeric/LAPACK/Scalar.hs +22/−0
• src/Numeric/LAPACK/Shape.hs +40/−0
• src/Numeric/LAPACK/Shape/Private.hs +15/−0
• src/Numeric/LAPACK/ShapeStatic.hs +3/−0
• src/Numeric/LAPACK/Singular.hs +94/−400
• src/Numeric/LAPACK/Singular/Plain.hs +475/−0
• src/Numeric/LAPACK/Split.hs +5/−9
• src/Numeric/LAPACK/Vector.hs +220/−162
• src/Numeric/LAPACK/Vector/Private.hs +55/−0
• test/Main.hs +4/−0
• test/Test/Banded.hs +64/−58
• test/Test/Banded/Utility.hs +16/−19
• test/Test/BandedHermitian.hs +91/−76
• test/Test/Divide.hs +130/−0
• test/Test/Example.hs +21/−0
• test/Test/Format.hs +52/−37
• test/Test/Generator.hs +370/−245
• test/Test/Generic.hs +52/−0
• test/Test/Hermitian.hs +229/−127
• test/Test/Indexed.hs +24/−28
• test/Test/Logic.hs +23/−6
• test/Test/LowerUpper.hs +221/−0
• test/Test/Matrix.hs +252/−202
• test/Test/Multiply.hs +50/−0
• test/Test/Orthogonal.hs +120/−103
• test/Test/Permutation.hs +47/−55
• test/Test/Shape.hs +37/−21
• test/Test/Singular.hs +42/−30
• test/Test/Square.hs +75/−47
• test/Test/Triangular.hs +230/−194
• test/Test/Utility.hs +109/−43
• test/Test/Vector.hs +153/−63

@@ -0,0 +1,15 @@+# Change log for the `lapack` package++## 0.3++ * Matrix data family++ * `Matrix`: `ZeroInt` -> `ShapeInt`, `zeroInt` -> `shapeInt`++ * `Hermitian`, `BandedHermitian`: `covariance` -> `gramian`++ * `Square.eigensystem`:+   Return left eigenvectors as rows of the last matrix.+   This is adjoint with respect to the definition in `lapack-0.2`+   but it is consistent+   with the other eigenvalue and singular value decompositions.

@@ -1,6 +1,10 @@+all:	ReadMe.html Changes.html+ run-test: 	runhaskell Setup configure --user --enable-tests -fbuildExamples 	runhaskell Setup build 	runhaskell Setup haddock 	./dist/build/lapack-test/lapack-test-	./dist/build/lapack-economic/lapack-economic++%.html: %.md+	pandoc --output=$@ $<

@@ -0,0 +1,225 @@+## Installation++Before installing the Haskell bindings+you need to install the BLAS and LAPACK packages.+Please note, that additionally to the reference implementation in FORTRAN 77+there are alternative optimized implementations+like OpenBLAS, ATLAS, Intel MKL.++### Debian, Ubuntu++~~~~+sudo apt-get install libblas-dev liblapack-dev+~~~~++### MacOS++You may install pkgconfig and LAPACK via MacPorts:++~~~~+sudo port install pkgconfig lapack+~~~~++However, the pkg-config files for LAPACK+seem to be installed in a non-standard location.+You must make them visible to pkg-config by++~~~~+export PKG_CONFIG_PATH=/opt/local/lib/lapack/pkgconfig:$PKG_CONFIG_PATH+~~~~++You may set the search PATH permanently by adding+the above command line to your `$HOME/.profile` file.++Alternatively, a solution for all users of your machine+would be to set symbolic links:++~~~~+sudo ln -s /opt/local/lib/lapack/pkgconfig/blas.pc /opt/local/lib/pkgconfig/blas.pc+sudo ln -s /opt/local/lib/lapack/pkgconfig/lapack.pc /opt/local/lib/pkgconfig/lapack.pc+~~~~+++## Introduction++Here is a small example for constructing and formatting matrices:++~~~~+Prelude> import qualified Numeric.LAPACK.Matrix as Matrix+Prelude Matrix> import Numeric.LAPACK.Format as Fmt ((##))+Prelude Matrix Fmt> let a = Matrix.fromList (Matrix.shapeInt 3) (Matrix.shapeInt 4) [(0::Float)..]+Prelude Matrix Fmt> a ## "%.4f"+ 0.0000 1.0000  2.0000  3.0000+ 4.0000 5.0000  6.0000  7.0000+ 8.0000 9.0000 10.0000 11.0000+Prelude Matrix Fmt> import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape+Prelude Matrix Fmt MatrixShape> import qualified Numeric.LAPACK.Matrix.Triangular as Triangular+Prelude Matrix Fmt MatrixShape Triangular> let u = Triangular.upperFromList MatrixShape.RowMajor (Matrix.shapeInt 4) [(0::Float)..]+Prelude Matrix Fmt MatrixShape Triangular> (u, Triangular.transpose u) ## "%.4f"+ 0.0000 1.0000 2.0000 3.0000+        4.0000 5.0000 6.0000+               7.0000 8.0000+                      9.0000++ 0.0000+ 1.0000 4.0000+ 2.0000 5.0000 7.0000+ 3.0000 6.0000 8.0000 9.0000+~~~~++You may find a more complex introductory example at:+<http://code.henning-thielemann.de/bob2019/main.pdf>+++## Formatting++We do not try to do fancy formatting for the `Show` instance.+The `Show` instances of matrices emit plain valid Haskell code in one line+where all numbers are printed in full precision.+If matrices are part of larger Haskell data structures+this kind of formatting works best.+For human-friendly formatting in GHCi you need to add something like `## "%.4f"`+after a matrix or vector expression.+It means: Print all numbers in fixed point representation+using four digits for the fractional part.+You can use the formatting placeholders provided by `printf`.+The matrices have `Hyper` instances for easy usage in+[HyperHaskell](https://hackage.haskell.org/package/hyper-haskell-server).++Formatting is based on the `Output` type class+that currently supports output as+[Text boxes](https://hackage.haskell.org/package/boxes)+for GHCi and+[HTML](https://hackage.haskell.org/package/blaze-html)+for HyperHaskell.+++## Matrix vs. Vector++Vectors are `Storable.Array`s from the+[comfort-array](https://hackage.haskell.org/package/comfort-array) package.+An array can have a fancy shape+like a shape defined by an Haskell enumeration type,+the shape of two appended arrays,+a shape that is compatible to a Haskell container type,+a rectangular or triangular shape.++All operations check dynamically+whether corresponding shapes match structurally.+E.g. `a|||b === c|||d` composes a matrix from four quadrants `a`, `b`, `c`, `d`.+It is not enough that `a|||b` and `c|||d` have the same width,+but the widths of `a` and `c` as well as `b` and `d` must match.+The type variables for shapes show which dimensions must be compatible.+We recommend to use type variables for the shapes as long as possible+because this increases type safety even+if you eventually use only `ShapeInt`.+If you use statically sized shapes you get static size checks.++A matrix can have any of these shapes as height and as width.+E.g. it is no problem to define a matrix+that maps a triangular shaped array to a rectangular shaped one.+There are actually practical applications to such matrices.+A matrix can be treated as vector, but there are limitations.+E.g. if you scale a Hermitian matrix by a complex factor+it will in general be no longer Hermitian.+Another problem: Two equally sized rectangular matrices+may differ in the element order (row major vs. column major).+You cannot simply add them by adding the flattened arrays element-wise.+Thus if you want to perform vector operations on a matrix+the package requires you to "unpack" a matrix to a vector+using `Matrix.Array.toVector`.+This conversion is almost a no-op and preserves most of the shape information.+The reverse operation is `Matrix.Array.fromVector`.++There are more matrix types that are not based on a single array.+E.g. we provide a symbolic inverse, a scaling matrix, a permutation matrix.+We also support arrays that represent factors of a matrix factorization.+You obtain these by LU and QR decompositions.+You can extract the matrix factors of it,+but you can also multiply the factors to other matrices+or use the decompositions for solving simultaneous linear equations.+++## Type tags for content constraints++Full matrices have additional type tags to distinguish+between four cases of the size relations between+the height and the width of a full matrix.+In a matrix of type `Matrix.Full vert horiz height width a`+the type variables mean:++~~~~+vert  height+Small Small - Square matrix   height==width+Big   Small - Tall matrix     size height >= size width+Small Big   - Wide matrix     size height <= size width+Big   Big   - General matrix  height and width arbitrary+~~~~++The relations are defined using two type tags+in order to support matrix transposition flawlessly.+Using `Small Small` for square matrices+and `Big Big` for general matrices appears to be arbitrary,+but is chosen such that altering `Small` to `Big`+generalizes the size relation.++Likewise we use the `Triangular` matrix type+also to represent diagonal and symmetric matrices.+For `Matrix.Triangular lo diag up size a` we have the cases:++~~~~+lo    up+Empty Empty - Diagonal matrix+Empty Full  - Upper triangular matrix+Full  Empty - Lower triangular matrix+Full  Full  - Symmetric matrix+~~~~++The `diag` type tag can be `NonUnit` or `Unit`.+`Unit` can be used for matrices+that always have a unit diagonal by construction.+The property of a unit diagonal is preserved by some operations+and enables some optimizations by LAPACK.+`NonUnit` is a bit of a misnomer.+A `NonUnit` matrix can still have a unit diagonal,+but in general it has not and no optimizations will take place.+++## Infix operators++The package provides fancy infix operators like `#*|` and `\*#`.+They symbolize both operands and operations.+E.g. in `#*|` the hash means Matrix, the star means Multiplication+and the bar means Column Vector.++Possible operations are:++  * `a_*_b` - `a` multiplied by `b`++  * `a_/_b` - `a` multiplied by `inverse b`++  * `a_\_b` - `inverse a` multiplied by `b`++Possible operands are:++  * `#` - a matrix that is generalized through a type class++  * `##` - a full matrix++  * `\` - a diagonal matrix represented by a `Vector`++  * `-` - a row vector++  * `|` - a column vector++  * `.` - a scalar++For multiplication of equally shaped matrices+we also provide instances of `Semigroup.<>`.++Precedence of the operators is chosen analogously to plain `*` and `/`.+Associativity is chosen such that the same operator can be applied+multiple times without parentheses.+But sometimes this may mean that you have to mix+left and right associative operators,+and thus you may still need parentheses.

@@ -1,71 +0,0 @@-{-# LANGUAGE TypeOperators #-}-module Main where--import qualified Numeric.LAPACK.Matrix.Square as Square-import qualified Numeric.LAPACK.Matrix as Matrix-import qualified Numeric.LAPACK.Vector as Vector-import Numeric.LAPACK.Matrix (ZeroInt, (#>), (|||))-import Numeric.LAPACK.Format ((##))--import qualified Data.Array.Comfort.Storable as Array-import qualified Data.Array.Comfort.Shape as Shape-import Data.Array.Comfort.Shape ((:+:)((:+:)))-import Data.Function.HT (nest)---type ZeroInt2 = ZeroInt:+:ZeroInt-type Vector sh = Vector.Vector sh Double-type Matrix height width = Matrix.General height width Double-type SquareMatrix size = Square.Square size Double---balances :: Vector ZeroInt2-balances =-   Vector.fromList (Matrix.zeroInt 2 :+: Matrix.zeroInt 2)-      [100000, 90000, -50000, -120000]--expenses :: Matrix ZeroInt ZeroInt2-expenses =-   Matrix.fromList (Matrix.zeroInt 2) (Matrix.zeroInt 2 :+: Matrix.zeroInt 2) $-   [16000,  4000,  8000, 12000,-    10000, 30000, 40000, 20000]--normalize ::-   (Eq height, Shape.C height, Shape.C width) =>-   Matrix height width -> Matrix height width-normalize x = Matrix.scaleRows (Array.map recip (Matrix.rowSums x)) x---subtractIdentity :: (Eq sh, Shape.C sh) => SquareMatrix sh -> SquareMatrix sh-subtractIdentity x = Matrix.sub x $ Square.identityFrom x--completeIdSquare :: Matrix ZeroInt2 ZeroInt -> SquareMatrix ZeroInt2-completeIdSquare x =-   Square.fromGeneral $-      (Matrix.takeLeftColumns $ Matrix.fromFull $ Square.identityFromHeight x)-      |||-      x--iterationMatrix :: SquareMatrix ZeroInt2-iterationMatrix =-   completeIdSquare $ Matrix.transpose $ normalize expenses--iterated :: Vector ZeroInt2-iterated = nest 30 (iterationMatrix #>) balances----compensated :: Vector ZeroInt-compensated =-   let a = Matrix.transpose $ normalize expenses-       p = Matrix.takeTopRows a-       k = Square.fromGeneral $ Matrix.takeBottomRows a-       x = Vector.takeLeft balances-       y = Vector.takeRight balances-   in Vector.sub x $ p #> Matrix.solveVector (subtractIdentity k) y---main :: IO ()-main = do-   Array.mapShape (Shape.ZeroBased . Shape.size) iterated ## "%10.2f"-   compensated ## "%10.2f"

@@ -1,5 +1,5 @@ Name:             lapack-Version:          0.2.4+Version:          0.3 License:          BSD3 License-File:     LICENSE Author:           Henning Thielemann <haskell@henning-thielemann.de>@@ -28,24 +28,27 @@   .   * Works with matrices and vectors with zero dimensions.   .+  * No automatic (and dangerous) implicit expansion+    of singleton vectors or matrices.+    Instead there are special operators for scaling of vectors and matrices.+  .   * Separate formatting operator @(##)@:     Works better for tuples of matrices and vectors than 'show'.     'Show' is used for code one-liners     that can be copied back into Haskell modules.+    Support for nice formatting in HyperHaskell.   .   See also: @hmatrix@. Tested-With:      GHC==7.4.2, GHC==7.8.4, GHC==8.2.2 Cabal-Version:    >=1.8 Build-Type:       Simple Extra-Source-Files:+  ReadMe.md+  Changes.md   Makefile -Flag buildExamples-  description: Build example executables-  default:     False- Source-Repository this-  Tag:         0.2.4+  Tag:         0.3   Type:        darcs   Location:    http://hub.darcs.net/thielema/lapack/ @@ -58,15 +61,18 @@     lapack-ffi >=0.0.1 && <0.1,     blas-ffi >=0.0 && <0.2,     netlib-ffi >=0.1.1 && <0.2,-    comfort-array >=0.3 && <0.4,+    comfort-array >=0.4 && <0.5,     guarded-allocation >=0.0 && <0.1,+    hyper >=0.1 && <0.2,+    blaze-html >=0.7 && <0.10,+    text >=1.2 && <1.3,     boxes >=0.1.5 && <0.2,     deepseq >=1.3 && <1.5,     lazyio >=0.1 && <0.2,-    transformers >=0.3 && <0.6,+    transformers >=0.4 && <0.6,     tfp >=1.0.1 && <1.1,     fixed-length >=0.2 && <0.3,-    semigroups >=0.18.3 && <0.19,+    semigroups >=0.18.3 && <1.0,     non-empty >=0.3 && <0.4,     utility-ht >=0.0.10 && <0.1,     base >=4.5 && <5@@ -75,6 +81,8 @@   Hs-Source-Dirs:   src   Exposed-Modules:     Numeric.LAPACK.Matrix+    Numeric.LAPACK.Matrix.Special+    Numeric.LAPACK.Matrix.Array     Numeric.LAPACK.Matrix.Extent     Numeric.LAPACK.Matrix.Shape.Box     Numeric.LAPACK.Matrix.Shape@@ -85,6 +93,7 @@     Numeric.LAPACK.Matrix.Banded     Numeric.LAPACK.Matrix.BandedHermitian     Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite+    Numeric.LAPACK.Matrix.Permutation     Numeric.LAPACK.Vector     Numeric.LAPACK.Scalar     Numeric.LAPACK.Orthogonal@@ -93,9 +102,15 @@     Numeric.LAPACK.Linear.LowerUpper     Numeric.LAPACK.Singular     Numeric.LAPACK.ShapeStatic+    Numeric.LAPACK.Shape     Numeric.LAPACK.Format+    Numeric.LAPACK.Output+    Numeric.LAPACK.Example.EconomicAllocation   Other-Modules:+    Numeric.LAPACK.Singular.Plain+    Numeric.LAPACK.Orthogonal.Plain     Numeric.LAPACK.Orthogonal.Private+    Numeric.LAPACK.Linear.Plain     Numeric.LAPACK.Linear.Private     Numeric.LAPACK.Split     Numeric.LAPACK.Permutation.Private@@ -120,12 +135,28 @@     Numeric.LAPACK.Matrix.Shape.Private     Numeric.LAPACK.Matrix.Extent.Private     Numeric.LAPACK.Matrix.Extent.Kind-    Numeric.LAPACK.Matrix.Multiply     Numeric.LAPACK.Matrix.Divide-    Numeric.LAPACK.Matrix.Indexed+    Numeric.LAPACK.Matrix.Multiply+    Numeric.LAPACK.Matrix.Class     Numeric.LAPACK.Matrix.Basic+    Numeric.LAPACK.Matrix.Plain.Class+    Numeric.LAPACK.Matrix.Plain.Divide+    Numeric.LAPACK.Matrix.Plain.Multiply+    Numeric.LAPACK.Matrix.Plain.Indexed+    Numeric.LAPACK.Matrix.Plain.Format+    Numeric.LAPACK.Matrix.Plain+    Numeric.LAPACK.Matrix.Indexed+    Numeric.LAPACK.Matrix.Type+    Numeric.LAPACK.Matrix.Inverse+    Numeric.LAPACK.Matrix.Array.Banded+    Numeric.LAPACK.Matrix.Array.Triangular+    Numeric.LAPACK.Matrix.Array.Basic+    Numeric.LAPACK.Matrix.Modifier+    Numeric.LAPACK.Matrix.RowMajor     Numeric.LAPACK.Matrix.Private+    Numeric.LAPACK.Vector.Private     Numeric.LAPACK.Private+    Numeric.LAPACK.Shape.Private     Numeric.LAPACK.Wrapper  Test-Suite lapack-test@@ -142,7 +173,7 @@     QuickCheck >=2.5 && <3,     ChasingBottoms >=1.2.2 && <1.4,     transformers,-    monoid-transformer >=0.0.3 && <0.1,+    monoid-transformer >=0.0.4 && <0.1,     semigroups,     non-empty >=0.3.1,     utility-ht,@@ -153,7 +184,11 @@   Main-Is:          Main.hs   Other-Modules:     Test.Shape+    Test.Example     Test.Indexed+    Test.Divide+    Test.Multiply+    Test.Generic     Test.Permutation     Test.Vector     Test.Matrix@@ -161,6 +196,7 @@     Test.Triangular     Test.Hermitian     Test.Orthogonal+    Test.LowerUpper     Test.Banded     Test.BandedHermitian     Test.Banded.Utility@@ -169,17 +205,3 @@     Test.Logic     Test.Format     Test.Utility--Executable lapack-economic-  If flag(buildExamples)-    Build-Depends:-      lapack,-      comfort-array,-      utility-ht,-      base-  Else-    Buildable: False--  GHC-Options:      -Wall-  Hs-Source-Dirs:   example-  Main-Is:          EconomicAllocation.hs

@@ -0,0 +1,81 @@+{-# LANGUAGE TypeOperators #-}+{- |+Do not import this module. It is only for demonstration purposes.+-}+module Numeric.LAPACK.Example.EconomicAllocation where++import qualified Numeric.LAPACK.Matrix.Square as Square+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix.Square ((|=|))+import Numeric.LAPACK.Matrix (ShapeInt, shapeInt, (#-#), (#*|), (#\|), (\\#))+import Numeric.LAPACK.Vector ((|-|))+import Numeric.LAPACK.Format ((##))++import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Shape ((:+:)((:+:)))+++type ZeroInt2 = ShapeInt:+:ShapeInt+type Vector sh = Vector.Vector sh Double+type Matrix height width = Matrix.General height width Double+type SquareMatrix size = Square.Square size Double+++balances0 :: Vector ZeroInt2+balances0 =+   Vector.fromList (shapeInt 2 :+: shapeInt 2)+      [100000, 90000, -50000, -120000]++expenses0 :: Matrix ShapeInt ZeroInt2+expenses0 =+   Matrix.fromList (shapeInt 2) (shapeInt 2 :+: shapeInt 2) $+   [16000,  4000,  8000, 12000,+    10000, 30000, 40000, 20000]++normalize ::+   (Eq height, Shape.C height, Shape.C width) =>+   Matrix height width -> Matrix height width+normalize x = Matrix.rowSums x \\# x++normalizeSplit ::+   (Shape.C sh0, Shape.C sh1, Eq sh1) =>+   Matrix sh1 (sh0:+:sh1) -> (Matrix sh0 sh1, SquareMatrix sh1)+normalizeSplit expenses =+   let a = Matrix.transpose $ normalize expenses+   in (Matrix.takeTop a, Square.fromGeneral $ Matrix.takeBottom a)+++completeIdSquare ::+   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1) =>+   Matrix sh1 (sh0:+:sh1) -> SquareMatrix (sh0:+:sh1)+completeIdSquare x =+   let (p,k) = normalizeSplit x+   in (Square.identityFromHeight p, p)+      |=|+      (Matrix.zero $ ArrMatrix.shape $ Matrix.transpose p, k)++iterated ::+   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1) =>+   Matrix sh1 (sh0:+:sh1) -> Vector (sh0:+:sh1) -> Vector (sh0:+:sh1)+iterated expenses =+   -- 'Stream.head' would be total+   head . dropWhile ((>=1e-5) . Vector.normInf . Vector.takeRight) .+   iterate (completeIdSquare expenses #*|)+++compensated ::+   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1) =>+   Matrix sh1 (sh0:+:sh1) -> Vector (sh0:+:sh1) -> Vector sh0+compensated expenses balances =+   let (p,k) = normalizeSplit expenses+       x = Vector.takeLeft balances+       y = Vector.takeRight balances+   in x |-| p #*| (k #-# Square.identityFrom k) #\| y+++main :: IO ()+main = do+   iterated expenses0 balances0 ## "%10.2f"+   compensated expenses0 balances0 ## "%10.2f"

@@ -2,43 +2,33 @@ {-# LANGUAGE ConstraintKinds #-} module Numeric.LAPACK.Format (    (##),+   hyper,    Format(format),    FormatArray(formatArray),-   deflt,+   Type.FormatMatrix(formatMatrix),+   ArrFormat.deflt,    ) where -import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape-import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent-import Numeric.LAPACK.Matrix.Shape.Private-         (Order(RowMajor, ColumnMajor), Filled(Filled), UnaryProxy)-import Numeric.LAPACK.Matrix.Private (Full)-import Numeric.LAPACK.Scalar (conjugate)-import Numeric.LAPACK.Wrapper (Flip(Flip, getFlip))+import qualified Numeric.LAPACK.Matrix.Plain.Format as ArrFormat+import qualified Numeric.LAPACK.Matrix.Type as Type+import qualified Numeric.LAPACK.Output as Output+import qualified Numeric.LAPACK.Permutation.Private as Perm+import Numeric.LAPACK.Matrix.Plain.Format+         (FormatArray, formatArray, printfComplex)+import Numeric.LAPACK.Output (Output, (/+/))  import qualified Numeric.Netlib.Class as Class -import qualified Type.Data.Num.Unary.Literal as TypeNum-import qualified Type.Data.Num.Unary as Unary-import Type.Data.Num (integralFromProxy)--import qualified Data.Array.Comfort.Storable.Unchecked as Array import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable (Array)  import qualified Text.PrettyPrint.Boxes as TextBox-import Text.PrettyPrint.Boxes (Box, (/+/))-import Text.Printf (PrintfArg, printf)+import qualified Hyper+import Text.Printf (printf)  import qualified Data.List.Reverse.StrictSpine as ListRev-import qualified Data.List.Match as Match-import qualified Data.List.HT as ListHT-import qualified Data.List as List-import Data.Functor.Compose (Compose(Compose, getCompose))-import Data.Foldable (foldMap)-import Data.List (mapAccumL, transpose)-import Data.Complex (Complex((:+)))-import Data.Maybe.HT (toMaybe)-import Data.Maybe (fromMaybe)+import Data.Foldable (fold)+import Data.Complex (Complex) import Data.Char (isSpace)  @@ -50,28 +40,27 @@ trim :: String -> String trim = unlines . map (ListRev.dropWhile isSpace) . lines --deflt :: String-deflt = "%.4g"+hyper :: (Format a) => String -> a -> Hyper.Graphic+hyper fmt = Output.hyper . format fmt   class Format a where-   format :: String -> a -> Box+   format :: (Output out) => String -> a -> out  instance Format Int where-   format _fmt = TextBox.text . show+   format _fmt = Output.text . show  instance Format Float where-   format fmt = TextBox.text . printf fmt+   format fmt = Output.text . printf fmt  instance Format Double where-   format fmt = TextBox.text . printf fmt+   format fmt = Output.text . printf fmt  instance (Class.Real a) => Format (Complex a) where-   format fmt = TextBox.text . concat . printfComplex fmt+   format fmt = Output.text . fold . printfComplex fmt  instance (Format a) => Format [a] where-   format fmt = TextBox.vsep 1 TextBox.right . map (format fmt)+   format fmt = Output.formatColumn . map (format fmt)  instance (Format a, Format b) => Format (a,b) where    format fmt (a,b) = format fmt a /+/ format fmt b@@ -79,296 +68,13 @@ instance (Format a, Format b, Format c) => Format (a,b,c) where    format fmt (a,b,c) = format fmt a /+/ format fmt b /+/ format fmt c +instance (Shape.C sh) => Format (Perm.Permutation sh) where+   format _fmt = Perm.format+ instance (FormatArray sh, Class.Floating a) => Format (Array sh a) where    format = formatArray --class (Shape.C sh) => FormatArray sh where-   {--   We use constraint @(Class.Floating a)@ and not @(Format a)@-   because it allows us to align the components of complex numbers.-   -}-   formatArray :: (Class.Floating a) => String -> Array sh a -> Box--instance (Integral i) => FormatArray (Shape.ZeroBased i) where-   formatArray = formatVector--instance (Integral i) => FormatArray (Shape.OneBased i) where-   formatArray = formatVector--formatVector :: (Shape.C sh, Class.Floating a) => String -> Array sh a -> Box-formatVector fmt =-   TextBox.hsep 1 TextBox.right .-   map (TextBox.text . concat . printfFloating fmt) . Array.toList- instance-   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>-      FormatArray (MatrixShape.Full vert horiz height width) where-   formatArray = formatFull--formatFull ::-   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,-    Class.Floating a) =>-   String -> Full vert horiz height width a -> Box-formatFull fmt m =-   let MatrixShape.Full order extent = Array.shape m-   in  formatAligned (printfFloating fmt) $-       splitRows order (Extent.dimensions extent) $ Array.toList m--instance-   (Eq lower, Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>-      FormatArray (MatrixShape.Split lower vert horiz height width) where-   formatArray = formatHouseholder--formatHouseholder ::-   (Eq lower, Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,-    Class.Floating a) =>-   String -> Array (MatrixShape.Split lower vert horiz height width) a -> Box-formatHouseholder fmt m =-   let MatrixShape.Split _ order extent = Array.shape m-   in formatSeparateTriangle (printfFloating fmt) $-      splitRows order (Extent.dimensions extent) $ Array.toList m--instance (Shape.C size) => FormatArray (MatrixShape.Hermitian size) where-   formatArray = formatHermitian--formatHermitian ::-   (Shape.C size, Class.Floating a) =>-   String -> Array (MatrixShape.Hermitian size) a -> Box-formatHermitian fmt m =-   let MatrixShape.Hermitian order size = Array.shape m-   in  formatSeparateTriangle (printfFloating fmt) $-       complementTriangle conjugate order (Shape.size size) $ Array.toList m--formatSymmetric ::-   (Shape.C size, Class.Floating a) =>-   String -> Array (MatrixShape.Symmetric size) a -> Box-formatSymmetric fmt m =-   let MatrixShape.Triangular _diag (Filled, Filled) order size = Array.shape m-   in  formatSeparateTriangle (printfFloating fmt) $-       complementTriangle id order (Shape.size size) $ Array.toList m--complementTriangle ::-   (Class.Floating a) => (a -> a) -> Order -> Int -> [a] -> [[a]]-complementTriangle adapt order n xs =-   let mergeTriangles lower upper =-         zipWith (++) (map (map adapt . init) lower) upper-   in case order of-         RowMajor ->-            let tri = slice (take n $ iterate pred n) xs-                trans = reverse $ transpose $ map reverse tri-            in  mergeTriangles trans tri-         ColumnMajor ->-            let tri = slice (take n [1..]) xs-            in  mergeTriangles tri (transpose tri)--instance-   (MatrixShape.Content lo, MatrixShape.Content up,-    MatrixShape.TriDiag diag, Shape.C size) =>-      FormatArray (MatrixShape.Triangular lo diag up size) where-   formatArray fmt =-      getFormatTriangular $-      MatrixShape.switchDiagUpLoSym-         (FormatTriangular $ \m ->-            let MatrixShape.Triangular _diag _uplo order size = Array.shape m-                n0 = Unary.unary TypeNum.u0-            in formatAligned (printfFloatingMaybe fmt) $-               formatBanded (n0,n0) order (size,size) $ Array.toList m)-         (FormatTriangular $ formatTriangular fmt)-         (FormatTriangular $ formatTriangular fmt)-         (FormatTriangular $-            formatSymmetric fmt .-            Array.mapShape MatrixShape.strictNonUnitDiagonal)--newtype FormatTriangular diag size a b lo up =-   FormatTriangular {-      getFormatTriangular ::-         Array (MatrixShape.Triangular lo diag up size) a -> b-   }--formatTriangular ::-   (MatrixShape.TriDiag diag, MatrixShape.UpLo lo up,-    Shape.C size, Class.Floating a) =>-   String -> Array (MatrixShape.Triangular lo diag up size) a -> Box-formatTriangular fmt m =-   let MatrixShape.Triangular _diag uplo order size = Array.shape m-   in  formatAligned (printfFloatingMaybe fmt) $-       MatrixShape.caseLoUp uplo-         padLowerTriangle padUpperTriangle order (Shape.size size) $-       Array.toList m--padUpperTriangle :: Order -> Int -> [a] -> [[Maybe a]]-padUpperTriangle order n xs =-   let mxs = map Just xs-       nothings = iterate (Nothing:) []-   in case order of-         RowMajor ->-            zipWith (++) nothings (slice (take n $ iterate pred n) mxs)-         ColumnMajor ->-            transpose $-            zipWith (++)-               (slice (take n [1..]) mxs)-               (reverse $ take n nothings)--padLowerTriangle :: Order -> Int -> [a] -> [[Maybe a]]-padLowerTriangle order n xs =-   map (map Just) $-   case order of-      RowMajor -> slice (take n [1..]) xs-      ColumnMajor ->-         foldr (\(y:ys) zs -> [y] : zipWith (:) ys zs) []-            (slice (take n $ iterate pred n) xs)--slice :: [Int] -> [a] -> [[a]]-slice ns xs =-   snd $ mapAccumL (\ys n -> let (vs,ws) = splitAt n ys in (ws,vs)) xs ns--formatSeparateTriangle :: (a -> [String]) -> [[a]] -> Box-formatSeparateTriangle printFmt =-   alignSeparated . map concat .-   zipWith-      (zipWith (\sep -> attachSeparators sep . printFmt))-      (ListHT.outerProduct-         (\row col -> if row==col then Bar else Space)-         [(0::Int)..] [0..])---instance-   (Unary.Natural sub, Unary.Natural super,-    Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>-      FormatArray (MatrixShape.Banded sub super vert horiz height width) where-   formatArray fmt m =-      let MatrixShape.Banded offDiag order extent = Array.shape m-      in  formatAligned (printfFloatingMaybe fmt) $-          formatBanded offDiag order (Extent.dimensions extent) $-          Array.toList m--formatBanded ::-   (Shape.C height, Shape.C width, Unary.Natural sub, Unary.Natural super) =>-   (UnaryProxy sub, UnaryProxy super) -> Order ->-   (height, width) -> [a] -> [[Maybe a]]-formatBanded (sub,super) order (height,width) xs =-   let slices =-         ListHT.sliceVertical (MatrixShape.bandedBreadth (sub,super)) xs-       m = Shape.size height-       n = Shape.size width-   in case order of-         RowMajor ->-            map (take n) $-            zipWith (shiftRow Nothing)-               (iterate (1+) (- integralFromProxy sub))-               (map (map Just) slices)-         ColumnMajor ->-            let ku = integralFromProxy super-            in take m $ drop ku $-               foldr-                  (\col mat ->-                     zipWith (:) (map Just col ++ repeat Nothing) ([]:mat))-                  (replicate (ku + m - n) [])-                  slices---instance-   (Unary.Natural offDiag, Shape.C size) =>-      FormatArray (MatrixShape.BandedHermitian offDiag size) where-   formatArray fmt m =-      let MatrixShape.BandedHermitian offDiag order size = Array.shape m-      in  formatSeparateTriangle (printfFloatingMaybe fmt) $-          formatBandedHermitian offDiag order size $ Array.toList m--formatBandedHermitian ::-   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>-   UnaryProxy offDiag -> Order -> size -> [a] -> [[Maybe a]]-formatBandedHermitian offDiag order _size xs =-   let k = integralFromProxy offDiag-       slices = ListHT.sliceVertical (k + 1) xs-   in case order of-         RowMajor ->-            foldr-               (\row square ->-                  Match.take ([]:square) (map Just row)-                  :-                  zipWith (:)-                     (tail $ map (Just . conjugate) row ++ repeat Nothing)-                     square)-               [] slices-         ColumnMajor ->-            zipWith (shiftRow Nothing) (iterate (1+) (-k)) $ map (map Just) $-            zipWith (++)-               (map (map conjugate . init) slices)-               (drop k $-                foldr-                  (\column band ->-                     zipWith (++) (map (:[]) column ++ repeat []) ([]:band))-                  (replicate k [])-                  slices)--shiftRow :: a -> Int -> [a] -> [a]-shiftRow pad k = if k<=0 then drop (-k) else (replicate k pad ++)--splitRows ::-   (Shape.C height, Shape.C width) =>-   Order -> (height, width) -> [a] -> [[a]]-splitRows order (height,width) =-   case order of-      RowMajor -> ListHT.sliceVertical (Shape.size width)-      ColumnMajor -> ListHT.sliceHorizontal (Shape.size height)--formatAligned :: (a -> [String]) -> [[a]] -> Box-formatAligned printFmt =-   alignSeparated . map (concatMap (attachSeparators Space . printFmt))---data Separator = Empty | Space | Bar-   deriving (Eq, Ord, Show)--alignSeparated :: [[(Separator, String)]] -> Box-alignSeparated =-   TextBox.hcat TextBox.top .-   map (TextBox.vcat TextBox.right . map TextBox.text) .-   concatMap ((\(seps,column) -> [map formatSeparator seps, column]) . unzip) .-   List.unfoldr (viewLAll (Empty,""))--viewLAll :: a -> [[a]] -> Maybe ([a], [[a]])-viewLAll x0 xs =-   toMaybe (any (not.null) xs)-      (unzip $ map (fromMaybe (x0,[]) . ListHT.viewL) xs)--formatSeparator :: Separator -> String-formatSeparator sep = case sep of Empty -> ""; Space -> " "; Bar -> "|"--attachSeparators :: Separator -> [str] -> [(Separator, str)]-attachSeparators sep = zip (sep:repeat Empty)---printfFloating :: (Class.Floating a) => String -> a -> [String]-printfFloating fmt =-   getFlip $-   Class.switchFloating-      (Flip $ (:[]) . printf fmt)-      (Flip $ (:[]) . printf fmt)-      (Flip $ printfComplex fmt)-      (Flip $ printfComplex fmt)--printfFloatingMaybe :: (Class.Floating a) => String -> Maybe a -> [String]-printfFloatingMaybe fmt =-   getFlip $ getCompose $-   Class.switchFloating-      (Compose $ Flip $ (:[]) . foldMap (printf fmt))-      (Compose $ Flip $ (:[]) . foldMap (printf fmt))-      (Compose $ Flip $ maybe ["",""] (printfComplex fmt))-      (Compose $ Flip $ maybe ["",""] (printfComplex fmt))--printfComplex :: (Class.Real a) => String -> Complex a -> [String]-printfComplex fmt =-   getFlip $ getCompose $-   Class.switchReal-      (Compose $ Flip $ printfComplexAux fmt)-      (Compose $ Flip $ printfComplexAux fmt)--printfComplexAux ::-   (PrintfArg a, Class.Real a) => String -> Complex a -> [String]-printfComplexAux fmt (r:+i) =-   if i<0 || isNegativeZero i-     then [printf (fmt ++ "-") r, printf (fmt ++ "i") (-i)]-     else [printf (fmt ++ "+") r, printf (fmt ++ "i") i]+   (Type.FormatMatrix typ, Class.Floating a) =>+      Format (Type.Matrix typ a) where+   format = Type.formatMatrix

@@ -1,16 +1,19 @@+{-# LANGUAGE TypeFamilies #-} module Numeric.LAPACK.Linear.LowerUpper (    LowerUpper,-   Square,-   Transposition(..),-   Conjugation(..),-   Inversion(..),-   mapExtent,+   Plain.Square,+   Plain.Tall,+   Plain.Wide,+   Plain.Transposition(..),+   Plain.Conjugation(..),+   Plain.Inversion(..),+   Plain.mapExtent,    fromMatrix,    toMatrix,    solve,-   multiplyFullRight,+   multiplyFull, -   determinant,+   Plain.determinant,     extractP,    multiplyP,@@ -25,154 +28,55 @@    tallMultiplyU,    tallSolveU, -   caseTallWide,+   Plain.caseTallWide,    ) where -import qualified Numeric.LAPACK.Matrix.Multiply as MM-import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Linear.Plain as Plain+import Numeric.LAPACK.Linear.Plain (LowerUpper)++import qualified Numeric.LAPACK.Matrix.Array.Triangular as Tri+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent-import qualified Numeric.LAPACK.Matrix.Basic as Basic-import qualified Numeric.LAPACK.Matrix.Private as Matrix-import qualified Numeric.LAPACK.Permutation.Private as Perm-import qualified Numeric.LAPACK.Split as Split-import Numeric.LAPACK.Matrix.Triangular.Basic (UnitLower, Upper)-import Numeric.LAPACK.Matrix.Shape.Private-         (Order(RowMajor, ColumnMajor), Triangle(Triangle))-import Numeric.LAPACK.Matrix.Private-         (Full, ZeroInt, zeroInt,-          Transposition(NonTransposed, Transposed),-          Conjugation(NonConjugated, Conjugated),-          Inversion(NonInverted, Inverted), flipInversion)-import Numeric.LAPACK.Linear.Private (solver, withInfo)-import Numeric.LAPACK.Vector (Vector)-import Numeric.LAPACK.Format (Format(format))-import Numeric.LAPACK.Private-         (pointerSeq, peekCInt,-          copyBlock, copyTransposed, copyToColumnMajor)+import qualified Numeric.LAPACK.Matrix.Permutation as PermMatrix+import qualified Numeric.LAPACK.Matrix as Matrix+import Numeric.LAPACK.Matrix.Array (Full)+import Numeric.LAPACK.Matrix.Modifier (Transposition, Conjugation, Inversion) -import qualified Numeric.LAPACK.FFI.Generic as LapackGen-import qualified Numeric.BLAS.FFI.Generic as BlasGen-import qualified Numeric.Netlib.Utility as Call import qualified Numeric.Netlib.Class as Class -import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO-import qualified Data.Array.Comfort.Storable.Unchecked as Array import qualified Data.Array.Comfort.Shape as Shape-import Data.Array.Comfort.Storable.Unchecked (Array(Array)) -import Foreign.Marshal.Array (advancePtr)-import Foreign.C.Types (CInt)-import Foreign.ForeignPtr (withForeignPtr)-import Foreign.Ptr (Ptr) -import Control.Monad.Trans.Cont (ContT(ContT), evalContT)-import Control.Monad.IO.Class (liftIO)-import Control.Monad (forM_)-import Control.Applicative ((<$>))---data LowerUpper vert horiz height width a =-   LowerUpper {-      _pivot :: Vector ZeroInt CInt,-      split_ ::-         Array-            (MatrixShape.Split MatrixShape.Triangle vert horiz height width) a-   } deriving (Show)--type Square sh = LowerUpper Extent.Small Extent.Small sh sh--instance-   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,-    Class.Floating a) =>-      Format (LowerUpper vert horiz height width a) where-   format fmt lu@(LowerUpper _ipiv m) = format fmt (extractP NonInverted lu, m)--mapExtent ::-   (Extent.C vertA, Extent.C horizA) =>-   (Extent.C vertB, Extent.C horizB) =>-   Extent.Map vertA horizA vertB horizB height width ->-   LowerUpper vertA horizA height width a ->-   LowerUpper vertB horizB height width a-mapExtent f (LowerUpper pivot split) =-   LowerUpper pivot $ Array.mapShape (MatrixShape.splitMapExtent f) split- {- | @LowerUpper.fromMatrix a@ computes the LU decomposition of matrix @a@ with row pivotisation. -You can reconstruct @a@ from @lu@ depending on wether @a@ is tall or wide.+You can reconstruct @a@ from @lu@ depending on whether @a@ is tall or wide. -> LU.multiplyP False lu $ LU.extractL lu <#> LU.tallExtractU lu-> LU.multiplyP False lu $ LU.wideExtractL lu <#> LU.extractU lu+> LU.multiplyP NonInverted lu $ LU.extractL lu ##*# LU.tallExtractU lu+> LU.multiplyP NonInverted lu $ LU.wideExtractL lu #*## LU.extractU lu -} fromMatrix ::    (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,     Class.Floating a) =>    Full vert horiz height width a ->    LowerUpper vert horiz height width a-fromMatrix (Array (MatrixShape.Full order extent) a) =-   let (height,width) = Extent.dimensions extent-       m = Shape.size height-       n = Shape.size width-   in uncurry LowerUpper $-      Array.unsafeCreateWithSizeAndResult (zeroInt $ min m n) $ \_ ipivPtr ->-      ArrayIO.unsafeCreate-         (MatrixShape.Split MatrixShape.Triangle ColumnMajor extent) $ \luPtr ->--   evalContT $ do-      mPtr <- Call.cint m-      nPtr <- Call.cint n-      aPtr <- ContT $ withForeignPtr a-      ldaPtr <- Call.leadingDim m-      liftIO $ do-         copyToColumnMajor order m n aPtr luPtr-         withInfo "getrf" $ LapackGen.getrf mPtr nPtr luPtr ldaPtr ipivPtr+fromMatrix = Plain.fromMatrix . ArrMatrix.toVector  solve ::    (Extent.C vert, Extent.C horiz, Eq height, Shape.C height, Shape.C width,     Class.Floating a) =>-   Square height a ->+   Plain.Square height a ->    Full vert horiz height width a ->    Full vert horiz height width a-solve-   (LowerUpper-      (Array _ ipiv)-      (Array (MatrixShape.Split MatrixShape.Triangle orderLU extentLU) lu)) =--   solver "LowerUpper.solve" (Extent.squareSize extentLU) $-         \n nPtr nrhsPtr xPtr ldxPtr -> do-      let lda = n-      transPtr <- Call.char 'N'-      aPtr <--         case orderLU of-            RowMajor -> do-               aPtr <- ContT $ withForeignPtr lu-               atmpPtr <- Call.allocaArray (n*n)-               liftIO $ copyToColumnMajor orderLU n n aPtr atmpPtr-               return atmpPtr-            ColumnMajor -> ContT $ withForeignPtr lu-      ldaPtr <- Call.leadingDim lda-      ipivPtr <- ContT $ withForeignPtr ipiv-      liftIO $-         withInfo "getrs" $-            LapackGen.getrs transPtr-               nPtr nrhsPtr aPtr ldaPtr ipivPtr xPtr ldxPtr--{- |-Caution:-@LU.determinant . LU.fromMatrix@ will fail for singular matrices.--}-determinant :: (Shape.C sh, Class.Floating a, Eq a) => Square sh a -> a-determinant (LowerUpper ipiv split) =-   let det = Split.determinantR split-   in if Split.oddPermutation $ Array.toList ipiv then -det else det+solve = ArrMatrix.lift1 . Plain.solve   extractP ::-   (Extent.C vert, Extent.C horiz, Shape.C height) =>-   Inversion -> LowerUpper vert horiz height width a -> Perm.Permutation height-extractP inverted (LowerUpper ipiv (Array shape _)) =-   Perm.fromPivots (flipInversion inverted) (MatrixShape.splitHeight shape) ipiv+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>+   Inversion -> LowerUpper vert horiz height width a ->+   Matrix.Permutation height a+extractP inverted = PermMatrix.fromPermutation . Plain.extractP inverted  multiplyP ::    (Extent.C vertA, Extent.C horizA, Extent.C vertB, Extent.C horizB,@@ -182,56 +86,7 @@    LowerUpper vertA horizA height widthA a ->    Full vertB horizB height widthB a ->    Full vertB horizB height widthB a-multiplyP inverted-      (LowerUpper (Array shapeIPiv ipiv)-         (Array (MatrixShape.Split _ _ extentLU) _lu))-      (Array shape@(MatrixShape.Full order extent) a) =-   Array.unsafeCreate shape $ \bPtr -> do--   Call.assert "multiplyP: heights mismatch"-      (Extent.height extentLU == Extent.height extent)--   let (height,width) = Extent.dimensions extent-   let m = Shape.size height-   let n = Shape.size width-   let k = Shape.size shapeIPiv--   evalContT $ do-      aPtr <- ContT $ withForeignPtr a-      ipivPtr <- ContT $ withForeignPtr ipiv-      liftIO $ copyBlock (n*m) aPtr bPtr-      case order of-         ColumnMajor -> do-            nPtr <- Call.cint n-            ldaPtr <- Call.leadingDim m-            k1Ptr <- Call.cint 1-            k2Ptr <- Call.cint k-            incxPtr <--               Call.cint $-               case inverted of-                  Inverted -> 1-                  NonInverted -> -1-            liftIO $-               LapackGen.laswp nPtr bPtr ldaPtr k1Ptr k2Ptr ipivPtr incxPtr-         RowMajor ->-            liftIO $ swapColumns m bPtr $ take k $-            case inverted of-               Inverted -> zip [0..] $ pointerSeq 1 ipivPtr-               NonInverted ->-                  zip (iterate (subtract 1) (k-1)) $-                  pointerSeq (-1) (advancePtr ipivPtr (k-1))--{-# INLINE swapColumns #-}-swapColumns ::-   (Class.Floating a) =>-   Int -> Ptr a -> [(Int, Ptr CInt)] -> IO ()-swapColumns m xPtr ptrs = evalContT $ do-   mPtr <- Call.cint m-   incPtr <- Call.cint 1-   let columnPtr k = advancePtr xPtr (m*k)-   liftIO $ forM_ ptrs $ \(i,ipivPtr) -> do-      j <- subtract 1 <$> peekCInt ipivPtr-      BlasGen.swap mPtr (columnPtr i) incPtr (columnPtr j) incPtr+multiplyP inverted = ArrMatrix.lift1 . Plain.multiplyP inverted   @@ -240,19 +95,19 @@     Class.Floating a) =>    LowerUpper vert horiz height width a ->    Full vert horiz height width a-extractL = Split.extractTriangle (Left Triangle) . split_+extractL = ArrMatrix.lift0 . Plain.extractL  wideExtractL ::    (Extent.C horiz, Shape.C height, Shape.C width, Class.Floating a) =>-   LowerUpper Extent.Small horiz height width a -> UnitLower height a-wideExtractL = Split.wideExtractL . split_+   LowerUpper Extent.Small horiz height width a -> Tri.UnitLower height a+wideExtractL = ArrMatrix.lift0 . Plain.wideExtractL  {- | @wideMultiplyL transposed lu a@ multiplies the square part of @lu@ containing the lower triangular matrix with @a@. -> wideMultiplyL False lu a == wideExtractL lu <#> a-> wideMultiplyL True lu a == wideExtractL (Tri.transposeUp lu) <#> a+> wideMultiplyL NonTransposed lu a == wideExtractL lu #*## a+> wideMultiplyL Transposed lu a == Tri.transpose (wideExtractL lu) #*## a -} wideMultiplyL ::    (Extent.C horizA, Extent.C vert, Extent.C horiz, Shape.C height, Eq height,@@ -261,7 +116,7 @@    LowerUpper Extent.Small horizA height widthA a ->    Full vert horiz height widthB a ->    Full vert horiz height widthB a-wideMultiplyL transposed = Split.wideMultiplyL transposed . split_+wideMultiplyL transposed = ArrMatrix.lift1 . Plain.wideMultiplyL transposed  wideSolveL ::    (Extent.C horizA, Extent.C vert, Extent.C horiz,@@ -270,7 +125,7 @@    LowerUpper Extent.Small horizA height width a ->    Full vert horiz height nrhs a -> Full vert horiz height nrhs a wideSolveL transposed conjugated =-   Split.wideSolveL transposed conjugated . split_+   ArrMatrix.lift1 . Plain.wideSolveL transposed conjugated   extractU ::@@ -278,19 +133,19 @@     Class.Floating a) =>    LowerUpper vert horiz height width a ->    Full vert horiz height width a-extractU = Split.extractTriangle (Right Triangle) . split_+extractU = ArrMatrix.lift0 . Plain.extractU  tallExtractU ::    (Extent.C vert, Shape.C height, Shape.C width, Class.Floating a) =>-   LowerUpper vert Extent.Small height width a -> Upper width a-tallExtractU = Split.tallExtractR . split_+   LowerUpper vert Extent.Small height width a -> Tri.Upper width a+tallExtractU = ArrMatrix.lift0 . Plain.tallExtractU  {- | @tallMultiplyU transposed lu a@ multiplies the square part of @lu@ containing the upper triangular matrix with @a@. -> tallMultiplyU False lu a == tallExtractU lu <#> a-> tallMultiplyU True lu a == tallExtractU (Tri.transposeDown lu) <#> a+> tallMultiplyU NonTransposed lu a == tallExtractU lu #*## a+> tallMultiplyU Transposed lu a == Tri.transpose (tallExtractU lu) #*## a -} tallMultiplyU ::    (Extent.C vertA, Extent.C vert, Extent.C horiz, Shape.C height, Eq height,@@ -299,7 +154,7 @@    LowerUpper vertA Extent.Small heightA height a ->    Full vert horiz height widthB a ->    Full vert horiz height widthB a-tallMultiplyU transposed = Split.tallMultiplyR transposed . split_+tallMultiplyU transposed = ArrMatrix.lift1 . Plain.tallMultiplyU transposed  tallSolveU ::    (Extent.C vertA, Extent.C vert, Extent.C horiz,@@ -308,7 +163,7 @@    LowerUpper vertA Extent.Small height width a ->    Full vert horiz width nrhs a -> Full vert horiz width nrhs a tallSolveU transposed conjugated =-   Split.tallSolveR transposed conjugated . split_+   ArrMatrix.lift1 . Plain.tallSolveU transposed conjugated   @@ -317,121 +172,14 @@     Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>    LowerUpper vert horiz height width a ->    Full vert horiz height width a-toMatrix =-   getToMatrix $-   Extent.switchTagPair-      (ToMatrix wideToMatrix)-      (ToMatrix wideToMatrix)-      (ToMatrix tallToMatrix)-      (ToMatrix $-         either-            (Matrix.fromFull . tallToMatrix)-            (Matrix.fromFull . wideToMatrix) .-         caseTallWide)--newtype ToMatrix height width a vert horiz =-   ToMatrix {-      getToMatrix ::-         LowerUpper vert horiz height width a ->-         Full vert horiz height width a-   }--tallToMatrix ::-   (Extent.C vert, Shape.C height, Shape.C width, Eq height, Eq width,-    Class.Floating a) =>-   LowerUpper vert Extent.Small height width a ->-   Full vert Extent.Small height width a-tallToMatrix a =-   multiplyP NonInverted a $ Basic.transpose $-   tallMultiplyU Transposed a $ Basic.transpose $ extractL a--wideToMatrix ::-   (Extent.C horiz, Shape.C height, Shape.C width, Eq height, Eq width,-    Class.Floating a) =>-   LowerUpper Extent.Small horiz height width a ->-   Full Extent.Small horiz height width a-wideToMatrix a =-   multiplyP NonInverted a $ wideMultiplyL NonTransposed a $ extractU a+toMatrix = ArrMatrix.lift0 . Plain.toMatrix  -multiplyFullRight ::+multiplyFull ::    (Extent.C vert, Extent.C horiz,     Shape.C height, Eq height, Shape.C width, Shape.C fuse, Eq fuse,     Class.Floating a) =>    LowerUpper vert horiz height fuse a ->    Full vert horiz fuse width a ->    Full vert horiz height width a-multiplyFullRight =-   getMultiplyFullRight $-   Extent.switchTagPair-      (MultiplyFullRight wideMultiplyFullRight)-      (MultiplyFullRight wideMultiplyFullRight)-      (MultiplyFullRight tallMultiplyFullRight)-      (MultiplyFullRight $-         either tallMultiplyFullRight wideMultiplyFullRight . caseTallWide)--newtype MultiplyFullRight height fuse width a vert horiz =-   MultiplyFullRight {-      getMultiplyFullRight ::-         LowerUpper vert horiz height fuse a ->-         Full vert horiz fuse width a ->-         Full vert horiz height width a-   }--tallMultiplyFullRight ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Shape.C fuse, Eq height, Eq fuse,-    Class.Floating a) =>-   LowerUpper vert Extent.Small height fuse a ->-   Full vert horiz fuse width a ->-   Full vert horiz height width a-tallMultiplyFullRight a =-   multiplyP NonInverted a .-   MM.multiply (Matrix.generalizeTall (extractL a)) .-   tallMultiplyU NonTransposed a--wideMultiplyFullRight ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Shape.C fuse, Eq height, Eq fuse,-    Class.Floating a) =>-   LowerUpper Extent.Small horiz height fuse a ->-   Full vert horiz fuse width a ->-   Full vert horiz height width a-wideMultiplyFullRight a =-   multiplyP NonInverted a . wideMultiplyL NonTransposed a .-   MM.multiply (Matrix.generalizeWide (extractU a))---type Tall = LowerUpper Extent.Big Extent.Small-type Wide = LowerUpper Extent.Small Extent.Big--caseTallWide ::-   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>-   LowerUpper vert horiz height width a ->-   Either (Tall height width a) (Wide height width a)-caseTallWide (LowerUpper ipiv (Array shape a)) =-   either-      (Left . LowerUpper ipiv . flip Array a)-      (Right . LowerUpper ipiv . flip Array a) $-   MatrixShape.caseTallWideSplit shape---_toRowMajor ::-   (Extent.C vert, Extent.C horiz, Eq height, Shape.C height, Shape.C width,-    Class.Floating a) =>-   LowerUpper vert horiz height width a ->-   LowerUpper vert horiz height width a-_toRowMajor-   (LowerUpper ipiv-      arr@(Array (MatrixShape.Split MatrixShape.Triangle order extent) a)) =-   LowerUpper ipiv $-   case order of-      RowMajor -> arr-      ColumnMajor ->-         Array.unsafeCreate-            (MatrixShape.Split MatrixShape.Triangle RowMajor extent) $ \bPtr ->-         withForeignPtr a $ \aPtr -> do-            let (height, width) = Extent.dimensions extent-            let n = Shape.size width-            let m = Shape.size height-            copyTransposed n m aPtr n bPtr+multiplyFull = ArrMatrix.lift1 . Plain.multiplyFull

@@ -0,0 +1,522 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE EmptyDataDecls #-}+module Numeric.LAPACK.Linear.Plain (+   LowerUpper,+   Square, Tall, Wide,+   Transposition(..),+   Conjugation(..),+   Inversion(..),+   mapExtent,+   fromMatrix,+   toMatrix,+   solve,+   multiplyFull,++   determinant,++   extractP,+   multiplyP,++   extractL,+   wideExtractL,+   wideMultiplyL,+   wideSolveL,++   extractU,+   tallExtractU,+   tallMultiplyU,+   tallSolveU,++   caseTallWide,+   ) where++import qualified Numeric.LAPACK.Matrix.Divide as Divide+import qualified Numeric.LAPACK.Matrix.Multiply as Multiply+import qualified Numeric.LAPACK.Matrix.Type as Type+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent+import qualified Numeric.LAPACK.Matrix.Extent as ExtentMap+import qualified Numeric.LAPACK.Matrix.Basic as Basic+import qualified Numeric.LAPACK.Matrix.Private as Matrix+import qualified Numeric.LAPACK.Permutation.Private as Perm+import qualified Numeric.LAPACK.Shape as ExtShape+import qualified Numeric.LAPACK.Split as Split+import Numeric.LAPACK.Output ((/+/))+import Numeric.LAPACK.Matrix.Plain.Format (formatArray)+import Numeric.LAPACK.Matrix.Type (FormatMatrix(formatMatrix))+import Numeric.LAPACK.Matrix.Triangular.Basic (UnitLower, Upper)+import Numeric.LAPACK.Matrix.Shape.Private+         (Order(RowMajor, ColumnMajor), Triangle(Triangle))+import Numeric.LAPACK.Matrix.Modifier+         (Transposition(NonTransposed, Transposed),+          Conjugation(NonConjugated, Conjugated),+          Inversion(NonInverted, Inverted))+import Numeric.LAPACK.Matrix.Private (Full)+import Numeric.LAPACK.Linear.Private (solver, withInfo)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Private (copyBlock, copyTransposed, copyToColumnMajor)++import qualified Numeric.LAPACK.FFI.Generic as LapackGen+import qualified Numeric.BLAS.FFI.Generic as BlasGen+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO+import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Unchecked (Array(Array), (!))++import Foreign.Marshal.Array (advancePtr)+import Foreign.ForeignPtr (withForeignPtr)+import Foreign.Ptr (Ptr)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Monad (forM_)++import Data.Monoid ((<>))+++data LU vert horiz height width++data instance Type.Matrix (LU vert horiz height width) a =+   LowerUpper {+      _pivot ::+         Vector (ExtShape.Min width (Perm.Shape height)) (Perm.Element height),+      split_ ::+         Array+            (MatrixShape.Split MatrixShape.Triangle vert horiz height width) a+   } deriving (Show)++type LowerUpper vert horiz height width =+         Type.Matrix (LU vert horiz height width)+type Square sh = LowerUpper Extent.Small Extent.Small sh sh+type Tall height width = LowerUpper Extent.Big Extent.Small height width+type Wide height width = LowerUpper Extent.Small Extent.Big height width++instance+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>+      FormatMatrix (LU vert horiz height width) where+   formatMatrix fmt lu@(LowerUpper _ipiv m) =+      Perm.format (extractP NonInverted lu)+      /+/+      formatArray fmt m++mapExtent ::+   (Extent.C vertA, Extent.C horizA) =>+   (Extent.C vertB, Extent.C horizB) =>+   Extent.Map vertA horizA vertB horizB height width ->+   LowerUpper vertA horizA height width a ->+   LowerUpper vertB horizB height width a+mapExtent f (LowerUpper pivot split) =+   LowerUpper pivot $ Array.mapShape (MatrixShape.splitMapExtent f) split++fromMatrix ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,+    Class.Floating a) =>+   Full vert horiz height width a ->+   LowerUpper vert horiz height width a+fromMatrix (Array (MatrixShape.Full order extent) a) =+   let (height,width) = Extent.dimensions extent+   in uncurry LowerUpper $+      Array.unsafeCreateWithSizeAndResult+         (ExtShape.Min width $ Perm.Shape height) $ \_ ipivPtr ->+      ArrayIO.unsafeCreate+         (MatrixShape.Split MatrixShape.Triangle ColumnMajor extent) $ \luPtr ->++   evalContT $ do+      let m = Shape.size height+      let n = Shape.size width+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      aPtr <- ContT $ withForeignPtr a+      ldaPtr <- Call.leadingDim m+      liftIO $ do+         copyToColumnMajor order m n aPtr luPtr+         withInfo "getrf" $+            LapackGen.getrf mPtr nPtr luPtr ldaPtr+               (Perm.deconsElementPtr ipivPtr)++solve ::+   (Extent.C vert, Extent.C horiz, Eq height, Shape.C height, Shape.C width,+    Class.Floating a) =>+   Square height a ->+   Full vert horiz height width a ->+   Full vert horiz height width a+solve = solveTrans NonTransposed++solveTrans ::+   (Extent.C vert, Extent.C horiz, Eq height, Shape.C height, Shape.C width,+    Class.Floating a) =>+   Transposition -> Square height a ->+   Full vert horiz height width a ->+   Full vert horiz height width a+solveTrans trans+   (LowerUpper+      (Array _ ipiv)+      (Array (MatrixShape.Split MatrixShape.Triangle orderLU extentLU) lu)) =++   solver "LowerUpper.solve" (Extent.squareSize extentLU) $+         \n nPtr nrhsPtr xPtr ldxPtr -> do+      let lda = n+      transPtr <- Call.char $+         case trans of+            NonTransposed -> 'N'+            Transposed -> 'T'+      aPtr <-+         case orderLU of+            RowMajor -> do+               aPtr <- ContT $ withForeignPtr lu+               atmpPtr <- Call.allocaArray (n*n)+               liftIO $ copyToColumnMajor orderLU n n aPtr atmpPtr+               return atmpPtr+            ColumnMajor -> ContT $ withForeignPtr lu+      ldaPtr <- Call.leadingDim lda+      ipivPtr <- fmap Perm.deconsElementPtr $ ContT $ withForeignPtr ipiv+      liftIO $+         withInfo "getrs" $+            LapackGen.getrs transPtr+               nPtr nrhsPtr aPtr ldaPtr ipivPtr xPtr ldxPtr++{- |+Caution:+@LU.determinant . LU.fromMatrix@ will fail for singular matrices.+-}+determinant :: (Shape.C sh, Class.Floating a) => Square sh a -> a+determinant (LowerUpper ipiv split) =+   Perm.condNegate (map Perm.deconsElement $ Array.toList ipiv) $+   Split.determinantR split+++extractP ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>+   Inversion -> LowerUpper vert horiz height width a -> Perm.Permutation height+extractP inverted (LowerUpper ipiv _) =+   Perm.fromTruncatedPivots (Inverted <> inverted) ipiv++multiplyP ::+   (Extent.C vertA, Extent.C horizA, Extent.C vertB, Extent.C horizB,+    Eq height, Shape.C height, Shape.C widthA, Shape.C widthB,+    Class.Floating a) =>+   Inversion ->+   LowerUpper vertA horizA height widthA a ->+   Full vertB horizB height widthB a ->+   Full vertB horizB height widthB a+multiplyP inverted+      (LowerUpper ipiv@(Array shapeIPiv ipivFPtr)+         (Array (MatrixShape.Split _ _ extentLU) _lu))+      (Array shape@(MatrixShape.Full order extent) a) =+   Array.unsafeCreate shape $ \bPtr -> do++   Call.assert "multiplyP: heights mismatch"+      (Extent.height extentLU == Extent.height extent)++   let (height,width) = Extent.dimensions extent+   let m = Shape.size height+   let n = Shape.size width+   let k = Shape.size shapeIPiv++   evalContT $ do+      aPtr <- ContT $ withForeignPtr a+      ipivPtr <- ContT $ withForeignPtr ipivFPtr+      liftIO $ copyBlock (n*m) aPtr bPtr+      case order of+         ColumnMajor -> do+            nPtr <- Call.cint n+            ldaPtr <- Call.leadingDim m+            k1Ptr <- Call.cint 1+            k2Ptr <- Call.cint k+            incxPtr <-+               Call.cint $+               case inverted of+                  Inverted -> 1+                  NonInverted -> -1+            liftIO $+               LapackGen.laswp nPtr bPtr ldaPtr k1Ptr k2Ptr+                  (Perm.deconsElementPtr ipivPtr) incxPtr+         RowMajor ->+            liftIO $ swapColumns height n bPtr ipiv $+            Perm.indices (Inverted <> inverted) shapeIPiv++{-# INLINE swapColumns #-}+swapColumns ::+   (Shape.C sh, Shape.C width, Class.Floating a) =>+   sh -> Int -> Ptr a ->+   Array (ExtShape.Min width (Perm.Shape sh)) (Perm.Element sh) ->+   [Perm.Element sh] -> IO ()+swapColumns sh n xPtr ipiv is = evalContT $ do+   nPtr <- Call.cint n+   incPtr <- Call.cint 1+   let columnPtr ix =+         advancePtr xPtr (n * Shape.uncheckedOffset (Perm.Shape sh) ix)+   liftIO $ forM_ is $ \i ->+      BlasGen.swap nPtr (columnPtr i) incPtr (columnPtr (ipiv!i)) incPtr++++extractL ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,+    Class.Floating a) =>+   LowerUpper vert horiz height width a ->+   Full vert horiz height width a+extractL = Split.extractTriangle (Left Triangle) . split_++wideExtractL ::+   (Extent.C horiz, Shape.C height, Shape.C width, Class.Floating a) =>+   LowerUpper Extent.Small horiz height width a -> UnitLower height a+wideExtractL = Split.wideExtractL . split_++wideMultiplyL ::+   (Extent.C horizA, Extent.C vert, Extent.C horiz, Shape.C height, Eq height,+    Shape.C widthA, Shape.C widthB, Class.Floating a) =>+   Transposition ->+   LowerUpper Extent.Small horizA height widthA a ->+   Full vert horiz height widthB a ->+   Full vert horiz height widthB a+wideMultiplyL transposed = Split.wideMultiplyL transposed . split_++wideSolveL ::+   (Extent.C horizA, Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>+   Transposition -> Conjugation ->+   LowerUpper Extent.Small horizA height width a ->+   Full vert horiz height nrhs a -> Full vert horiz height nrhs a+wideSolveL transposed conjugated =+   Split.wideSolveL transposed conjugated . split_+++extractU ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,+    Class.Floating a) =>+   LowerUpper vert horiz height width a ->+   Full vert horiz height width a+extractU = Split.extractTriangle (Right Triangle) . split_++tallExtractU ::+   (Extent.C vert, Shape.C height, Shape.C width, Class.Floating a) =>+   LowerUpper vert Extent.Small height width a -> Upper width a+tallExtractU = Split.tallExtractR . split_++tallMultiplyU ::+   (Extent.C vertA, Extent.C vert, Extent.C horiz, Shape.C height, Eq height,+    Shape.C heightA, Shape.C widthB, Class.Floating a) =>+   Transposition ->+   LowerUpper vertA Extent.Small heightA height a ->+   Full vert horiz height widthB a ->+   Full vert horiz height widthB a+tallMultiplyU transposed = Split.tallMultiplyR transposed . split_++tallSolveU ::+   (Extent.C vertA, Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Eq width, Shape.C nrhs, Class.Floating a) =>+   Transposition -> Conjugation ->+   LowerUpper vertA Extent.Small height width a ->+   Full vert horiz width nrhs a -> Full vert horiz width nrhs a+tallSolveU transposed conjugated =+   Split.tallSolveR transposed conjugated . split_++++toMatrix ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+   LowerUpper vert horiz height width a ->+   Full vert horiz height width a+toMatrix =+   getToMatrix $+   Extent.switchTagPair+      (ToMatrix wideToMatrix)+      (ToMatrix wideToMatrix)+      (ToMatrix tallToMatrix)+      (ToMatrix $+         either+            (Matrix.fromFull . tallToMatrix)+            (Matrix.fromFull . wideToMatrix) .+         caseTallWide)++newtype ToMatrix height width a vert horiz =+   ToMatrix {+      getToMatrix ::+         LowerUpper vert horiz height width a ->+         Full vert horiz height width a+   }++tallToMatrix ::+   (Extent.C vert, Shape.C height, Shape.C width, Eq height, Eq width,+    Class.Floating a) =>+   LowerUpper vert Extent.Small height width a ->+   Full vert Extent.Small height width a+tallToMatrix a =+   multiplyP NonInverted a $ Basic.transpose $+   tallMultiplyU Transposed a $ Basic.transpose $ extractL a++wideToMatrix ::+   (Extent.C horiz, Shape.C height, Shape.C width, Eq height, Eq width,+    Class.Floating a) =>+   LowerUpper Extent.Small horiz height width a ->+   Full Extent.Small horiz height width a+wideToMatrix a =+   multiplyP NonInverted a $ wideMultiplyL NonTransposed a $ extractU a+++multiplyFull ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Shape.C fuse, Eq fuse,+    Class.Floating a) =>+   LowerUpper vert horiz height fuse a ->+   Full vert horiz fuse width a ->+   Full vert horiz height width a+multiplyFull =+   getMultiplyFullRight $+   Extent.switchTagPair+      {-+      We cannot simply use squareFull here,+      because this requires height~width.+      -}+      (MultiplyFullRight wideMultiplyFullRight)+      (MultiplyFullRight wideMultiplyFullRight)+      (MultiplyFullRight tallMultiplyFullRight)+      (MultiplyFullRight $+         either tallMultiplyFullRight wideMultiplyFullRight . caseTallWide)++newtype MultiplyFullRight height fuse width a vert horiz =+   MultiplyFullRight {+      getMultiplyFullRight ::+         LowerUpper vert horiz height fuse a ->+         Full vert horiz fuse width a ->+         Full vert horiz height width a+   }++tallMultiplyFullRight ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Shape.C fuse, Eq height, Eq fuse,+    Class.Floating a) =>+   LowerUpper vert Extent.Small height fuse a ->+   Full vert horiz fuse width a ->+   Full vert horiz height width a+tallMultiplyFullRight a =+   multiplyP NonInverted a .+   Basic.multiply (Matrix.generalizeTall (extractL a)) .+   tallMultiplyU NonTransposed a++wideMultiplyFullRight ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Shape.C fuse, Eq height, Eq fuse,+    Class.Floating a) =>+   LowerUpper Extent.Small horiz height fuse a ->+   Full vert horiz fuse width a ->+   Full vert horiz height width a+wideMultiplyFullRight a =+   multiplyP NonInverted a . wideMultiplyL NonTransposed a .+   Basic.multiply (Matrix.generalizeWide (extractU a))+++tallTransMultiplyFullRight ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Shape.C fuse, Eq height, Eq fuse,+    Class.Floating a) =>+   LowerUpper horiz Extent.Small fuse height a ->+   Full vert horiz fuse width a ->+   Full vert horiz height width a+tallTransMultiplyFullRight a =+   tallMultiplyU Transposed a .+   Basic.multiply (Basic.transpose $ Matrix.generalizeTall $ extractL a) .+   multiplyP Inverted a++wideTransMultiplyFullRight ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Shape.C fuse, Eq height, Eq fuse,+    Class.Floating a) =>+   LowerUpper Extent.Small vert fuse height a ->+   Full vert horiz fuse width a ->+   Full vert horiz height width a+wideTransMultiplyFullRight a =+   Basic.multiply (Basic.transpose $ Matrix.generalizeWide $ extractU a) .+   wideMultiplyL Transposed a .+   multiplyP Inverted a+++caseTallWide ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>+   LowerUpper vert horiz height width a ->+   Either (Tall height width a) (Wide height width a)+caseTallWide (LowerUpper ipiv (Array shape a)) =+   either+      (Left . LowerUpper ipiv . flip Array a)+      (Right . LowerUpper ipiv . flip Array a) $+   MatrixShape.caseTallWideSplit shape+++_toRowMajor ::+   (Extent.C vert, Extent.C horiz, Eq height, Shape.C height, Shape.C width,+    Class.Floating a) =>+   LowerUpper vert horiz height width a ->+   LowerUpper vert horiz height width a+_toRowMajor+   (LowerUpper ipiv+      arr@(Array (MatrixShape.Split MatrixShape.Triangle order extent) a)) =+   LowerUpper ipiv $+   case order of+      RowMajor -> arr+      ColumnMajor ->+         Array.unsafeCreate+            (MatrixShape.Split MatrixShape.Triangle RowMajor extent) $ \bPtr ->+         withForeignPtr a $ \aPtr -> do+            let (height, width) = Extent.dimensions extent+            let n = Shape.size width+            let m = Shape.size height+            copyTransposed n m aPtr n bPtr+++instance+   (Extent.C vert, Extent.C horiz) =>+      Type.Box (LU vert horiz height width) where+   type HeightOf (LU vert horiz height width) = height+   type WidthOf (LU vert horiz height width) = width+   height = MatrixShape.splitHeight . Array.shape . split_+   width = MatrixShape.splitWidth . Array.shape . split_++instance+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Eq width) =>+      Multiply.MultiplyVector (LU vert horiz height width) where+   matrixVector lu x =+      Basic.unliftColumn MatrixShape.ColumnMajor+         (multiplyFull (mapExtent ExtentMap.toGeneral lu)) x+   vectorMatrix x lu =+      Basic.unliftColumn MatrixShape.ColumnMajor+         (either tallTransMultiplyFullRight wideTransMultiplyFullRight $+          caseTallWide lu)+         x++instance+   (vert ~ Extent.Small, horiz ~ Extent.Small,+    Shape.C height, height ~ width) =>+      Multiply.MultiplySquare (LU vert horiz height width) where++   squareFull lu =+      ArrMatrix.lift1 $+         multiplyP NonInverted lu .+         wideMultiplyL NonTransposed lu .+         tallMultiplyU NonTransposed lu++   fullSquare = flip $ \lu ->+      ArrMatrix.lift1 $+         Basic.transpose .+         tallMultiplyU Transposed lu .+         wideMultiplyL Transposed lu .+         multiplyP Inverted lu .+         Basic.transpose++instance+   (vert ~ Extent.Small, horiz ~ Extent.Small,+    Shape.C height, height ~ width) =>+      Divide.Determinant (LU vert horiz height width) where+   determinant = determinant++instance+   (vert ~ Extent.Small, horiz ~ Extent.Small,+    Shape.C height, height ~ width) =>+      Divide.Solve (LU vert horiz height width) where+   solve trans = ArrMatrix.lift1 . solveTrans trans

@@ -1,926 +1,768 @@ {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE TypeOperators #-} module Numeric.LAPACK.Matrix (-   Full,-   General, Tall, Wide,-   ZeroInt, zeroInt,-   transpose, adjoint,-   Matrix.height, Matrix.width,-   caseTallWide,-   fromScalar, toScalar,-   fromList,-   mapExtent, fromFull,-   generalizeTall, generalizeWide,-   mapHeight, mapWidth,-   identity,-   diagonal,-   fromRowsNonEmpty,    fromRowArray,    fromRows,-   fromColumnsNonEmpty, fromColumnArray, fromColumns,-   Basic.singleRow,   Basic.singleColumn,-   Basic.flattenRow,  Basic.flattenColumn,-   Basic.liftRow,     Basic.liftColumn,-   Basic.unliftRow,   Basic.unliftColumn,-   toRows, toColumns,-   toRowArray, toColumnArray,-   takeRow, takeColumn,-   takeRows, takeColumns, takeEqually,-   dropRows, dropColumns, dropEqually,-   takeTopRows, takeBottomRows,-   takeLeftColumns, takeRightColumns,-   takeRowArray, takeColumnArray,-   swapRows, swapColumns,-   reverseRows, reverseColumns,-   fromRowMajor, toRowMajor, flatten,-   forceOrder, adaptOrder,-   (|||),-   (===),--   tensorProduct,-   outer,-   sumRank1,--   RealOf,-   add, sub,-   rowSums, columnSums,-   scaleRows, scaleColumns,-   scaleRowsComplex, scaleColumnsComplex,-   scaleRowsReal, scaleColumnsReal,-   multiply,-   multiplyVector,--   Multiply, (<#>),-   MultiplyLeft, (<#),-   MultiplyRight, (#>),-   Indexed, (#!),--   Solve, solve, solveVector,-   Inverse, inverse,-   ) where--import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape-import qualified Numeric.LAPACK.Matrix.Square.Basic as Square-import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent-import qualified Numeric.LAPACK.Matrix.Basic as Basic-import qualified Numeric.LAPACK.Matrix.Private as Matrix-import qualified Numeric.LAPACK.Vector as Vector-import qualified Numeric.LAPACK.Private as Private-import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor, ColumnMajor))-import Numeric.LAPACK.Matrix.Multiply-         (Multiply((<#>)), MultiplyLeft((<#)), MultiplyRight((#>)),-          multiplyVector, multiply, multiplyVectorUnchecked)-import Numeric.LAPACK.Matrix.Indexed (Indexed((#!)))-import Numeric.LAPACK.Matrix.Divide-         (Solve(solve), solveVector, Inverse(inverse))-import Numeric.LAPACK.Matrix.Basic-         (transpose, forceOrder, forceRowMajor, scaleRows, scaleColumns)-import Numeric.LAPACK.Matrix.Private-         (Full, Tall, Wide, General, argGeneral, ZeroInt, zeroInt,-          mapExtent, fromFull, generalizeTall, generalizeWide)-import Numeric.LAPACK.Vector (Vector)-import Numeric.LAPACK.Scalar (RealOf, zero, one)-import Numeric.LAPACK.Private-         (pointerSeq, fill, copyTransposed, copySubMatrix, copyBlock)--import qualified Numeric.LAPACK.FFI.Generic as LapackGen-import qualified Numeric.BLAS.FFI.Generic as BlasGen-import qualified Numeric.Netlib.Utility as Call-import qualified Numeric.Netlib.Class as Class--import qualified Data.Array.Comfort.Boxed as BoxedArray-import qualified Data.Array.Comfort.Storable.Unchecked as Array-import qualified Data.Array.Comfort.Storable as CheckedArray-import qualified Data.Array.Comfort.Shape as Shape-import Data.Array.Comfort.Storable.Unchecked (Array(Array))-import Data.Array.Comfort.Shape ((:+:)((:+:)))--import Foreign.Marshal.Array (copyArray, advancePtr, pokeArray)-import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)-import Foreign.Ptr (Ptr, castPtr)-import Foreign.Storable (Storable, poke, peek)--import System.IO.Unsafe (unsafePerformIO)--import Control.Monad.Trans.Cont (ContT(ContT), evalContT)-import Control.Monad.IO.Class (liftIO)-import Control.Monad (when, mfilter)--import qualified Data.NonEmpty.Mixed as NonEmptyM-import qualified Data.NonEmpty as NonEmpty-import Data.Complex (Complex)-import Data.Foldable (forM_)-import Data.Maybe (listToMaybe)-import Data.Bool.HT (if')---{- |-conjugate transpose--Problem: @adjoint a <#> a@ is always square,-but how to convince the type checker to choose the Square type?-Anser: Use @Hermitian.toSquare $ Hermitian.covariance a@ instead.--}-adjoint ::-   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,-    Class.Floating a) =>-   Full vert horiz height width a -> Full horiz vert width height a-adjoint = transpose . Vector.conjugate---{- |-Square matrices will be classified as 'Tall'.--}-caseTallWide ::-   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>-   Full vert horiz height width a ->-   Either (Tall height width a) (Wide height width a)-caseTallWide (Array shape a) =-   either (Left . flip Array a) (Right . flip Array a) $-   MatrixShape.caseTallWide shape---fromScalar :: (Storable a) => a -> General () () a-fromScalar = Square.toGeneral . Square.fromScalar--toScalar :: (Storable a) => General () () a -> a-toScalar = argGeneral $ \_ () () a ->-   unsafePerformIO $ withForeignPtr a peek--fromList ::-   (Shape.C height, Shape.C width, Storable a) =>-   height -> width -> [a] -> General height width a-fromList height width =-   CheckedArray.fromList (MatrixShape.general RowMajor height width)---identity ::-   (Shape.C sh, Class.Floating a) =>-   sh -> General sh sh a-identity = Square.toGeneral . Square.identity--diagonal ::-   (Shape.C sh, Class.Floating a) =>-   Vector sh a -> General sh sh a-diagonal = Square.toGeneral . Square.diagonal---{- |-The number of rows must be maintained by the height mapping function.--}-mapHeight ::-   (Shape.C heightA, Shape.C heightB,-    Extent.GeneralTallWide vert horiz,-    Extent.GeneralTallWide horiz vert) =>-   (heightA -> heightB) ->-   Full vert horiz heightA width a ->-   Full vert horiz heightB width a-mapHeight f = Basic.mapHeight $ withSizeCheck "mapHeight" f--{- |-The number of columns must be maintained by the width mapping function.--}-mapWidth ::-   (Shape.C widthA, Shape.C widthB,-    Extent.GeneralTallWide vert horiz,-    Extent.GeneralTallWide horiz vert) =>-   (widthA -> widthB) ->-   Full vert horiz height widthA a ->-   Full vert horiz height widthB a-mapWidth f = Basic.mapWidth $ withSizeCheck "mapWidth" f--withSizeCheck ::-   (Shape.C sha, Shape.C shb) =>-   String -> (sha -> shb) -> sha -> shb-withSizeCheck name f sizeA =-   let sizeB = f sizeA-   in if Shape.size sizeA == Shape.size sizeB-         then sizeB-         else error $ name ++ ": sizes mismatch"---fromRowsNonEmpty ::-   (Shape.C width, Eq width, Storable a) =>-   NonEmpty.T [] (Vector width a) -> General ZeroInt width a-fromRowsNonEmpty (NonEmpty.Cons row rows) =-   fromRows (Array.shape row) (row:rows)--fromRowArray ::-   (Shape.C height, Shape.C width, Eq width, Storable a) =>-   width -> BoxedArray.Array height (Vector width a) -> General height width a-fromRowArray width rows =-   Basic.mapHeight (const $ BoxedArray.shape rows) $-   fromRows width $ BoxedArray.toList rows--fromRows ::-   (Shape.C width, Eq width, Storable a) =>-   width -> [Vector width a] -> General ZeroInt width a-fromRows width rows =-   Array.unsafeCreate-      (MatrixShape.general RowMajor (zeroInt $ length rows) width)-      (gather width rows)--fromColumnsNonEmpty ::-   (Shape.C height, Eq height, Storable a) =>-   NonEmpty.T [] (Vector height a) -> General height ZeroInt a-fromColumnsNonEmpty (NonEmpty.Cons column columns) =-   fromColumns (Array.shape column) (column:columns)--fromColumnArray ::-   (Shape.C height, Eq height, Shape.C width, Storable a) =>-   height -> BoxedArray.Array width (Vector height a) -> General height width a-fromColumnArray height columns =-   Basic.mapWidth (const $ BoxedArray.shape columns) $-   fromColumns height $ BoxedArray.toList columns--fromColumns ::-   (Shape.C height, Eq height, Storable a) =>-   height -> [Vector height a] -> General height ZeroInt a-fromColumns height columns =-   Array.unsafeCreate-      (MatrixShape.general ColumnMajor height (zeroInt $ length columns))-      (gather height columns)--gather ::-   (Shape.C width, Eq width, Storable a) =>-   width -> [Array width a] -> Ptr a -> IO ()-gather width rows dstPtr =-   let widthSize = Shape.size width-   in forM_ (zip (pointerSeq widthSize dstPtr) rows) $-         \(dstRowPtr, Array.Array rowWidth srcFPtr) ->-         withForeignPtr srcFPtr $ \srcPtr -> do-            Call.assert-               "Matrix.fromRows/fromColumnsNonEmpty: non-matching vector size"-               (width == rowWidth)-            copyArray dstRowPtr srcPtr widthSize---toRows ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Class.Floating a) =>-   Full vert horiz height width a -> [Vector width a]-toRows a =-   let ad = Basic.mapHeight Shape.Deferred $ fromFull a-   in map (takeRow ad) $ Shape.indices $ Matrix.height ad--toColumns ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Class.Floating a) =>-   Full vert horiz height width a -> [Vector height a]-toColumns a =-   let ad = Basic.mapWidth Shape.Deferred $ fromFull a-   in map (takeColumn ad) $ Shape.indices $ Matrix.width ad--toRowArray ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Class.Floating a) =>-   Full vert horiz height width a -> BoxedArray.Array height (Vector width a)-toRowArray a = BoxedArray.fromList (Matrix.height a) (toRows a)--toColumnArray ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Class.Floating a) =>-   Full vert horiz height width a -> BoxedArray.Array width (Vector height a)-toColumnArray a = BoxedArray.fromList (Matrix.width a) (toColumns a)---takeRow ::-   (Extent.C vert, Extent.C horiz,-    Shape.Indexed height, Shape.C width, Shape.Index height ~ ix,-    Class.Floating a) =>-   Full vert horiz height width a -> ix -> Vector width a-takeRow (Array (MatrixShape.Full order extent) x) ix =-   let (height,width) = Extent.dimensions extent-   in case order of-         RowMajor -> pickConsecutive height width x ix-         ColumnMajor -> pickScattered width height x ix--takeColumn ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.Indexed width, Shape.Index width ~ ix,-    Class.Floating a) =>-   Full vert horiz height width a -> ix -> Vector height a-takeColumn (Array (MatrixShape.Full order extent) x) ix =-   let (height,width) = Extent.dimensions extent-   in case order of-         RowMajor -> pickScattered height width x ix-         ColumnMajor -> pickConsecutive width height x ix--pickConsecutive ::-   (Shape.Indexed height, Shape.C width, Shape.Index height ~ ix,-    Class.Floating a) =>-   height -> width -> ForeignPtr a -> ix -> Vector width a-pickConsecutive height width x ix =-   Array.unsafeCreateWithSize width $ \n yPtr -> evalContT $ do-      let offset = Shape.offset height ix-      nPtr <- Call.cint n-      xPtr <- ContT $ withForeignPtr x-      incxPtr <- Call.cint 1-      incyPtr <- Call.cint 1-      liftIO $-         BlasGen.copy nPtr (advancePtr xPtr (n*offset)) incxPtr yPtr incyPtr--pickScattered ::-   (Shape.C height, Shape.Indexed width, Shape.Index width ~ ix,-    Class.Floating a) =>-   height -> width -> ForeignPtr a -> ix -> Vector height a-pickScattered height width x ix =-   Array.unsafeCreateWithSize height $ \n yPtr -> evalContT $ do-      let offset = Shape.offset width ix-      nPtr <- Call.cint n-      xPtr <- ContT $ withForeignPtr x-      incxPtr <- Call.cint $ Shape.size width-      incyPtr <- Call.cint 1-      liftIO $-         BlasGen.copy nPtr (advancePtr xPtr offset) incxPtr yPtr incyPtr---takeTopRows ::-   (Extent.C vert, Shape.C height0, Shape.C height1, Shape.C width,-    Class.Floating a) =>-   Full vert Extent.Big (height0:+:height1) width a ->-   Full vert Extent.Big height0 width a-takeTopRows (Array (MatrixShape.Full order extentA) a) =-   let (heightA@(heightB:+:_), width) = Extent.dimensions extentA-       extentB = Extent.reduceWideHeight heightB extentA-       ma = Shape.size heightA-       mb = Shape.size heightB-       n = Shape.size width-   in Array.unsafeCreateWithSize (MatrixShape.Full order extentB) $-            \blockSize bPtr ->-      withForeignPtr a $ \aPtr ->-      case order of-         RowMajor -> copyBlock blockSize aPtr bPtr-         ColumnMajor -> copySubMatrix mb n ma aPtr mb bPtr--takeBottomRows ::-   (Extent.C vert, Shape.C height0, Shape.C height1, Shape.C width,-    Class.Floating a) =>-   Full vert Extent.Big (height0:+:height1) width a ->-   Full vert Extent.Big height1 width a-takeBottomRows (Array (MatrixShape.Full order extentA) a) =-   let (heightA@(height0:+:heightB), width) = Extent.dimensions extentA-       extentB = Extent.reduceWideHeight heightB extentA-       k = Shape.size height0-       ma = Shape.size heightA-       mb = Shape.size heightB-       n = Shape.size width-   in Array.unsafeCreateWithSize (MatrixShape.Full order extentB) $-            \blockSize bPtr ->-      withForeignPtr a $ \aPtr ->-      case order of-         RowMajor -> copyBlock blockSize (advancePtr aPtr (k*n)) bPtr-         ColumnMajor -> copySubMatrix mb n ma (advancePtr aPtr k) mb bPtr--takeLeftColumns ::-   (Extent.C vert, Shape.C height, Shape.C width0, Shape.C width1,-    Class.Floating a) =>-   Full Extent.Big vert height (width0:+:width1) a ->-   Full Extent.Big vert height width0 a-takeLeftColumns = transpose . takeTopRows . transpose--takeRightColumns ::-   (Extent.C vert, Shape.C height, Shape.C width0, Shape.C width1,-    Class.Floating a) =>-   Full Extent.Big vert height (width0:+:width1) a ->-   Full Extent.Big vert height width1 a-takeRightColumns = transpose . takeBottomRows . transpose---splitRows ::-   (Extent.C vert, Shape.C width, Class.Floating a) =>-   Int ->-   Full vert Extent.Big ZeroInt width a ->-   Full vert Extent.Big (ZeroInt:+:ZeroInt) width a-splitRows k =-   Array.mapShape-      (\(MatrixShape.Full order extentA) ->-         let (Shape.ZeroBased heightA) = Extent.height extentA-             heightB = min k heightA-         in if' (k<0) (error "split: negative number") $-            MatrixShape.Full order $-            Extent.reduceWideHeight-               (Shape.ZeroBased heightB :+: Shape.ZeroBased (heightA-heightB))-               extentA)--takeRows, dropRows ::-   (Extent.C vert, Shape.C width, Class.Floating a) =>-   Int ->-   Full vert Extent.Big ZeroInt width a ->-   Full vert Extent.Big ZeroInt width a-takeRows k = takeTopRows . splitRows k-dropRows k = takeBottomRows . splitRows k--takeColumns, dropColumns ::-   (Extent.C horiz, Shape.C height, Class.Floating a) =>-   Int ->-   Full Extent.Big horiz height ZeroInt a ->-   Full Extent.Big horiz height ZeroInt a-takeColumns k = transpose . takeRows k . transpose-dropColumns k = transpose . dropRows k . transpose---{- |-Take a left-top aligned square or as much as possible of it.-The advantange of this function is that it maintains the matrix size relation,-e.g. Square remains Square, Tall remains Tall.--}-takeEqually ::-   (Extent.C vert, Extent.C horiz, Class.Floating a) =>-   Int ->-   Full vert horiz ZeroInt ZeroInt a ->-   Full vert horiz ZeroInt ZeroInt a-takeEqually k (Array (MatrixShape.Full order extentA) a) =-   let (Shape.ZeroBased heightA, Shape.ZeroBased widthA) =-         Extent.dimensions extentA-       heightB = min k heightA-       widthB  = min k widthA-       extentB =-         Extent.reduceConsistent-            (Shape.ZeroBased heightB) (Shape.ZeroBased widthB) extentA-   in if' (k<0) (error "take: negative number") $-      Array.unsafeCreate (MatrixShape.Full order extentB) $ \bPtr ->-      withForeignPtr a $ \aPtr ->-      case order of-         RowMajor -> copySubMatrix widthB heightB widthA aPtr widthB bPtr-         ColumnMajor -> copySubMatrix heightB widthB heightA aPtr heightB bPtr--{- |-Drop the same number of top-most rows and left-most columns.-The advantange of this function is that it maintains the matrix size relation,-e.g. Square remains Square, Tall remains Tall.--}-dropEqually ::-   (Extent.C vert, Extent.C horiz, Class.Floating a) =>-   Int ->-   Full vert horiz ZeroInt ZeroInt a ->-   Full vert horiz ZeroInt ZeroInt a-dropEqually k (Array (MatrixShape.Full order extentA) a) =-   let (Shape.ZeroBased heightA, Shape.ZeroBased widthA) =-         Extent.dimensions extentA-       heightB = heightA - top; top  = min k heightA-       widthB  = widthA - left; left = min k widthA-       extentB =-         Extent.reduceConsistent-            (Shape.ZeroBased heightB) (Shape.ZeroBased widthB) extentA-   in if' (k<0) (error "drop: negative number") $-      Array.unsafeCreate (MatrixShape.Full order extentB) $ \bPtr ->-      withForeignPtr a $ \aPtr ->-      case order of-         RowMajor ->-            copySubMatrix widthB heightB-               widthA (advancePtr aPtr (top*widthA+left)) widthB bPtr-         ColumnMajor ->-            copySubMatrix heightB widthB-               heightA (advancePtr aPtr (left*heightA+top)) heightB bPtr---swapRows ::-   (Extent.C vert, Extent.C horiz,-    Shape.Indexed height, Shape.C width, Class.Floating a) =>-   Shape.Index height -> Shape.Index height ->-   Full vert horiz height width a -> Full vert horiz height width a-swapRows i j (Array shape@(MatrixShape.Full order extent) a) =-   Array.unsafeCreateWithSize shape $ \blockSize bPtr -> evalContT $ do-      let (height,width) = Extent.dimensions extent-      let m = Shape.size height-      let n = Shape.size width-      nPtr <- Call.cint n-      aPtr <- ContT $ withForeignPtr a-      let offsetI = Shape.offset height i-      let offsetJ = Shape.offset height j-      let (incVert,incHoriz) =-            case order of-               RowMajor -> (n,1)-               ColumnMajor -> (1,m)-      incPtr <- Call.cint incHoriz-      liftIO $ do-         copyBlock blockSize aPtr bPtr-         when (offsetI/=offsetJ) $-            BlasGen.swap nPtr-               (advancePtr bPtr (incVert*offsetI)) incPtr-               (advancePtr bPtr (incVert*offsetJ)) incPtr--swapColumns ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.Indexed width, Class.Floating a) =>-   Shape.Index width -> Shape.Index width ->-   Full vert horiz height width a -> Full vert horiz height width a-swapColumns i j = transpose . swapRows i j . transpose----- alternative: laswp-reverseRows ::-   (Extent.C vert, Extent.C horiz, Shape.C width, Class.Floating a) =>-   Full vert horiz ZeroInt width a -> Full vert horiz ZeroInt width a-reverseRows (Array shape@(MatrixShape.Full order extent) a) =-   Array.unsafeCreateWithSize shape $ \blockSize bPtr -> evalContT $ do-      let (height,width) = Extent.dimensions extent-      let n = Shape.size height-      let m = Shape.size width-      fwdPtr <- Call.bool True-      nPtr <- Call.cint n-      mPtr <- Call.cint m-      kPtr <- Call.allocaArray n-      aPtr <- ContT $ withForeignPtr a-      liftIO $ do-         copyBlock blockSize aPtr bPtr-         pokeArray kPtr $ take n $ iterate (subtract 1) $ fromIntegral n-         case order of-            RowMajor -> LapackGen.lapmt fwdPtr mPtr nPtr bPtr mPtr kPtr-            ColumnMajor -> LapackGen.lapmr fwdPtr nPtr mPtr bPtr nPtr kPtr--reverseColumns ::-   (Extent.C vert, Extent.C horiz, Shape.C height, Class.Floating a) =>-   Full vert horiz height ZeroInt a -> Full vert horiz height ZeroInt a-reverseColumns = transpose . reverseRows . transpose---{--The function is optimized for blocks of consecutive rows.-For scattered rows in column major order-the function has quite ugly memory access patterns.--}-takeRowArray ::-   (Shape.Indexed height, Shape.C width, Shape.C sh, Class.Floating a) =>-   BoxedArray.Array sh (Shape.Index height) ->-   General height width a -> General sh width a-takeRowArray ixs (Array (MatrixShape.Full order extent) a) =-   let (heightA, width) = Extent.dimensions extent-       heightB = BoxedArray.shape ixs-       offsets = map (Shape.offset heightA) $ BoxedArray.toList ixs-       startBlocks blocks = zip (scanl (+) 0 $ map fst blocks) blocks-       ma = Shape.size heightA-       mb = Shape.size heightB-       n = Shape.size width-   in Array.unsafeCreate (MatrixShape.general order heightB width) $ \bPtr ->-      withForeignPtr a $ \aPtr ->-      case order of-         RowMajor -> do-            forM_ (startBlocks $ chopRowBlocks offsets) $-               \(dest, (numRows, (start,step))) ->-                  copySubMatrix n numRows-                     (step*n) (advancePtr aPtr (start*n))-                     n (advancePtr bPtr (dest*n))-         ColumnMajor -> do-            forM_ (startBlocks $ chopColumnBlocks offsets) $-               \(dest, (numRows, start)) ->-                  copySubMatrix numRows n-                     ma (advancePtr aPtr start)-                     mb (advancePtr bPtr dest)--chopRowBlocks :: (Integral i) => [i] -> [(Int,(i,i))]-chopRowBlocks =-   let go [] = []-       go is@(i0:is0) =-         case mfilter (i0<) $ listToMaybe is0 of-            Nothing -> (1,(i0,0)) : go is0-            Just i1 ->-               let (consecutive,remainder) =-                     span (uncurry (==)) $ zip [i0,i1..] is-               in (length consecutive, (i0,i1-i0)) : go (map snd remainder)-   in go--chopColumnBlocks :: (Integral i) => [i] -> [(Int,i)]-chopColumnBlocks =-   map (\is -> (length $ NonEmpty.flatten is, NonEmpty.head is)) .-   NonEmptyM.groupBy (\i j -> i+1 == j)---takeColumnArray ::-   (Shape.C height, Shape.Indexed width, Shape.C sh, Class.Floating a) =>-   BoxedArray.Array sh (Shape.Index width) ->-   General height width a -> General height sh a-takeColumnArray ixs = transpose . takeRowArray ixs . transpose----fromRowMajor ::-   (Shape.C height, Shape.C width, Class.Floating a) =>-   Array (height,width) a -> General height width a-fromRowMajor = Array.mapShape (uncurry $ MatrixShape.general RowMajor)--toRowMajor ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Class.Floating a) =>-   Full vert horiz height width a -> Array (height,width) a-toRowMajor =-   Array.mapShape-      (\shape -> (MatrixShape.fullHeight shape, MatrixShape.fullWidth shape)) .-   forceRowMajor--{- |-@adaptOrder x y@ contains the data of @y@ with the layout of @x@.--}-adaptOrder ::-   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,-    Class.Floating a) =>-   Full vert horiz height width a ->-   Full vert horiz height width a ->-   Full vert horiz height width a-adaptOrder x = forceOrder (MatrixShape.fullOrder $ Array.shape x)--flatten ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Class.Floating a) =>-   Full vert horiz height width a -> Vector ZeroInt a-flatten = Array.mapShape (zeroInt . Shape.size) . toRowMajor---infixl 3 |||-infixl 2 ===--(|||) ::-   (Extent.C vert, Shape.C height, Eq height, Shape.C widtha, Shape.C widthb,-    Class.Floating a) =>-   Full vert Extent.Big height widtha a ->-   Full vert Extent.Big height widthb a ->-   Full vert Extent.Big height (widtha:+:widthb) a-(|||)-      (Array (MatrixShape.Full orderA extentA) a)-      (Array (MatrixShape.Full orderB extentB) b) =-   let (heightA,widthA) = Extent.dimensions extentA-       (heightB,widthB) = Extent.dimensions extentB-       extent = Extent.widen (widthA:+:widthB) extentA-       shape order = MatrixShape.Full order extent-   in-    if heightA /= heightB-      then error "(|||): mismatching heights"-      else-         case (orderA,orderB) of-            (RowMajor,RowMajor) ->-               Array.unsafeCreate (shape RowMajor) $-               \cPtr -> evalContT $ do-                  let n = Shape.size heightA-                  let ma = Shape.size widthA-                  let mb = Shape.size widthB-                  let m = ma+mb-                  maPtr <- Call.cint ma-                  mbPtr <- Call.cint mb-                  aPtr <- ContT $ withForeignPtr a-                  bPtr <- ContT $ withForeignPtr b-                  incxPtr <- Call.cint 1-                  incyPtr <- Call.cint 1-                  liftIO $-                     sequence_ $ take n $-                     zipWith3-                        (\akPtr bkPtr ckPtr -> do-                           BlasGen.copy maPtr akPtr incxPtr ckPtr incyPtr-                           BlasGen.copy mbPtr bkPtr incxPtr-                              (ckPtr `advancePtr` ma) incyPtr)-                        (pointerSeq ma aPtr)-                        (pointerSeq mb bPtr)-                        (pointerSeq m cPtr)-            (RowMajor,ColumnMajor) ->-               Array.unsafeCreate (shape ColumnMajor) $-               \cPtr -> evalContT $ do-                  let n = Shape.size heightA-                  let ma = Shape.size widthA-                  let mb = Shape.size widthB-                  aPtr <- ContT $ withForeignPtr a-                  bPtr <- ContT $ withForeignPtr b-                  liftIO $ do-                     copyTransposed n ma aPtr n cPtr-                     copyBlock (n*mb) bPtr (advancePtr cPtr (n*ma))-            (ColumnMajor,RowMajor) ->-               Array.unsafeCreate (shape ColumnMajor) $-               \cPtr -> evalContT $ do-                  let n = Shape.size heightA-                  let ma = Shape.size widthA-                  let mb = Shape.size widthB-                  let volA = n*ma-                  aPtr <- ContT $ withForeignPtr a-                  bPtr <- ContT $ withForeignPtr b-                  liftIO $ do-                     copyBlock volA aPtr cPtr-                     copyTransposed n mb bPtr n (advancePtr cPtr volA)-            (ColumnMajor,ColumnMajor) ->-               Array.unsafeCreate (shape ColumnMajor) $-               \cPtr -> evalContT $ do-                  let n = Shape.size heightA-                  let na = n * Shape.size widthA-                  let nb = n * Shape.size widthB-                  naPtr <- Call.cint na-                  nbPtr <- Call.cint nb-                  aPtr <- ContT $ withForeignPtr a-                  bPtr <- ContT $ withForeignPtr b-                  incxPtr <- Call.cint 1-                  incyPtr <- Call.cint 1-                  liftIO $ do-                     BlasGen.copy naPtr aPtr incxPtr cPtr incyPtr-                     BlasGen.copy nbPtr bPtr incxPtr-                        (cPtr `advancePtr` na) incyPtr--(===) ::-   (Extent.C horiz, Shape.C width, Eq width, Shape.C heighta, Shape.C heightb,-    Class.Floating a) =>-   Full Extent.Big horiz heighta width a ->-   Full Extent.Big horiz heightb width a ->-   Full Extent.Big horiz (heighta:+:heightb) width a-(===) a b = transpose (transpose a ||| transpose b)---add, sub ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Eq height, Eq width,-    Class.Floating a) =>-   Full vert horiz height width a ->-   Full vert horiz height width a ->-   Full vert horiz height width a-add x y = Vector.add (adaptOrder y x) y-sub x y = Vector.sub (adaptOrder y x) y---rowSums ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Class.Floating a) =>-   Full vert horiz height width a -> Vector height a-rowSums m =-   let width = MatrixShape.fullWidth $ Array.shape m-   in  multiplyVectorUnchecked m (Vector.constant width one)--columnSums ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Class.Floating a) =>-   Full vert horiz height width a -> Vector width a-columnSums m =-   let height = MatrixShape.fullHeight $ Array.shape m-   in  multiplyVectorUnchecked (transpose m) (Vector.constant height one)---scaleRowsComplex ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Eq height, Shape.C width, Class.Real a) =>-   Vector height a ->-   Full vert horiz height width (Complex a) ->-   Full vert horiz height width (Complex a)-scaleRowsComplex-   (Array heightX x) (Array shape@(MatrixShape.Full order extent) a) =-      Array.unsafeCreate shape $ \bComplexPtr -> do-   let (height,width) = Extent.dimensions extent-   Call.assert "scaleRowsComplex: sizes mismatch" (heightX == height)-   let bPtr = castPtr bComplexPtr-   case order of-      RowMajor -> evalContT $ do-         let m = Shape.size height-         let n = Shape.size width * 2-         alphaPtr <- Call.alloca-         nPtr <- Call.cint n-         xPtr <- ContT $ withForeignPtr x-         aPtr <- fmap castPtr $ ContT $ withForeignPtr a-         incaPtr <- Call.cint 1-         incbPtr <- Call.cint 1-         liftIO $ sequence_ $ take m $-            zipWith3-               (\xkPtr akPtr bkPtr -> do-                  poke alphaPtr =<< peek xkPtr-                  BlasGen.copy nPtr akPtr incaPtr bkPtr incbPtr-                  BlasGen.scal nPtr alphaPtr bkPtr incbPtr)-               (pointerSeq 1 xPtr)-               (pointerSeq n aPtr)-               (pointerSeq n bPtr)-      ColumnMajor -> evalContT $ do-         let m = Shape.size width-         let nr = Shape.size height-         let n = 2*nr-         transPtr <- Call.char 'N'-         nrPtr <- Call.cint nr-         nPtr <- Call.cint n-         klPtr <- Call.cint 0-         kuPtr <- Call.cint 0-         alphaPtr <- Call.number one-         xrPtr <- ContT $ withForeignPtr x-         xPtr <- Call.allocaArray n-         incxrPtr <- Call.cint 1-         incxPtr <- Call.cint 2-         ldxPtr <- Call.leadingDim 1-         aPtr <- fmap castPtr $ ContT $ withForeignPtr a-         incaPtr <- Call.cint 1-         betaPtr <- Call.number zero-         incbPtr <- Call.cint 1-         liftIO $ do-            BlasGen.copy nrPtr xrPtr incxrPtr xPtr incxPtr-            BlasGen.copy nrPtr xrPtr incxrPtr (advancePtr xPtr 1) incxPtr-            sequence_ $ take m $-               zipWith-                  (\akPtr bkPtr ->-                     Private.gbmv transPtr-                        nPtr nPtr klPtr kuPtr alphaPtr xPtr ldxPtr-                        akPtr incaPtr betaPtr bkPtr incbPtr)-                  (pointerSeq n aPtr)-                  (pointerSeq n bPtr)--scaleColumnsComplex ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Eq width, Class.Real a) =>-   Vector width a ->-   Full vert horiz height width (Complex a) ->-   Full vert horiz height width (Complex a)-scaleColumnsComplex x = transpose . scaleRowsComplex x . transpose---scaleRowsReal ::-   (Extent.C vert, Extent.C horiz, Shape.C height, Eq height, Shape.C width,-    Class.Floating a) =>-   Vector height (RealOf a) ->-   Full vert horiz height width a ->-   Full vert horiz height width a-scaleRowsReal =-   getScaleRowsReal $-   Class.switchFloating-      (ScaleRowsReal scaleRows)-      (ScaleRowsReal scaleRows)-      (ScaleRowsReal scaleRowsComplex)-      (ScaleRowsReal scaleRowsComplex)--newtype ScaleRowsReal f g a =-   ScaleRowsReal {getScaleRowsReal :: f (RealOf a) -> g a -> g a}--scaleColumnsReal ::-   (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.C width, Eq width, Class.Floating a) =>-   Vector width (RealOf a) ->-   Full vert horiz height width a ->-   Full vert horiz height width a-scaleColumnsReal x = transpose . scaleRowsReal x . transpose---{- |-> tensorProduct order x y = singleColumn order x <#> singleRow order y--}-tensorProduct ::-   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>-   Order -> Vector height a -> Vector width a -> General height width a-tensorProduct order x y =-   case order of-      ColumnMajor -> tensorProd 'T' order x y-      RowMajor -> transpose $ tensorProd 'T' order y x--{- |-> outer order x y = tensorProduct order x (Vector.conjugate y)--}-outer ::-   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>-   Order -> Vector height a -> Vector width a -> General height width a-outer order x y =-   case order of-      ColumnMajor -> tensorProd 'C' ColumnMajor x y-      RowMajor -> transpose $ tensorProd 'C' RowMajor y x--{-# INLINE tensorProd #-}-tensorProd ::-   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>-   Char -> Order ->-   Vector height a -> Vector width a -> General height width a-tensorProd trans order (Array shX x) (Array shY y) =-   Array.unsafeCreate (MatrixShape.general MatrixShape.ColumnMajor shX shY) $-      \cPtr -> do-   let m = Shape.size shX-   let n = Shape.size shY-   let ((transa,transb),(lda,ldb)) =-         case order of-            ColumnMajor -> (('N',trans),(m,n))-            RowMajor -> ((trans,'N'),(1,1))-   evalContT $ do-      transaPtr <- Call.char transa-      transbPtr <- Call.char transb-      mPtr <- Call.cint m-      nPtr <- Call.cint n-      kPtr <- Call.cint 1-      alphaPtr <- Call.number one-      aPtr <- ContT $ withForeignPtr x-      ldaPtr <- Call.leadingDim lda-      bPtr <- ContT $ withForeignPtr y-      ldbPtr <- Call.leadingDim ldb-      betaPtr <- Call.number zero-      ldcPtr <- Call.leadingDim m-      liftIO $-         BlasGen.gemm-            transaPtr transbPtr mPtr nPtr kPtr alphaPtr-            aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr---sumRank1 ::-   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>-   (height,width) ->-   [(a, (Vector height a, Vector width a))] -> General height width a-sumRank1 (height,width) xys =-   Array.unsafeCreateWithSize (MatrixShape.general ColumnMajor height width) $-      \size aPtr ->-   evalContT $ do-      let m = Shape.size height-      let n = Shape.size width-      mPtr <- Call.cint m-      nPtr <- Call.cint n-      alphaPtr <- Call.alloca-      incxPtr <- Call.cint 1-      incyPtr <- Call.cint 1-      ldaPtr <- Call.leadingDim m-      liftIO $ do-         fill zero size aPtr-         forM_ xys $ \(alpha, (Array shX x, Array shY y)) ->-            withForeignPtr x $ \xPtr ->-            withForeignPtr y $ \yPtr -> do-               Call.assert "Matrix.sumRank1: non-matching height" (height==shX)-               Call.assert "Matrix.sumRank1: non-matching width" (width==shY)-               poke alphaPtr alpha-               BlasGen.gerc mPtr nPtr-                  alphaPtr xPtr incxPtr yPtr incyPtr aPtr ldaPtr+   Type.Matrix,+   Full,+   General, Tall, Wide, Square.Square,+   Triangular.Triangular, Triangular.Upper, Triangular.Lower,+   Triangular.Diagonal, Triangular.Symmetric,+   Hermitian.Hermitian,+   Permutation,++   ShapeInt, shapeInt,+   transpose, adjoint,+   Type.height, Type.width,+   Type.HeightOf, Type.WidthOf,+   Type.Box, Type.indices,+   ArrMatrix.reshape,+   ArrMatrix.mapShape,+   caseTallWide,+   fromScalar, toScalar,+   fromList,+   mapExtent, fromFull,+   asGeneral, asTall, asWide,+   tallFromGeneral, wideFromGeneral,+   generalizeTall, generalizeWide,+   mapHeight, mapWidth,+   identity,+   diagonal,+   fromRowsNonEmpty,    fromRowArray,    fromRows,+   fromRowsNonEmptyContainer,    fromRowContainer,+   fromColumnsNonEmpty, fromColumnArray, fromColumns,+   fromColumnsNonEmptyContainer, fromColumnContainer,+   singleRow,   singleColumn,+   flattenRow,  flattenColumn,+   liftRow,     liftColumn,+   unliftRow,   unliftColumn,+   toRows, toColumns,+   toRowArray, toColumnArray,+   toRowContainer, toColumnContainer,+   takeRow, takeColumn,+   takeRows, takeColumns, takeEqually,+   dropRows, dropColumns, dropEqually,+   takeTop, takeBottom,+   takeLeft, takeRight,+   takeRowArray, takeColumnArray,+   swapRows, swapColumns,+   reverseRows, reverseColumns,+   fromRowMajor, toRowMajor,+   ArrMatrix.forceOrder, ArrMatrix.adaptOrder,+   Basic.OrderBias, leftBias, rightBias, contiguousBias,+   (|||), beside,+   (===), above,++   (|*-),+   tensorProduct,+   outer,+   kronecker,+   sumRank1,++   map,+   MatrixClass.Complex, MatrixClass.conjugate,+   MatrixClass.fromReal, MatrixClass.toComplex,+   MatrixClass.SquareShape, MatrixClass.toSquare,+   MatrixClass.identityFrom,+   MatrixClass.identityFromHeight, MatrixClass.identityFromWidth,+   MatrixClass.takeDiagonal, MatrixClass.trace,++   RealOf,+   rowSums, columnSums,+   rowArgAbsMaximums, columnArgAbsMaximums,+   scaleRows, scaleColumns,+   scaleRowsReal, scaleColumnsReal,+   (\*#), (#*\),+   (\\#), (#/\),+   multiply,+   multiplyVector,++   ArrMatrix.zero, ArrMatrix.negate,+   ArrMatrix.scale, ArrMatrix.scaleReal, ArrMatrix.scaleRealReal,+   (ArrMatrix..*#),+   ArrMatrix.add, ArrMatrix.sub,+   (ArrMatrix.#+#), (ArrMatrix.#-#),+   Multiply.Multiply, (Multiply.#*#),+   Multiply.MultiplyVector, (Multiply.#*|), (Multiply.-*#),+   Multiply.MultiplySquare, multiplySquare,+   Multiply.Power, Multiply.square, Multiply.power,+   (Multiply.##*#), (Multiply.#*##),+   Indexed.Indexed, (Indexed.#!),++   Divide.Determinant, Divide.determinant,+   Divide.Solve, Divide.solve, Divide.solveLeft, Divide.solveRight,+   (Divide.##/#), (Divide.#\##),+   Divide.solveVector, (Divide.-/#), (Divide.#\|),+   Divide.Inverse, Divide.inverse,+   Mod.Transposition(..),+   ) where++import qualified Numeric.LAPACK.Matrix.Permutation as PermMatrix+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular+import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian+import qualified Numeric.LAPACK.Matrix.Square as Square++import qualified Numeric.LAPACK.Matrix.Extent as Extent+import qualified Numeric.LAPACK.Matrix.Basic as Basic+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix.Type as Type+import qualified Numeric.LAPACK.Matrix.Plain as Plain+import qualified Numeric.LAPACK.Matrix.Modifier as Mod+import qualified Numeric.LAPACK.Matrix.Divide as Divide+import qualified Numeric.LAPACK.Matrix.Multiply as Multiply+import qualified Numeric.LAPACK.Matrix.Indexed as Indexed+import qualified Numeric.LAPACK.Matrix.Class as MatrixClass+import qualified Numeric.LAPACK.Matrix.Private as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix.Shape.Private (Order)+import Numeric.LAPACK.Matrix.Array (Full, General, Tall, Wide)+import Numeric.LAPACK.Matrix.Private (ShapeInt, shapeInt)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Scalar (RealOf)++import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Boxed as BoxedArray+import qualified Data.Array.Comfort.Container as Container+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Unchecked (Array, (!))+import Data.Array.Comfort.Shape ((:+:))++import Foreign.Storable (Storable)++import qualified Data.NonEmpty as NonEmpty+import qualified Data.Either.HT as EitherHT++import Prelude hiding (map)+++type Permutation sh = Type.Matrix (PermMatrix.Permutation sh)++mapExtent ::+   (Extent.C vertA, Extent.C horizA) =>+   (Extent.C vertB, Extent.C horizB) =>+   Extent.Map vertA horizA vertB horizB height width ->+   Full vertA horizA height width a -> Full vertB horizB height width a+mapExtent = ArrMatrix.lift1 . Plain.mapExtent++fromFull ::+   (Extent.C vert, Extent.C horiz) =>+   Full vert horiz height width a -> General height width a+fromFull = ArrMatrix.lift1 Plain.fromFull++tallFromGeneral ::+   (Shape.C height, Shape.C width, Storable a) =>+   General height width a -> Tall height width a+tallFromGeneral = ArrMatrix.lift1 Plain.tallFromGeneral++wideFromGeneral ::+   (Shape.C height, Shape.C width, Storable a) =>+   General height width a -> Wide height width a+wideFromGeneral = ArrMatrix.lift1 Plain.wideFromGeneral++generalizeTall ::+   (Extent.C vert, Extent.C horiz) =>+   Full vert Extent.Small height width a -> Full vert horiz height width a+generalizeTall = mapExtent Extent.generalizeTall++generalizeWide ::+   (Extent.C vert, Extent.C horiz) =>+   Full Extent.Small horiz height width a -> Full vert horiz height width a+generalizeWide = mapExtent Extent.generalizeWide+++asGeneral :: General height width a -> General height width a+asGeneral = id++asTall :: Tall height width a -> Tall height width a+asTall = id++asWide :: Wide height width a -> Wide height width a+asWide = id++++fromScalar :: (Storable a) => a -> General () () a+fromScalar = fromFull . Square.fromScalar++toScalar :: (Storable a) => General () () a -> a+toScalar a = either id id (Matrix.revealOrder (ArrMatrix.toVector a)) ! ((),())++fromList ::+   (Shape.C height, Shape.C width, Storable a) =>+   height -> width -> [a] -> General height width a+fromList height width = ArrMatrix.lift0 . Plain.fromList height width+++identity ::+   (Shape.C sh, Class.Floating a) =>+   sh -> General sh sh a+identity = ArrMatrix.lift0 . Plain.identity++diagonal ::+   (Shape.C sh, Class.Floating a) =>+   Vector sh a -> General sh sh a+diagonal = ArrMatrix.lift0 . Plain.diagonal+++{- |+Square matrices will be classified as 'Tall'.+-}+caseTallWide ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>+   Full vert horiz height width a ->+   Either (Tall height width a) (Wide height width a)+caseTallWide =+   EitherHT.mapBoth ArrMatrix.lift0 ArrMatrix.lift0 .+   Basic.caseTallWide . ArrMatrix.toVector+++transpose ::+   (Extent.C vert, Extent.C horiz) =>+   Full vert horiz height width a -> Full horiz vert width height a+transpose = ArrMatrix.lift1 Basic.transpose++{- |+conjugate transpose++Problem: @adjoint a \<\> a@ is always square,+but how to convince the type checker to choose the Square type?+Anser: Use @Hermitian.toSquare $ Hermitian.gramian a@ instead.+-}+adjoint ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,+    Class.Floating a) =>+   Full vert horiz height width a -> Full horiz vert width height a+adjoint = ArrMatrix.lift1 Basic.adjoint+++{- |+The number of rows must be maintained by the height mapping function.+-}+mapHeight ::+   (Shape.C heightA, Shape.C heightB,+    Extent.GeneralTallWide vert horiz,+    Extent.GeneralTallWide horiz vert) =>+   (heightA -> heightB) ->+   Full vert horiz heightA width a ->+   Full vert horiz heightB width a+mapHeight = ArrMatrix.lift1 . Plain.mapHeight++{- |+The number of columns must be maintained by the width mapping function.+-}+mapWidth ::+   (Shape.C widthA, Shape.C widthB,+    Extent.GeneralTallWide vert horiz,+    Extent.GeneralTallWide horiz vert) =>+   (widthA -> widthB) ->+   Full vert horiz height widthA a ->+   Full vert horiz height widthB a+mapWidth = ArrMatrix.lift1 . Plain.mapWidth+++singleRow :: Order -> Vector width a -> General () width a+singleRow order = ArrMatrix.lift0 . Basic.singleRow order++singleColumn :: Order -> Vector height a -> General height () a+singleColumn order = ArrMatrix.lift0 . Basic.singleColumn order++flattenRow :: General () width a -> Vector width a+flattenRow = Basic.flattenRow . ArrMatrix.toVector++flattenColumn :: General height () a -> Vector height a+flattenColumn = Basic.flattenColumn . ArrMatrix.toVector++liftRow ::+   Order ->+   (Vector height0 a -> Vector height1 b) ->+   General () height0 a -> General () height1 b+liftRow order = ArrMatrix.lift1 . Basic.liftRow order++liftColumn ::+   Order ->+   (Vector height0 a -> Vector height1 b) ->+   General height0 () a -> General height1 () b+liftColumn order = ArrMatrix.lift1 . Basic.liftColumn order++unliftRow ::+   Order ->+   (General () height0 a -> General () height1 b) ->+   Vector height0 a -> Vector height1 b+unliftRow order = Basic.unliftRow order .  ArrMatrix.unlift1++unliftColumn ::+   Order ->+   (General height0 () a -> General height1 () b) ->+   Vector height0 a -> Vector height1 b+unliftColumn order = Basic.unliftColumn order . ArrMatrix.unlift1+++fromRowsNonEmpty ::+   (Shape.C width, Eq width, Storable a) =>+   NonEmpty.T [] (Vector width a) -> General ShapeInt width a+fromRowsNonEmpty = ArrMatrix.lift0 . Plain.fromRowsNonEmpty++fromRowArray ::+   (Shape.C height, Shape.C width, Eq width, Storable a) =>+   width -> BoxedArray.Array height (Vector width a) -> General height width a+fromRowArray width = ArrMatrix.lift0 . Plain.fromRowArray width++-- ToDo: generalize to a new NonEmpty.Head class+fromRowsNonEmptyContainer ::+   (f ~ NonEmpty.T g, Container.C g,+    Shape.C width, Eq width, Storable a) =>+   f (Vector width a) -> General (Container.Shape f) width a+fromRowsNonEmptyContainer = ArrMatrix.lift0 . Plain.fromRowsNonEmptyContainer++fromRowContainer ::+   (Container.C f, Shape.C width, Eq width, Storable a) =>+   width -> f (Vector width a) -> General (Container.Shape f) width a+fromRowContainer width = ArrMatrix.lift0 . Plain.fromRowContainer width++fromRows ::+   (Shape.C width, Eq width, Storable a) =>+   width -> [Vector width a] -> General ShapeInt width a+fromRows width = ArrMatrix.lift0 . Plain.fromRows width+++fromColumnsNonEmpty ::+   (Shape.C height, Eq height, Storable a) =>+   NonEmpty.T [] (Vector height a) -> General height ShapeInt a+fromColumnsNonEmpty = ArrMatrix.lift0 . Plain.fromColumnsNonEmpty++fromColumnArray ::+   (Shape.C height, Eq height, Shape.C width, Storable a) =>+   height -> BoxedArray.Array width (Vector height a) -> General height width a+fromColumnArray height = ArrMatrix.lift0 . Plain.fromColumnArray height++fromColumnsNonEmptyContainer ::+   (f ~ NonEmpty.T g, Container.C g,+    Shape.C height, Eq height, Storable a) =>+   f (Vector height a) -> General height (Container.Shape f) a+fromColumnsNonEmptyContainer =+   ArrMatrix.lift0 . Plain.fromColumnsNonEmptyContainer++fromColumnContainer ::+   (Container.C f, Shape.C height, Eq height, Storable a) =>+   height -> f (Vector height a) -> General height (Container.Shape f) a+fromColumnContainer height = ArrMatrix.lift0 . Plain.fromColumnContainer height++fromColumns ::+   (Shape.C height, Eq height, Storable a) =>+   height -> [Vector height a] -> General height ShapeInt a+fromColumns height = ArrMatrix.lift0 . Plain.fromColumns height+++toRows ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Class.Floating a) =>+   Full vert horiz height width a -> [Vector width a]+toRows = Plain.toRows . ArrMatrix.toVector++toColumns ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Class.Floating a) =>+   Full vert horiz height width a -> [Vector height a]+toColumns = Plain.toColumns . ArrMatrix.toVector++toRowArray ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Class.Floating a) =>+   Full vert horiz height width a -> BoxedArray.Array height (Vector width a)+toRowArray = Plain.toRowArray . ArrMatrix.toVector++toColumnArray ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Class.Floating a) =>+   Full vert horiz height width a -> BoxedArray.Array width (Vector height a)+toColumnArray = Plain.toColumnArray . ArrMatrix.toVector++toRowContainer ::+   (Extent.C vert, Extent.C horiz,+    Container.C f, Shape.C width, Class.Floating a) =>+   Full vert horiz (Container.Shape f) width a -> f (Vector width a)+toRowContainer = Plain.toRowContainer . ArrMatrix.toVector++toColumnContainer ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Container.C f, Class.Floating a) =>+   Full vert horiz height (Container.Shape f) a -> f (Vector height a)+toColumnContainer = Plain.toColumnContainer . ArrMatrix.toVector++++{-+The parameter order is swapped with respect to 'takeRowArray'+but it is the order that is used most oftenly.+-}+takeRow ::+   (Extent.C vert, Extent.C horiz,+    Shape.Indexed height, Shape.C width, Shape.Index height ~ ix,+    Class.Floating a) =>+   Full vert horiz height width a -> ix -> Vector width a+takeRow = Plain.takeRow . ArrMatrix.toVector++takeColumn ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.Indexed width, Shape.Index width ~ ix,+    Class.Floating a) =>+   Full vert horiz height width a -> ix -> Vector height a+takeColumn = Plain.takeColumn . ArrMatrix.toVector+++takeTop ::+   (Extent.C vert, Shape.C height0, Shape.C height1, Shape.C width,+    Class.Floating a) =>+   Full vert Extent.Big (height0:+:height1) width a ->+   Full vert Extent.Big height0 width a+takeTop = ArrMatrix.lift1 Basic.takeTop++takeBottom ::+   (Extent.C vert, Shape.C height0, Shape.C height1, Shape.C width,+    Class.Floating a) =>+   Full vert Extent.Big (height0:+:height1) width a ->+   Full vert Extent.Big height1 width a+takeBottom = ArrMatrix.lift1 Basic.takeBottom++takeLeft ::+   (Extent.C vert, Shape.C height, Shape.C width0, Shape.C width1,+    Class.Floating a) =>+   Full Extent.Big vert height (width0:+:width1) a ->+   Full Extent.Big vert height width0 a+takeLeft = ArrMatrix.lift1 Basic.takeLeft++takeRight ::+   (Extent.C vert, Shape.C height, Shape.C width0, Shape.C width1,+    Class.Floating a) =>+   Full Extent.Big vert height (width0:+:width1) a ->+   Full Extent.Big vert height width1 a+takeRight = ArrMatrix.lift1 Basic.takeRight++takeRows, dropRows ::+   (Extent.C vert, Shape.C width, Class.Floating a) =>+   Int ->+   Full vert Extent.Big ShapeInt width a ->+   Full vert Extent.Big ShapeInt width a+takeRows = ArrMatrix.lift1 . Basic.takeRows+dropRows = ArrMatrix.lift1 . Basic.dropRows++takeColumns, dropColumns ::+   (Extent.C horiz, Shape.C height, Class.Floating a) =>+   Int ->+   Full Extent.Big horiz height ShapeInt a ->+   Full Extent.Big horiz height ShapeInt a+takeColumns = ArrMatrix.lift1 . Basic.takeColumns+dropColumns = ArrMatrix.lift1 . Basic.dropColumns+++{- |+Take a left-top aligned square or as much as possible of it.+The advantange of this function is that it maintains the matrix size relation,+e.g. Square remains Square, Tall remains Tall.+-}+takeEqually ::+   (Extent.C vert, Extent.C horiz, Class.Floating a) =>+   Int ->+   Full vert horiz ShapeInt ShapeInt a ->+   Full vert horiz ShapeInt ShapeInt a+takeEqually = ArrMatrix.lift1 . Plain.takeEqually++{- |+Drop the same number of top-most rows and left-most columns.+The advantange of this function is that it maintains the matrix size relation,+e.g. Square remains Square, Tall remains Tall.+-}+dropEqually ::+   (Extent.C vert, Extent.C horiz, Class.Floating a) =>+   Int ->+   Full vert horiz ShapeInt ShapeInt a ->+   Full vert horiz ShapeInt ShapeInt a+dropEqually = ArrMatrix.lift1 . Plain.dropEqually+++swapRows ::+   (Extent.C vert, Extent.C horiz,+    Shape.Indexed height, Shape.C width, Class.Floating a) =>+   Shape.Index height -> Shape.Index height ->+   Full vert horiz height width a -> Full vert horiz height width a+swapRows i j = ArrMatrix.lift1 $ Plain.swapRows i j++swapColumns ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.Indexed width, Class.Floating a) =>+   Shape.Index width -> Shape.Index width ->+   Full vert horiz height width a -> Full vert horiz height width a+swapColumns i j =  ArrMatrix.lift1 $ Plain.swapColumns i j+++reverseRows ::+   (Extent.C vert, Extent.C horiz, Shape.C width, Class.Floating a) =>+   Full vert horiz ShapeInt width a -> Full vert horiz ShapeInt width a+reverseRows = ArrMatrix.lift1 Plain.reverseRows++reverseColumns ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Class.Floating a) =>+   Full vert horiz height ShapeInt a -> Full vert horiz height ShapeInt a+reverseColumns = ArrMatrix.lift1 Plain.reverseColumns+++{- |+The function is optimized for blocks of consecutive rows.+For scattered rows in column major order+the function has quite ugly memory access patterns.+-}+takeRowArray ::+   (Shape.Indexed height, Shape.C width, Shape.C sh, Class.Floating a) =>+   BoxedArray.Array sh (Shape.Index height) ->+   General height width a -> General sh width a+takeRowArray = ArrMatrix.lift1 . Plain.takeRowArray++takeColumnArray ::+   (Shape.C height, Shape.Indexed width, Shape.C sh, Class.Floating a) =>+   BoxedArray.Array sh (Shape.Index width) ->+   General height width a -> General height sh a+takeColumnArray = ArrMatrix.lift1 . Plain.takeColumnArray++++fromRowMajor ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   Array (height,width) a -> General height width a+fromRowMajor = ArrMatrix.lift0 . Plain.fromRowMajor++toRowMajor ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Class.Floating a) =>+   Full vert horiz height width a -> Array (height,width) a+toRowMajor = Plain.toRowMajor . ArrMatrix.toVector+++infixr 3 |||+infixr 2 ===++(|||) ::+   (Extent.C vertA, Extent.C vertB, Extent.C vertC,+    Extent.Append vertA vertB ~ vertC,+    Shape.C height, Eq height, Shape.C widthA, Shape.C widthB,+    Class.Floating a) =>+   Full vertA Extent.Big height widthA a ->+   Full vertB Extent.Big height widthB a ->+   Full vertC Extent.Big height (widthA:+:widthB) a+(|||) = beside rightBias Extent.appendAny++(===) ::+   (Extent.C horizA, Extent.C horizB, Extent.C horizC,+    Extent.Append horizA horizB ~ horizC,+    Shape.C width, Eq width, Shape.C heightA, Shape.C heightB,+    Class.Floating a) =>+   Full Extent.Big horizA heightA width a ->+   Full Extent.Big horizB heightB width a ->+   Full Extent.Big horizC (heightA:+:heightB) width a+(===) = above rightBias Extent.appendAny++beside ::+   (Extent.C vertA, Extent.C vertB, Extent.C vertC,+    Shape.C height, Eq height, Shape.C widthA, Shape.C widthB,+    Class.Floating a) =>+   Basic.OrderBias ->+   Extent.AppendMode vertA vertB vertC height widthA widthB ->+   Full vertA Extent.Big height widthA a ->+   Full vertB Extent.Big height widthB a ->+   Full vertC Extent.Big height (widthA:+:widthB) a+beside orderBias = ArrMatrix.lift2 . Basic.beside orderBias++above ::+   (Extent.C horizA, Extent.C horizB, Extent.C horizC,+    Shape.C width, Eq width, Shape.C heightA, Shape.C heightB,+    Class.Floating a) =>+   Basic.OrderBias ->+   Extent.AppendMode horizA horizB horizC width heightA heightB ->+   Full Extent.Big horizA heightA width a ->+   Full Extent.Big horizB heightB width a ->+   Full Extent.Big horizC (heightA:+:heightB) width a+above orderBias = ArrMatrix.lift2 . Basic.above orderBias++{- |+Use the element order of the first operand.+-}+leftBias :: Basic.OrderBias+leftBias = Basic.LeftBias++{- |+Use the element order of the second operand.+-}+rightBias :: Basic.OrderBias+rightBias = Basic.RightBias++{- |+Choose element order such that, if possible,+one part can be copied as one block.+For 'above' this means that 'RowMajor' is chosen+whenever at least one operand is 'RowMajor'+and 'ColumnMajor' is chosen when both operands are 'ColumnMajor'.+-}+contiguousBias :: Basic.OrderBias+contiguousBias = Basic.ContiguousBias+++rowSums ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Class.Floating a) =>+   Full vert horiz height width a -> Vector height a+rowSums = Plain.rowSums . ArrMatrix.toVector++columnSums ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Class.Floating a) =>+   Full vert horiz height width a -> Vector width a+columnSums = Plain.columnSums . ArrMatrix.toVector+++rowArgAbsMaximums ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.InvIndexed width, Shape.Index width ~ ix, Storable ix,+    Class.Floating a) =>+   Full vert horiz height width a -> (Vector height ix, Vector height a)+rowArgAbsMaximums = Plain.rowArgAbsMaximums . ArrMatrix.toVector++columnArgAbsMaximums ::+   (Extent.C vert, Extent.C horiz,+    Shape.InvIndexed height, Shape.C width,+    Shape.Index height ~ ix, Storable ix,+    Class.Floating a) =>+   Full vert horiz height width a -> (Vector width ix, Vector width a)+columnArgAbsMaximums = Plain.columnArgAbsMaximums . ArrMatrix.toVector+++map ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,+    Storable a, Storable b) =>+   (a -> b) ->+   Full vert horiz height width a -> Full vert horiz height width b+map = ArrMatrix.lift1 . Array.map+++infixl 7 |*-++(|*-) ::+   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+   Vector height a -> Vector width a -> General height width a+x|*-y = ArrMatrix.lift0 $ x Plain.|*- y++{- |+> tensorProduct order x y = singleColumn order x #*# singleRow order y+-}+tensorProduct ::+   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+   Order -> Vector height a -> Vector width a -> General height width a+tensorProduct order x y = ArrMatrix.lift0 $ Plain.tensorProduct order x y++{- |+> outer order x y = tensorProduct order x (Vector.conjugate y)+-}+outer ::+   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+   Order -> Vector height a -> Vector width a -> General height width a+outer order x y = ArrMatrix.lift0 $ Plain.outer order x y++kronecker ::+   (Extent.C vert, Extent.C horiz,+    Shape.C heightA, Shape.C widthA, Shape.C heightB, Shape.C widthB,+    Class.Floating a) =>+   Full vert horiz heightA widthA a ->+   Full vert horiz heightB widthB a ->+   Full vert horiz (heightA,heightB) (widthA,widthB) a+kronecker = ArrMatrix.lift2 Plain.kronecker++sumRank1 ::+   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+   (height,width) ->+   [(a, (Vector height a, Vector width a))] -> General height width a+sumRank1 dims = ArrMatrix.lift0 . Plain.sumRank1 dims+++infixl 7 #*\, #/\+infixr 7 \*#, \\#++scaleRows, (\*#) ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+   Vector height a ->+   Full vert horiz height width a ->+   Full vert horiz height width a+scaleRows = ArrMatrix.lift1 . Basic.scaleRows+(\*#) = scaleRows++scaleColumns ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   Vector width a ->+   Full vert horiz height width a ->+   Full vert horiz height width a+scaleColumns = ArrMatrix.lift1 . Basic.scaleColumns++(#*\) ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   Full vert horiz height width a ->+   Vector width a ->+   Full vert horiz height width a+(#*\) = flip scaleColumns++scaleRowsReal ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Eq height, Shape.C width,+    Class.Floating a) =>+   Vector height (RealOf a) ->+   Full vert horiz height width a ->+   Full vert horiz height width a+scaleRowsReal = ArrMatrix.lift1 . Basic.scaleRowsReal++scaleColumnsReal ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   Vector width (RealOf a) ->+   Full vert horiz height width a ->+   Full vert horiz height width a+scaleColumnsReal = ArrMatrix.lift1 . Basic.scaleColumnsReal++(\\#) ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+   Vector height a ->+   Full vert horiz height width a ->+   Full vert horiz height width a+(\\#) = scaleRows . Vector.recip++(#/\) ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   Full vert horiz height width a ->+   Vector width a ->+   Full vert horiz height width a+(#/\) a x = scaleColumns (Vector.recip x) a+++multiplySquare ::+   (Multiply.MultiplySquare typ,+    Type.HeightOf typ ~ height, Eq height, Shape.C width,+    Extent.C horiz, Extent.C vert, Class.Floating a) =>+   Mod.Transposition -> Type.Matrix typ a ->+   Full vert horiz height width a -> Full vert horiz height width a+multiplySquare = Multiply.transposableSquare++multiplyVector ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width, Eq width,+    Class.Floating a) =>+   Full vert horiz height width a -> Vector width a -> Vector height a+multiplyVector = Basic.multiplyVector . ArrMatrix.toVector++multiply ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height,+    Shape.C fuse, Eq fuse,+    Shape.C width,+    Class.Floating a) =>+   Full vert horiz height fuse a ->+   Full vert horiz fuse width a ->+   Full vert horiz height width a+multiply = ArrMatrix.lift2 Basic.multiply

@@ -0,0 +1,237 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE EmptyDataDecls #-}+module Numeric.LAPACK.Matrix.Array (+   Matrix(Array),+   ArrayMatrix,+   Array,++   Full,+   General,+   Tall,+   Wide,+   Square,++   shape,+   reshape,+   mapShape,+   toVector,+   fromVector,+   lift0,+   lift1,+   lift2,+   lift3,+   lift4,+   unlift1,+   unlift2,+   unliftRow,+   unliftColumn,++   Plain.Homogeneous, zero, negate, scaleReal, scale, scaleRealReal, (.*#),+   Plain.ShapeOrder, forceOrder, Plain.shapeOrder, adaptOrder,+   Plain.Additive, add, sub, (#+#), (#-#),+   Plain.Complex,+   Plain.SquareShape,+   Multiply.MultiplyVector,+   Multiply.MultiplySquare,+   Multiply.Power,+   Multiply.Multiply,+   Divide.Determinant,+   Divide.Solve,+   Divide.Inverse,+   ) where++import qualified Numeric.LAPACK.Matrix.Plain.Divide as Divide+import qualified Numeric.LAPACK.Matrix.Plain.Multiply as Multiply+import qualified Numeric.LAPACK.Matrix.Plain.Class as Plain+import qualified Numeric.LAPACK.Matrix.Type as Type+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Shape.Box as Box+import qualified Numeric.LAPACK.Matrix.Basic as Basic+import Numeric.LAPACK.Matrix.Plain.Format (FormatArray, formatArray)+import Numeric.LAPACK.Matrix.Type (Matrix)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Scalar (RealOf)++import qualified Numeric.Netlib.Class as Class++import qualified Control.DeepSeq as DeepSeq++import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Storable as CheckedArray+import qualified Data.Array.Comfort.Shape as Shape++import Prelude hiding (negate)+++data Array shape+newtype instance Matrix (Array shape) a = Array (Array.Array shape a)+   deriving (Show)++type ArrayMatrix shape = Matrix (Array shape)+++type Full vert horiz height width =+         ArrayMatrix (MatrixShape.Full vert horiz height width)+type General height width = ArrayMatrix (MatrixShape.General height width)+type Tall height width = ArrayMatrix (MatrixShape.Tall height width)+type Wide height width = ArrayMatrix (MatrixShape.Wide height width)+type Square sh = ArrayMatrix (MatrixShape.Square sh)+++instance (DeepSeq.NFData shape) => Type.NFData (Array shape) where+   rnf (Array arr) = DeepSeq.rnf arr++instance (Box.Box sh) => Type.Box (Array sh) where+   type HeightOf (Array sh) = Box.HeightOf sh+   type WidthOf (Array sh) = Box.WidthOf sh+   height (Array arr) = Box.height $ Array.shape arr+   width (Array arr) = Box.width $ Array.shape arr+++shape :: ArrayMatrix sh a -> sh+shape (Array a) = Array.shape a++reshape ::+   (Shape.C sh0, Shape.C sh1) =>+   sh1 -> ArrayMatrix sh0 a -> ArrayMatrix sh1 a+reshape = lift1 . CheckedArray.reshape++mapShape ::+   (Shape.C sh0, Shape.C sh1) =>+   (sh0 -> sh1) -> ArrayMatrix sh0 a -> ArrayMatrix sh1 a+mapShape = lift1 . CheckedArray.mapShape+++toVector :: ArrayMatrix sh a -> Array.Array sh a+toVector (Array a) = a++fromVector ::+   (Plain.Admissible sh, Class.Floating a) =>+   Array.Array sh a -> ArrayMatrix sh a+fromVector arr =+   Array $+   case Plain.check arr of+      Nothing -> arr+      Just msg -> error $ "Matrix.Array.fromVector: " ++ msg++{- |+'lift0' is a synonym for 'fromVector' but lacks the admissibility check.+You may thus fool the type tags.+This applies to the other lift functions, too.+-}+lift0 :: Array.Array shA a -> ArrayMatrix shA a+lift0 = Array++lift1 ::+   (Array.Array shA a -> Array.Array shB b) ->+   ArrayMatrix shA a -> ArrayMatrix shB b+lift1 f (Array a) = Array $ f a++lift2 ::+   (Array.Array shA a -> Array.Array shB b -> Array.Array shC c) ->+   ArrayMatrix shA a -> ArrayMatrix shB b -> ArrayMatrix shC c+lift2 f (Array a) (Array b) = Array $ f a b++lift3 ::+   (Array.Array shA a -> Array.Array shB b ->+    Array.Array shC c -> Array.Array shD d) ->+   ArrayMatrix shA a -> ArrayMatrix shB b ->+   ArrayMatrix shC c -> ArrayMatrix shD d+lift3 f (Array a) (Array b) (Array c) = Array $ f a b c++lift4 ::+   (Array.Array shA a -> Array.Array shB b ->+    Array.Array shC c -> Array.Array shD d ->+    Array.Array shE e) ->+   ArrayMatrix shA a -> ArrayMatrix shB b ->+   ArrayMatrix shC c -> ArrayMatrix shD d ->+   ArrayMatrix shE e+lift4 f (Array a) (Array b) (Array c) (Array d) = Array $ f a b c d+++unlift1 ::+   (ArrayMatrix shA a -> ArrayMatrix shB b) ->+   Array.Array shA a -> Array.Array shB b+unlift1 f a = toVector $ f $ Array a++unlift2 ::+   (ArrayMatrix shA a -> ArrayMatrix shB b -> ArrayMatrix shC c) ->+   Array.Array shA a -> Array.Array shB b -> Array.Array shC c+unlift2 f a b = toVector $ f (Array a) (Array b)+++unliftRow ::+   MatrixShape.Order ->+   (General () height0 a -> General () height1 b) ->+   Vector height0 a -> Vector height1 b+unliftRow order = Basic.unliftRow order . unlift1++unliftColumn ::+   MatrixShape.Order ->+   (General height0 () a -> General height1 () b) ->+   Vector height0 a -> Vector height1 b+unliftColumn order = Basic.unliftColumn order . unlift1+++instance (FormatArray sh) => Type.FormatMatrix (Array sh) where+   formatMatrix fmt (Array a) = formatArray fmt a++instance (Multiply.MultiplySame sh) => Type.MultiplySame (Array sh) where+   multiplySame = lift2 Multiply.same+++zero ::+   (Plain.Homogeneous shape, Class.Floating a) => shape -> ArrayMatrix shape a+zero = lift0 . Plain.zero++negate ::+   (Plain.Homogeneous shape, Class.Floating a) =>+   ArrayMatrix shape a -> ArrayMatrix shape a+negate = lift1 Plain.negate++scaleReal ::+   (Plain.Homogeneous shape, Class.Floating a) =>+   RealOf a -> ArrayMatrix shape a -> ArrayMatrix shape a+scaleReal = lift1 . Plain.scaleReal++newtype ScaleReal f a = ScaleReal {getScaleReal :: a -> f a -> f a}++scaleRealReal ::+   (Plain.Homogeneous shape, Class.Real a) =>+   a -> ArrayMatrix shape a -> ArrayMatrix shape a+scaleRealReal =+   getScaleReal $ Class.switchReal (ScaleReal scaleReal) (ScaleReal scaleReal)+++scale, (.*#) ::+   (Multiply.Scale shape, Class.Floating a) =>+   a -> ArrayMatrix shape a -> ArrayMatrix shape a+scale = lift1 . Multiply.scale+(.*#) = scale++infixl 7 .*#+++forceOrder ::+   (Plain.ShapeOrder shape, Class.Floating a) =>+   MatrixShape.Order -> ArrayMatrix shape a -> ArrayMatrix shape a+forceOrder = lift1 . Plain.forceOrder++{- |+@adaptOrder x y@ contains the data of @y@ with the layout of @x@.+-}+adaptOrder ::+   (Plain.ShapeOrder shape, Class.Floating a) =>+   ArrayMatrix shape a -> ArrayMatrix shape a -> ArrayMatrix shape a+adaptOrder = lift2 Plain.adaptOrder+++infixl 6 #+#, #-#, `add`, `sub`++add, sub, (#+#), (#-#) ::+   (Plain.Additive shape, Class.Floating a) =>+   ArrayMatrix shape a -> ArrayMatrix shape a -> ArrayMatrix shape a+add = lift2 Plain.add+sub = lift2 Plain.sub+(#+#) = add+(#-#) = sub

@@ -0,0 +1,24 @@+module Numeric.LAPACK.Matrix.Array.Banded where++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Numeric.LAPACK.Matrix.Array (ArrayMatrix)++import qualified Type.Data.Num.Unary.Literal as TypeNum+++type Banded sub super vert horiz height width =+      ArrayMatrix (MatrixShape.Banded sub super vert horiz height width)++type General sub super height width =+      ArrayMatrix (MatrixShape.BandedGeneral sub super height width)++type Square sub super size =+      ArrayMatrix (MatrixShape.BandedSquare sub super size)++type Lower sub size = Square sub TypeNum.U0 size+type Upper super size = Square TypeNum.U0 super size++type Diagonal size = Square TypeNum.U0 TypeNum.U0 size++type Hermitian offDiag sh =+         ArrayMatrix (MatrixShape.BandedHermitian offDiag sh)

@@ -0,0 +1,22 @@+module Numeric.LAPACK.Matrix.Array.Basic where++import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix.Basic as Basic+import qualified Numeric.LAPACK.Matrix.Extent as Extent+import Numeric.LAPACK.Matrix.Array (Full)+++transpose ::+   (Extent.C vert, Extent.C horiz) =>+   Full vert horiz height width a -> Full horiz vert width height a+transpose = ArrMatrix.lift1 Basic.transpose++swapMultiply ::+   (Extent.C vertA, Extent.C vertB, Extent.C horizA, Extent.C horizB) =>+   (matrix ->+    Full horizA vertA widthA heightA a ->+    Full horizB vertB widthB heightB a) ->+   Full vertA horizA heightA widthA a ->+   matrix ->+   Full vertB horizB heightB widthB a+swapMultiply multiplyTrans a b = transpose $ multiplyTrans b $ transpose a

@@ -0,0 +1,26 @@+module Numeric.LAPACK.Matrix.Array.Triangular where++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Numeric.LAPACK.Matrix.Array (ArrayMatrix)+import Numeric.LAPACK.Matrix.Shape.Private (NonUnit, Unit)+++type Hermitian sh = ArrayMatrix (MatrixShape.Hermitian sh)++type Lower sh = FlexLower NonUnit sh+type Upper sh = FlexUpper NonUnit sh+type Symmetric sh = FlexSymmetric NonUnit sh+type Diagonal sh = FlexDiagonal NonUnit sh++type UnitLower sh = FlexLower Unit sh+type UnitUpper sh = FlexUpper Unit sh++type FlexDiagonal diag sh =+      ArrayMatrix+         (MatrixShape.Triangular MatrixShape.Empty diag MatrixShape.Empty sh)+type FlexLower diag sh = ArrayMatrix (MatrixShape.LowerTriangular diag sh)+type FlexUpper diag sh = ArrayMatrix (MatrixShape.UpperTriangular diag sh)+type FlexSymmetric diag sh = ArrayMatrix (MatrixShape.FlexSymmetric diag sh)++type Triangular lo diag up sh =+         ArrayMatrix (MatrixShape.Triangular lo diag up sh)

@@ -1,28 +1,184 @@ {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE TypeOperators #-} module Numeric.LAPACK.Matrix.Banded (-   module Numeric.LAPACK.Matrix.Banded.Basic,+   Banded,+   General,+   Square,+   Upper,+   Lower,+   Diagonal,+   Hermitian,    height, width,+   fromList,+   squareFromList,+   lowerFromList,+   upperFromList,+   mapExtent,+   diagonal,+   takeDiagonal,+   toFull,+   toLowerTriangular,+   toUpperTriangular,+   transpose,+   adjoint,+   multiplyVector,+   multiply,+   multiplyFull,     solve,    determinant,    ) where -import Numeric.LAPACK.Matrix.Banded.Basic-import Numeric.LAPACK.Matrix.Banded.Linear+import qualified Numeric.LAPACK.Matrix.Banded.Linear as Linear+import qualified Numeric.LAPACK.Matrix.Banded.Basic as Basic +import qualified Numeric.LAPACK.Matrix.Array.Triangular as Tri+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent+import Numeric.LAPACK.Matrix.Array.Banded+         (Banded, General, Square, Lower, Upper, Diagonal, Hermitian)+import Numeric.LAPACK.Matrix.Array (Full)+import Numeric.LAPACK.Matrix.Shape.Private (Order, UnaryProxy)+import Numeric.LAPACK.Vector (Vector) -import qualified Data.Array.Comfort.Storable as Array+import qualified Numeric.Netlib.Class as Class +import qualified Type.Data.Num.Unary as Unary+import Type.Data.Num.Unary ((:+:)) +import qualified Data.Array.Comfort.Shape as Shape++import Foreign.Storable (Storable)++ height ::    (Extent.C vert, Extent.C horiz) =>    Banded sub super vert horiz height width a -> height-height = MatrixShape.bandedHeight . Array.shape+height = MatrixShape.bandedHeight . ArrMatrix.shape  width ::    (Extent.C vert, Extent.C horiz) =>    Banded sub super vert horiz height width a -> width-width = MatrixShape.bandedWidth . Array.shape+width = MatrixShape.bandedWidth . ArrMatrix.shape++++fromList ::+   (Unary.Natural sub, Unary.Natural super,+    Shape.C height, Shape.C width, Storable a) =>+   (UnaryProxy sub, UnaryProxy super) -> Order -> height -> width -> [a] ->+   General sub super height width a+fromList offDiag order height_ width_ =+   ArrMatrix.lift0 . Basic.fromList offDiag order height_ width_++squareFromList ::+   (Unary.Natural sub, Unary.Natural super, Shape.C size, Storable a) =>+   (UnaryProxy sub, UnaryProxy super) -> Order -> size -> [a] ->+   Square sub super size a+squareFromList offDiag order size =+   ArrMatrix.lift0 . Basic.squareFromList offDiag order size++lowerFromList ::+   (Unary.Natural sub, Shape.C size, Storable a) =>+   UnaryProxy sub -> Order -> size -> [a] -> Lower sub size a+lowerFromList numOff order size =+   ArrMatrix.lift0 . Basic.lowerFromList numOff order size++upperFromList ::+   (Unary.Natural super, Shape.C size, Storable a) =>+   UnaryProxy super -> Order -> size -> [a] -> Upper super size a+upperFromList numOff order size =+   ArrMatrix.lift0 . Basic.upperFromList numOff order size++mapExtent ::+   (Extent.C vertA, Extent.C horizA) =>+   (Extent.C vertB, Extent.C horizB) =>+   Extent.Map vertA horizA vertB horizB height width ->+   Banded super sub vertA horizA height width a ->+   Banded super sub vertB horizB height width a+mapExtent = ArrMatrix.lift1 . Basic.mapExtent++transpose ::+   (Extent.C vert, Extent.C horiz) =>+   Banded sub super vert horiz height width a ->+   Banded super sub horiz vert width height a+transpose = ArrMatrix.lift1 Basic.transpose++adjoint ::+   (Unary.Natural sub, Unary.Natural super, Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Class.Floating a) =>+   Banded sub super vert horiz height width a ->+   Banded super sub horiz vert width height a+adjoint = ArrMatrix.lift1 Basic.adjoint++diagonal ::+   (Shape.C sh, Class.Floating a) => Order -> Vector sh a -> Diagonal sh a+diagonal order = ArrMatrix.lift0 . Basic.diagonal order++takeDiagonal ::+   (Unary.Natural sub, Unary.Natural super, Shape.C sh, Class.Floating a) =>+   Square sub super sh a -> Vector sh a+takeDiagonal = Basic.takeDiagonal . ArrMatrix.toVector++multiplyVector ::+   (Unary.Natural sub, Unary.Natural super,+    Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width, Eq width,+    Class.Floating a) =>+   Banded sub super vert horiz height width a ->+   Vector width a -> Vector height a+multiplyVector = Basic.multiplyVector . ArrMatrix.toVector++multiply ::+   (Unary.Natural subA, Unary.Natural superA,+    Unary.Natural subB, Unary.Natural superB,+    (subA :+: subB) ~ subC,+    (superA :+: superB) ~ superC,+    Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Shape.C fuse, Eq fuse,+    Class.Floating a) =>+   Banded subA superA vert horiz height fuse a ->+   Banded subB superB vert horiz fuse width a ->+   Banded subC superC vert horiz height width a+multiply = ArrMatrix.lift2 Basic.multiply++multiplyFull ::+   (Unary.Natural sub, Unary.Natural super,+    Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Shape.C fuse, Eq fuse,+    Class.Floating a) =>+   Banded sub super vert horiz height fuse a ->+   Full vert horiz fuse width a -> Full vert horiz height width a+multiplyFull = ArrMatrix.lift2 Basic.multiplyFull++toLowerTriangular ::+   (Unary.Natural sub, Shape.C sh, Class.Floating a) =>+   Lower sub sh a -> Tri.Lower sh a+toLowerTriangular = ArrMatrix.lift1 Basic.toLowerTriangular++toUpperTriangular ::+   (Unary.Natural super, Shape.C sh, Class.Floating a) =>+   Upper super sh a -> Tri.Upper sh a+toUpperTriangular = ArrMatrix.lift1 Basic.toUpperTriangular++toFull ::+   (Unary.Natural sub, Unary.Natural super,+    Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,+    Class.Floating a) =>+   Banded sub super vert horiz height width a ->+   Full vert horiz height width a+toFull = ArrMatrix.lift1 Basic.toFull++++solve ::+   (Unary.Natural sub, Unary.Natural super, Extent.C vert, Extent.C horiz,+    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>+   Square sub super sh a ->+   Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a+solve = ArrMatrix.lift2 Linear.solve++determinant ::+   (Unary.Natural sub, Unary.Natural super, Shape.C sh, Class.Floating a) =>+   Square sub super sh a -> a+determinant = Linear.determinant . ArrMatrix.toVector

@@ -12,7 +12,11 @@    lowerFromList,    upperFromList,    mapExtent,+   mapHeight,+   mapWidth,+   identityFatOrder,    diagonal,+   fromDiagonal,    takeDiagonal,    toFull,    toLowerTriangular,@@ -28,12 +32,15 @@ import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent import qualified Numeric.LAPACK.Matrix.Triangular.Private as TriangularPriv import qualified Numeric.LAPACK.Matrix.Triangular.Basic as Triangular+import qualified Numeric.LAPACK.Matrix.RowMajor as RowMajor import qualified Numeric.LAPACK.Matrix.Private as Matrix import qualified Numeric.LAPACK.Vector as Vector import qualified Numeric.LAPACK.Private as Private+import qualified Numeric.LAPACK.ShapeStatic as ShapeStatic import Numeric.LAPACK.Matrix.Shape.Private          (Order(RowMajor,ColumnMajor), transposeFromOrder, swapOnRowMajor,           UnaryProxy, addOffDiagonals)+import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated)) import Numeric.LAPACK.Vector (Vector) import Numeric.LAPACK.Scalar (zero, one) import Numeric.LAPACK.Private@@ -132,6 +139,28 @@    Banded super sub vertB horizB height width a mapExtent f = Array.mapShape $ MatrixShape.bandedMapExtent f +mapHeight ::+   (Extent.GeneralTallWide vert horiz,+    Extent.GeneralTallWide horiz vert) =>+   (heightA -> heightB) ->+   Banded super sub vert horiz heightA width a ->+   Banded super sub vert horiz heightB width a+mapHeight f =+   Array.mapShape+      (\(MatrixShape.Banded offDiag order extent) ->+         MatrixShape.Banded offDiag order $ Extent.mapHeight f extent)++mapWidth ::+   (Extent.GeneralTallWide vert horiz,+    Extent.GeneralTallWide horiz vert) =>+   (widthA -> widthB) ->+   Banded super sub vert horiz height widthA a ->+   Banded super sub vert horiz height widthB a+mapWidth f =+   Array.mapShape+      (\(MatrixShape.Banded offDiag order extent) ->+         MatrixShape.Banded offDiag order $ Extent.mapWidth f extent)+ transpose ::    (Extent.C vert, Extent.C horiz) =>    Banded sub super vert horiz height width a ->@@ -139,16 +168,67 @@ transpose = Array.mapShape MatrixShape.bandedTranspose  adjoint ::-   (Unary.Natural super, Unary.Natural sub, Extent.C vert, Extent.C horiz,-    Shape.C width, Shape.C height, Class.Floating a) =>+   (Unary.Natural sub, Unary.Natural super, Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Class.Floating a) =>    Banded sub super vert horiz height width a ->    Banded super sub horiz vert width height a adjoint = Vector.conjugate . transpose  -diagonal :: (Shape.C sh, Class.Floating a) => Vector sh a -> Diagonal sh a-diagonal (Array sh x) =-   Array (MatrixShape.bandedSquare (Proxy,Proxy) ColumnMajor sh) x+identityFatOrder ::+   (Unary.Natural sub, Unary.Natural super, Shape.C size, Class.Floating a) =>+   Order -> size -> Square sub super size a+identityFatOrder order =+   case order of+      RowMajor -> fromRowMajor Extent.square . identityStripe+      ColumnMajor -> fromColumnMajor Extent.square . identityStripe++type Section sub super =+      ShapeStatic.ZeroBased sub Shape.:+: ()+         Shape.:+: ShapeStatic.ZeroBased super++identityStripe ::+   (Unary.Natural sub, Unary.Natural super, Shape.C height, Class.Floating a) =>+   height -> RowMajor.Matrix height (Section sub super) a+identityStripe sh =+   RowMajor.tensorProduct (Left NonConjugated)+      (Vector.one sh) (Vector.unit Shape.static $ Right (Left ()))++fromRowMajor ::+   (Unary.Natural sub, Unary.Natural super, Shape.C height) =>+   (height -> Extent.Extent vert horiz height width) ->+   RowMajor.Matrix height (Section sub super) a ->+   Banded sub super vert horiz height width a+fromRowMajor extent =+   Array.mapShape+      (\(size,+         ShapeStatic.ZeroBased kl Shape.:+: ()+            Shape.:+: ShapeStatic.ZeroBased ku) ->+         MatrixShape.Banded (kl,ku) RowMajor $ extent size)++fromColumnMajor ::+   (Unary.Natural sub, Unary.Natural super, Shape.C width) =>+   (width -> Extent.Extent vert horiz height width) ->+   RowMajor.Matrix width (Section super sub) a ->+   Banded sub super vert horiz height width a+fromColumnMajor extent =+   Array.mapShape+      (\(size,+         ShapeStatic.ZeroBased ku Shape.:+: ()+            Shape.:+: ShapeStatic.ZeroBased kl) ->+         MatrixShape.Banded (kl,ku) ColumnMajor $ extent size)+++diagonal ::+   (Shape.C sh, Class.Floating a) => Order -> Vector sh a -> Diagonal sh a+diagonal order (Array sh x) =+   Array (MatrixShape.bandedSquare (Proxy,Proxy) order sh) x++fromDiagonal ::+   (Shape.C sh, Class.Floating a) =>+   TriangularPriv.FlexDiagonal diag sh a -> Diagonal sh a+fromDiagonal (Array (MatrixShape.Triangular _diag _uplo order sh) x) =+   Array (MatrixShape.bandedSquare (Proxy,Proxy) order sh) x  takeDiagonal ::    (Unary.Natural sub, Unary.Natural super, Shape.C sh, Class.Floating a) =>

@@ -8,7 +8,7 @@ import qualified Numeric.LAPACK.Matrix.Banded.Basic as Banded import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent-import qualified Numeric.LAPACK.Split as Split+import qualified Numeric.LAPACK.Permutation.Private as Perm import qualified Numeric.LAPACK.Private as Private import Numeric.LAPACK.Linear.Private (solver, withDeterminantInfo, withInfo) import Numeric.LAPACK.Matrix.Shape.Private@@ -112,4 +112,4 @@             (do                det <- Private.product n (advancePtr abPtr (kl+ku)) ldab                ipiv <- peekArray n ipivPtr-               return $ if Split.oddPermutation ipiv then -det else det)+               return $ Perm.condNegate ipiv det)

@@ -1,34 +1,131 @@ {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-} module Numeric.LAPACK.Matrix.BandedHermitian (-   module Numeric.LAPACK.Matrix.BandedHermitian.Basic,+   BandedHermitian,+   Transposition(..),+   size,+   fromList,+   identity,+   diagonal,+   takeDiagonal,+   toHermitian,+   toBanded,+   multiplyVector,+   multiplyFull,+   gramian,+   sumRank1,     eigenvalues,    eigensystem,    ) where  import qualified Numeric.LAPACK.Matrix.BandedHermitian.Eigen as Eigen-import Numeric.LAPACK.Matrix.BandedHermitian.Basic+import qualified Numeric.LAPACK.Matrix.BandedHermitian.Basic as Basic -import qualified Numeric.LAPACK.Matrix.Private as Matrix+import qualified Numeric.LAPACK.Matrix.Array.Banded as Banded+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Extent as Extent+import Numeric.LAPACK.Matrix.Array.Banded (Square)+import Numeric.LAPACK.Matrix.Array.Triangular (Hermitian)+import Numeric.LAPACK.Matrix.Array (Full)+import Numeric.LAPACK.Matrix.Shape.Private (Order, UnaryProxy)+import Numeric.LAPACK.Matrix.Modifier (Transposition(NonTransposed, Transposed)) import Numeric.LAPACK.Vector (Vector) import Numeric.LAPACK.Scalar (RealOf)  import qualified Numeric.Netlib.Class as Class +import qualified Type.Data.Num.Unary.Literal as TypeNum import qualified Type.Data.Num.Unary as Unary+import Type.Data.Num.Unary ((:+:))  import qualified Data.Array.Comfort.Shape as Shape +import Foreign.Storable (Storable) +import Data.Tuple.HT (mapFst)+++type BandedHermitian offDiag sh = Banded.Hermitian offDiag sh+type Diagonal size = BandedHermitian TypeNum.U0 size+++size :: BandedHermitian offDiag sh a -> sh+size = MatrixShape.bandedHermitianSize . ArrMatrix.shape+++fromList ::+   (Unary.Natural offDiag, Shape.C size, Storable a) =>+   UnaryProxy offDiag -> Order -> size -> [a] ->+   BandedHermitian offDiag size a+fromList numOff order size_ =+   ArrMatrix.lift0 . Basic.fromList numOff order size_++identity :: (Shape.C sh, Class.Floating a) => sh -> Diagonal sh a+identity = ArrMatrix.lift0 . Basic.identity++diagonal ::+   (Shape.C sh, Class.Floating a) => Vector sh (RealOf a) -> Diagonal sh a+diagonal = ArrMatrix.lift0 . Basic.diagonal++takeDiagonal ::+   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>+   BandedHermitian offDiag size a -> Vector size (RealOf a)+takeDiagonal = Basic.takeDiagonal . ArrMatrix.toVector++toHermitian ::+   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>+   BandedHermitian offDiag size a -> Hermitian size a+toHermitian = ArrMatrix.lift1 Basic.toHermitian++toBanded ::+   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>+   BandedHermitian offDiag size a ->+   Square offDiag offDiag size a+toBanded = ArrMatrix.lift1 Basic.toBanded++multiplyVector ::+   (Unary.Natural offDiag, Shape.C size, Eq size, Class.Floating a) =>+   Transposition -> BandedHermitian offDiag size a ->+   Vector size a -> Vector size a+multiplyVector transposed = Basic.multiplyVector transposed .  ArrMatrix.toVector++gramian ::+   (Shape.C size, Eq size, Class.Floating a,+    Unary.Natural sub, Unary.Natural super) =>+   Square sub super size a ->+   BandedHermitian (sub :+: super) size a+gramian = ArrMatrix.lift1 Basic.gramian++multiplyFull ::+   (Unary.Natural offDiag, Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+   Transposition -> BandedHermitian offDiag height a ->+   Full vert horiz height width a ->+   Full vert horiz height width a+multiplyFull = ArrMatrix.lift2 . Basic.multiplyFull++{- |+The list represents ragged rows of a sparse matrix.+-}+sumRank1 ::+   (Unary.Natural k, Shape.Indexed sh, Class.Floating a) =>+   Order -> sh ->+   [(RealOf a, (Shape.Index sh, Basic.StaticVector (Unary.Succ k) a))] ->+   BandedHermitian k sh a+sumRank1 order sh = ArrMatrix.lift0 . Basic.sumRank1 order sh++ eigenvalues ::    (Unary.Natural offDiag, Shape.C sh, Class.Floating a) =>    BandedHermitian offDiag sh a -> Vector sh (RealOf a)-eigenvalues = Eigen.values+eigenvalues = Eigen.values . ArrMatrix.toVector  {- | For symmetric eigenvalue problems, @eigensystem@ and @schur@ coincide. -} eigensystem ::    (Unary.Natural offDiag, Shape.C sh, Class.Floating a) =>-   BandedHermitian offDiag sh a -> (Matrix.Square sh a, Vector sh (RealOf a))-eigensystem = Eigen.decompose+   BandedHermitian offDiag sh a -> (ArrMatrix.Square sh a, Vector sh (RealOf a))+eigensystem = mapFst ArrMatrix.lift0 . Eigen.decompose . ArrMatrix.toVector

@@ -3,37 +3,44 @@ {-# LANGUAGE GADTs #-} module Numeric.LAPACK.Matrix.BandedHermitian.Basic (    BandedHermitian,+   StaticVector,    Transposition(..),    fromList,    identity,+   identityFatOrder,    diagonal,    takeDiagonal,    toHermitian,    toBanded,    multiplyVector,    multiplyFull,-   covariance,+   gramian,    sumRank1,+   takeUpper,    ) where -import qualified Numeric.LAPACK.ShapeStatic as ShapeStatic import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape import qualified Numeric.LAPACK.Matrix.Extent as Extent import qualified Numeric.LAPACK.Matrix.Banded.Basic as Banded import qualified Numeric.LAPACK.Matrix.Triangular.Private as TriangularPriv+import qualified Numeric.LAPACK.Matrix.RowMajor as RowMajor import qualified Numeric.LAPACK.Matrix.Private as Matrix+import qualified Numeric.LAPACK.Vector.Private as VectorPriv import qualified Numeric.LAPACK.Vector as Vector+import qualified Numeric.LAPACK.ShapeStatic as ShapeStatic import Numeric.LAPACK.Matrix.Hermitian.Private (TakeDiagonal(..)) import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian) import Numeric.LAPACK.Matrix.Shape.Private          (Order(RowMajor,ColumnMajor), flipOrder, uploFromOrder,           UnaryProxy, natFromProxy)-import Numeric.LAPACK.Matrix.Private-         (Transposition(NonTransposed, Transposed), transposeOrder)+import Numeric.LAPACK.Matrix.Modifier+         (Transposition(NonTransposed, Transposed), transposeOrder,+          Conjugation(NonConjugated), conjugatedOnRowMajor) import Numeric.LAPACK.Vector (Vector) import Numeric.LAPACK.Scalar (RealOf, zero, one) import Numeric.LAPACK.Private-         (fill, lacgv, copyConjugate, condConjugateToTemp,+         (fill, lacgv, caseRealComplexFunc, realPtr,+          copyConjugate, condConjugate, condConjugateToTemp,           pointerSeq, pokeCInt, copySubMatrix)  import qualified Numeric.BLAS.FFI.Generic as BlasGen@@ -57,12 +64,11 @@ import Foreign.Marshal.Array (advancePtr) import Foreign.C.Types (CInt, CChar) import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)-import Foreign.Ptr (Ptr, castPtr)+import Foreign.Ptr (Ptr) import Foreign.Storable (Storable, poke, peek, peekElemOff)  import Control.Monad.Trans.Cont (ContT(ContT), evalContT) import Control.Monad.IO.Class (liftIO)-import Control.Monad (when)  import Data.Foldable (for_) import Data.Tuple.HT (mapPair)@@ -83,11 +89,42 @@ fromList numOff order size =    CheckedArray.fromList (MatrixShape.BandedHermitian numOff order size) +identityFatOrder ::+   (Unary.Natural offDiag, Shape.C sh, Class.Floating a) =>+   Order -> sh -> BandedHermitian offDiag sh a+identityFatOrder order sh =+   case order of+      RowMajor ->+         fromRowMajor $+         RowMajor.tensorProduct (Left NonConjugated)+            (Vector.one sh) (Vector.unit Shape.static $ Left ())+      ColumnMajor ->+         fromColumnMajor $+         RowMajor.tensorProduct (Left NonConjugated)+            (Vector.one sh) (Vector.unit Shape.static $ Right ())++fromRowMajor ::+   (Unary.Natural offDiag, Shape.C size) =>+   RowMajor.Matrix size (() Shape.:+: ShapeStatic.ZeroBased offDiag) a ->+   BandedHermitian offDiag size a+fromRowMajor =+   Array.mapShape+      (\(size, () Shape.:+: ShapeStatic.ZeroBased k) ->+         MatrixShape.BandedHermitian k RowMajor size)++fromColumnMajor ::+   (Unary.Natural offDiag, Shape.C size) =>+   RowMajor.Matrix size (ShapeStatic.ZeroBased offDiag Shape.:+: ()) a ->+   BandedHermitian offDiag size a+fromColumnMajor =+   Array.mapShape+      (\(size, ShapeStatic.ZeroBased k Shape.:+: ()) ->+         MatrixShape.BandedHermitian k ColumnMajor size)+ identity ::    (Shape.C sh, Class.Floating a) => sh -> Diagonal sh a-identity sh =-   Array.mapShape (MatrixShape.BandedHermitian Proxy ColumnMajor) $-   Vector.constant sh one+identity =+   Array.mapShape (MatrixShape.BandedHermitian Proxy ColumnMajor) . Vector.one  diagonal ::    (Shape.C sh, Class.Floating a) => Vector sh (RealOf a) -> Diagonal sh a@@ -101,24 +138,24 @@ takeDiagonal =    runTakeDiagonal $    Class.switchFloating-      (TakeDiagonal $ takeDiagonalAux 1) (TakeDiagonal $ takeDiagonalAux 1)-      (TakeDiagonal $ takeDiagonalAux 2) (TakeDiagonal $ takeDiagonalAux 2)+      (TakeDiagonal takeDiagonalAux) (TakeDiagonal takeDiagonalAux)+      (TakeDiagonal takeDiagonalAux) (TakeDiagonal takeDiagonalAux)  takeDiagonalAux ::    (Unary.Natural offDiag, Shape.C size,     Class.Floating a, RealOf a ~ ar, Class.Real ar) =>-   Int -> BandedHermitian offDiag size a -> Vector size ar-takeDiagonalAux dim (Array (MatrixShape.BandedHermitian numOff order size) a) =+   BandedHermitian offDiag size a -> Vector size ar+takeDiagonalAux (Array (MatrixShape.BandedHermitian numOff order size) a) =    let k = integralFromProxy numOff    in Array.unsafeCreateWithSize size $ \n yPtr -> evalContT $ do          nPtr <- Call.cint n          aPtr <- ContT $ withForeignPtr a          let xPtr =-               castPtr $ advancePtr aPtr $+               realPtr $ advancePtr aPtr $                case order of                   RowMajor -> 0                   ColumnMajor -> k-         incxPtr <- Call.cint (dim * (k+1))+         incxPtr <- Call.cint (caseRealComplexFunc aPtr 1 2 * (k+1))          incyPtr <- Call.cint 1          liftIO $ BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr @@ -210,7 +247,7 @@       (size == sizeX)    let k = integralFromProxy numOff    evalContT $ do-      let conj = transposeOrder transposed order == RowMajor+      let conj = conjugatedOnRowMajor $ transposeOrder transposed order       uploPtr <- Call.char $ uploFromOrder order       nPtr <- Call.cint n       kPtr <- Call.cint k@@ -224,15 +261,15 @@       liftIO $ do          BlasGen.hbmv uploPtr nPtr kPtr             alphaPtr aPtr ldaPtr xPtr incxPtr betaPtr yPtr incyPtr-         when conj $ lacgv nPtr yPtr incyPtr+         condConjugate conj nPtr yPtr incyPtr  -covariance ::+gramian ::    (Shape.C size, Eq size, Class.Floating a,     Unary.Natural sub, Unary.Natural super) =>    Banded.Square sub super size a ->    BandedHermitian (sub :+: super) size a-covariance a =+gramian a =    case mapPair (natFromProxy,natFromProxy) $         MatrixShape.bandedOffDiagonals $ Array.shape a of       (sub,super) ->@@ -259,7 +296,7 @@  multiplyFull ::    (Unary.Natural offDiag, Extent.C vert, Extent.C horiz,-    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>    Transposition -> BandedHermitian offDiag height a ->    Matrix.Full vert horiz height width a ->    Matrix.Full vert horiz height width a@@ -332,7 +369,7 @@  multiplyFullGeneric ::    (Unary.Natural offDiag, Extent.C vert, Extent.C horiz,-    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>    Transposition -> BandedHermitian offDiag height a ->    Matrix.Full vert horiz height width a ->    Matrix.Full vert horiz height width a@@ -342,9 +379,9 @@          case transposed of             Transposed -> (Banded.transpose upper, Banded.transpose lower)             NonTransposed -> (lower,upper)-   in Banded.multiplyFull (Banded.mapExtent Extent.fromSquare lowerT) b-      `Vector.add`-      Banded.multiplyFull (Banded.mapExtent Extent.fromSquare upperT) b+   in VectorPriv.mac one+         (Banded.multiplyFull (Banded.mapExtent Extent.fromSquare lowerT) b)+         (Banded.multiplyFull (Banded.mapExtent Extent.fromSquare upperT) b)  takeUpper ::    (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>@@ -380,9 +417,6 @@  type StaticVector n = Vector (ShapeStatic.ZeroBased n) -{--The list represents ragged rows of a sparse matrix.--} sumRank1 ::    (Unary.Natural k, Shape.Indexed sh, Class.Floating a) =>    Order -> sh ->@@ -430,9 +464,7 @@             let bPtr = advancePtr aPtr (lda*i)             hbr order k alpha                uploPtr mPtr kPtr alphaPtr xPtr incxPtr bPtr incxPtr ldbPtr-      case order of-         RowMajor -> lacgv bSizePtr aPtr incxPtr-         ColumnMajor -> return ()+      condConjugate (conjugatedOnRowMajor order) bSizePtr aPtr incxPtr   type HBR_ ar a =@@ -466,14 +498,14 @@    case order of       ColumnMajor -> do          let aPtr = advancePtr a0Ptr k-         let alphaRealPtr = castPtr alphaPtr+         let alphaRealPtr = realPtr alphaPtr          poke alphaRealPtr alpha          BlasComplex.her uploPtr kPtr alphaRealPtr xPtr incxPtr aPtr ldaPtr          poke alphaPtr . fmap (alpha*) . conjugate =<< peekElemOff xPtr k          BlasGen.axpy nPtr alphaPtr xPtr incxPtr (advancePtr aPtr (k*k)) incaPtr       RowMajor -> do          let aPtr = a0Ptr-         let alphaRealPtr = castPtr alphaPtr+         let alphaRealPtr = realPtr alphaPtr          poke alphaPtr . fmap (alpha*) . conjugate =<< peek xPtr          BlasGen.axpy nPtr alphaPtr xPtr incxPtr aPtr incaPtr          poke alphaRealPtr alpha

@@ -9,8 +9,8 @@ import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape import qualified Numeric.LAPACK.Matrix.Private as Matrix import Numeric.LAPACK.Matrix.Hermitian.Private (TakeDiagonal(..))-import Numeric.LAPACK.Matrix.Shape.Private-         (Order(RowMajor,ColumnMajor), uploFromOrder)+import Numeric.LAPACK.Matrix.Shape.Private (Order(ColumnMajor), uploFromOrder)+import Numeric.LAPACK.Matrix.Modifier (conjugatedOnRowMajor) import Numeric.LAPACK.Vector (Vector) import Numeric.LAPACK.Scalar (RealOf) import Numeric.LAPACK.Private@@ -96,7 +96,7 @@       jobzPtr <- Call.char 'V'       uploPtr <- Call.char $ uploFromOrder order       kPtr <- Call.cint k-      aPtr <- copyCondConjugateToTemp (order==RowMajor) (n*lda) a+      aPtr <- copyCondConjugateToTemp (conjugatedOnRowMajor order) (n*lda) a       ldaPtr <- Call.leadingDim lda       ldzPtr <- Call.leadingDim n       liftIO $ withInfo eigenMsg "hbev" $

@@ -1,5 +1,52 @@ module Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite (-   module Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite.Linear,+   solve,+   solveDecomposed,+   decompose,+   determinant,    ) where -import Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite.Linear+import qualified Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite.Linear+                                                                  as Linear+import qualified Numeric.LAPACK.Matrix.Array.Banded as Banded+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix.Extent as Extent+import Numeric.LAPACK.Matrix.Array (Full)+import Numeric.LAPACK.Scalar (RealOf)++import qualified Numeric.Netlib.Class as Class++import qualified Type.Data.Num.Unary as Unary++import qualified Data.Array.Comfort.Shape as Shape+++solve ::+   (Unary.Natural offDiag, Shape.C size, Eq size,+    Extent.C vert, Extent.C horiz, Shape.C nrhs, Class.Floating a) =>+   Banded.Hermitian offDiag size a ->+   Full vert horiz size nrhs a -> Full vert horiz size nrhs a+solve = ArrMatrix.lift2 Linear.solve++{- |+> solve a b == solveDecomposed (decompose a) b+> solve (gramian u) b == solveDecomposed u b+-}+solveDecomposed ::+   (Unary.Natural offDiag, Shape.C size, Eq size,+    Extent.C vert, Extent.C horiz, Shape.C nrhs, Class.Floating a) =>+   Banded.Upper offDiag size a ->+   Full vert horiz size nrhs a -> Full vert horiz size nrhs a+solveDecomposed = ArrMatrix.lift2 Linear.solveDecomposed++{- |+Cholesky decomposition+-}+decompose ::+   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>+   Banded.Hermitian offDiag size a -> Banded.Upper offDiag size a+decompose = ArrMatrix.lift1 Linear.decompose++determinant ::+   (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>+   Banded.Hermitian offDiag size a -> RealOf a+determinant = Linear.determinant . ArrMatrix.toVector

@@ -14,7 +14,8 @@ import Numeric.LAPACK.Matrix.Hermitian.Private (Determinant(..)) import Numeric.LAPACK.Matrix.Triangular.Private (copyTriangleToTemp) import Numeric.LAPACK.Matrix.Shape.Private (uploFromOrder)-import Numeric.LAPACK.Matrix.Private (Full, Conjugation(Conjugated))+import Numeric.LAPACK.Matrix.Modifier (Conjugation(Conjugated))+import Numeric.LAPACK.Matrix.Private (Full) import Numeric.LAPACK.Scalar (RealOf, realPart) import Numeric.LAPACK.Private (copyBlock, withInfo, rankMsg, definiteMsg) @@ -53,10 +54,6 @@          withInfo definiteMsg "pbsv" $             LapackGen.pbsv uploPtr nPtr kPtr nrhsPtr aPtr ldaPtr xPtr ldxPtr -{- |-> solve a b == solveDecomposed (decompose a) b-> solve (covariance u) b == solveDecomposed u b--} solveDecomposed ::    (Unary.Natural offDiag, Shape.C size, Eq size,     Extent.C vert, Extent.C horiz, Shape.C nrhs, Class.Floating a) =>@@ -76,9 +73,6 @@             LapackGen.pbtrs uploPtr nPtr kPtr nrhsPtr aPtr ldaPtr xPtr ldxPtr  -{- |-Cholesky decomposition--} decompose ::    (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>    BandedHermitian offDiag size a -> Banded.Upper offDiag size a

@@ -4,12 +4,19 @@  import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent+import qualified Numeric.LAPACK.Matrix.RowMajor as RowMajor+import qualified Numeric.LAPACK.Vector as Vector import qualified Numeric.LAPACK.Private as Private-import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor, ColumnMajor))-import Numeric.LAPACK.Matrix.Private (Full, General)+import Numeric.LAPACK.Matrix.Shape.Private+         (Order(RowMajor, ColumnMajor), transposeFromOrder, flipOrder)+import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated))+import Numeric.LAPACK.Matrix.Private+         (Full, Tall, Wide, General, ShapeInt, revealOrder) import Numeric.LAPACK.Vector (Vector)-import Numeric.LAPACK.Scalar (zero, one)-import Numeric.LAPACK.Private (pointerSeq)+import Numeric.LAPACK.Scalar (RealOf, zero, one)+import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))+import Numeric.LAPACK.Private+         (pointerSeq, copyTransposed, copySubMatrix, copyBlock)  import qualified Numeric.BLAS.FFI.Generic as BlasGen import qualified Numeric.Netlib.Utility as Call@@ -18,19 +25,49 @@ import qualified Data.Array.Comfort.Storable.Unchecked as Array import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable.Unchecked (Array(Array))+import Data.Array.Comfort.Shape ((:+:)((:+:))) -import Foreign.ForeignPtr (withForeignPtr)-import Foreign.Storable (poke, peek)+import Foreign.Marshal.Array (advancePtr)+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)  import Control.Monad.Trans.Cont (ContT(ContT), evalContT) import Control.Monad.IO.Class (liftIO) +import Data.Complex (Complex) ++caseTallWide ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>+   Full vert horiz height width a ->+   Either (Tall height width a) (Wide height width a)+caseTallWide (Array shape a) =+   either (Left . flip Array a) (Right . flip Array a) $+   MatrixShape.caseTallWide shape++ transpose ::    (Extent.C vert, Extent.C horiz) =>    Full vert horiz height width a -> Full horiz vert width height a transpose = Array.mapShape MatrixShape.transpose +adjoint ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,+    Class.Floating a) =>+   Full vert horiz height width a -> Full horiz vert width height a+adjoint = transpose . Vector.conjugate+++swapMultiply ::+   (Extent.C vertA, Extent.C vertB, Extent.C horizA, Extent.C horizB) =>+   (matrix ->+    Full horizA vertA widthA heightA a ->+    Full horizB vertB widthB heightB a) ->+   Full vertA horizA heightA widthA a ->+   matrix ->+   Full vertB horizB heightB widthB a+swapMultiply multiplyTrans a b = transpose $ multiplyTrans b $ transpose a++ mapHeight ::    (Extent.GeneralTallWide vert horiz,     Extent.GeneralTallWide horiz vert) =>@@ -53,7 +90,25 @@       (\(MatrixShape.Full order extent) ->          MatrixShape.Full order $ Extent.mapWidth f extent) +uncheck ::+   (Extent.C vert, Extent.C horiz) =>+   Full vert horiz height width a ->+   Full vert horiz (Unchecked height) (Unchecked width) a+uncheck =+   Array.mapShape+      (\(MatrixShape.Full order extent) ->+         MatrixShape.Full order $ Extent.mapWrap Unchecked Unchecked extent) +recheck ::+   (Extent.C vert, Extent.C horiz) =>+   Full vert horiz (Unchecked height) (Unchecked width) a ->+   Full vert horiz height width a+recheck =+   Array.mapShape+      (\(MatrixShape.Full order extent) ->+         MatrixShape.Full order $ Extent.recheck extent)++ singleRow :: Order -> Vector width a -> General () width a singleRow order = Array.mapShape (MatrixShape.general order ()) @@ -119,58 +174,200 @@       ColumnMajor -> transpose . forceRowMajor . transpose  +takeSub ::+   (Extent.C vert, Extent.C horiz,+    Shape.C heightA, Shape.C height, Shape.C width, Class.Floating a) =>+   heightA -> Int -> ForeignPtr a ->+   MatrixShape.Full vert horiz height width ->+   Full vert horiz height width a+takeSub heightA k a shape@(MatrixShape.Full order extentB) =+   Array.unsafeCreateWithSize shape $ \blockSize bPtr ->+   withForeignPtr a $ \aPtr ->+   let ma = Shape.size heightA+       mb = Shape.size $ Extent.height extentB+       n  = Shape.size $ Extent.width  extentB+   in case order of+         RowMajor -> copyBlock blockSize (advancePtr aPtr (k*n)) bPtr+         ColumnMajor -> copySubMatrix mb n ma (advancePtr aPtr k) mb bPtr++takeTop ::+   (Extent.C vert, Shape.C height0, Shape.C height1, Shape.C width,+    Class.Floating a) =>+   Full vert Extent.Big (height0:+:height1) width a ->+   Full vert Extent.Big height0 width a+takeTop (Array (MatrixShape.Full order extentA) a) =+   let heightA@(heightB:+:_) = Extent.height extentA+       extentB = Extent.reduceWideHeight heightB extentA+   in takeSub heightA 0 a $ MatrixShape.Full order extentB++takeBottom ::+   (Extent.C vert, Shape.C height0, Shape.C height1, Shape.C width,+    Class.Floating a) =>+   Full vert Extent.Big (height0:+:height1) width a ->+   Full vert Extent.Big height1 width a+takeBottom (Array (MatrixShape.Full order extentA) a) =+   let heightA@(height0:+:heightB) = Extent.height extentA+       extentB = Extent.reduceWideHeight heightB extentA+   in takeSub heightA (Shape.size height0) a $ MatrixShape.Full order extentB++takeLeft ::+   (Extent.C vert, Shape.C height, Shape.C width0, Shape.C width1,+    Class.Floating a) =>+   Full Extent.Big vert height (width0:+:width1) a ->+   Full Extent.Big vert height width0 a+takeLeft = transpose . takeTop . transpose++takeRight ::+   (Extent.C vert, Shape.C height, Shape.C width0, Shape.C width1,+    Class.Floating a) =>+   Full Extent.Big vert height (width0:+:width1) a ->+   Full Extent.Big vert height width1 a+takeRight = transpose . takeBottom . transpose+++splitRows ::+   (Extent.C vert, Shape.C width, Class.Floating a) =>+   Int ->+   Full vert Extent.Big ShapeInt width a ->+   Full vert Extent.Big (ShapeInt:+:ShapeInt) width a+splitRows k =+   Array.mapShape+      (\(MatrixShape.Full order extent) ->+         MatrixShape.Full order $+         Extent.reduceWideHeight+            (Shape.zeroBasedSplit k $ Extent.height extent)+            extent)++takeRows, dropRows ::+   (Extent.C vert, Shape.C width, Class.Floating a) =>+   Int ->+   Full vert Extent.Big ShapeInt width a ->+   Full vert Extent.Big ShapeInt width a+takeRows k = takeTop . splitRows k+dropRows k = takeBottom . splitRows k++takeColumns, dropColumns ::+   (Extent.C horiz, Shape.C height, Class.Floating a) =>+   Int ->+   Full Extent.Big horiz height ShapeInt a ->+   Full Extent.Big horiz height ShapeInt a+takeColumns k = transpose . takeRows k . transpose+dropColumns k = transpose . dropRows k . transpose+++data OrderBias = LeftBias | RightBias | ContiguousBias+   deriving (Eq, Ord, Enum, Show)++beside ::+   (Extent.C vertA, Extent.C vertB, Extent.C vertC,+    Shape.C height, Eq height, Shape.C widthA, Shape.C widthB,+    Class.Floating a) =>+   OrderBias ->+   Extent.AppendMode vertA vertB vertC height widthA widthB ->+   Full vertA Extent.Big height widthA a ->+   Full vertB Extent.Big height widthB a ->+   Full vertC Extent.Big height (widthA:+:widthB) a+beside orderBias (Extent.AppendMode appendMode)+      (Array (MatrixShape.Full orderA extentA) a)+      (Array (MatrixShape.Full orderB extentB) b) =+   let (heightA,widthA) = Extent.dimensions extentA+       (heightB,widthB) = Extent.dimensions extentB+       n = Shape.size heightA+       ma = Shape.size widthA; volA = n*ma+       mb = Shape.size widthB; volB = n*mb+       m = ma+mb+       create order act =+          Array.unsafeCreate+             (MatrixShape.Full order $ appendMode extentA extentB) $ \cPtr ->+          withForeignPtr a $ \aPtr ->+          withForeignPtr b $ \bPtr ->+          act aPtr bPtr cPtr $ advancePtr cPtr $+          case order of+             RowMajor -> ma+             ColumnMajor -> volA+   in+    if heightA /= heightB+      then error "beside: mismatching heights"+      else+         case (orderA,orderB) of+            (RowMajor,RowMajor) ->+               create RowMajor $ \aPtr bPtr cPtr _ -> evalContT $ do+                  maPtr <- Call.cint ma+                  mbPtr <- Call.cint mb+                  incxPtr <- Call.cint 1+                  incyPtr <- Call.cint 1+                  liftIO $+                     sequence_ $ take n $+                     zipWith3+                        (\akPtr bkPtr ckPtr -> do+                           BlasGen.copy maPtr akPtr incxPtr ckPtr incyPtr+                           BlasGen.copy mbPtr bkPtr incxPtr+                              (ckPtr `advancePtr` ma) incyPtr)+                        (pointerSeq ma aPtr)+                        (pointerSeq mb bPtr)+                        (pointerSeq m cPtr)+            (RowMajor,ColumnMajor) ->+               case orderBias of+                  LeftBias ->+                     create RowMajor $ \aPtr bPtr clPtr crPtr -> do+                        copySubMatrix ma n ma aPtr m clPtr+                        copyTransposed mb n bPtr m crPtr+                  _ ->+                     create ColumnMajor $ \aPtr bPtr clPtr crPtr -> do+                        copyTransposed n ma aPtr n clPtr+                        copyBlock volB bPtr crPtr+            (ColumnMajor,RowMajor) ->+               case orderBias of+                  RightBias ->+                     create RowMajor $ \aPtr bPtr clPtr crPtr -> do+                        copyTransposed ma n aPtr m clPtr+                        copySubMatrix mb n mb bPtr m crPtr+                  _ ->+                     create ColumnMajor $ \aPtr bPtr clPtr crPtr -> do+                        copyBlock volA aPtr clPtr+                        copyTransposed n mb bPtr n crPtr+            (ColumnMajor,ColumnMajor) ->+               create ColumnMajor $ \aPtr bPtr clPtr crPtr -> evalContT $ do+                  naPtr <- Call.cint volA+                  nbPtr <- Call.cint volB+                  incxPtr <- Call.cint 1+                  incyPtr <- Call.cint 1+                  liftIO $ do+                     BlasGen.copy naPtr aPtr incxPtr clPtr incyPtr+                     BlasGen.copy nbPtr bPtr incxPtr crPtr incyPtr++above ::+   (Extent.C horizA, Extent.C horizB, Extent.C horizC,+    Shape.C width, Eq width, Shape.C heightA, Shape.C heightB,+    Class.Floating a) =>+   OrderBias ->+   Extent.AppendMode horizA horizB horizC width heightA heightB ->+   Full Extent.Big horizA heightA width a ->+   Full Extent.Big horizB heightB width a ->+   Full Extent.Big horizC (heightA:+:heightB) width a+above orderBias appendMode a b =+   transpose $ beside orderBias appendMode (transpose a) (transpose b)+++liftRowMajor ::+   (Extent.C vert, Extent.C horiz) =>+   (Array (height, width) a -> Array (height, width) b) ->+   (Array (width, height) a -> Array (width, height) b) ->+   Full vert horiz height width a ->+   Full vert horiz height width b+liftRowMajor fr fc a =+   either+      (Array.reshape (Array.shape a) . fr)+      (Array.reshape (Array.shape a) . fc) $+   revealOrder a+ scaleRows ::    (Extent.C vert, Extent.C horiz,     Shape.C height, Eq height, Shape.C width, Class.Floating a) =>    Vector height a ->    Full vert horiz height width a ->    Full vert horiz height width a-scaleRows-   (Array heightX x) (Array shape@(MatrixShape.Full order extent) a) =-      Array.unsafeCreate shape $ \bPtr -> do-   let (height,width) = Extent.dimensions extent-   Call.assert "scaleRows: sizes mismatch" (heightX == height)-   case order of-      RowMajor -> evalContT $ do-         let m = Shape.size height-         let n = Shape.size width-         alphaPtr <- Call.alloca-         nPtr <- Call.cint n-         xPtr <- ContT $ withForeignPtr x-         aPtr <- ContT $ withForeignPtr a-         incaPtr <- Call.cint 1-         incbPtr <- Call.cint 1-         liftIO $ sequence_ $ take m $-            zipWith3-               (\xkPtr akPtr bkPtr -> do-                  poke alphaPtr =<< peek xkPtr-                  BlasGen.copy nPtr akPtr incaPtr bkPtr incbPtr-                  BlasGen.scal nPtr alphaPtr bkPtr incbPtr)-               (pointerSeq 1 xPtr)-               (pointerSeq n aPtr)-               (pointerSeq n bPtr)-      ColumnMajor -> evalContT $ do-         let m = Shape.size width-         let n = Shape.size height-         transPtr <- Call.char 'N'-         nPtr <- Call.cint n-         klPtr <- Call.cint 0-         kuPtr <- Call.cint 0-         alphaPtr <- Call.number one-         xPtr <- ContT $ withForeignPtr x-         ldxPtr <- Call.leadingDim 1-         aPtr <- ContT $ withForeignPtr a-         incaPtr <- Call.cint 1-         betaPtr <- Call.number zero-         incbPtr <- Call.cint 1-         liftIO $ sequence_ $ take m $-            zipWith-               (\akPtr bkPtr ->-                  Private.gbmv transPtr-                     nPtr nPtr klPtr kuPtr alphaPtr xPtr ldxPtr-                     akPtr incaPtr betaPtr bkPtr incbPtr)-               (pointerSeq n aPtr)-               (pointerSeq n bPtr)+scaleRows x = liftRowMajor (RowMajor.scaleRows x) (RowMajor.scaleColumns x)  scaleColumns ::    (Extent.C vert, Extent.C horiz,@@ -179,3 +376,145 @@    Full vert horiz height width a ->    Full vert horiz height width a scaleColumns x = transpose . scaleRows x . transpose+++scaleRowsComplex ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Class.Real a) =>+   Vector height a ->+   Full vert horiz height width (Complex a) ->+   Full vert horiz height width (Complex a)+scaleRowsComplex x =+   liftRowMajor+      (RowMajor.recomplex . RowMajor.scaleRows x . RowMajor.decomplex)+      (RowMajor.recomplex .+       RowMajor.scaleColumns+         (RowMajor.tensorProduct (Left NonConjugated) x+            (Vector.one Shape.Enumeration)) .+       RowMajor.decomplex)++scaleColumnsComplex ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Eq width, Class.Real a) =>+   Vector width a ->+   Full vert horiz height width (Complex a) ->+   Full vert horiz height width (Complex a)+scaleColumnsComplex x = transpose . scaleRowsComplex x . transpose+++scaleRowsReal ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Eq height, Shape.C width,+    Class.Floating a) =>+   Vector height (RealOf a) ->+   Full vert horiz height width a ->+   Full vert horiz height width a+scaleRowsReal =+   getScaleRowsReal $+   Class.switchFloating+      (ScaleRowsReal scaleRows)+      (ScaleRowsReal scaleRows)+      (ScaleRowsReal scaleRowsComplex)+      (ScaleRowsReal scaleRowsComplex)++newtype ScaleRowsReal f g a =+   ScaleRowsReal {getScaleRowsReal :: f (RealOf a) -> g a -> g a}++scaleColumnsReal ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   Vector width (RealOf a) ->+   Full vert horiz height width a ->+   Full vert horiz height width a+scaleColumnsReal x = transpose . scaleRowsReal x . transpose++++multiplyVector ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width, Eq width,+    Class.Floating a) =>+   Full vert horiz height width a -> Vector width a -> Vector height a+multiplyVector a x =+   let width = MatrixShape.fullWidth $ Array.shape a+   in if width == Array.shape x+         then multiplyVectorUnchecked a x+         else error "multiplyVector: width shapes mismatch"++multiplyVectorUnchecked ::+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,+    Class.Floating a) =>+   Full vert horiz height width a -> Vector width a -> Vector height a+multiplyVectorUnchecked+   (Array shape@(MatrixShape.Full order extent) a) (Array _ x) =+      Array.unsafeCreate (Extent.height extent) $ \yPtr -> do+   let (m,n) = MatrixShape.dimensions shape+   let lda = m+   evalContT $ do+      transPtr <- Call.char $ transposeFromOrder order+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      alphaPtr <- Call.number one+      aPtr <- ContT $ withForeignPtr a+      ldaPtr <- Call.leadingDim lda+      xPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint 1+      betaPtr <- Call.number zero+      incyPtr <- Call.cint 1+      liftIO $+         Private.gemv+            transPtr mPtr nPtr alphaPtr aPtr ldaPtr+            xPtr incxPtr betaPtr yPtr incyPtr++{- |+Multiply two matrices with the same dimension constraints.+E.g. you can multiply 'General' and 'General' matrices,+or 'Square' and 'Square' matrices.+It may seem to be overly strict in this respect,+but that design supports type inference the best.+You can lift the restrictions by generalizing operands+with 'Square.toFull', 'Matrix.fromFull',+'Matrix.generalizeTall' or 'Matrix.generalizeWide'.+-}+multiply, multiplyColumnMajor ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height,+    Shape.C fuse, Eq fuse,+    Shape.C width,+    Class.Floating a) =>+   Full vert horiz height fuse a ->+   Full vert horiz fuse width a ->+   Full vert horiz height width a+-- preserve order of the right factor+multiply+   (Array (MatrixShape.Full orderA extentA) a)+   (Array (MatrixShape.Full orderB extentB) b) =+   case Extent.fuse extentA extentB of+      Nothing -> error "multiply: fuse shapes mismatch"+      Just extent ->+         Array.unsafeCreate (MatrixShape.Full orderB extent) $ \cPtr -> do++      let (height,fuse) = Extent.dimensions extentA+      let width = Extent.width extentB+      let m = Shape.size height+      let n = Shape.size width+      let k = Shape.size fuse+      case orderB of+         RowMajor ->+            Private.multiplyMatrix (flipOrder orderB) (flipOrder orderA)+               n k m b a cPtr+         ColumnMajor -> Private.multiplyMatrix orderA orderB m k n a b cPtr++-- always return ColumnMajor+multiplyColumnMajor+   (Array (MatrixShape.Full orderA extentA) a)+   (Array (MatrixShape.Full orderB extentB) b) =+   case Extent.fuse extentA extentB of+      Nothing -> error "multiply: fuse shapes mismatch"+      Just extent ->+         Array.unsafeCreate (MatrixShape.Full ColumnMajor extent) $ \cPtr -> do++      let (height,fuse) = Extent.dimensions extentA+      let width = Extent.width extentB+      let m = Shape.size height+      let n = Shape.size width+      let k = Shape.size fuse+      Private.multiplyMatrix orderA orderB m k n a b cPtr

@@ -0,0 +1,107 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Matrix.Class (+   SquareShape(toSquare, identityOrder, takeDiagonal),+   identityFrom,+   identityFromHeight,+   identityFromWidth,+   trace,+   Complex(conjugate, fromReal, toComplex),+   ) where++import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix.Type as Type+import qualified Numeric.LAPACK.Matrix.Plain.Class as Plain+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape+import qualified Numeric.LAPACK.Matrix.Permutation as Permutation+import qualified Numeric.LAPACK.Permutation.Private as Perm+import qualified Numeric.LAPACK.Vector as Vector+import qualified Numeric.LAPACK.Scalar as Scalar+import Numeric.LAPACK.Matrix.Array (ArrayMatrix)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Scalar (RealOf, ComplexOf)++import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Shape as Shape+++class Complex typ where+   conjugate ::+      (Class.Floating a) => Type.Matrix typ a -> Type.Matrix typ a+   fromReal ::+      (Class.Floating a) => Type.Matrix typ (RealOf a) -> Type.Matrix typ a+   toComplex ::+      (Class.Floating a) => Type.Matrix typ a -> Type.Matrix typ (ComplexOf a)++instance (Plain.Complex sh) => Complex (ArrMatrix.Array sh) where+   conjugate = ArrMatrix.lift1 Plain.conjugate+   fromReal  = ArrMatrix.lift1 Plain.fromReal+   toComplex = ArrMatrix.lift1 Plain.toComplex++instance (Shape.C shape) => Complex (Type.Scale shape) where+   conjugate (Type.Scale sh m) = Type.Scale sh $ Scalar.conjugate m+   fromReal (Type.Scale sh m) = Type.Scale sh $ Scalar.fromReal m+   toComplex (Type.Scale sh m) = Type.Scale sh $ Scalar.toComplex m++instance (Shape.C shape) => Complex (Perm.Permutation shape) where+   conjugate = id+   fromReal (Type.Permutation p) = Type.Permutation p+   toComplex (Type.Permutation p) = Type.Permutation p+++class SquareShape typ where+   toSquare ::+      (Type.HeightOf typ ~ sh, Class.Floating a) =>+      Type.Matrix typ a -> ArrMatrix.Square sh a+   identityOrder ::+      (Type.HeightOf typ ~ sh, Class.Floating a) =>+      MatrixShape.Order -> sh -> Type.Matrix typ a+   takeDiagonal ::+      (Type.HeightOf typ ~ sh, Class.Floating a) =>+      Type.Matrix typ a -> Vector sh a++instance (ArrMatrix.SquareShape sh) => SquareShape (ArrMatrix.Array sh) where+   toSquare = ArrMatrix.lift1 Plain.toSquare+   identityOrder order = ArrMatrix.lift0 . Plain.identityOrder order+   takeDiagonal = Plain.takeDiagonal . ArrMatrix.toVector++instance (Shape.C sh) => SquareShape (Type.Scale sh) where+   toSquare (Type.Scale sh a) =+      Triangular.toSquare $ Triangular.asDiagonal $+      Triangular.diagonal MatrixShape.RowMajor $ Vector.constant sh a+   identityOrder _ sh = Type.Scale sh Scalar.one+   takeDiagonal (Type.Scale sh a) = Vector.constant sh a++instance (Shape.C sh) => SquareShape (Perm.Permutation sh) where+   toSquare = Permutation.toMatrix+   identityOrder _ = Permutation.identity+   takeDiagonal = Perm.takeDiagonal . Permutation.toPermutation+++identityFrom ::+   (Plain.SquareShape shape, ArrMatrix.ShapeOrder shape, Class.Floating a) =>+   ArrayMatrix shape a -> ArrayMatrix shape a+identityFrom m =+   identityOrder (ArrMatrix.shapeOrder $ ArrMatrix.shape m) (Type.height m)++identityFromHeight ::+   (ArrMatrix.ShapeOrder shape, MatrixShape.Box shape,+    MatrixShape.HeightOf shape ~ Type.HeightOf typ, SquareShape typ,+    Class.Floating a) =>+   ArrayMatrix shape a -> Type.Matrix typ a+identityFromHeight m =+   identityOrder (ArrMatrix.shapeOrder $ ArrMatrix.shape m) (Type.height m)++identityFromWidth ::+   (ArrMatrix.ShapeOrder shape, MatrixShape.Box shape,+    MatrixShape.WidthOf shape ~ Type.HeightOf typ, SquareShape typ,+    Class.Floating a) =>+   ArrayMatrix shape a -> Type.Matrix typ a+identityFromWidth m =+   identityOrder (ArrMatrix.shapeOrder $ ArrMatrix.shape m) (Type.width m)++trace ::+   (SquareShape typ, Type.HeightOf typ ~ sh, Shape.C sh, Class.Floating a) =>+   Type.Matrix typ a -> a+trace = Vector.sum . takeDiagonal

@@ -1,97 +1,131 @@ {-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE ConstraintKinds #-}-{-# LANGUAGE UndecidableInstances #-} module Numeric.LAPACK.Matrix.Divide where -import qualified Numeric.LAPACK.Matrix.Square.Linear-                                           as Square-import qualified Numeric.LAPACK.Matrix.Triangular.Linear-                                           as Triangular-import qualified Numeric.LAPACK.Matrix.Hermitian.Linear-                                           as Hermitian-import qualified Numeric.LAPACK.Matrix.Banded.Linear-                                           as Banded-import qualified Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite.Linear-                                           as BandedHermitianPositiveDefinite--import qualified Numeric.LAPACK.Matrix.Basic as Basic+import qualified Numeric.LAPACK.Matrix.Plain.Divide as ArrDivide+import qualified Numeric.LAPACK.Matrix.Array.Basic as Basic+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix.Permutation as PermMatrix+import qualified Numeric.LAPACK.Matrix.Multiply as Multiply+import qualified Numeric.LAPACK.Matrix.Type as Type import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape-import qualified Numeric.LAPACK.Matrix.Shape.Box as Box import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent-import Numeric.LAPACK.Matrix.Extent.Private (Small)-import Numeric.LAPACK.Matrix.Private (Full)+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix.Array (Full)+import Numeric.LAPACK.Matrix.Type (Matrix, scaleWithCheck)+import Numeric.LAPACK.Matrix.Modifier+         (Transposition(NonTransposed,Transposed),+          Inversion(Inverted)) import Numeric.LAPACK.Vector (Vector)  import qualified Numeric.Netlib.Class as Class -import qualified Type.Data.Num.Unary as Unary- import qualified Data.Array.Comfort.Shape as Shape-import Data.Array.Comfort.Storable (Array) +import Data.Semigroup ((<>)) -class (Shape.C shape) => Solve shape where++class+   (Type.Box typ, Type.HeightOf typ ~ Type.WidthOf typ) =>+      Determinant typ where+   determinant :: (Class.Floating a) => Matrix typ a -> a++class (Type.Box typ, Type.HeightOf typ ~ Type.WidthOf typ) => Solve typ where+   {-# MINIMAL solve | solveLeft,solveRight #-}    solve ::-      (Class.Floating a, Box.HeightOf shape ~ height, Eq height,-       Extent.C horiz, Extent.C vert, Shape.C nrhs) =>-      Array shape a ->-      Full vert horiz height nrhs a -> Full vert horiz height nrhs a+      (Type.HeightOf typ ~ height, Eq height, Shape.C width,+       Extent.C horiz, Extent.C vert, Class.Floating a) =>+      Transposition -> Matrix typ a ->+      Full vert horiz height width a -> Full vert horiz height width a+   solve NonTransposed a b = solveRight a b+   solve Transposed a b = Basic.transpose $ solveLeft (Basic.transpose b) a -class (Solve shape) => Inverse shape where-   inverse :: (Class.Floating a) => Array shape a -> Array shape a+   solveRight ::+      (Type.HeightOf typ ~ height, Eq height, Shape.C width,+       Extent.C horiz, Extent.C vert, Class.Floating a) =>+      Matrix typ a ->+      Full vert horiz height width a -> Full vert horiz height width a+   solveRight = solve NonTransposed +   solveLeft ::+      (Type.WidthOf typ ~ width, Eq width, Shape.C height,+       Extent.C horiz, Extent.C vert, Class.Floating a) =>+      Full vert horiz height width a ->+      Matrix typ a -> Full vert horiz height width a+   solveLeft = Basic.swapMultiply $ solve Transposed++class (Solve typ, Multiply.Power typ) => Inverse typ where+   inverse :: (Class.Floating a) => Matrix typ a -> Matrix typ a++infixl 7 ##/#+infixr 7 #\##++(#\##) ::+   (Solve typ, Type.HeightOf typ ~ height, Eq height, Shape.C width,+    Extent.C horiz, Extent.C vert, Class.Floating a) =>+   Matrix typ a ->+   Full vert horiz height width a -> Full vert horiz height width a+(#\##) = solveRight++(##/#) ::+   (Solve typ, Type.WidthOf typ ~ width, Eq width, Shape.C height,+    Extent.C horiz, Extent.C vert, Class.Floating a) =>+   Full vert horiz height width a ->+   Matrix typ a -> Full vert horiz height width a+(##/#) = solveLeft++ solveVector ::-   (Solve shape, Box.HeightOf shape ~ height, Eq height, Class.Floating a) =>-   Array shape a -> Vector height a -> Vector height a-solveVector m =-   Basic.flattenColumn . solve m . Basic.singleColumn MatrixShape.ColumnMajor+   (Solve typ, Type.HeightOf typ ~ height, Eq height, Class.Floating a) =>+   Transposition -> Matrix typ a -> Vector height a -> Vector height a+solveVector trans =+   ArrMatrix.unliftColumn MatrixShape.ColumnMajor . solve trans +infixl 7 -/#+infixr 7 #\| -instance-   (vert ~ Small, horiz ~ Small,-    Shape.C width, Shape.C height, height ~ width) =>-      Solve (MatrixShape.Full vert horiz height width) where-   solve = Square.solve+(#\|) ::+   (Solve typ, Type.HeightOf typ ~ height, Eq height, Class.Floating a) =>+   Matrix typ a -> Vector height a -> Vector height a+(#\|) = solveVector NonTransposed -instance-   (vert ~ Small, horiz ~ Small,-    Shape.C width, Shape.C height, height ~ width) =>-      Inverse (MatrixShape.Full vert horiz height width) where-   inverse = Square.inverse+(-/#) ::+   (Solve typ, Type.HeightOf typ ~ height, Eq height, Class.Floating a) =>+   Vector height a -> Matrix typ a -> Vector height a+(-/#) = flip $ solveVector Transposed  -instance (Shape.C shape) => Solve (MatrixShape.Hermitian shape) where-   solve = Hermitian.solve+instance (Shape.C shape, Eq shape) => Determinant (Type.Scale shape) where+   determinant (Type.Scale sh a) = a ^ Shape.size sh -instance (Shape.C shape) => Inverse (MatrixShape.Hermitian shape) where-   inverse = Hermitian.inverse+instance (Shape.C shape, Eq shape) => Solve (Type.Scale shape) where+   solve _trans =+      scaleWithCheck "Matrix.Scale.solve" Type.height $+         ArrMatrix.lift1 . Vector.scale . recip +instance (Shape.C shape, Eq shape) => Inverse (Type.Scale shape) where+   inverse (Type.Scale shape a) = Type.Scale shape $ recip a -instance-   (MatrixShape.Content lo, MatrixShape.Content up,-    MatrixShape.TriDiag diag, Shape.C shape) =>-      Solve (MatrixShape.Triangular lo diag up shape) where-   solve = Triangular.solve -instance-   (MatrixShape.DiagUpLo lo up,-    MatrixShape.TriDiag diag, Shape.C shape) =>-      Inverse (MatrixShape.Triangular lo diag up shape) where-   inverse = Triangular.inverse+instance (Shape.C shape) => Determinant (PermMatrix.Permutation shape) where+   determinant = PermMatrix.determinant +instance (Shape.C shape) => Solve (PermMatrix.Permutation shape) where+   solve trans =+      PermMatrix.multiplyFull+         (Inverted <> PermMatrix.inversionFromTransposition trans) -instance-   (Unary.Natural sub, Unary.Natural super, vert ~ Small, horiz ~ Small,-    Shape.C width, Shape.C height, width ~ height) =>-      Solve (MatrixShape.Banded sub super vert horiz height width) where-   solve = Banded.solve+instance (Shape.C shape) => Inverse (PermMatrix.Permutation shape) where+   inverse = PermMatrix.transpose  -{- |-There is no solver for indefinite matrices.-Thus the instance will fail for indefinite but solvable systems.--} instance-   (Unary.Natural offDiag, Shape.C size) =>-      Solve (MatrixShape.BandedHermitian offDiag size) where-   solve = BandedHermitianPositiveDefinite.solve+      (ArrDivide.Determinant shape) => Determinant (ArrMatrix.Array shape) where+   determinant = ArrDivide.determinant . ArrMatrix.toVector++instance (ArrDivide.Solve shape) => Solve (ArrMatrix.Array shape) where+   solve = ArrMatrix.lift2 . ArrDivide.solve+   solveLeft = ArrMatrix.lift2 ArrDivide.solveLeft+   solveRight = ArrMatrix.lift2 ArrDivide.solveRight++instance (ArrDivide.Inverse shape) => Inverse (ArrMatrix.Array shape) where+   inverse = ArrMatrix.lift1 ArrDivide.inverse

@@ -17,14 +17,21 @@    fromSquareLiberal,    generalizeTall,    generalizeWide,+   Extent.GeneralTallWide,++   Extent.AppendMode,+   Extent.appendSame,+   Extent.appendLeft,+   Extent.appendRight,+   Extent.Append,+   Extent.appendAny,    ) where  import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent import Numeric.LAPACK.Matrix.Extent.Private (C, Small, Big, Map(Map))  -toGeneral ::-   (C vert, C horiz) => Map vert horiz Big Big height width+toGeneral :: (C vert, C horiz) => Map vert horiz Big Big height width toGeneral = Map Extent.toGeneral  fromSquare :: (C vert, C horiz) => Map Small Small vert horiz size size

@@ -1,11 +1,15 @@ {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE GADTs #-} module Numeric.LAPACK.Matrix.Extent.Private where  import qualified Numeric.LAPACK.Matrix.Extent.Kind as EK+import Numeric.LAPACK.Shape.Private (Unchecked(deconsUnchecked)) import Numeric.LAPACK.Wrapper (Flip(Flip, getFlip)) +import Data.Array.Comfort.Shape ((:+:)((:+:)))+ import Control.DeepSeq (NFData, rnf)  import Data.Maybe.HT (toMaybe)@@ -382,8 +386,60 @@       (Adapt $ \(Extent o (EK.Tall h w)) -> Extent o (EK.Tall h (f w)))       (Adapt $ \(Extent o (EK.General h w)) -> Extent o (EK.General h (f w))) +mapSquareSize :: (shA -> shB) -> Square shA -> Square shB+mapSquareSize f (Extent o (EK.Square s)) = Extent o (EK.Square (f s))  +mapWrap ::+   (C vert, C horiz) =>+   (height -> f height) ->+   (width -> f width) ->+   Extent vert horiz height width ->+   Extent vert horiz (f height) (f width)+mapWrap fh fw =+   getAdapt $+   switchTagPair+      (Adapt $ \(Extent o (EK.Square h)) -> Extent o (EK.Square (fh h)))+      (Adapt $ \(Extent o (EK.Wide h w)) -> Extent o (EK.Wide (fh h) (fw w)))+      (Adapt $ \(Extent o (EK.Tall h w)) -> Extent o (EK.Tall (fh h) (fw w)))+      (Adapt $ \(Extent o (EK.General h w)) ->+                  Extent o (EK.General (fh h) (fw w)))++{- only admissible since GHC-7.8+mapUnwrap ::+   (C vert, C horiz) =>+   (f height -> height) ->+   (f width -> width) ->+   Extent vert horiz (f height) (f width) ->+   Extent vert horiz height width+mapUnwrap fh fw =+   getAdapt $+   switchTagPair+      (Adapt $ \(Extent o (EK.Square h)) -> Extent o (EK.Square (fh h)))+      (Adapt $ \(Extent o (EK.Wide h w)) -> Extent o (EK.Wide (fh h) (fw w)))+      (Adapt $ \(Extent o (EK.Tall h w)) -> Extent o (EK.Tall (fh h) (fw w)))+      (Adapt $ \(Extent o (EK.General h w)) ->+                  Extent o (EK.General (fh h) (fw w)))+-}++recheck ::+   (C vert, C horiz) =>+   Extent vert horiz (Unchecked height) (Unchecked width) ->+   Extent vert horiz height width+recheck =+   getAdapt $+   switchTagPair+      (Adapt $ \(Extent o (EK.Square h)) ->+         Extent o (EK.Square (deconsUnchecked h)))+      (Adapt $ \(Extent o (EK.Wide h w)) ->+         Extent o (EK.Wide (deconsUnchecked h) (deconsUnchecked w)))+      (Adapt $ \(Extent o (EK.Tall h w)) ->+         Extent o (EK.Tall (deconsUnchecked h) (deconsUnchecked w)))+      (Adapt $ \(Extent o (EK.General h w)) ->+         Extent o (EK.General (deconsUnchecked h) (deconsUnchecked w)))+++ newtype Fuse height fuse width vert horiz =    Fuse {       getFuse ::@@ -412,6 +468,100 @@       (Fuse $        \(Extent o (EK.General h f0)) (Extent _ (EK.General f1 w)) ->          toMaybe (f0==f1) $ Extent o (EK.General h w))+++newtype Kronecker heightA widthA heightB widthB vert horiz =+   Kronecker {+      getKronecker ::+         Extent vert horiz heightA widthA ->+         Extent vert horiz heightB widthB ->+         Extent vert horiz (heightA,heightB) (widthA,widthB)+   }++kronecker ::+   (C vert, C horiz) =>+   Extent vert horiz heightA widthA ->+   Extent vert horiz heightB widthB ->+   Extent vert horiz (heightA,heightB) (widthA,widthB)+kronecker =+   getKronecker $+   switchTagPair+      (Kronecker $+       \(Extent o (EK.Square s0)) (Extent _ (EK.Square s1)) ->+         Extent o (EK.Square (s0,s1)))+      (Kronecker $+       \(Extent o (EK.Wide h0 w0)) (Extent _ (EK.Wide h1 w1)) ->+         Extent o (EK.Wide (h0,h1) (w0,w1)))+      (Kronecker $+       \(Extent o (EK.Tall h0 w0)) (Extent _ (EK.Tall h1 w1)) ->+         Extent o (EK.Tall (h0,h1) (w0,w1)))+      (Kronecker $+       \(Extent o (EK.General h0 w0)) (Extent _ (EK.General h1 w1)) ->+         Extent o (EK.General (h0,h1) (w0,w1)))++++{-+Tag table for 'beside'.++Small Small  Small Small -> Small Big+Small Small  Small Big   -> Small Big+Small Small  Big   Small -> Small Big+Small Small  Big   Big   -> Small Big+Small Big    Small Small -> Small Big+Small Big    Small Big   -> Small Big+Small Big    Big   Small -> Small Big+Small Big    Big   Big   -> Small Big+Big   Small  Small Small -> Small Big+Big   Small  Small Big   -> Small Big+Big   Small  Big   Small -> Big   Big+Big   Small  Big   Big   -> Big   Big+Big   Big    Small Small -> Small Big+Big   Big    Small Big   -> Small Big+Big   Big    Big   Small -> Big   Big+Big   Big    Big   Big   -> Big   Big+-}+newtype AppendMode vertA vertB vertC height widthA widthB =+   AppendMode (+      Extent vertA Big height widthA ->+      Extent vertB Big height widthB ->+      Extent vertC Big height (widthA:+:widthB)+   )++appendLeftAux ::+   (C vertA, C vertB) => AppendMode vertA vertB vertA height widthA widthB+appendLeftAux =+   AppendMode $ \extentA extentB ->+      widen (width extentA :+: width extentB) extentA++appendSame :: (C vert) => AppendMode vert vert vert height widthA widthB+appendSame = appendLeftAux++appendLeft :: (C vert) => AppendMode vert Big vert height widthA widthB+appendLeft = appendLeftAux++appendRight :: (C vert) => AppendMode Big vert vert height widthA widthB+appendRight =+   AppendMode $ \extentA extentB ->+      widen (width extentA :+: width extentB) extentB++type family Append a b+type instance Append Small b = Small+type instance Append Big   b = b++newtype+   AppendAny vertB height widthA widthB vertA =+      AppendAny {+         getAppendAny ::+            AppendMode vertA vertB (Append vertA vertB) height widthA widthB+      }++appendAny ::+   (C vertA, C vertB) =>+   AppendMode vertA vertB (Append vertA vertB) height widthA widthB+appendAny =+   getAppendAny $ switchTag (AppendAny appendLeftAux) (AppendAny appendRight)+   type family Multiply a b

@@ -1,29 +1,310 @@ {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-} module Numeric.LAPACK.Matrix.Hermitian (-   module Numeric.LAPACK.Matrix.Hermitian.Basic,-   module Numeric.LAPACK.Matrix.Hermitian.Linear,+   Hermitian,+   Transposition(..), +   size,+   fromList,+   autoFromList,+   identity,+   diagonal,+   takeDiagonal,+   forceOrder,++   stack, (*%%%#),+   split,+   takeTopLeft,+   takeTopRight,+   takeBottomRight,++   multiplyVector,+   square,+   multiplyFull,+   outer,+   sumRank1, sumRank1NonEmpty,+   sumRank2, sumRank2NonEmpty,++   toSquare,+   gramian,            gramianAdjoint,+   congruenceDiagonal, congruenceDiagonalAdjoint,+   congruence,         congruenceAdjoint,+   anticommutator,     anticommutatorAdjoint,+   addAdjoint,++   solve,+   inverse,+   determinant,+    eigenvalues,    eigensystem,    ) where  import qualified Numeric.LAPACK.Matrix.Hermitian.Eigen as Eigen-import Numeric.LAPACK.Matrix.Hermitian.Basic-import Numeric.LAPACK.Matrix.Hermitian.Linear--import Numeric.LAPACK.Matrix.Private (Square)+import qualified Numeric.LAPACK.Matrix.Hermitian.Linear as Linear+import qualified Numeric.LAPACK.Matrix.Hermitian.Basic as Basic+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Extent as Extent+import Numeric.LAPACK.Matrix.Array.Triangular (Hermitian)+import Numeric.LAPACK.Matrix.Array (Full, General, Square)+import Numeric.LAPACK.Matrix.Shape.Private (Order)+import Numeric.LAPACK.Matrix.Modifier (Transposition(NonTransposed, Transposed))+import Numeric.LAPACK.Matrix.Private (ShapeInt) import Numeric.LAPACK.Vector (Vector)-import Numeric.LAPACK.Scalar (RealOf)+import Numeric.LAPACK.Scalar (RealOf, one)  import qualified Numeric.Netlib.Class as Class +import qualified Data.Array.Comfort.Storable.Unchecked as Array import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Shape ((:+:)) +import Foreign.Storable (Storable) +import qualified Data.NonEmpty as NonEmpty+import Data.Tuple.HT (mapFst)+++size :: Hermitian sh a -> sh+size = MatrixShape.hermitianSize . ArrMatrix.shape++fromList :: (Shape.C sh, Storable a) => Order -> sh -> [a] -> Hermitian sh a+fromList order sh = ArrMatrix.lift0 . Basic.fromList order sh++autoFromList :: (Storable a) => Order -> [a] -> Hermitian ShapeInt a+autoFromList order = ArrMatrix.lift0 . Basic.autoFromList order++identity :: (Shape.C sh, Class.Floating a) => Order -> sh -> Hermitian sh a+identity order = ArrMatrix.lift0 . Basic.identity order++diagonal ::+   (Shape.C sh, Class.Floating a) =>+   Order -> Vector sh (RealOf a) -> Hermitian sh a+diagonal order = ArrMatrix.lift0 . Basic.diagonal order++takeDiagonal ::+   (Shape.C sh, Class.Floating a) =>+   Hermitian sh a -> Vector sh (RealOf a)+takeDiagonal = Basic.takeDiagonal . ArrMatrix.toVector++forceOrder ::+   (Shape.C sh, Class.Floating a) =>+   Order -> Hermitian sh a -> Hermitian sh a+forceOrder = ArrMatrix.lift1 . Basic.forceOrder++{- |+> toSquare (stack a b c)+>+> =+>+> toSquare a ||| b+> ===+> adjoint b ||| toSquare c++It holds @order (stack a b c) = order b@.+The function is most efficient when the order of all blocks match.+-}+stack ::+   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>+   Hermitian sh0 a -> General sh0 sh1 a -> Hermitian sh1 a ->+   Hermitian (sh0:+:sh1) a+stack = ArrMatrix.lift3 Basic.stack++infixr 2 *%%%#++(*%%%#) ::+   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>+   (Hermitian sh0 a, General sh0 sh1 a) -> Hermitian sh1 a ->+   Hermitian (sh0:+:sh1) a+(*%%%#) = uncurry stack+++split ::+   (Shape.C sh0, Shape.C sh1, Class.Floating a) =>+   Hermitian (sh0:+:sh1) a ->+   (Hermitian sh0 a, General sh0 sh1 a, Hermitian sh1 a)+split a = (takeTopLeft a, takeTopRight a, takeBottomRight a)++takeTopLeft ::+   (Shape.C sh0, Shape.C sh1, Class.Floating a) =>+   Hermitian (sh0:+:sh1) a -> Hermitian sh0 a+takeTopLeft = ArrMatrix.lift1 Basic.takeTopLeft++takeTopRight ::+   (Shape.C sh0, Shape.C sh1, Class.Floating a) =>+   Hermitian (sh0:+:sh1) a -> General sh0 sh1 a+takeTopRight = ArrMatrix.lift1 Basic.takeTopRight++takeBottomRight ::+   (Shape.C sh0, Shape.C sh1, Class.Floating a) =>+   Hermitian (sh0:+:sh1) a -> Hermitian sh1 a+takeBottomRight = ArrMatrix.lift1 Basic.takeBottomRight+++multiplyVector ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   Transposition -> Hermitian sh a -> Vector sh a -> Vector sh a+multiplyVector order = Basic.multiplyVector order . ArrMatrix.toVector++square ::+   (Shape.C sh, Eq sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a+square = ArrMatrix.lift1 Basic.square++multiplyFull ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width,+    Class.Floating a) =>+   Transposition -> Hermitian height a ->+   Full vert horiz height width a ->+   Full vert horiz height width a+multiplyFull = ArrMatrix.lift2 . Basic.multiplyFull++outer ::+   (Shape.C sh, Class.Floating a) => Order -> Vector sh a -> Hermitian sh a+outer order = ArrMatrix.lift0 . Basic.outer order++sumRank1 ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   Order -> sh -> [(RealOf a, Vector sh a)] -> Hermitian sh a+sumRank1 order sh = ArrMatrix.lift0 . Basic.sumRank1 order sh++sumRank1NonEmpty ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   Order -> NonEmpty.T [] (RealOf a, Vector sh a) -> Hermitian sh a+sumRank1NonEmpty order (NonEmpty.Cons x xs) =+   sumRank1 order (Array.shape $ snd x) (x:xs)++sumRank2 ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   Order -> sh -> [(a, (Vector sh a, Vector sh a))] -> Hermitian sh a+sumRank2 order sh = ArrMatrix.lift0 . Basic.sumRank2 order sh++sumRank2NonEmpty ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   Order -> NonEmpty.T [] (a, (Vector sh a, Vector sh a)) -> Hermitian sh a+sumRank2NonEmpty order (NonEmpty.Cons xy xys) =+   sumRank2 order (Array.shape $ fst $ snd xy) (xy:xys)++toSquare ::+   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Square sh a+toSquare = ArrMatrix.lift1 Basic.toSquare++{- |+gramian A = A^H * A+-}+gramian ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   General height width a -> Hermitian width a+gramian = ArrMatrix.lift1 Basic.gramian++{- |+gramianAdjoint A = A * A^H = gramian (A^H)+-}+gramianAdjoint ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   General height width a -> Hermitian height a+gramianAdjoint = ArrMatrix.lift1 Basic.gramianAdjoint++{- |+congruenceDiagonal D A = A^H * D * A+-}+congruenceDiagonal ::+   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+   Vector height (RealOf a) -> General height width a -> Hermitian width a+congruenceDiagonal = ArrMatrix.lift1 . Basic.congruenceDiagonal++{- |+congruenceDiagonalAdjoint A D = A * D * A^H+-}+congruenceDiagonalAdjoint ::+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   General height width a -> Vector width (RealOf a) -> Hermitian height a+congruenceDiagonalAdjoint a =+   ArrMatrix.lift0 . Basic.congruenceDiagonalAdjoint (ArrMatrix.toVector a)++{- |+congruence B A = A^H * B * A+-}+congruence ::+   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+   Hermitian height a -> General height width a -> Hermitian width a+congruence = ArrMatrix.lift2 Basic.congruence++{- |+congruenceAdjoint B A = A * B * A^H+-}+congruenceAdjoint ::+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   General height width a -> Hermitian width a -> Hermitian height a+congruenceAdjoint = ArrMatrix.lift2 Basic.congruenceAdjoint+++{- |+anticommutator A B  =  A^H * B + B^H * A++Not exactly a matrix anticommutator,+thus I like to call it Hermitian anticommutator.+-}+anticommutator ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Eq width,+    Class.Floating a) =>+   Full vert horiz height width a ->+   Full vert horiz height width a -> Hermitian width a+anticommutator = ArrMatrix.lift2 $ Basic.scaledAnticommutator one++{- |+anticommutatorAdjoint A B+   = A * B^H + B * A^H+   = anticommutator (adjoint A) (adjoint B)+-}+anticommutatorAdjoint ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Eq width,+    Class.Floating a) =>+   Full vert horiz height width a ->+   Full vert horiz height width a -> Hermitian height a+anticommutatorAdjoint = ArrMatrix.lift2 $ Basic.scaledAnticommutatorAdjoint one++{- |+scaledAnticommutator alpha A B  =  alpha * A^H * B + conj alpha * B^H * A+-}+_scaledAnticommutator ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Eq width,+    Class.Floating a) =>+   a ->+   Full vert horiz height width a ->+   Full vert horiz height width a -> Hermitian width a+_scaledAnticommutator = ArrMatrix.lift2 . Basic.scaledAnticommutator++{- |+addAdjoint A = A^H + A+-}+addAdjoint :: (Shape.C sh, Class.Floating a) => Square sh a -> Hermitian sh a+addAdjoint = ArrMatrix.lift1 Basic.addAdjoint++++solve ::+   (Extent.C vert, Extent.C horiz,+    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>+   Hermitian sh a -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a+solve = ArrMatrix.lift2 Linear.solve++inverse :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a+inverse = ArrMatrix.lift1 Linear.inverse++determinant :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> RealOf a+determinant = Linear.determinant . ArrMatrix.toVector+++ eigenvalues ::    (Shape.C sh, Class.Floating a) =>    Hermitian sh a -> Vector sh (RealOf a)-eigenvalues = Eigen.values+eigenvalues = Eigen.values . ArrMatrix.toVector  {- | For symmetric eigenvalue problems, @eigensystem@ and @schur@ coincide.@@ -31,4 +312,4 @@ eigensystem ::    (Shape.C sh, Class.Floating a) =>    Hermitian sh a -> (Square sh a, Vector sh (RealOf a))-eigensystem = Eigen.decompose+eigensystem = mapFst ArrMatrix.lift0 . Eigen.decompose . ArrMatrix.toVector

@@ -5,22 +5,32 @@    Transposition(..),    fromList,    autoFromList,+   recheck,    identity,    diagonal,    takeDiagonal,    forceOrder,+    stack,+   takeTopLeft,+   takeTopRight,+   takeBottomRight,     multiplyVector,-   square,    multiplyFull,+   square, power,    outer,-   sumRank1, sumRank1NonEmpty,-   sumRank2, sumRank2NonEmpty,+   sumRank1,+   sumRank2,     toSquare,-   covariance,+   gramian,              gramianAdjoint,+   congruenceDiagonal,   congruenceDiagonalAdjoint,+   congruence,           congruenceAdjoint,+   scaledAnticommutator, scaledAnticommutatorAdjoint,    addAdjoint,++   takeUpper,    ) where  import qualified Numeric.LAPACK.Matrix.Symmetric.Private as Symmetric@@ -28,22 +38,26 @@ import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent import qualified Numeric.LAPACK.Matrix.Basic as Basic+import qualified Numeric.LAPACK.Split as Split import Numeric.LAPACK.Matrix.Hermitian.Private (Diagonal(..), TakeDiagonal(..)) import Numeric.LAPACK.Matrix.Triangular.Private-         (forPointers, pack, packRect, unpack, unpackToTemp,+         (forPointers, pack, unpack, unpackToTemp,           diagonalPointers, diagonalPointerPairs,           rowMajorPointers, columnMajorPointers) import Numeric.LAPACK.Matrix.Shape.Private          (Order(RowMajor,ColumnMajor), flipOrder, sideSwapFromOrder,           uploFromOrder)+import Numeric.LAPACK.Matrix.Modifier+         (Transposition(NonTransposed, Transposed), transposeOrder,+          Conjugation(Conjugated), conjugatedOnRowMajor) import Numeric.LAPACK.Matrix.Private-         (Full, General, argGeneral, Square, argSquare, ZeroInt, zeroInt,-          Transposition(NonTransposed, Transposed), transposeOrder,-          Conjugation(Conjugated))+         (Full, General, Square, argSquare, ShapeInt, shapeInt) import Numeric.LAPACK.Vector (Vector)-import Numeric.LAPACK.Scalar (RealOf, zero, one, fromReal, realPart)+import Numeric.LAPACK.Scalar (RealOf, zero, one)+import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked)) import Numeric.LAPACK.Private-         (fill, lacgv, copyConjugate, conjugateToTemp, condConjugateToTemp)+         (fill, lacgv, realPtr,+          copyConjugate, condConjugate, conjugateToTemp, condConjugateToTemp)  import qualified Numeric.BLAS.FFI.Generic as BlasGen import qualified Numeric.BLAS.FFI.Complex as BlasComplex@@ -55,10 +69,10 @@ import qualified Data.Array.Comfort.Storable as CheckedArray import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable.Unchecked (Array(Array))-import Data.Array.Comfort.Shape ((:+:))+import Data.Array.Comfort.Shape ((:+:)((:+:)))  import Foreign.C.Types (CInt, CChar)-import Foreign.ForeignPtr (withForeignPtr)+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr) import Foreign.Ptr (Ptr) import Foreign.Storable (Storable, poke, peek) @@ -66,10 +80,11 @@ import Control.Monad.IO.Class (liftIO) import Control.Monad (when) -import qualified Data.NonEmpty as NonEmpty import Data.Foldable (forM_) +import Data.Function.HT (powerAssociative) + type Hermitian sh = Array (MatrixShape.Hermitian sh)  @@ -77,14 +92,26 @@ fromList order sh =    CheckedArray.fromList (MatrixShape.Hermitian order sh) -autoFromList :: (Storable a) => Order -> [a] -> Hermitian ZeroInt a+autoFromList :: (Storable a) => Order -> [a] -> Hermitian ShapeInt a autoFromList order xs =    fromList order-      (zeroInt $ MatrixShape.triangleExtent "Hermitian.autoFromList" $+      (shapeInt $ MatrixShape.triangleExtent "Hermitian.autoFromList" $        length xs)       xs +uncheck :: Hermitian sh a -> Hermitian (Unchecked sh) a+uncheck =+   Array.mapShape $+      \(MatrixShape.Hermitian order sh) ->+         MatrixShape.Hermitian order (Unchecked sh) +recheck :: Hermitian (Unchecked sh) a -> Hermitian sh a+recheck =+   Array.mapShape $+      \(MatrixShape.Hermitian order (Unchecked sh)) ->+         MatrixShape.Hermitian order sh++ identity :: (Shape.C sh, Class.Floating a) => Order -> sh -> Hermitian sh a identity order sh =    Array.unsafeCreateWithSize (MatrixShape.Hermitian order sh) $@@ -110,7 +137,7 @@    fill zero triSize aPtr    withForeignPtr x $ \xPtr ->       forM_ (diagonalPointerPairs order (Shape.size sh) xPtr aPtr) $-         \(srcPtr,dstPtr) -> poke dstPtr . fromReal =<< peek srcPtr+         \(srcPtr,dstPtr) -> poke (realPtr dstPtr) =<< peek srcPtr   takeDiagonal ::@@ -123,13 +150,13 @@       (TakeDiagonal takeDiagonalAux) (TakeDiagonal takeDiagonalAux)  takeDiagonalAux ::-   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>+   (Shape.C sh, Storable a, RealOf a ~ ar, Storable ar) =>    Hermitian sh a -> Vector sh ar takeDiagonalAux (Array (MatrixShape.Hermitian order sh) a) =    Array.unsafeCreateWithSize sh $ \n xPtr ->    withForeignPtr a $ \aPtr ->       forM_ (diagonalPointerPairs order n xPtr aPtr) $-         \(dstPtr,srcPtr) -> poke dstPtr . realPart =<< peek srcPtr+         \(dstPtr,srcPtr) -> poke dstPtr =<< peek (realPtr srcPtr)   {-@@ -149,27 +176,23 @@ fromUpperPart ::    (Extent.C vert, Shape.C height, Shape.C width, Class.Floating a) =>    Full vert Extent.Small height width a -> Hermitian width a-fromUpperPart (Array (MatrixShape.Full order extent) a) =-   let (height,width) = Extent.dimensions extent-       m = Shape.size height-       n = Shape.size width-       k = case order of RowMajor -> n; ColumnMajor -> m-   in Array.unsafeCreate (MatrixShape.Hermitian order width) $ \bPtr ->-      withForeignPtr a $ \aPtr -> packRect order n k aPtr bPtr+fromUpperPart = Triangular.fromUpperPart MatrixShape.Hermitian +{-+Naming is inconsistent to Triangular.takeUpper,+because here Hermitian is the input+and in Triangular.takeUpper, Triangular is the output.+-}+takeUpper ::+   (Shape.C sh, Class.Floating a) =>+   Hermitian sh a ->+   Array (MatrixShape.UpperTriangular MatrixShape.NonUnit sh) a+takeUpper =+   Array.mapShape+      (\(MatrixShape.Hermitian order sh) ->+         MatrixShape.Triangular MatrixShape.NonUnit MatrixShape.upper order sh) -{- |-> toSquare (stack a b c)->-> =->-> toSquare a ||| b-> ===-> adjoint b ||| toSquare c -It holds @order (stack a b c) = order b@.-The function is most efficient when the order of all blocks match.--} stack ::    (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>    Hermitian sh0 a -> General sh0 sh1 a -> Hermitian sh1 a ->@@ -179,7 +202,29 @@    in Triangular.stack "Hermitian" (MatrixShape.Hermitian order)          (forceOrder order a) b (forceOrder order c) +takeTopLeft ::+   (Shape.C sh0, Shape.C sh1, Class.Floating a) =>+   Hermitian (sh0:+:sh1) a -> Hermitian sh0 a+takeTopLeft =+   Triangular.takeTopLeft+      (\(MatrixShape.Hermitian order sh@(sh0:+:_sh1)) ->+         (MatrixShape.Hermitian order sh0, (order,sh))) +takeTopRight ::+   (Shape.C sh0, Shape.C sh1, Class.Floating a) =>+   Hermitian (sh0:+:sh1) a -> General sh0 sh1 a+takeTopRight =+   Triangular.takeTopRight (\(MatrixShape.Hermitian order sh) -> (order,sh))++takeBottomRight ::+   (Shape.C sh0, Shape.C sh1, Class.Floating a) =>+   Hermitian (sh0:+:sh1) a -> Hermitian sh1 a+takeBottomRight =+   Triangular.takeBottomRight+      (\(MatrixShape.Hermitian order sh@(_sh0:+:sh1)) ->+         (MatrixShape.Hermitian order sh1, (order,sh)))++ multiplyVector ::    (Shape.C sh, Eq sh, Class.Floating a) =>    Transposition -> Hermitian sh a -> Vector sh a -> Vector sh a@@ -188,7 +233,7 @@       Array.unsafeCreateWithSize shX $ \n yPtr -> do    Call.assert "Hermitian.multiplyVector: width shapes mismatch" (shA == shX)    evalContT $ do-      let conj = transposeOrder transposed order == RowMajor+      let conj = conjugatedOnRowMajor $ transposeOrder transposed order       uploPtr <- Call.char $ uploFromOrder order       nPtr <- Call.cint n       alphaPtr <- Call.number one@@ -200,17 +245,28 @@       liftIO $ do          BlasGen.hpmv             uploPtr nPtr alphaPtr aPtr xPtr incxPtr betaPtr yPtr incyPtr-         when conj $ lacgv nPtr yPtr incyPtr+         condConjugate conj nPtr yPtr incyPtr  -square ::-   (Shape.C sh, Eq sh, Class.Floating a) =>-   Hermitian sh a -> Hermitian sh a+square :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a square (Array shape@(MatrixShape.Hermitian order sh) a) =    Array.unsafeCreate shape $       Symmetric.square Conjugated order (Shape.size sh) a +{-+Requires frequent unpacking and packing of triangles.+-}+power ::+   (Shape.C sh, Class.Floating a) => Integer -> Hermitian sh a -> Hermitian sh a+power n a0@(Array (MatrixShape.Hermitian order sh) _) =+   recheck $+   powerAssociative+      (\a b -> fromUpperPart $ multiplyFull NonTransposed a $ toSquare b)+      (identity order $ Unchecked sh)+      (uncheck a0)+      n + multiplyFull ::    (Extent.C vert, Extent.C horiz,     Shape.C height, Eq height, Shape.C width,@@ -263,26 +319,15 @@    liftIO $ do       fill zero triSize aPtr       act uploPtr nPtr incxPtr-      case order of-         RowMajor -> lacgv sizePtr aPtr incxPtr-         ColumnMajor -> return ()+      condConjugate (conjugatedOnRowMajor order) sizePtr aPtr incxPtr  outer ::    (Shape.C sh, Class.Floating a) => Order -> Vector sh a -> Hermitian sh a-outer order =-   getMap $-   Class.switchFloating-      (Map $ outerAux order) (Map $ outerAux order)-      (Map $ outerAux order) (Map $ outerAux order)--outerAux ::-   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>-   Order -> Vector sh a -> Hermitian sh a-outerAux order (Array sh x) =+outer order (Array sh x) =    Array.unsafeCreateWithSize (MatrixShape.Hermitian order sh) $       \triSize aPtr ->    evalContT $ do-      alphaPtr <- Call.real one+      alphaPtr <- realOneArg aPtr       xPtr <- ContT $ withForeignPtr x       withConjBuffer order sh triSize aPtr $ \uploPtr nPtr incxPtr ->          hpr uploPtr nPtr alphaPtr xPtr incxPtr aPtr@@ -302,7 +347,7 @@ newtype SumRank1 sh a = SumRank1 {getSumRank1 :: SumRank1_ sh (RealOf a) a}  sumRank1Aux ::-   (Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   (Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>    SumRank1_ sh ar a sumRank1Aux order sh xs =    Array.unsafeCreateWithSize (MatrixShape.Hermitian order sh) $@@ -318,13 +363,6 @@                hpr uploPtr nPtr alphaPtr xPtr incxPtr aPtr  -sumRank1NonEmpty ::-   (Shape.C sh, Eq sh, Class.Floating a) =>-   Order -> NonEmpty.T [] (RealOf a, Vector sh a) -> Hermitian sh a-sumRank1NonEmpty order (NonEmpty.Cons x xs) =-   sumRank1 order (Array.shape $ snd x) (x:xs)-- type HPR_ a =    Ptr CChar -> Ptr CInt ->    Ptr (RealOf a) -> Ptr a -> Ptr CInt -> Ptr a -> IO ()@@ -358,13 +396,7 @@                poke alphaPtr alpha                BlasGen.hpr2 uploPtr nPtr alphaPtr xPtr incPtr yPtr incPtr aPtr -sumRank2NonEmpty ::-   (Shape.C sh, Eq sh, Class.Floating a) =>-   Order -> NonEmpty.T [] (a, (Vector sh a, Vector sh a)) -> Hermitian sh a-sumRank2NonEmpty order (NonEmpty.Cons xy xys) =-   sumRank2 order (Array.shape $ fst $ snd xy) (xy:xys) - {- It is not strictly necessary to keep the 'order'. It would be neither more complicated nor less efficient@@ -388,56 +420,85 @@       Symmetric.unpack Conjugated order (Shape.size sh) aPtr bPtr  -{- |-A^H * A--}-covariance ::-   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+gramian ::+   (Shape.C height, Shape.C width, Class.Floating a) =>    General height width a -> Hermitian width a-covariance =-   getMap $-   Class.switchFloating-      (Map covarianceAux) (Map covarianceAux)-      (Map covarianceAux) (Map covarianceAux)+gramian (Array (MatrixShape.Full order extent) a) =+   Array.unsafeCreate (MatrixShape.Hermitian order $ Extent.width extent) $+   \bPtr -> gramianIO order a bPtr $ gramianParameters order extent -newtype Map f g a = Map {getMap :: f a -> g a}+gramianParameters ::+   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>+   Order ->+   Extent.Extent vert horiz height width ->+   ((Int, Int), (Char, Char, Int))+gramianParameters order extent =+   let (height, width) = Extent.dimensions extent+       n = Shape.size width+       k = Shape.size height+    in ((n,k),+         case order of+            ColumnMajor -> ('U', 'C', k)+            RowMajor -> ('L', 'N', n)) -covarianceAux ::-   (Shape.C height, Shape.C width, Eq width,-    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>-   General height width a -> Hermitian width a-covarianceAux = argGeneral $ \order height width a ->-   Array.unsafeCreate (MatrixShape.Hermitian order width) $ \bPtr -> do -   let n = Shape.size width-   let k = Shape.size height+gramianAdjoint ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   General height width a -> Hermitian height a+gramianAdjoint (Array (MatrixShape.Full order extent) a) =+   Array.unsafeCreate (MatrixShape.Hermitian order $ Extent.height extent) $+   \bPtr -> gramianIO order a bPtr $ gramianAdjointParameters order extent++gramianAdjointParameters ::+   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>+   Order ->+   Extent.Extent vert horiz height width ->+   ((Int, Int), (Char, Char, Int))+gramianAdjointParameters order extent =+   let (height, width) = Extent.dimensions extent+       n = Shape.size height+       k = Shape.size width+   in ((n,k),+         case order of+            ColumnMajor -> ('U', 'N', n)+            RowMajor -> ('L', 'C', k))++{-+Another way to unify 'gramian' and 'gramianAdjoint'+would have been this function:++> gramianConjugation ::+>    Conjugation -> General height width a -> Hermitian width a++with++> gramianAdjoint a = gramianConjugation (transpose a)++but I would like to have++> order (gramianAdjoint a) = order a+-}+gramianIO ::+   (Class.Floating a) =>+   Order ->+   ForeignPtr a -> Ptr a ->+   ((Int, Int), (Char, Char, Int)) -> IO ()+gramianIO order a bPtr ((n,k), (uplo,trans,lda)) =    evalContT $ do+      uploPtr <- Call.char uplo+      transPtr <- Call.char trans       nPtr <- Call.cint n       kPtr <- Call.cint k-      alphaPtr <- Call.number one+      alphaPtr <- realOneArg a       aPtr <- ContT $ withForeignPtr a-      betaPtr <- Call.number zero+      ldaPtr <- Call.leadingDim lda+      betaPtr <- realZeroArg a       cPtr <- Call.allocaArray (n*n)       ldcPtr <- Call.leadingDim n--      case order of-         ColumnMajor -> do-            uploPtr <- Call.char 'U'-            transPtr <- Call.char 'C'-            ldaPtr <- Call.leadingDim k-            liftIO $ do-               herk uploPtr transPtr-                  nPtr kPtr alphaPtr aPtr ldaPtr betaPtr cPtr ldcPtr-               pack ColumnMajor n cPtr bPtr--         RowMajor -> do-            uploPtr <- Call.char 'L'-            transPtr <- Call.char 'N'-            ldaPtr <- Call.leadingDim n-            liftIO $ do-               herk uploPtr transPtr-                  nPtr kPtr alphaPtr aPtr ldaPtr betaPtr cPtr ldcPtr-               pack RowMajor n cPtr bPtr+      liftIO $ do+         herk uploPtr transPtr+            nPtr kPtr alphaPtr aPtr ldaPtr betaPtr cPtr ldcPtr+         pack order n cPtr bPtr   type HERK_ a =@@ -456,9 +517,142 @@       (HERK BlasComplex.herk)  +skipCheckCongruence ::+   ((sh -> Unchecked sh) -> matrix0 -> matrix1) ->+   (matrix1 -> Hermitian (Unchecked sh) a) -> matrix0 -> Hermitian sh a+skipCheckCongruence mapSize f a =+   recheck $ f $ mapSize Unchecked a+++congruenceDiagonal ::+   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+   Vector height (RealOf a) -> General height width a -> Hermitian width a+congruenceDiagonal d =+   skipCheckCongruence Basic.mapWidth $ \a ->+      scaledAnticommutator 0.5 a $ Basic.scaleRowsReal d a++congruenceDiagonalAdjoint ::+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   General height width a -> Vector width (RealOf a) -> Hermitian height a+congruenceDiagonalAdjoint =+   flip $ \d -> skipCheckCongruence Basic.mapHeight $ \a ->+      scaledAnticommutatorAdjoint 0.5 a $ Basic.scaleColumnsReal d a+++congruence ::+   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+   Hermitian height a -> General height width a -> Hermitian width a+congruence b =+   skipCheckCongruence Basic.mapWidth $ \a ->+      scaledAnticommutator one a $+      Split.tallMultiplyR NonTransposed (takeHalf b) a++congruenceAdjoint ::+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   General height width a -> Hermitian width a -> Hermitian height a+congruenceAdjoint =+   flip $ \b -> skipCheckCongruence Basic.mapHeight $ \a ->+      scaledAnticommutatorAdjoint one a $+      Basic.swapMultiply (Split.tallMultiplyR Transposed) a (takeHalf b)+++data Corrupt = Corrupt+   deriving (Eq)+ {- |-A^H + A+> let b = takeHalf a+> ==>+> isTriangular b && a == addAdjoint b -}+takeHalf ::+   (Shape.C sh, Class.Floating a) =>+   Hermitian sh a -> Split.Square Corrupt sh a+takeHalf (Array (MatrixShape.Hermitian order sh) a) =+   Array.unsafeCreate (MatrixShape.Split Corrupt order (Extent.square sh)) $+      \bPtr -> evalContT $ do+   let n = Shape.size sh+   aPtr <- ContT $ withForeignPtr a+   nPtr <- Call.cint n+   alphaPtr <- Call.number 0.5+   incxPtr <- Call.cint (n+1)+   liftIO $ do+      unpack order n aPtr bPtr+      BlasGen.scal nPtr alphaPtr bPtr incxPtr+++scaledAnticommutator ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+   a ->+   Full vert horiz height width a ->+   Full vert horiz height width a -> Hermitian width a+scaledAnticommutator alpha arr (Array (MatrixShape.Full order extentB) b) = do+   let (Array (MatrixShape.Full _ extentA) a) = Basic.forceOrder order arr+   Array.unsafeCreate (MatrixShape.Hermitian order $ Extent.width extentB) $+         \cpPtr -> do+      Call.assert "Hermitian.anticommutator: extents mismatch"+         (extentA==extentB)+      scaledAnticommutatorIO alpha order a b cpPtr $+         gramianParameters order extentB++scaledAnticommutatorAdjoint ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+   a ->+   Full vert horiz height width a ->+   Full vert horiz height width a -> Hermitian height a+scaledAnticommutatorAdjoint+      alpha arr (Array (MatrixShape.Full order extentB) b) = do+   let (Array (MatrixShape.Full _ extentA) a) = Basic.forceOrder order arr+   Array.unsafeCreate (MatrixShape.Hermitian order $ Extent.height extentB) $+         \cpPtr -> do+      Call.assert "Hermitian.anticommutatorAdjoint: extents mismatch"+         (extentA==extentB)+      scaledAnticommutatorIO alpha order a b cpPtr $+         gramianAdjointParameters order extentB++scaledAnticommutatorIO ::+   (Class.Floating a) =>+   a ->+   Order -> ForeignPtr a -> ForeignPtr a -> Ptr a ->+   ((Int, Int), (Char, Char, Int)) -> IO ()+scaledAnticommutatorIO alpha order a b cpPtr ((n,k), (uplo,trans,lda)) =+   evalContT $ do+      uploPtr <- Call.char uplo+      transPtr <- Call.char trans+      nPtr <- Call.cint n+      kPtr <- Call.cint k+      alphaPtr <- Call.number alpha+      aPtr <- ContT $ withForeignPtr a+      ldaPtr <- Call.leadingDim lda+      bPtr <- ContT $ withForeignPtr b+      let ldbPtr = ldaPtr+      betaPtr <- realZeroArg aPtr+      cPtr <- Call.allocaArray (n*n)+      ldcPtr <- Call.leadingDim n+      liftIO $ do+         her2k uploPtr transPtr nPtr kPtr alphaPtr+            aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr+         pack order n cPtr cpPtr+++type HER2K_ a =+   Ptr CChar -> Ptr CChar -> Ptr CInt -> Ptr CInt -> Ptr a ->+   Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->+   Ptr (RealOf a) -> Ptr a -> Ptr CInt -> IO ()++newtype HER2K a = HER2K {getHER2K :: HER2K_ a}++her2k :: Class.Floating a => HER2K_ a+her2k =+   getHER2K $+   Class.switchFloating+      (HER2K BlasReal.syr2k)+      (HER2K BlasReal.syr2k)+      (HER2K BlasComplex.her2k)+      (HER2K BlasComplex.her2k)++ addAdjoint, _addAdjoint ::    (Shape.C sh, Class.Floating a) => Square sh a -> Hermitian sh a _addAdjoint =@@ -513,3 +707,23 @@                forPointers (rowMajorPointers n fullPtr packedPtr) $                   \nPtr (srcPtr,dstPtr) ->                      BlasGen.copy nPtr srcPtr incxPtr dstPtr incxPtr+++realZeroArg, realOneArg ::+   (Class.Floating a) => f a -> ContT r IO (Ptr (RealOf a))+realZeroArg =+   runRealArg $+   Class.switchFloating+      (RealArg $ const $ Call.number zero)+      (RealArg $ const $ Call.number zero)+      (RealArg $ const $ Call.number zero)+      (RealArg $ const $ Call.number zero)+realOneArg =+   runRealArg $+   Class.switchFloating+      (RealArg $ const $ Call.number one)+      (RealArg $ const $ Call.number one)+      (RealArg $ const $ Call.number one)+      (RealArg $ const $ Call.number one)++newtype RealArg f g a = RealArg {runRealArg :: f a -> g (Ptr (RealOf a))}

@@ -4,13 +4,12 @@    decompose,    ) where -import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian)-import Numeric.LAPACK.Matrix.Square.Basic (Square)- import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian) import Numeric.LAPACK.Matrix.Hermitian.Private (TakeDiagonal(..))-import Numeric.LAPACK.Matrix.Shape.Private-         (Order(RowMajor,ColumnMajor), uploFromOrder)+import Numeric.LAPACK.Matrix.Square.Basic (Square)+import Numeric.LAPACK.Matrix.Shape.Private (Order(ColumnMajor), uploFromOrder)+import Numeric.LAPACK.Matrix.Modifier (conjugatedOnRowMajor) import Numeric.LAPACK.Vector (Vector) import Numeric.LAPACK.Scalar (RealOf) import Numeric.LAPACK.Private@@ -84,7 +83,8 @@    evalContT $ do       jobzPtr <- Call.char 'V'       uploPtr <- Call.char $ uploFromOrder order-      aPtr <- copyCondConjugateToTemp (order==RowMajor) (triangleSize n) a+      aPtr <-+         copyCondConjugateToTemp (conjugatedOnRowMajor order) (triangleSize n) a       ldzPtr <- Call.leadingDim n       liftIO $ withInfo eigenMsg "hpev" $          hpev jobzPtr uploPtr n aPtr wPtr zPtr ldzPtr

@@ -5,15 +5,15 @@    determinant,    ) where -import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian)-import Numeric.LAPACK.Matrix.Private (Full)- import qualified Numeric.LAPACK.Matrix.Symmetric.Private as Symmetric import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent+import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian) import Numeric.LAPACK.Matrix.Hermitian.Private (Determinant(..))-import Numeric.LAPACK.Matrix.Private (Conjugation(Conjugated))+import Numeric.LAPACK.Matrix.Modifier (Conjugation(Conjugated))+import Numeric.LAPACK.Matrix.Private (Full) import Numeric.LAPACK.Scalar (RealOf, absoluteSquared)+import Numeric.LAPACK.Private (realPtr)  import qualified Numeric.Netlib.Class as Class @@ -23,28 +23,25 @@  import System.IO.Unsafe (unsafePerformIO) -import Foreign.Ptr (Ptr, castPtr)+import Foreign.Ptr (Ptr) import Foreign.Storable (peek)   solve ::    (Extent.C vert, Extent.C horiz,     Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>-   Hermitian sh a ->-   Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a+   Hermitian sh a -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a solve (Array (MatrixShape.Hermitian orderA shA) a) =    Symmetric.solve "Hermitian.solve" Conjugated orderA shA a  -inverse ::-   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a+inverse :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a inverse (Array shape@(MatrixShape.Hermitian order sh) a) =    Array.unsafeCreateWithSize shape $       Symmetric.inverse Conjugated order (Shape.size sh) a  -determinant ::-   (Shape.C sh, Class.Floating a) => Hermitian sh a -> RealOf a+determinant :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> RealOf a determinant =    getDeterminant $    Class.switchFloating@@ -63,7 +60,7 @@    (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>    (Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar peekBlockDeterminant (a0Ptr,ext) = do-   let peekReal = peek . castPtr+   let peekReal = peek . realPtr    a0 <- peekReal a0Ptr    case ext of       Nothing -> return a0

@@ -1,14 +1,12 @@ {-# LANGUAGE TypeFamilies #-} module Numeric.LAPACK.Matrix.Hermitian.Private where -import Numeric.LAPACK.Vector (Vector) import Numeric.LAPACK.Scalar (RealOf)  -newtype Diagonal f sh a =-   Diagonal {runDiagonal :: Vector sh (RealOf a) -> f a}+newtype Diagonal f g a = Diagonal {runDiagonal :: g (RealOf a) -> f a} -newtype TakeDiagonal f sh a =-   TakeDiagonal {runTakeDiagonal :: f a -> Vector sh (RealOf a)}+newtype TakeDiagonal f g a =+   TakeDiagonal {runTakeDiagonal :: f a -> g (RealOf a)}  newtype Determinant f a = Determinant {getDeterminant :: f a -> RealOf a}

@@ -1,5 +1,48 @@ module Numeric.LAPACK.Matrix.HermitianPositiveDefinite (-   module Numeric.LAPACK.Matrix.HermitianPositiveDefinite.Linear,+   solve,+   solveDecomposed,+   inverse,+   decompose,+   determinant,    ) where -import Numeric.LAPACK.Matrix.HermitianPositiveDefinite.Linear+import qualified Numeric.LAPACK.Matrix.HermitianPositiveDefinite.Linear+                                                                  as Linear+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix.Extent as Extent+import Numeric.LAPACK.Matrix.Array.Triangular (Hermitian, Upper)+import Numeric.LAPACK.Matrix.Array (Full)+import Numeric.LAPACK.Scalar (RealOf)++import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Shape as Shape+++solve ::+   (Extent.C vert, Extent.C horiz,+    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>+   Hermitian sh a -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a+solve = ArrMatrix.lift2 Linear.solve++{- |+> solve a b == solveDecomposed (decompose a) b+> solve (gramian u) b == solveDecomposed u b+-}+solveDecomposed ::+   (Extent.C vert, Extent.C horiz,+    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>+   Upper sh a -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a+solveDecomposed = ArrMatrix.lift2 Linear.solveDecomposed++inverse :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a+inverse = ArrMatrix.lift1 Linear.inverse++{- |+Cholesky decomposition+-}+decompose :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> Upper sh a+decompose = ArrMatrix.lift1 Linear.decompose++determinant :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> RealOf a+determinant = Linear.determinant . ArrMatrix.toVector

@@ -7,16 +7,15 @@    determinant,    ) where -import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian)-import Numeric.LAPACK.Matrix.Triangular.Basic (Upper, takeDiagonal)-import Numeric.LAPACK.Matrix.Private (Full)- import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent+import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian) import Numeric.LAPACK.Matrix.Hermitian.Private (Determinant(..))+import Numeric.LAPACK.Matrix.Triangular.Basic (Upper, takeDiagonal) import Numeric.LAPACK.Matrix.Triangular.Private (copyTriangleToTemp) import Numeric.LAPACK.Matrix.Shape.Private (NonUnit(NonUnit), uploFromOrder)-import Numeric.LAPACK.Matrix.Private (Conjugation(Conjugated))+import Numeric.LAPACK.Matrix.Modifier (Conjugation(Conjugated))+import Numeric.LAPACK.Matrix.Private (Full) import Numeric.LAPACK.Linear.Private (solver) import Numeric.LAPACK.Scalar (RealOf, realPart) import Numeric.LAPACK.Private (copyBlock, withInfo, rankMsg, definiteMsg)@@ -49,10 +48,6 @@          withInfo definiteMsg "ppsv" $             LapackGen.ppsv uploPtr nPtr nrhsPtr apPtr xPtr ldxPtr -{- |-> solve a b == solveDecomposed (decompose a) b-> solve (covariance u) b == solveDecomposed u b--} solveDecomposed ::    (Extent.C vert, Extent.C horiz,     Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>@@ -80,9 +75,6 @@          withInfo definiteMsg "pptrf" $ LapackGen.pptrf uploPtr nPtr bPtr          withInfo rankMsg "pptri" $ LapackGen.pptri uploPtr nPtr bPtr -{- |-Cholesky decomposition--} decompose ::    (Shape.C sh, Class.Floating a) => Hermitian sh a -> Upper sh a decompose@@ -111,4 +103,5 @@    (Shape.C sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>    Hermitian sh a -> ar determinantAux =-   (^(2::Int)) . product . map realPart . Array.toList . takeDiagonal . decompose+   (^(2::Int)) . product . map realPart .+   Array.toList . takeDiagonal . decompose

@@ -1,86 +1,38 @@ {-# LANGUAGE TypeFamilies #-} module Numeric.LAPACK.Matrix.Indexed where -import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape-import qualified Numeric.LAPACK.Matrix.Shape.Box as Box-import qualified Numeric.LAPACK.Matrix.Extent as Extent-import Numeric.LAPACK.Scalar (conjugate, zero)+import qualified Numeric.LAPACK.Matrix.Plain.Indexed as ArrIndexed+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix.Type as Type+import qualified Numeric.LAPACK.Permutation as Perm+import Numeric.LAPACK.Matrix.Type (Matrix)+import Numeric.LAPACK.Scalar (one, zero)+import Numeric.LAPACK.Permutation.Private (Permutation(Permutation))  import qualified Numeric.Netlib.Class as Class -import qualified Type.Data.Num.Unary as Unary- import qualified Data.Array.Comfort.Storable.Unchecked as UArray import qualified Data.Array.Comfort.Shape as Shape-import Data.Array.Comfort.Storable (Array, (!))--import Data.Tuple.HT (swap)-import Data.Bool.HT (if')+import Data.Array.Comfort.Storable ((!))   infixl 9 #! -class (Box.Box sh) => Indexed sh where+class (Type.Box typ) => Indexed typ where    (#!) ::       (Class.Floating a) =>-      Array sh a ->-      (Shape.Index (Box.HeightOf sh), Shape.Index (Box.WidthOf sh)) -> a--instance-   (Extent.C vert, Extent.C horiz, Shape.Indexed height, Shape.Indexed width) =>-      Indexed (MatrixShape.Full vert horiz height width) where-   (#!) = (!)--instance (Shape.Indexed size) => Indexed (MatrixShape.Hermitian size) where-   arr#!ij =-      if Shape.inBounds (UArray.shape arr) ij-         then arr UArray.! ij-         else conjugate $ arr ! swap ij--instance-   (MatrixShape.Content lo, MatrixShape.TriDiag diag, MatrixShape.Content up,-    Shape.Indexed size) =>-      Indexed (MatrixShape.Triangular lo diag up size) where-   arr#!ij =-      let sh = UArray.shape arr-      in if Shape.inBounds sh ij-            then arr UArray.! ij-            else-               MatrixShape.caseDiagUpLoSym (MatrixShape.triangularUplo sh)-                  (checkedZero "Diagonal" sh ij)-                  (checkedZero "UpperTriangular" sh ij)-                  (checkedZero "LowerTriangular" sh ij)-                  (arr ! swap ij)+      Matrix typ a ->+      (Shape.Index (Type.HeightOf typ), Shape.Index (Type.WidthOf typ)) -> a -instance-   (Unary.Natural sub, Unary.Natural super,-    Extent.C vert, Extent.C horiz, Shape.Indexed height, Shape.Indexed width) =>-      Indexed (MatrixShape.Banded sub super vert horiz height width) where-   arr#!ij =-      let boxIx = uncurry MatrixShape.InsideBox-          sh = UArray.shape arr-      in if Shape.inBounds sh $ boxIx ij-            then arr UArray.! boxIx ij-            else checkedZero "Banded" sh ij+instance (ArrIndexed.Indexed sh) => Indexed (ArrMatrix.Array sh) where+   ArrMatrix.Array arr #! ij = arr ArrIndexed.#! ij -instance-   (Unary.Natural off, Shape.Indexed size) =>-      Indexed (MatrixShape.BandedHermitian off size) where-   arr#!ij =-      let boxIx = uncurry MatrixShape.InsideBox-          sh = UArray.shape arr-      in if' (Shape.inBounds sh $ boxIx ij)-            (arr UArray.! boxIx ij) $-         if' (Shape.inBounds sh $ boxIx $ swap ij)-            (conjugate $ arr UArray.! boxIx (swap ij))-         (checkedZero "BandedHermitian" sh ij)+instance (Shape.Indexed size) => Indexed (Type.Scale size) where+   Type.Scale sh a #! (i,j) =+      if Shape.offset sh i == Shape.offset sh j then a else zero -checkedZero ::-   (Box.Box shape, Class.Floating a,-    Box.HeightOf shape ~ height, Shape.Indexed height,-    Box.WidthOf shape ~ width, Shape.Indexed width) =>-   String -> shape -> (Shape.Index height, Shape.Index width) -> a-checkedZero name sh ij =-   if Shape.inBounds (Box.height sh, Box.width sh) ij-      then zero-      else error $ "Matrix.Indexed." ++ name ++ ": index out of range"+instance (Shape.Indexed size) => Indexed (Permutation size) where+   Type.Permutation (Permutation perm) #! (i,j) =+      let psh@(Perm.Shape sh) = UArray.shape perm+          reindex = Shape.indexFromOffset psh . Shape.offset sh+      in if perm ! reindex i == reindex j then one else zero