pardiso

Calculates the solution of a set of sparse linear equations with multiple right-hand sides.

Syntax

Fortran:

call pardiso (pt, maxfct, mnum, type, phase, n, a, ia, ja, perm, nrhs, iparm, msglvl, b, x, error)

pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);

Interface:

SUBROUTINE pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, perm, nrhs, iparm, msglvl, b, x, error)

PARDISO_POINTER_TYPE pt (64)

INTEGER maxfct, mnum, mtype, phase, n, nrhs, error, ia(*), ja(*), perm(*), iparm(*)

PARDISO_DATA_TYPE a(*), b(n,nrhs), x(n,nrhs)

PARDISO_POINTER_TYPE can be INTEGER for 32-bit architectures, or INTEGER*8 for 64-bit architectures.

PARDISO_DATA_TYPE can be one of the following types:

DOUBLE PRECISION - for real types of matrices (mtype=1, 2, -2 and 11) and for double precision pardiso (iparm(28)=0)

REAL - for real types of matrices (mtype=1, 2, -2 and 11) and for single precision pardiso (iparm(28)=1)

DOUBLE COMPLEX - for complex types of matrices (mtype=3, 6, 13, 14 and -4) and for double precision pardiso (iparm(28)=0)

COMPLEX - for complex types of matrices (mtype=3, 6, 13, 14 and -4) and for single precision pardiso *(iparm(28)=1)

Include Files

The FORTRAN 77 interfaces are specified in the mkl_pardiso.f77 include file, the Fortran 90 interfaces are specified in the mkl_pardiso.f90 include file, and the C interfaces are specified in the mkl_pardiso.h include file.

Description

The routine pardiso calculates the solution of a set of sparse linear equations
A*X = B
with multiple right-hand sides, using a parallel LU, LDL or LL^T factorization, where A is an n-by-n matrix, and X and B are n-by-nrhs matrices. pardiso performs the following analysis steps depending on the structure of the input matrix A.

Symmetric Matrices: The solver first computes a symmetric fill-in reducing permutation P based on either the minimum degree algorithm [Liu85] or the nested dissection algorithm from the METIS package [Karypis98] (included with Intel MKL), followed by the parallel left-right looking numerical Cholesky factorization [Schenk00-2] of PAP^T = LL^T for symmetric positive-definite matrices, or PAP^T = LDL^T for symmetric indefinite matrices. The solver uses diagonal pivoting, or 1x1 and 2x2 Bunch and Kaufman pivoting for symmetric indefinite matrices, and an approximation of X is found by forward and backward substitution and iterative refinements.

The coefficient matrix is perturbed whenever numerically acceptable 1x1 and 2x2 pivots cannot be found within the diagonal supernode block. One or two passes of iterative refinements may be required to correct the effect of the perturbations. This restricting notion of pivoting with iterative refinements is effective for highly indefinite symmetric systems. Furthermore the accuracy of this method is for a large set of matrices from different applications areas as accurate as a direct factorization method that uses complete sparse pivoting techniques [Schenk04].

Another possibility to improve the pivoting accuracy is to use symmetric weighted matching algorithms. These methods identify large entries in the coefficient matrix A that, if permuted close to the diagonal, permit the factorization process to identify more acceptable pivots and proceed with fewer pivot perturbations. The methods are based on maximum weighted matchings and improve the quality of the factor in a complementary way to the alternative idea of using more complete pivoting techniques.

The inertia is also computed for real symmetric indefinite matrices.

Structurally Symmetric Matrices: The solver first computes a symmetric fill-in reducing permutation P followed by the parallel numerical factorization of PAP^T = QLU^T. The solver uses partial pivoting in the supernodes and an approximation of X is found by forward and backward substitution and iterative refinements.

Unsymmetric Matrices: The solver first computes a non-symmetric permutation P_MPS and scaling matrices D_r and D_c with the aim to place large entries on the diagonal that enhances reliability of the numerical factorization process [Duff99]. In the next step the solver computes a fill-in reducing permutation P based on the matrix P_MPSA + (P_MPSA)_T followed by the parallel numerical factorization

QLUR = PP_MPSD_rAD_cP

with supernode pivoting matrices Q and R. When the factorization algorithm reaches a point where it cannot factorize the supernodes with this pivoting strategy, it uses a pivoting perturbation strategy similar to [Li99]. The magnitude of the potential pivot is tested against a constant threshold of alpha = eps*||A2||_inf , where eps is the machine precision, A2 = P*P_MPS*D_r*A*D_c*P, and ||A2||_inf is the infinity norm of the scaled and permuted matrix A. Any tiny pivots encountered during elimination are set to the sign (l_II)*eps*||A2||_inf, which trades off some numerical stability for the ability to keep pivots from getting too small. Although many failures could render the factorization well-defined but essentially useless, in practice the diagonal elements are rarely modified for a large class of matrices. The result of this pivoting approach is that the factorization is, in general, not exact and iterative refinement may be needed.

Direct-Iterative Preconditioning for Unsymmetric Linear Systems. The solver also allows for a combination of direct and iterative methods [Sonn89] to accelerate the linear solution process for transient simulation. Most of applications of sparse solvers require solutions of systems with gradually changing values of the nonzero coefficient matrix, but the same identical sparsity pattern. In these applications, the analysis phase of the solvers has to be performed only once and the numerical factorizations are the important time-consuming steps during the simulation. PARDISO uses a numerical factorization A = LU for the first system and applies the factors L and U for the next steps in a preconditioned Krylow-Subspace iteration. If the iteration does not converge, the solver automatically switches back to the numerical factorization. This method can be applied to unsymmetric matrices in PARDISO and the user can select the method using only one input parameter. For further details see the parameter description (iparm(4), iparm(20)).

Single and Double Precision Computations

PARDISO solves tasks using single or double precision. Each precision has its own pros and cons. Double precision variables have more digits to store value, so solver uses more memory for keeping data. But this mode allows to solve matrices with better accuracy, and input matrices can have big condition numbers.

Single precision variables have fewer digits to store values, so solver uses less memory than in the double precision mode. Additionally this mode usually takes less time. But as computations are made more roughly, only numerically stable process can be made using single precision.

Separate Forward and Backward Substitution.

The solver execution step ( see parameter phase =33 below) can be divided into two or three separate substitutions: forward, backward and diagonal (if possible). This separation can be explained on the examples of solving system with different matrix types.

Real symmetric positive definite matrix A (mtype = 2) is factorized by PARDISO as A = L*L^T . In this case solution of the system A*x=b can be found as sequence of substitutions: L*y=b (forward substitution, phase =331) and L^T*x=y (backward substitution, phase =333).

Real unsymmetric matrix A (mtype = 11) is factorized by PARDISO as A = L*U . In this case solution of the system A*x=b can be found by the following sequence: L*y=b (forward substitution, phase =331) andU*x=y (backward substitution, phase =333).

Solving system with a real symmetric indefinite matrix A (mtype = -2) is slightly different from above cases. PARDISO factorizes this matrix as A=LDL^T, and solution of the system A*x=b can be calculated as the following sequence of substitutions: L*y=b (forward substitution, phase =331) s: D*v=y (diagonal substitution, phase =332) and, finally L^T*x=v (backward substitution, phase =333). Diagonal substitution makes sense only for indefinite matrices (mtype = -2, -4, 6). For matrices of other types solution can be found as described in the first two examples.

Note

Number of refinement steps (iparm(8)) must be set to zero if solution is calculated step by step (phase = 331, 332, 333), otherwise PARDISO produces wrong result.

Note

Different pivoting (iparm(21)) produces different LDL^T factorization. Therefore results of forward, diagonal and backward substitutions with diagonal pivoting can differ from results of the same steps with Bunch and Kaufman pivoting. Of course the final results of sequential execution of forward, diagonal and backward substitution are equal to the results of the full solving step (phase=33) regardless of the pivoting used.

The sparse data storage in PARDISO follows the scheme described in Sparse Matrix Storage Format section with ja standing for columns, ia for rowIndex, and a for values. The algorithms in PARDISO require column indices ja to be increasingly ordered per row and the presence of the diagonal element per row for any symmetric or structurally symmetric matrix. For unsymmetric matrices the diagonal elements are not necessary: they may be present or not.

PARDISO performs four tasks:

analysis and symbolic factorization
numerical factorization
forward and backward substitution including iterative refinement
termination to release all internal solver memory.

When an input data structure is not accessed in a call, a NULL pointer or any valid address can be passed as a place holder for that argument.

Note

This routine supports the Progress Routine feature (for the numeric factorization phase only). See Progress Function for details.

Input Parameters

pt

INTEGER

Array, DIMENSION (64)

On entry, solver internal data address pointer. These addresses are passed to the solver and all related internal memory management is organized through this pointer

Note

pt is an integer array with 64 entries. It is very important that the pointer is initialized with zero at the first call of PARDISO. After that first call you should never modify the pointer, as a serious memory leak can occur. The integer length should be 4-byte on 32-bit operating systems and 8-byte on 64-bit operating systems.

maxfct

INTEGER

Maximal number of factors with identical nonzero sparsity structure that the user would like to keep at the same time in memory. It is possible to store several different factorizations with the same nonzero structure at the same time in the internal data management of the solver. In most of the applications this value is equal to 1.

PARDISO can process several matrices with identical matrix sparsity pattern and store the factors of these matrices at the same time. Matrices with a different sparsity structure can be kept in memory with different memory address pointers pt.

mnum

INTEGER

Actual matrix for the solution phase. With this scalar you can define the matrix that you would like to factorize. The value must be: 1 ≤ mnum ≤ maxfct.

In most applications this value is equal to 1.

mtype

INTEGER

This scalar value defines the matrix type. The PARDISO solver supports the following matrices:

= 1: real and structurally symmetric matrix
= 2: real and symmetric positive definite matrix
= -2: real and symmetric indefinite matrix
= 3: complex and structurally symmetric matrix
= 4: complex and Hermitian positive definite matrix
= -4: complex and Hermitian indefinite matrix
= 6: complex and symmetric matrix
= 11: real and unsymmetric matrix
= 13: complex and unsymmetric matrix

This parameter influences the pivoting method.

phase

INTEGER

Controls the execution of the solver. Usually it is a two-digit integer ij (10I + j, 1≤I≤3, I≤j≤3 for normal execution modes). The I digit indicates the starting phase of execution, and j indicates the ending phase. PARDISO has the following phases of execution:

Phase 1: Fill-reduction analysis and symbolic factorization
Phase 2: Numerical factorization
Phase 3: Forward and Backward solve including iterative refinements

This phase can be divided into two or three separate substitutions: forward, backward, and diagonal (see above).
Termination and Memory Release Phase (phase≤ 0)

If a previous call to the routine has computed information from previous phases, execution may start at any phase. The phase parameter can have the following values:

phase: Solver Execution Steps
11: Analysis
12: Analysis, numerical factorization
13: Analysis, numerical factorization, solve, iterative refinement
22: Numerical factorization
23: Numerical factorization, solve, iterative refinement
33: Solve, iterative refinement
331: like phase=33, but only forward substitution
332: like phase=33, but only diagonal substitution
333: like phase=33, but only backward substitution
0: Release internal memory for L and U matrix number mnum
-1: Release all internal memory for all matrices

n

INTEGER

Number of equations in the sparse linear systems of equations A*X = B. Constraint: n > 0.

a

DOUBLE PRECISION - for real types of matrices (mtype=1, 2, -2 and 11) and for double precision PARDISO (iparm(28)=0)

REAL - for real types of matrices (mtype=1, 2, -2 and 11) and for single precision PARDISO (iparm(28)=1)

DOUBLE COMPLEX - for complex types of matrices (mtype=3, 6, 13, 14 and -4) and for double precision PARDISO (iparm(28)=0)

COMPLEX - for complex types of matrices (mtype=3, 6, 13, 14 and -4) and for single precision PARDISO (iparm(28)=1)

Array. Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja. The size of a is the same as that of ja and the coefficient matrix can be either real or complex. The matrix must be stored in compressed sparse row format with increasing values of ja for each row. Refer to values array description in Sparse Matrix Storage Format for more details.

Note

The non-zero elements of each row of the matrix A must be stored in increasing order. For symmetric or structural symmetric matrices, it is also important that the diagonal elements are available and stored in the matrix. If the matrix is symmetric, the array a is only accessed in the factorization phase, in the triangular solution and iterative refinement phase. Unsymmetric matrices are accessed in all phases of the solution process.

ia

INTEGER

Array, dimension (n+1). For I≤n, ia(I) points to the first column index of row I in the array ja in compressed sparse row format. That is, ia(I) gives the index of the element in array a that contains the first non-zero element from row I of A. The last element ia(n+1) is taken to be equal to the number of non-zero elements in A, plus one. Refer to rowIndex array description in Sparse Matrix Storage Format for more details. The array ia is also accessed in all phases of the solution process. Note that the row and columns numbers start from 1 by default, but indexing base for input matrices can be changed to C style indexing by iparm(35).

by default, but indexing base for input matrices can be changed to C style indexing by iparm(35)

ja

INTEGER

Array. ja(*) contains column indices of the sparse matrix A stored in compressed sparse row format. The indices in each row must be sorted in increasing order. The array ja is also accessed in all phases of the solution process. For symmetric and structurally symmetric matrices it is assumed that zero diagonal elements are also stored in the list of non-zero elements in a and ja. For symmetric matrices, the solver needs only the upper triangular part of the system as is shown for columns array in Sparse Matrix Storage Format.

perm

INTEGER

Array, dimension (n). Holds the permutation vector of size n. The array perm is defined as follows. Let A be the original matrix and B = P*A*P^T be the permuted matrix. Row (column) I of A is the perm(I) row (column) of B. The numbering of the array must start with 1 and must describe a permutation.

On entry, you can apply your own fill-in reducing ordering to the solver. The permutation vector perm is used by the solver if iparm(5) = 1. The second purpose of the array perm is to return the permutation vector calculated during fill-in reducing ordering stage of PARDISO to the user. The permutation vector is returned into the perm array if iparm(5) = 2.

nrhs

INTEGER

Number of right-hand sides that need to be solved for.

iparm

INTEGER

Array, dimension (64). This array is used to pass various parameters to PARDISO and to return some useful information after execution of the solver. If iparm(1) = 0, PARDISO fills iparm(2) through iparm(64)with default values and uses them.

Individual components of the iparm array are described below (some of them are described in the Output Parameters section).

iparm(1)- use default values.

If iparm(1) = 0, iparm(2) through iparm(64) are filled with default values, otherwise the user has to supply all values in iparm from iparm(2) to iparm(64).

iparm(2) - fill-in reducing ordering.

iparm(2) controls the fill-in reducing ordering for the input matrix.

If iparm(2) = 0, the minimum degree algorithm is applied [Li99].

If iparm(2) = 2, the solver uses the nested dissection algorithm from the METIS package [Karypis98].

If iparm(2) = 3, the parallel (OpenMP) version of the nested dissection algorithm is used. It can decrease the time of computations on multi-core computers, especially when the time of the PARDISO Phase 1 is significant for your task.

The default value of iparm(2) is 2.

Caution

You can control the parallel execution of the solver by explicitly setting MKL_NUM_THREADS. If less processors are available than specified, the execution may slow down instead of speeding up. If the variable MKL_NUM_THREADS is not defined, then the solver uses all available processors.

iparm(3)- currently is not used.

iparm(4) - preconditioned CGS.

This parameter controls preconditioned CGS [Sonn89] for unsymmetric or structural symmetric matrices and Conjugate-Gradients for symmetric matrices. iparm(4) has the form iparm(4)= 10*L+K

The K and L values have the meanings as follow.

Value of K: Description
0: The factorization is always computed as required by phase.
1: CGS iteration replaces the computation of LU. The preconditioner is LU that was computed at a previous step (the first step or last step with a failure) in a sequence of solutions needed for identical sparsity patterns.
2: CGS iteration for symmetric matrices replaces the computation of LU. The preconditioner is LU that was computed at a previous step (the first step or last step with a failure) in a sequence of solutions needed for identical sparsity patterns.

Value L:

The value L controls the stopping criterion of the Krylow-Subspace iteration:

eps_CGS = 10^-L is used in the stopping criterion

||dx_I|| / ||dx_I|| < eps_CGS

with ||dx_I|| = ||inv(L*U)*r_I|| and r_I is the residuum at iteration I of the preconditioned Krylow-Subspace iteration.

Strategy: A maximum number of 150 iterations is fixed by expecting that the iteration will converge before consuming half the factorization time. Intermediate convergence rates and residuum excursions are checked and can terminate the iteration process. If phase =23, then the factorization for a given A is automatically recomputed in these cases where the Krylow-Subspace iteration failed, and the corresponding direct solution is returned. Otherwise the solution from the preconditioned Krylow-Subspace iteration is returned. Using phase =33 results in an error message (error =4) if the stopping criteria for the Krylow-Subspace iteration can not be reached. More information on the failure can be obtained from iparm(20).

The default is iparm(4)=0, and other values are only recommended for an advanced user. iparm(4) must be greater or equal to zero.

Examples:

iparm(4): Description
31: LU-preconditioned CGS iteration with a stopping criterion of 1.0E-3 for unsymmetric matrices
61: LU-preconditioned CGS iteration with a stopping criterion of 1.0E-6 for unsymmetric matrices
62: LU-preconditioned CGS iteration with a stopping criterion of 1.0E-6 for symmetric matrices

iparm(5)- user permutation.

This parameter controls whether user supplied fill-in reducing permutation is used instead of the integrated multiple-minimum degree or nested dissection algorithms. Another possible use of this parameter is to control obtaining the fill-in reducing permutation vector calculated during the reordering stage of PARDISO.

This option may be useful for testing reordering algorithms, adapting the code to special applications problems (for instance, to move zero diagonal elements to the end P*A*P^T), or for using the permutation vector more than once for equal or similar matrices. For definition of the permutation, see description of the perm parameter.

If parm(5)=0 (default value), then the array perm is not used by PARDISO;

if parm(5)=1, then the user supplied fill-in reducing permutation is used;

if parm(5)=2, then PARDISO returns the permutation vector into the array perm.

iparm(6)- write solution on x.

If iparm(6) = 0 (default value), then the array x contains the solution and the value of b is not changed.

If iparm(6) = 1, then the solver stores the solution on the right-hand side b.

Note that the array x is always used. The default value of iparm(6) is 0.

iparm(8) - iterative refinement step.

On entry to the solve and iterative refinement step, iparm(8)must be set to the maximum number of iterative refinement steps that the solver performs. The solver does not perform more than the absolute value of iparm(8)steps of iterative refinement and stops the process if a satisfactory level of accuracy of the solution in terms of backward error is achieved.

If iparm(8)< 0, the accumulation of the residuum uses extended precision real and complex data types. Perturbed pivots result in iterative refinement (independent of iparm(8)=0) and the number of executed iterations is reported in iparm(7).

The solver automatically performs two steps of iterative refinements when perturbed pivots are obtained during the numerical factorization and iparm(8) = 0.

The number of performed iterative refinement steps is reported in iparm(7).

The default value for iparm(8) is 0.

iparm(9)

This parameter is reserved for future use. Its value must be set to 0.

iparm(10)- pivoting perturbation.

This parameter instructs PARDISO how to handle small pivots or zero pivots for unsymmetric matrices (mtype =11 or mtype =13) and symmetric matrices (mtype =-2, mtype =-4, or mtype =6). For these matrices the solver uses a complete supernode pivoting approach. When the factorization algorithm reaches a point where it cannot factorize the supernodes with this pivoting strategy, it uses a pivoting perturbation strategy similar to [Li99], [Schenk04].

The magnitude of the potential pivot is tested against a constant threshold of

alpha = eps*||A2||_inf,

where eps = 10^(-iparm(10)), A2 = P*P_MPS*D_r*A*D_c*P, and ||A2||_inf is the infinity norm of the scaled and permuted matrix A. Any tiny pivots encountered during elimination are set to the sign (l_II)*eps*||A2||_inf - this trades off some numerical stability for the ability to keep pivots from getting too small. Small pivots are therefore perturbed with eps = 10^(-iparm(10)).

The default value of iparm(10) is 13 and therefore eps = 1.0E-13 for unsymmetric matrices (mtype =11 or mtype =13).

The default value of iparm(10) is 8, and therefore eps = 1.0E-8 for symmetric indefinite matrices (mtype =-2, mtype =-4, or mtype =6).

iparm(11)- scaling vectors.

PARDISO uses a maximum weight matching algorithm to permute large elements on the diagonal and to scale the matrix so that the diagonal elements are equal to 1 and the absolute values of the off-diagonal entries are less or equal to 1. This scaling method is applied only to unsymmetric matrices (mtype =11 or mtype =13). The scaling can also be used for symmetric indefinite matrices (mtype =-2, mtype =-4, or mtype =6) when the symmetric weighted matchings are applied (iparm(13)= 1).

Use iparm(11) = 1 (scaling) and iparm(13) = 1 (matching) for highly indefinite symmetric matrices, for example, from interior point optimizations or saddle point problems. Note that in the analysis phase (phase=11) you must provide the numerical values of the matrix A in case of scaling and symmetric weighted matching.

The default value of iparm(11) is 1 for unsymmetric matrices (mtype =11 or mtype =13). The default value of iparm(11) is 0 for symmetric indefinite matrices (mtype =-2, mtype =-4, or mtype =6).

iparm(12)

This parameter is reserved for future use. Its value must be set to 0.

iparm(13) - improved accuracy using (non-)symmetric weighted matchings.

PARDISO can use a maximum weighted matching algorithm to permute large elements close the diagonal. This strategy adds an additional level of reliability to our factorization methods and can be seen as a complement to the alternative idea of using more complete pivoting techniques during the numerical factorization.

It is recommended to use iparm(11)=1 (scalings) and iparm(13)=1 (matchings) for highly indefinite symmetric matrices, for example from interior point optimizations or saddle point problems. It is also very important to note that in the analysis phase (phase =11) you must provide the numerical values of the matrix A in the case of scalings and symmetric weighted matchings.

The default value of iparm(13) is 1 for unsymmetric matrices (mtype =11 or mtype =13). The default value of iparm(13) is 0 for symmetric matrices (mtype =-2, mtype =-4, or mtype =6).

iparm(18)

If iparm(18)< 0 on entry, the solver reports the numbers of non-zero elements in the factors .

The default value of iparm(18)is -1.

iparm(19)- MFlops of factorization.

If iparm(19)< 0 on entry, the solver reports MFlop (1.0E6) that are necessary to factor the matrix A. This increases the reordering time.

The default value of iparm(19) is 0.

iparm(21) - pivoting for symmetric indefinite matrices.

iparm(21) controls the pivoting method for sparse symmetric indefinite matrices.

If iparm(21) = 0, then 1x1 diagonal pivoting is used.

If iparm(21) = 1, then 1x1 and 2x2 Bunch and Kaufman pivoting is used in the factorization process.

Note

Use iparm(11) = 1 (scaling) and iparm(13) = 1 (matchings) for highly indefinite symmetric matrices, for example from interior point optimizations or saddle point problems.

The default value of iparm(21) is 1. Bunch and Kaufman pivoting is available for matrices: mtype=-2, mtype=-4, or mtype=6.

iparm(24) - parallel factorization control.

If iparm(24) = 0 (default value), then PARDISO uses the previous algorithm for factorization.

If iparm(24) = 1, then PARDISO uses new two-level factorization algorithm. This algorithm generally improves scalability in case of parallel factorization on many threads (>= 8).

Note

Two-level factorization algorithm is enabled by default in the previous MKL releases for matrices mtype=11. If you see performance degradation for such matrix with the default value, set manually iparm(24)=1.

iparm(27) - matrix checker.

If iparm(27)=0 (default value), PARDISO does not check the sparse matrix representation.

If iparm(27)=1, then PARDISO check integer arrays ia and ja. In particular, PARDISO checks whether column indices are sorted in increasing order within each row.

iparm(28) - single or double precision of PARDISO.

iparm(28) switches PARDISO precision between single and double modes.

If iparm(28)=0, then the input arrays (matrix a, vectors x and b) and all internal arrays must be presented in double precision.

If iparm(28)=1, then the input arrays must be presented in single precision. In this case all internal computations are performed in single precision.

Depending on the sign of iparm(8), refinement steps can be calculated in quad/double precision for double precision accuracy, and in double/single precision for single precision accuracy.

Default value of iparm(28) is 0 - double precision.

Important

iparm(28) value is stored in the PARDISO handle between PARDISO calls, so the precision mode can be changed only on the solver's phase 1.

iparm(31) - partial solution for sparse right-hand sides and sparse solution.

This parameter controls the solution method if the right hand side contains a few nonzero components. It can be used if only few components of the solution vector are needed, or if you want to reduce computation cost at solver step. To use this option define the input permutation vector perm so that perm(i) = 1 means that the i-the component in the right-hand side is nonzero.

On entry, if iparm(31) =1, and perm(i) = 1, the i-th component in the solution vector is computed.

If iparm(31) =2, all components of the solution vector are computed.

In the last case the computation cost at solver step is reduced due to reduced forward solver step.

To use iparm(31) =2, the i-th component of the right hand side must be set to zero explicitly if perm(i) is not equal to 1.

The permutation vector perm must be present in all phases of Intel MKL PARDISO software. At the reordering step, the software overwrites the input vector perm by a permutation vector used by the software at the factorization and solver step. If m is the number of components such that perm(i) = 1, then the last m components of the output vector perm are a set of the indices i satisfying the condition perm(i) = 1 on input.

The default value of iparm(31) is 0.

Note

Turning on this option often increases time used by PARDISO for factorization and reordering steps, but it enables time reducing at solver step.

Important

This feature is available only for in-core version, that is you must set iparm(60) =0. Set the parameters iparm(8) (iterative refinement steps) and iparm(4) (preconditioned CGS) to 0 as well.

iparm(35) - C or Fortran style array indexing.

iparm(35) determines the indexing base for input matrices.

If iparm(35)=0, (default value) then PARDISO uses Fortran style indexing: first value is referenced as array element 1, otherwise PARDISO uses C style indexing: first value is referenced as array element 0.

iparm(60) - version of PARDISO.

iparm(60) controls what version of PARDISO - out-of-core (OOC) version or in-core version - is used. The OOC PARDISO can solve very large problems by holding the matrix factors in files on the disk. Because of that the amount of main memory required by OOC PARDISO is significantly reduced.

If iparm(60) is set to 0, then the in-core PARDISO is used. If iparm(60) is set to 2 - the OOC PARDISO is used. If iparm(60) is set to 1 - the in-core PARDISO is used if the total memory (in MBytes) needed for storing the matrix factors is less than the value of the environment variable MKL_PARDISO_OOC_MAX_CORE_SIZE (its default value is 2000), and OOC PARDISO is used otherwise.

The default value of iparm(60) is 0.

If iparm(60) is equal to 1 or 2, and the total peak memory needed for strong local arrays is more than MKL_PARDISO_OOC_MAX_CORE_SIZE, the program stops with error -9. In this case, increase MKL_PARDISO_OOC_MAX_CORE_SIZE.

OOC parameters can be set in the configuration file. You can set the path to this file and its name using environmental variable MKL_PARDISO_OOC_CFG_PATH and MKL_PARDISO_OOC_CFG_FILE_NAME.

Path and name are as follows:

<MKL_PARDISO_OOC_CFG_PATH >/< MKL_PARDISO_OOC_CFG_FILE_NAME> for Linux* OS, and

<MKL_PARDISO_OOC_CFG_PATH >\< MKL_PARDISO_OOC_CFG_FILE_NAME> for Windows* OS.

By default, the name of the file is pardiso_ooc.cfg and it is placed to the current directory.

All temporary data files can be deleted or stored when the calculations are completed in accordance with the value of the environmental variable MKL_PARDISO_OOC_KEEP_FILE. If it is set to 1 (default value) - all files are deleted, if it is set to 0 - all files are stored.

By default, the OOC PARDISO uses the current directory for storing data, and all work arrays associated with the matrix factors are stored in files named ooc_temp with different extensions. These default values can be changed by using the environmental variable MKL_PARDISO_OOC_PATH.

To set the environmental variables MKL_PARDISO_OOC_MAX_CORE_SIZE, MKL_PARDISO_OOC_KEEP_FILE, and MKL_PARDISO_OOC_PATH, create the configuration file with the following lines:

MKL_PARDISO_OOC_PATH = <path>\ooc_file

MKL_PARDISO_OOC_MAX_CORE_SIZE = N

MKL_PARDISO_OOC_KEEP_FILE = 0 (or 1)

where <path> - the directory for storing data, ooc_file - file name without extension, N - memory size in MBytes, do not set this value greater than the size of the available RAM (default value - 2000 MBytes).

Important

Maximum length of the path lines in the configuration files is 1000 characters.

Alternatively the environment variables can be set via command line:

export MKL_PARDISO_OOC_PATH = <path>/ooc_file

export MKL_PARDISO_OOC_MAX_CORE_SIZE = N

export MKL_PARDISO_OOC_KEEP_FILE = 0 (or 1)

for Linux* OS, and

set MKL_PARDISO_OOC_PATH = <path>\ooc_file

set MKL_PARDISO_OOC_MAX_CORE_SIZE = N

set MKL_PARDISO_OOC_KEEP_FILE = 0 (or 1)

for Windows* OS.

Note

The values specified in a command line have higher priorities - it means that if a variable is changed in the configuration file and in the command line, OOC PARDISO uses only value defined in the command line.

msglvl

INTEGER

Message level information. If msglvl = 0 then PARDISO generates no output, if msglvl = 1 the solver prints statistical information to the screen.

b

DOUBLE PRECISION - for real types of matrices (mtype=1, 2, -2 and 11) and for double precision PARDISO (iparm(28)=0)

REAL - for real types of matrices (mtype=1, 2, -2 and 11) and for single precision PARDISO (iparm(28)=1)

DOUBLE COMPLEX - for complex types of matrices (mtype=3, 6, 13, 14 and -4) and for double precision PARDISO (iparm(28)=0)

COMPLEX - for complex types of matrices (mtype=3, 6, 13, 14 and -4) and for single precision PARDISO (iparm(28)=1)

Array, dimension (n, nrhs). On entry, contains the right-hand side vectors, which are placed in memory contiguously, namely b(I+(k-1)×nrhs) must hold the i-th component of k-th right-hand side vector. Note that b is only accessed in the solution phase.

Output Parameters

pt

This parameter contains internal address pointers.

iparm

On output, some iparm values report useful information, for example, numbers of non-zero elements in the factors, and so on.

iparm(7)- number of performed iterative refinement steps.

The number of iterative refinement steps that are actually performed during the solve step.

iparm(14)- number of perturbed pivots.

After factorization, iparm(14) contains the number of perturbed pivots during the elimination process for mtype =11, mtype =13, mtype =-2, mtype =-4, or mtype =-6.

iparm(15) - peak memory symbolic factorization.

The parameter iparm(15) reports the total peak memory in KBytes that the solver needs during the analysis and symbolic factorization phase. This value is only computed in phase 1.

iparm(16) - permanent memory symbolic factorization.

The parameter iparm(16) reports the permanent memory in KBytes from the analysis and symbolic factorization phase that the solver needs in the factorization and solve phases. This value is only computed in phase 1.

iparm(17) - size of factors /memory numerical factorization and solution.

This parameter is computed in phase 1 and 2 with different meanings.

The parameter iparm(17) computed in phase 1 provides the user with an estimate of size of factors (KBytes). In OOC mode, this information gives an estimate of disk space for storing data required.

The parameter iparm(17) computed in phase 2 provides the user with the total double precision memory consumption (KBytes) of the solver for the factorization and solve phases.

Note that the total peak memory solver consumption for all phases is max(iparm(15), iparm(16)+iparm(17))

iparm(18) - number of non-zero elements in factors.

The solver reports the numbers of non-zero elements on the factors if iparm(18) < 0 on entry.

iparm(19) - MFlops of factorization.

The solver reports the number of operations in MFlop (1.0E6 operations) that are necessary to factor the matrix A if iparm(19) < 0 on entry.

iparm(20) - CG/CGS diagnostics.

The value is used to give CG/CGS diagnostics (for example, the number of iterations and cause of failure):

If iparm(20)> 0, CGS succeeded, and the number of iterations executed are reported in iparm(20).

If iparm(20 )< 0, iterations executed, but CG/CGS failed. The error report details in iparm(20) are of the form: iparm(20)= - it_cgs*10 - cgs_error.

If phase= 23, then the factors L, U are recomputed for the matrix A and the error flag error=0 in case of a successful factorization. If phase =33, then error = -4 signals the failure.

Description of cgs_error is given in the table below:

cgs_error: Description
1: - fluctuations of the residuum are too large
2: - ||dx_{max_it_cgs/2}|| is too large (slow convergence)
3: - stopping criterion is not reached at max_it_cgs
4: - perturbed pivots causes iterative refinement
5: - factorization is too fast for this matrix. It is better to use the factorization method with iparm(4)=0

iparm(22) - inertia: number of positive eigenvalues.

The parameter iparm(22) reports the number of positive eigenvalues for symmetric indefinite matrices.

iparm(23) - inertia: number of negative eigenvalues.

The parameter iparm(23) reports the number of negative eigenvalues for symmetric indefinite matrices.

iparm(30) - the number of equation where PARDISO detects zero or negative pivot for MTYPE=2 (real positive definite matrix) and MTYPE=4 (complex and Hermitian positive definite matrices). If the solver detects a zero or negative pivot for these matrix types, the factorization is stopped, PARDISO returns immediately with an error ( error = -4) and iparm(30) contains the number of the equation where the first zero or negative pivot is detected.

iparm(31) - iparm(34), iparm(36) - iparm(59), iparm(61) - iparm(64)

These parameters are reserved for future use. Their values must be set to 0.

b

On output, the array is replaced with the solution if iparm(6) = 1.

x

DOUBLE PRECISION - for real types of matrices (mtype=1, 2, -2 and 11) and for double precision PARDISO (iparm(28)=0)

REAL - for real types of matrices (mtype=1, 2, -2 and 11) and for single precision PARDISO (iparm(28)=1)

DOUBLE COMPLEX - for complex types of matrices (mtype=3, 6, 13, 14 and -4) and for double precision PARDISO (iparm(28)=0)

COMPLEX - for complex types of matrices (mtype=3, 6, 13, 14 and -4) and for single precision PARDISO (iparm(28)=1)

Array, dimension (n,nrhs). On output, contains solution if iparm(6)=0. Like vector b, x(i+(k-1)× nrhs) holds the i-th component of the k-th solution vector. Note that x is only accessed in the solution phase.

x(i+(k-1)× nrhs) holds the i-th component of the k-th solution vector.

error

INTEGER

The error indicator according to the below table:

error: Information
0: no error
-1: input inconsistent
-2: not enough memory
-3: reordering problem
-4: zero pivot, numerical factorization or iterative refinement problem
-5: unclassified (internal) error
-6: preordering failed (matrix types 11, 13 only)
-7: diagonal matrix is singular
-8: 32-bit integer overflow problem
-9: not enough memory for OOC
-10: problems with opening OOC temporary files
-11: read/write problems with the OOC data file

pardiso

Syntax

Include Files

Description

Note

Note

Note

Input Parameters

Note

Note

Caution

Note

Note

Important

Note

Important

Important

Note

Output Parameters

See Also