HSL_MI35: Sparse least squares: incomplete factorization preconditioner

HSL_MI35 computes an incomplete factorization preconditioner that may be used in the solution of weighted sparse linear least squares problems. Given an \(m \times n\) (\(m \ge n\)) sparse matrix \({A} = \{ a_{ij} \}\) and an \(m \times m\) diagonal matrix of weights \(W\), HSL_MI35 computes an incomplete factorization of the matrix \[C = A^TW^2A.\] The matrix \(C\) may be optionally reordered, scaled and, if necessary, shifted to avoid breakdown of the factorization. The computed lower triangular matrix \(L\) is such that \(LL^T\) is an incomplete factorization of \[\overline{C} = SQ^T C QS + \alpha I,\] where \(Q\) is a permutation matrix, \(S\) is a diagonal scaling matrix and \(\alpha\) is a non-negative shift. The incomplete factorization may be used for preconditioning when solving the normal equations \[A^TW^2Ax = A^TW^2b.\] A separate entry performs the preconditioning operation \[{y = Pz}\] where \({P} = (\overline{L} {\overline{L}}^T)^{-1}\) is the incomplete factorization preconditioner and \(\overline{L} = Q S^{-1} L\).

The incomplete factorization is based on a matrix decomposition of the form \[\overline{C} = LL^T + LR^T + RL^T - E,\] where \(L\) is lower triangular with positive diagonal entries, \(R\) is a strictly lower triangular matrix with small entries that is used to stabilize the factorization process, and \(E\) has the structure \[E = RR^T + F + F^T,\] where \(F\) is strictly lower triangle. \(F\) is discarded while \(R\) is used in the computation of \(L\) but is then discarded. The user controls the dropping of small entries from \(L\) and \(R\) and the maximum number of entries within each column of \(L\) and \(R\) (and thus the amount of memory for \(L\) and the intermediate work and memory used in computing the incomplete factorization).

The incomplete factorization may be computed with or without forming \(C\) explicitly; the advantage of the latter option being that less memory is required.