API

This data type is used to store a sparse matrix in the COO (or coordinate) format through the irn, jcn and val fields containing the row indices, column indices and values, respectively and the m, n and nz containing the number of rows, columns and nonzeros, respectively. qr_mumps uses a Fortran-style 1-based numbering and thus all row indices are expected to be between 1 and m and all the column indices between 1 and n. Duplicate entries are summed during the factorization, out-of-bound entries are ignored. The sym field is used to specify if the matrix is symmetric and positive definite (true) or not (false).

This type is used to set the parameters that control the behavior of a sparse factorization, to collect information about its execution (number of flops, memory consumpnion etc) and store the result of the factorization, namely, the factors with all the symbolic information needed to use them in the solve phase.

qrm_init(ncpu, ngpu)

This routine initializes qr_mumps and should be called prior to any other qr_mumps routine. This function is automatically called if you use qr_mumps precompiled with Yggdrasil.

qrm_init()

ncpu and ngpu are optional arguments and their default value are, respectively, 1 and 0.

Input Arguments :

• ncpu: number of working threads on CPU.
• ngpu: number of working threads on GPU.

qrm_finalize()

This routine finalizes qr_mumps and no other qr_mumps routine should be called afterwards.

spmat = qrm_spmat_init(A; sym=false)
spmat = qrm_spmat_init(m, n, rows, cols, vals; sym=false)

qrm_spmat_init!(spmat, A; sym=false)
spmat = qrm_spmat_init!(spmat, m, n, rows, cols, vals; sym=false)

This routine initializes a qrm_spmat type data structure from a sparseMatrixCSC.

Input Arguments :

In the first form,

• spmat: the qrm_spmat sparse matrix to be initialized.
• A : a Julia sparse matrix stored in either SparseMatrixCSC or SparseMatrixCOO format (see SparseMatricesCOO.jl for the second case).
• sym : a boolean to specify if the matrix is symmetric / hermitian (true) or unsymmetric (false).

In the second form, the matrix A is specified using

• m: the number of rows.
• n: the number of columns.
• rows: the array of row indices of nonzero elements.
• cols: the array of column indices of nonzero elements.
• vals: the array of values of nonzero elements.

qrm_spmat_destroy!(spmat)

This routine cleans up a qrm_spmat type data structure.

Input Argument :

• spfct: the sparse matrix to be destroyed.

spfct = qrm_spfct_init(spmat)

qrm_spfct_init!(spmat, spfct)

This routine initializes a qrm_spfct type data structure. This amounts to setting all the control parameters to the default values.

Input Arguments :

• spmat: the input matrix of type qrm_spmat.
• spfct: the sparse factorization object to be initialized.

qrm_spfct_destroy!(spfct)

This routine cleans up a qrm_spfct type data structure by deleting the result of a sparse factorization.

Input Argument :

• spfct: the sparse factorization object to be destroyed.

spfct = qrm_analyse(spmat; transp='n')

qrm_analyse!(spmat, spfct; transp='n')

This routine performs the analysis phase and updates spfct.

Input Arguments :

• spmat: the input matrix of type qrm_spmat.
• spfct: the sparse factorization object of type qrm_spfct.
• transp: whether the input matrix should be transposed or not. Can be either 't', 'c' or 'n'.

qrm_update!(spmat, A)
qrm_update!(spmat, vals)

This routine updates a qrm_spmat type data structure from a sparseMatrixCSC. spmat and A must have the same sparsity pattern.

Input Arguments :

In the first form,

• spmat: the qrm_spmat sparse matrix to be updated.
• A : a Julia sparse matrix stored in SparseMatrixCSC format.

In the second form,

• vals: the array of values of nonzero elements of A.

qrm_factorize!(spmat, spfct; transp='n')

This routine performs the factorization phase. It can only be executed if the analysis is already done.

Input Arguments :

• spmat: the input matrix of type qrm_spmat.
• spfct: the sparse factorization object of type qrm_spfct.
• transp: whether the input matrix should be transposed or not. Can be either 't', 'c' or 'n'.

x = qrm_solve(spfct, b; transp='n')

qrm_solve!(spfct, b, x; transp='n')

This routine solves the triangular system Rx = b or Rᵀx = b. It can only be executed once the factorization is done.

Input Arguments :

• spfct: the sparse factorization object resulting from the qrm_factorize! function.
• b: the right-hand side(s).
• x: the solution vector(s).
• transp: whether to solve for R, Rᵀ or Rᴴ. Can be either 't', 'c' or 'n'.

z = qrm_apply(spfct, b; transp='n')

qrm_apply!(spfct, b; transp='n')

This routine computes z = Qb, z = Qᵀb or z = Qᴴb in place and overwrites b. It can only be executed once the factorization is done.

Input Arguments :

• spfct: the sparse factorization object resulting from the qrm_factorize! function.
• b: the vector(s) to which Q, Qᵀ or Qᴴ is applied.
• transp: whether to apply Q, Qᵀ or Qᴴ. Can be either 't', 'c' or 'n'.

cp = qrm_spfct_get_cp(spfct)

Returns the column permutation.

Input Arguments :

• spfct: a sparse factorization object of type qrm_spfct.

rp = qrm_spfct_get_rp(spfct)

Returns the row permutation.

Input Arguments :

• spfct: a sparse factorization object of type qrm_spfct.

R = qrm_spfct_get_r(spfct)

Returns the R factor as a SparseMatrixCSC matrix.

Input Arguments :

• spfct: a sparse factorization object of type qrm_spfct.

x = qrm_spbackslash(spmat, b; transp='n')

qrm_spbackslash!(spmat, b, x; transp='n')

x = qrm_spposv(spmat, b)

qrm_spposv!(spmat, b, x)

This function can be used to solve a linear symmetric, positive definite problem. It is a shortcut for the sequence

x = b
qrm_analyse!(spmat, spfct; transp='n')
qrm_factorize!(spmat, spfct; transp='n')
qrm_solve!(spfct, x, x; transp='t')
qrm_solve!(spfct, x, x; transp='t')

Input Arguments :

• spmat: the input matrix.
• b: the right-hand side(s).
• x: the solution vector(s).

x = qrm_least_squares(spmat, b; transp='n')

qrm_least_squares!(spmat, b, x; transp='n')

This function can be used to solve a linear least squares problem

\[\min \|Ax − b\|_2\]

in the case where the input matrix is square or overdetermined. It is a shortcut for the sequence

qrm_analyse!(spmat, spfct; transp='n')
qrm_factorize!(spmat, spfct; transp='n')
qrm_apply!(spfct, b; transp='t')
qrm_solve!(spfct, b, x; transp='n')

Input Arguments :

• spmat: the input matrix.
• b: the ight-hand side(s).
• x: the solution vector(s).
• transp: whether to use A, Aᵀ or Aᴴ. Can be either 't', 'c' or 'n'.

x = qrm_least_squares_semi_normal(spmat, b)

qrm_least_squares_semi_normal!(spmat, b, x, z, Δx, y; transp='n')

This function can be used to solve a linear least squares problem

\[\min \|Ax − b\|_2\]

in the case where A is square or overdetermined. Contrary to qrm_least_squares!, this function allows to solve the problem without storing the Q-factor of the QR factorization of A.

It is a shortcut for the sequence

qrm_set(spfct, "qrm_keeph", 0)
qrm_analyse!(spmat, spfct, transp  = 'n')
qrm_factorize!(spmat, spfct, transp = 'n')
qrm_spmat_mv!(spmat, T(1), b, T(0), z, transp = 't')
qrm_solve!(spfct, z, y, transp = 't')
qrm_solve!(spfct, y, x, transp = 'n')
qrm_refine!(spmat, spfct, x, z, Δx, y)

Note that the Q-factor is not used in this sequence; only A and R.

Input Arguments

• spmat: the input matrix.
• spfct: a sparse factorization object of type qrm_spfct.
• b: the right-hand side(s).
• x: the solution vector(s).
• Δx: an auxiliary vector (or matrix if b and x are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as x and can be uninitialized when the function is called.
• z: an auxiliary vector (or matrix if b and x are matrices) used to store the value z = Aᵀb, the size of this vector (resp. matrix) is the same as x and can be uninitialized when the function is called.
• y: an auxiliary vector (or matrix if b and x are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as b and can be uninitialized when the function is called.
• transp: whether to use A, Aᵀ or Aᴴ. Can be either 't', 'c' or 'n'.

x = qrm_min_norm(spmat, b; transp='n')

qrm_min_norm!(spmat, b, x; transp='n')

This function can be used to solve a linear minimum norm problem

\[\min \|x\|_2 \quad s.t. \quad Ax = b\]

in the case where the input matrix is square or underdetermined. It is a shortcut for the sequence

qrm_analyse!(spmat, spfct; transp='t')
qrm_factorize!(spmat, spfct; transp='t')
qrm_solve!(spfct, b, x; transp='t')
qrm_apply!(spfct, x; transp='n')

Input Arguments :

• spmat: the input matrix.
• b: the right-hand side(s).
• x: the solution vector(s).
• transp: whether to use A, Aᵀ or Aᴴ. Can be either 't', 'c' or 'n'.

x = qrm_min_norm_semi_normal(spmat, b)

qrm_min_norm_semi_normal!(spmat, spfct, b, x, Δx, y; transp='n')

This function can be used to solve a linear minimum-norm problem

\[\min \|x\|_2 \quad s.t. \quad Ax = b\]

in the case where A is square or underdetermined. Contrary to qrm_min_norm!, this function allows to solve the problem without storing the Q-factor in the QR factorization of Aᵀ. It is a shortcut for the sequence

qrm_set(spfct, "qrm_keeph", 0)
qrm_analyse!(spmat, spfct, transp = 't')
qrm_factorize!(spmat, spfct, transp = 't')
qrm_solve!(spfct, b, Δx, transp = 't')
qrm_solve!(spfct, Δx, y, transp = 'n')
qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = 't')

Remark that the Q-factor is not used in this sequence but rather A and R.

Input Arguments :

• spmat: the input matrix.
• spfct: a sparse factorization object of type qrm_spfct.
• b: the right-hand side(s).
• x: the solution vector(s).
• Δx: an auxiliary vector (or matrix if x and b are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as x and can be uninitialized when the function is called.
• y: an auxiliary vector (or matrix if x and b are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as b and can be uninitialized when the function is called.
• transp: whether to use A, Aᵀ or Aᴴ. Can be either 't', 'c' or 'n'.

x_refined = qrm_refine(spmat, spfct, x, z)

qrm_refine!(spmat, spfct, x, z, Δx, y)

Given an approximate solution x of the linear system RᵀRx ≈ z where R is the R-factor of some QR factorization of size (m, n), compute a refined solution.

Input Arguments :

• spmat: the input matrix.
• spfct: a sparse factorization object of type qrm_spfct.
• x: the approximate solution vector or matrix, the size of this vector (resp. matrix) is n (resp. n×k).
• z: the RHS vector or matrix of the linear system, the size of this vector (resp. matrix) is n (resp. n×k).
• Δx: an auxiliary vector (or matrix if x and z are matrices) used to compute the refinement, the size of this vector (resp. matrix) is n (resp. n×k) and can be uninitialized when the function is called.
• y: an auxiliary vector (or matrix if x and z are matrices) used to compute the refinement, the size of this vector (resp. matrix) is m (resp. m×k) and can be uninitialized when the function is called.

qrm_update_shift_spmat!(shifted_spmat, α)

Given a shifted block matrix of the form (A √α), update the parameter α in the matrix.

Input Arguments

• shifted_spmat: the input matrix of the qrmshiftedspmat type. See qrm_shifted_spmat for more information.
• α ≥ 0: the regularization parameter (note the square root in the block matrix representation).

qrm_shifted_spmat = qrm_shift_spmat(spmat, α = 0)

Given a spmat structure representing some sparse matrix A, return the block matrix (A √α) as a qrm_shifted_spmat type. See qrm_shifted_spmat for more information. This can be especially useful when A is rank deficient, as choosing α > 0 acts as a regularization.

Input Arguments

• spmat: the input matrix
• α ≥ 0: the regularization parameter (note the square root in the block matrix representation).

This data type represents a "shifted" matrix. When one wants to solve a "regularized" problem of the form (AᵀA + αI)x = b, one can use a QR factorization of the block matrix QR = (A √α)ᵀ and note that RᵀR = AᵀA + αI. This type of problem is especially useful when A is poorly conditioned or rank deficient. It only contains a qrm_spmat matrix representing the above block matrix and the regularization parameter α.

qrm_golub_riley!(shifted_spmat, spfct, x, b, Δx, y; α = ϵm, max_iter = 50, tol = ϵm, transp = 'n')

This method implements the Golub-Riley iteration. Given a (possibly ill-conditionned or rank deficient) system Ax = b where A can have any shape m×n, compute x = A†b = Aᵀ(AAᵀ)†b where A† is the Moore-Penrose pseudoinverse of A.

Input Arguments :

• shifted_spmat: a qrm_shifted_spmat type representing the matrix (A √α) where α is a regularization parameter for the Golub-Riley iteration. See qrm_shift_spmat for more information.
• spfct: a sparse factorization object of type qrm_spfct.
• x: the approximate solution vector x := A†b = Aᵀy, the size of this vector is n+m, its values are overwritten by this function and it can be uninitialized when the function is called.
• b: the RHS vector of the linear system, the size of this vector is m.
• Δx: an auxiliary vector used to compute the solution, the size of this vector is n+m and can be uninitialized when the function is called.
• y: the vector used to compute y := (AAᵀ)†b, the size of this vector is n and can be uninitialized when the function is called.
• Δy: an auxiliary vector used to compute the solution, the size of this vector is n and can be uninitialized when the function is called.

The definition of the shifted_spmat may look odd at first sight. We use such a block matrix as an argument of qrm_golub_riley! because the QR factorization of this matrix has the property that RᵀR = AᵀA + αI which is the system we need to repeatedly solve from for the Golub-Riley iteration.

The Golub-Riley method works as follows:

• Given A and b, choose a fixed regularization parameter α > 0 and let x := 0.
• At each iteration, compute Δx := Aᵀ(AAᵀ + αI)⁻¹ (b - Ax), and update x := x + Δx.

This method has the advantage of only costing one QR-factorization of the qrm_shifted_spmat matrix. For more details, see A. Dax and L. Eldén, Approximating minimum norm solutions of rank-deficient least squares problems, Numerical Linear Algebra with Applications, 5, pp. 79-99, 1998, DOI 10.1002/(SICI)1099-1506(199803/04)5:2<79::AID-NLA126>3.0.CO;2-4

x = qrm_golub_riley(spmat, b)

Input Arguments :

• spmat: the input matrix of the ill-conditionned system Ax = b.
• b: the RHS of the ill-conditionned system.

qrm_spmat_mv!(spmat, alpha, x, beta, y; transp='n')

This subroutine performs a matrix-vector product of the type y = αAx + βy, y = αAᵀx + βy or y = αAᴴx + βy.

Input Arguments :

• spmat: the input matrix.
• alpha, beta : the α and β scalars
• x: the x vector(s).
• y: the y vector(s).
• transp: whether to multiply by A, Aᵀ or Aᴴ. Can be either 't', 'c' or 'n'.

qrm_spmat_nrm(spmat; ntype='f')

This routine computes the one-norm $\|A\|_1$, the infinity-norm $\|x\|_{\infty}$ or the two-norm $\|x\|_2$ of a matrix.

Input Arguments :

• spmat: the input matrix.
• ntype: the type of norm to be computed. It can be either 'i', '1' or 'f' for the infinity, one and Frobenius norms, respectively.

qrm_vecnrm!(x, nrm; ntype='2')

This routine computes the one-norm $\|x\|_1$, the infinity-norm $\|x\|_{\infty}$ or the two-norm $\|x\|_2$ of a vector.

Input Arguments :

• x: the x vector(s).
• nrm: the computed norm(s). If x is a matrix this argument has to be a vector and each of its elements will contain the norm of the corresponding column of x.
• ntype: the type of norm to be computed. It can be either 'i', '1' or '2' for the infinity, one and two norms, respectively.

nrm = qrm_vecnrm(x; ntype='2')

qrm_residual_norm!(spmat, b, x, nrm; transp='n')

This function computes the scaled norm of the residual $\frac{\|b - Ax\|_{\infty}}{\|b\|_{\infty} + \|x\|_{\infty} \|A\|_{\infty}}$, i.e., the normwise backward error.

Input Arguments :

• spmat: the input matrix.
• b: the right-hand side(s).
• x: the solution vector(s).
• nrm: the computed norm(s).
• transp: whether to use A, Aᵀ or Aᴴ. Can be either 't', 'c' or 'n'.

nrm = qrm_residual_norm(spmat, b, x; transp='n')

qrm_residual_orth!(spmat, r, nrm; transp='n')

Computes the quantity $\frac{\|A^T r\|_2}{\|r\|_2}$ which can be used to evaluate the quality of the solution of a least squares problem.

Input Arguments :

• spmat: the input matrix.
• r: the residual(s).
• nrm: the computed norm(s).
• transp: whether to use A, Aᵀ or Aᴴ. Can be either 't', 'c' or 'n'.

nrm = qrm_residual_orth(spmat, r; transp='n')

qrm_set(str, val)
qrm_set(spfct, str, val)

Set control parameters that define the behavior of qr_mumps.

Input Arguments :

• spfct: a sparse factorization object of type qrm_spfct.
• str: either a string or a symbol describing the parameter to set.
• val: the parameter value.

val = qrm_get(str)
val = qrm_get(spfct, str)

Returns the value of a control parameter or an information parameter.

Input Arguments :

• spfct: a sparse factorization object of type qrm_spfct.
• str: either a string or a symbol describing the parameter to get.