(m      !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~ Safe<=DR9union binary lift : apply function on _union_ of two SetsGintersection binary lift : apply function on _intersection_ of two Sets inner productmultiplication by a scalarRing zero elementRing +negate the values in a functorsubtract two Additive objectslinear interpolation#`hilbertDistSq x y = || x - y ||^2`Squared 2-norm L1 norm!Euclidean norm"Lp norm (p > 0)# Infinity-norm$"Normalize w.r.t. p-norm (p finite)%Lp inner product (p > 0)& Reciprocal'Scale( unary dimension-checking bracket)!binary dimension-checking bracket*  !"#$%&'()*  !"#$%&'()* !"#$%&' ()     !"#$%&'(),(C) 2016 Marco Zocca, 2012-2015 Edward KmettGPL-3 (see LICENSE)zocca.marco gmail provisionalportableSafe 2Provides a test to see if a quantity is near zero.nearZero (1e-11 :: Double)FalsenearZero (1e-17 :: Double)TruenearZero (1e-5 :: Float)FalsenearZero (1e-7 :: Float)True%Determine if a quantity is near zero. Rounding rule Rounding rule Rounding ruleRound to respectively 0 or 1  a  1e-12  a  1e-6  a  1e-12  a  1e-6            None5Wrap a function with a null check, returning in MaybeIComponentwise tuple operations TODO : use semilattice properties insteadIComponentwise tuple operations TODO : use semilattice properties insteadinteger-indexed ziplist ", 2d arrays foldr over the results of a fmapstrict left foldindexed right foldSafe !"# !"# !"# Safe$Insert an element% Lookup a key&Lookup with default 0'APopulate an IM2 from a list of (row index, column index, value) (7Indexed left fold over an IM2, with general accumulator)Indexed left fold over an IM2*/Left fold over an IM2, with general accumulator+wInner indices become outer ones and vice versa. No loss of information because both inner and outer IntMaps are nubbed.,+Map over outer IM and filter all inner IM's-7Specialized filtering : keep only sub-diagonal elements.<List of (row, col) indices of (nonzero) subdiagonal elements/ Map over IM20Indexed map over IM21 Mapping keys$%&'()*+,-2.3/01456789$%&'()*+,-2.3/014$%&'()*+,-2.3/01456789(C) 2016 Marco ZoccaGPL-3 (see LICENSE)zocca.marco gmail provisionalportableNone9;<=DR63,`zeroSM m n` : Empty SpMatrix of size (m, n)5"`eye n` : identity matrix of rank n6Permutation matrix from a (possibly incomplete) list of row swaps starting from row 0 e.g. `permutationSM 5 [1,3]` first swaps rows (0, 1) and then rows (1, 3) :  0,1,0,0,0 0,0,0,1,0 0,0,1,0,0 1,0,0,0,0 0,0,0,0,17yPermutation matrix from a (possibly incomplete) list of row pair swaps e.g. `permutPairs 5 [(2,4)]` swaps rows 2 and 4 :  1,0,0,0,0 0,1,0,0,0 0,0,0,0,1 0,0,0,1,0 0,0,1,0,089`mkSubDiagonal n o xx` creates a square SpMatrix of size n with xx on the oth subdiagonal9DInsert an element in a preexisting Spmatrix at the specified indices:?Add to existing SpMatrix using data from list (row, col, value);:Create new SpMatrix using data from list (row, col, value)<ECreate new SpMatrix assuming contiguous, 0-based indexing of elements=HPopulate list with SpMatrix contents and populate missing entries with 0?`Looks up an element in the matrix with a default (if the element is not found, zero is returned)@3Zero-default lookup, infix form (no bound checking)`Looks up an element in the matrix with a default (if the element is not found, zero is returned)AIndexed filtering functionBbDiagonal, subdiagonal, superdiagonal partitions of a SpMatrix (useful for writing preconditioners)CbDiagonal, subdiagonal, superdiagonal partitions of a SpMatrix (useful for writing preconditioners)DbDiagonal, subdiagonal, superdiagonal partitions of a SpMatrix (useful for writing preconditioners)EfExtract a submatrix given the specified index bounds, rebalancing keys with the two supplied functionsF^Extract a submatrix given the specified index bounds NB : subtracts (i1, j1) from the indicesG\Extract a submatrix given the specified index bounds NB : submatrix indices are _preserved_JExtract whole columnK!Extract column within a row rangeL1Extract column within a row range, rebalance keysM.Are the supplied indices within matrix bounds?NIs the matrix square?OIs the matrix diagonal?P%Is the matrix lower/upper triangular?Q%Is the matrix lower/upper triangular?R,Is the matrix orthogonal? i.e. Q^t ## Q == IS/Data in internal representation (do not export)T#(Number of rows, Number of columns)U&Number of rows times number of columnsVNumber of rowsWNumber of columns_Vertical stacking`Vertical stackingaHorizontal stackingbHorizontal stackingdLeft fold over SpMatrixeIndexed left fold over SpMatrixfCount sub-diagonal nonzerosg`Filter the index subset that lies below the diagonal (used in the QR decomposition, for example)iSparsify an SpMatrixj3Round almost-0 and almost-1 to 0 and 1 respectivelyk0Swap two rows of a SpMatrix (bounds not checked)l.Swap two rows of a SpMatrix (bounds checked) mtransposeSM : Matrix transposerRemoves all elements x for which `| x | <= eps`)sRemoves all elements x for which `| x | <= eps`)tA^T BuA B^Tv Contract row i of A with column j of B up to an index nD, i.e. summing over repeated indices: Aij Bjk , for j in [0 .. n] w.ces are sparse containers too, i.e. any specific component may be missing (so it is assumed to be 0)z.-es are maps between finite-dimensional spaces{.Fes form a ring, in that they can be added and possess a zero element U*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~M*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvU./012~}|{zyxw3456789:;<=>?@ADCBEFGHIJKLMNOPQRSTUVW*+,-XYZ[\]^_`abcdefghijklmnopqrstuvO*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~(C) 2016 Marco ZoccaGPL-3 (see LICENSE)zocca.marco gmail provisionalportableNone9;<=DR#SpVector sparsityNumber of nonzeros*Empty sparse vector (length n, no entries)"Singleton sparse vector (length 1)Canonical basis vector in R^nQcreate a sparse vector from an association list while discarding all zero entries3", from logically dense array (consecutive indices)>Create new sparse vector, assumin 0-based, contiguous indexing^Map a function over a range of indices and filter the result (indices and values) to fit in a n-long SpVector0", using just the integer bounds of the interval_one-hot encoding : `oneHotSV n k` produces a SpVector of length n having 1 at the k-th positionDENSE vector of `1`sDENSE vector of `0`s3Populate a SpVector with the contents of a Vector. jPopulate a Vector with the entries of a SpVector, discarding the indices (NB: loses sparsity information).Z- | Populate a Vector with the entries of a SpVector, replacing the missing entries with 0insert element x at index i in a preexisting SpVectorTo dense list (default = 0)Indexed fold over SpVectorLookup an index in a SpVector7Lookup an index, return a default value if lookup fails8Lookup an index in a SpVector, returns 0 if lookup fails Tail elements Head element2Keep the first n components of the SpVector (like : for lists)LDiscard the first n components of the SpVector and rebalance the keys (like ; for lists)2Keep the first n components of the SpVector (like : for lists)Concatenate two sparse vectorsFilterIndexed filter*Generate an arbitrary (not random) vector u such that `v dot u = 0`Since 6s form a Hilbert space, we can define a norm for them Hs form a Hilbert space, in that we can define an inner product over thembs are sparse containers too, i.e. any specific component may be missing (so it is assumed to be 0) s are finite-dimensional vectors@s form a vector space because they can be multiplied by a scalar3'30 None<&Compressed Row Storage specification : Ahttp://netlib.org/utk/people/JackDongarra/etemplates/node373.html$The compressed row storage (CRS) format puts the subsequent nonzeros of the matrix rows in contiguous memory locations. Assuming we have a nonsymmetric sparse matrix $A$, we create three vectors: one for floating point numbers (val) and the other two for integers (col_ind, row_ptr).zThe val vector stores the values of the nonzero elements of the matrix $A$ as they are traversed in a row-wise fashion.The col_ind vector stores the column indexes of the elements in the val vector, that is, if val(k)=a_{i,j}, then col_ind(k)=j$.The row_ptr vector stores the locations in the val vector that start a row; that is, if val(k)=a_{i,j}, then row_ptr(i) <= k < row_ptr(i+1)<=>?@ABC<=>?@ABC<=>?@ABC(C) 2016 Marco ZoccaGPL-3 (see LICENSE)zocca.marco gmail provisionalportableNoneD2modify the size of a SpVector. Do not use directlyAInsert row , using the provided row index transformation functionInsert row CInsert column, using the provided row index transformation function Insert column#Outer product (all-with-all matrix)#Outer product (all-with-all matrix)AFill the diagonal of a SpMatrix with the components of a SpVectorpromote a SV to SM.Demote (n x 1) or (1 x n) SpMatrix to SpVectorRLookup a row in a SpMatrix; returns an SpVector with the row, if this is non-emptyExtract ith rowExtract jth columnGeneric extraction functionExtract ith row (dense)Extract jth columnExtract the diagonal&extract row interval (all entries between columns j1 and j2, INCLUDED, are returned) extractSubRow :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpVector a extractSubRow m i (j1, j2) = case lookupRowSM m i of Nothing -> zeroSV (ncols m) Just rv -> ifilterSV (j _ -> j >= j1 && j <= j2) rv", returning in Maybe extractSubRow :: SpMatrix a -> IxRow -> (Int, Int) -> Maybe (SpVector a) extractSubRow m i (j1, j2) = resizeSV (j2 - j1) . ifilterSV (j _ -> j >= j1 && j  =j2) <$ lookupRowSM m iExtract an interval of SpVector components, changing accordingly the resulting SpVector size. Keys are _not_ rebalanced, i.e. components are still labeled according with respect to the source matrix.>extract row interval, rebalance keys by subtracting lowest oneextract column intervalAextract column interval, rebalance keys by subtracting lowest oneMatrix-on-vectorMatrix-on-vector>Vector-on-matrix (FIXME : transposes matrix: more costly than  , I think)>Vector-on-matrix (FIXME : transposes matrix: more costly than  , I think)6Pack a V.Vector of SpVectors as columns of an SpMatrix&EDFGHIJKLM  !"#$%&'()  !"#*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuv%EDFGHIJKLMNone 9:;<=DR(Sparsify an SpVector+uses the R matrix from the QR factorization a vector xS uniquely defines an orthogonal plane; the Householder operator reflects any point v3 with respect to this plane: v' = (I - 2 x >< x) vNGivens coefficients (using stable algorithm shown in Anderson, Edward (4 December 2000). "Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem". LAPACK Working Note)zGivens method, row version: choose other row index i' s.t. i' is : * below the diagonal * corresponding element is nonzeroQR.C1 ) To zero out entry A(i, j) we must find row k such that A(k, j) is non-zero but A has zeros in row k for all columns less than j.LNB: the current version is quite inefficient in that: 1. the Givens' matrix G_iC is different from Identity only in 4 entries 2. at each iteration i we multiply G_i by the previous partial result M0. Since this corresponds to a rotation, and the NQ function already computes the value of the resulting non-zero component (output r@), `G_i ## M` can be simplified by just changing two entries of M3 (i.e. zeroing one out and changing the other into r).O,Returns a set of rows {k} that satisfy QR.C1PIs the k'th the first nonzero column in the row?Given a matrix A, returns a pair of matrices (Q, R) such that Q R = A, where Q is orthogonal and R is upper triangular. Applies Givens rotation iteratively to zero out sub-diagonal elements.`eigsQR n mm` performs n* iterations of the QR algorithm on matrix mm3, and returns a SpVector containing all eigenvalues`eigsRayleigh n mm` performs n0 iterations of the Rayleigh algorithm on matrix mm and returns the eigenpair closest to the initialization. It displays cubic-order convergence, but it also requires an educated guess on the initial eigenpair eigRayleigh :: Int -- max # iterations -> SpMatrix Double -- matrix -> (SpVector Double, Double) -- initial guess of (eigenvector, eigenvalue) -> (SpVector Double, Double) -- final estimate of (eigenvector, eigenvalue)_Given a positive semidefinite matrix A, returns a lower-triangular matrix L such that L L^T = A{Given a matrix A, returns a pair of matrices (L, U) where L is lower triangular and U is upper triangular such that L U = AQoApply a function over a range of integer indices, zip the result with it and filter out the almost-zero entries4Given a matrix A, a vector b and a positive integer n-, this procedure finds the basis of an order n Krylov subspace (as the columns of matrix Q), along with an upper Hessenberg matrix H, such that A = Q^T H Q. At the i`th iteration, it finds (i + 1) coefficients (the i`th column of the Hessenberg matrix H) and the (i + 1)`th Krylov vector.WPartition a matrix into strictly subdiagonal, diagonal and strictly superdiagonal parts+Used for Incomplete LU : remove entries in m" corresponding to zero entries in m2V`mSsor aa omega` : if `omega = 1` it returns the symmetric Gauss-Seidel preconditionerGDirect solver based on a triangular factorization of the system matrix.RNB in the computation of xi. we must rebalance the subrow indices because N does that too, in order to take the inner product with consistent index pairsS&Given a linear system `A x = b` where AR is an (m x m) real-valued matrix, the GMRES method finds an approximate solution xhatI such that the Euclidean norm of the residual `A xhat - b` is minimized. xhat is spanned by the order-n Krylov subspace of (A, b). In this implementation: 1) the Arnoldi iteration is carried out until numerical breakdown (therefore yielding _at_most_ mu Krylov basis vectors) 2) the resulting Hessenberg matrix is factorized in QR form 3) the Krylov-subspace solution yhatQ is found by backsubstitution 4) the approximate solution in the original space xhat6 is computed using the Krylov basis, `xhat = Q_n yhat`Many optimizations are possible, for example interleaving the QR factorization with the Arnoldi process (and employing an updating QR factorization which only requires one Givens' rotation at every update). Tone step of TFQMRUone step of BCGVone step of CGSWHiterate solver until convergence or until max # of iterations is reachedXone step of BiCGSTABYHiterate solver until convergence or until max # of iterations is reachedGLeast-squares approximation of a rectangular system of equaitons. Uses  \ for the linear solveKLinear solve with _deterministic_ starting vector (every component at 0.1)  \0 : linSolve using the BiCGSTAB method as defaultZ(transform state until a condition is met[0Keep a moving window buffer (length 2) of state xw to assess convergence, stop when either a condition on that list is satisfied or when max # of iterations is reached 0Keep a moving window buffer (length 2) of state x to assess convergence, stop when either a condition on that list is satisfied or when max # of iterations is reached (i.e. same thing as [& but this one runs in the State monad)Oiterate until convergence is verified or we run out of a fixed iteration budget\convergence check (FIXME)]run niter! iterations and append the state x to a list xs, stop when either the xs satisfies a predicate q or when the counter reaches 0", NO convergence check Dense SpMatrixDense SpVector Sparse SpMatrixSparse SpVectore^_`abcdefghijklmnopqrstuvwxyz{|}~NOPQRSTUVWXYZ[\]""@^_`abcdefghijklmnop qrstuvwxyz{|}~NOPQRSTUVWXYZ[\]    !"#$%&'()*+,-./01234556789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~      !"#$%&'()*+,-./01 2 3 4 5 6 7 8 9 : ; < = > ? @ A B C D E F GHIHJ K K L M N O P QRSTUVWXYZ[\]^_`abcdefghijkllmnoppqrstuuvwxyz{{|}~4sparse-linear-algebra-0.2.2.0-6BtTCE6r5KKHGNUOcJVHN2Numeric.LinearAlgebra.ClassData.Sparse.SpMatrixData.Sparse.SpVectorData.Sparse.CommonNumeric.LinearAlgebra.Sparse Numeric.EpsData.Sparse.UtilsData.Sparse.TypesData.Sparse.Internal.IntMap2Data.Sparse.Internal.CSR SpContainerScIxscInsertscLookup@@SetliftU2liftI2SparsespyHasDataHDDatadat FiniteDimFDSizedimNormednormHilbertdot VectorSpace.*Additivezero^+^one^*^negated^-^lerp hilbertDistSqnormSqnorm1norm2normP normInfty normalizedotLp reciprocalscalewithDimwithDim2SMInfosmNzsmSpySpMatrixSMsmDimsmDatasizeStrzeroSM mkDiagonaleye permutationSM permutPairsSM mkSubDiagonalinsertSpMatrix fromListSM' fromListSMfromListDenseSM toDenseListSMlookupSM lookupWD_SM@@!filterSMextractSubDiagextractSuperDiag extractDiagextractSubmatrixSMextractSubmatrixRebalanceKeysextractSubmatrixtakeRowstakeCols extractColSMextractSubColSMextractSubColSM_RK isValidIxSM isSquareSM isDiagonalSM isLowerTriSM isUpperTriSMisOrthogonalSMimmSMdimSMnelSMnrowsncolsinfoSMnzSMspySMnzRowbwMinSMbwMaxSM bwBoundsSM vertStackSM-=- horizStackSM-||- ifilterSMfoldlSMifoldlSMcountSubdiagonalNZSMsubdiagIndicesSM sparsifyIM2 sparsifySMroundZeroOneSMswapRows swapRowsSafe transposeSMmatScale normFrobeniusmatMat##matMatSparsified#~##^###^ contractSub$fSpContainerSpMatrixa$fSparseSpMatrixa$fHasDataSpMatrixa$fFiniteDimSpMatrix$fAdditiveSpMatrix $fSetSpMatrix$fFunctorSpMatrix$fShowSpMatrix $fEqSpMatrix $fEqSMInfo $fShowSMInfoSpVectorSVsvDimsvDataspySVnzSV sizeStrSVzeroSV singletonSVei mkSpVector mkSpVectorD mkSpVector1fromListDenseSVspVectorDenseIxspVectorDenseLoHi oneHotSVUoneHotSVonesSVzerosSV fromVectortoVector toVectorDenseinsertSpVector fromListSVtoListSV toDenseListSVifoldSVlookupSVlookupDefaultSV lookupDenseSVtailSVheadSVtakeSVdropSVconcatSVfilterSV ifilterSV orthogonalSV$fShowSpVector$fNormedSpVector$fHilbertSpVector$fSpContainerSpVectora$fSparseSpVectora$fHasDataSpVectora$fFiniteDimSpVector$fVectorSpaceSpVector$fAdditiveSpVector$fFoldableSpVector $fSetSpVector$fFunctorSpVector $fEqSpVectorprd insertRowWith insertRow insertColWith insertCol outerProdSV>< diagonalSMsvToSMtoSV lookupRowSM extractRow extractColextractVectorDenseWithextractRowDenseextractColDenseextractDiagDense extractSubRowextractSubRow_RK extractSubColextractSubCol_RKmatVec#>vecMat<#fromCols$fPrintDenseSpMatrix$fPrintDenseSpVectorLinSolveMethodGMRES_CGNE_TFQMR_BCG_CGS_ BICGSTAB_ sparsifySVconditionNumberSMhhMathhReflgivensqreigsQR eigRayleighcholluarnoldidiagPartitionsilu0mSsorluSolvepinvlinSolve<\>modifyInspectNdiffSqLuntilConverged runAppendN' randArrayrandMatrandVec randSpMat randSpVec$fShowBICGSTAB $fShowCGS $fShowBCG $fShowTFQMR $fShowCGNE$fEqCGNE $fEqTFQMR$fEqBCG$fEqCGS $fEqBICGSTAB$fEqLinSolveMethod$fShowLinSolveMethodEpsilonnearZero almostZero almostOneisNz roundZeroOne$fEpsilonCDoublebaseGHC.Numabsghc-prim GHC.Classes<=$fEpsilonCFloat$fEpsilonDouble$fEpsilonFloat withDefault roundZeroroundOne with2DefaultsharnessmaxTupminTup denseIxArray denseIxArray2foldrMap foldlStrictifoldrUBLBinBounds inBounds2 inBounds0 inBounds02head'tail'IxColIxRowColsRows insertIM2 lookupIM2 lookupWD_IM fromListIM2 ifoldlIM2' ifoldlIM2foldlIM2 transposeIM2 ifilterIM2 filterSubdiagsubdiagIndicesmapIM2imapIM2 mapKeysIM2countSubdiagonalNZrpairs mapColumnIM2$fNormedIntMap$fHilbertIntMap$fVectorSpaceIntMap$fAdditiveIntMap $fSetIntMapGHC.Listtakedrop CsrMatrixcsrVal csrColInd csrRowPtrcsrNnzcsrNrowscsrNcolsresizeSV PrintDense mapKeysSV showNonZero toDenseRowtoDenseRowClipnewline printDenseSMtoDenseListClip printDenseSV givensCoef candidateRowsfirstNonZeroColumn onRangeSparse triUpperSolvegmres tfqmrStepbcgStepcgsStepcgs bicgstabStepbicgstab modifyUntil loopUntilAccnormDiffConverged runAppendNBICGSTAB _xBicgstab _rBicgstab _pBicgstabCGS_x_r_p_uBCG_xBcg_rBcg_rHatBcg_pBcg_pHatBcgTFQMR_xTfq_wTfq_uTfq_vTfq_dTfq_mTfq_tauTfq _thetaTfq_etaTfq_rhoTfq _alphaTfqCGNE_xCgne_rCgne_pCgnehypotsignhhV triLowerSolvecgneStepcgnetfqmrbcgmeanlnorm2l