pardiso

Syntax

Fortran:

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

C:

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

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.

Supported Matrix Types.

The analysis steps performed by pardiso depends 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] (both 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.

Whenever numerically acceptable 1x1 and 2x2 pivots cannot be found within the diagonal supernode block, the coefficient matrix is perturbed. 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, for a large set of matrices from different applications areas, this method is as accurate as a direct factorization method that uses complete sparse pivoting techniques [Schenk04].

Another method of improving the pivoting accuracy is to use symmetric weighted matching algorithms. These algorithms 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. These algorithms 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 of placing large entries on the diagonal to enhance 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 factor 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 enables to use 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. You 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 benefits and drawbacks. Double precision variables have more digits to store value, so the solver uses more memory for keeping data. But this mode solves matrices with better accuracy, and input matrices can have large condition numbers.

Single precision variables have fewer digits to store values, so the 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 use 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 possible diagonal . This separation can be explained by the examples of solving systems with different matrix types.

A real symmetric positive definite matrix A (mtype = 2) is factored by PARDISO as A = L*L^T . In this case the 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).

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

Solving a system with a real symmetric indefinite matrix A (mtype = -2) is slightly different from the cases above. PARDISO factors this matrix as A=LDL^T, and the 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 a solution can be found as described in the first two examples.

Note

The number of refinement steps (iparm(8)) must be set to zero if a solution is calculated with separate substitutions (phase = 331, 332, 333), otherwise PARDISO produces the 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.

Sparse Data Storage.

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

Note

The presence of diagonal elements for symmetric matrices is not mandatory starting from the Intel MKL 10.3 beta release.

Caution

It's recommended to set explicitly zero diagonal elements for symmetric matrices because in the opposite case PARDISO creates internal copies of arrays ia, ja and a full of diagonal elements that requires additional memory and computational time. However, in general, memory and time overheads are not significant comparing to the memory and the time needed to factor and solve the matrix.

Input Parameters

Note

Parameters types in this section are specified in FORTRAN 77 notation. See PARDISO parameters in tabular form for detailed description of types of PARDISO parameters in C/Fortran 90 notations.

pt

INTEGER

Array, DIMENSION (64)

Pointer to the address of solver internal data. 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 do not modify the pointer, as a serious memory leak can occur. The integer length must be 4 bytes on 32-bit operating systems and 8 bytes on 64-bit operating systems.

maxfct

INTEGER

Maximum number of factors with identical nonzero sparsity structure that must be keep at the same time in memory. In most applications this value is equal to 1. 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.

pardiso can process several matrices with an identical matrix sparsity pattern and it can 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

Indicates the actual matrix for the solution phase. With this scalar you can define which matrix to factorize. The value must be: 1 ≤ mnum ≤ maxfct.

In most applications this value is 1.

mtype

INTEGER

Defines the matrix type, which influences the pivoting method. The PARDISO solver supports the following matrices:

1

real and structurally symmetric

2

real and symmetric positive definite

-2

real and symmetric indefinite

3

complex and structurally symmetric

4

complex and Hermitian positive definite

-4

complex and Hermitian indefinite

6

complex and symmetric

11

real and unsymmetric

13

complex and unsymmetric

phase

INTEGER

Controls the execution of the solver. Usually it is a two- or three-digit integer ij (10i + j, 1≤i≤3, i≤j≤3 for normal execution modes). The i digit indicates the starting phase of execution, 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 Storage Formats for the Direct Sparse Solvers 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 Storage Formats for the Direct Sparse Solvers for more details. The array ia is also accessed in all phases of the solution process.

Indexing of ia is one-based by default, but it can be changed to zero-based by setting the appropriate value to the parameter 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 structurally symmetric matrices it is assumed that diagonal elements, which are zero, 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 Storage Formats for the Direct Sparse Solvers.

Indexing of ja is one-based by default, but it can be changed to zero-based by setting the appropriate value to the parameter iparm(35).

perm

INTEGER

Array, dimension (n). Holds the permutation vector of size n. You can use it to apply your own fill-in reducing ordering to the solver. 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 permutation vector perm is used by the solver if iparm(5) = 1.

The array perm is also used to return the permutation vector calculated during fill-in reducing ordering stage. The permutation vector is returned into the perm array if iparm(5) = 2.

Indexing of perm is one-based by default, but it can be changed to zero-based by setting the appropriate value to the parameter iparm(35).

Note

The first elements of row, column and permutation are numbered as array elements 1 by default (Fortran style, or one based array indexing), but these first elements can be numbered as array elements 0 (C style, or zero based array indexing) by setting the appropriate value to the parameter iparm(35).

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 uses default values for iparm(2) through iparm(64).

The 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 you must set 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 PARDISO Phase 1 takes significant time.

The default value of iparm(2) is 2.

Caution

You can control the parallel execution of the solver by explicitly setting the environment variable MKL_NUM_THREADS. If fewer 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 structurally 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 follows.

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 residue at iteration i of the preconditioned Krylow-Subspace iteration.

Strategy: A maximum number of 150 iterations is fixed with the assumption that the iteration will converge before consuming half the factorization time. Intermediate convergence rates and residue excursions are checked and can terminate the iteration process. If phase =23, then the factorization for a given A is automatically recomputed in 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 is useful for testing reordering algorithms, adapting the code to special applications problems (for instance, to move zero diagonal elements to the end of P*A*P^T), or for using the permutation vector more than once for equal or similar matrices. For definition of the permutation, see the 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 in the array perm 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 in 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 residue 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 factor 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)).

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

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

    

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) - solving with conjugate or transposed matrix.

If iparm(12)= 0, PARDISO solves a linear system Ax = b (default value).

If iparm(12)= 1, PARDISO solves a conjugate system A^Hx = b based on the factorization of the matrix A.

If iparm(12)= 2, PARDISO solves a transposed system A^Tx = b based on the factorization of the matrix A.

Note

For real matrices conjugate and transposed terms are equivalent.

  

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.

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. 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) - numbers of non-zero elements in the factors.

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 the number of MFLOPS (1.0E6) that are necessary to factor the matrix A. Reporting this number 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.

This parameter selects the scheduling method for the parallel numerical factorization.

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

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

Note

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

  

iparm(25) - parallel forward/backward solve control.

If iparm(25) = 0 (default value), then PARDISO uses a parallel algorithm for the solve step.

If iparm(25) = 1, then PARDISO uses sequential forward and backward solve.

Important

This feature is available only for in-core version in the case of one right hand side.

  

iparm(27) - matrix checker.

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

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

  

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

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 or double precision for double precision accuracy, and in double or 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 during 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 also 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 or the i-th component in the solution vector is computed.

If iparm(31) =0 (default value), this option is disabled.

If iparm(31) =1, the right hand side must be sparse, and the i-th component in the solution vector is computed if perm(i) = 1. You can set perm(i) = 1 only if the i-th component of the right hand side is nonzero.

If iparm(31) =2, the right hand side must be sparse, all components of the solution vector are computed. perm(i) = 1 means that the i-th component of the right hand side is nonzero.

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

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

If iparm(31) =3, the right hand side can be of any type and you must set perm(i) = 1 to compute the i-th component in the solution vector.

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.

Note

Turning on this option often increases time used by PARDISO for factorization and reordering steps, but it enables time to be reduced for the solver step.

Important

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

  

iparm(32) - iparm(34) - these parameters are reserved for future use. Their values must be set to 0.

  

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: the first value is referenced as array element 0.

  

iparm(35) - iparm(59) - these parameters are reserved for future use. Their values must be set to 0.

  

iparm(60) - version of PARDISO.

iparm(60) controls what version of PARDISO - out-of-core (OC) version or in-core (IC) version - is used. The OC 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 OC PARDISO is significantly reduced.

If iparm(60) = 0 (default value), then IC PARDISO is used.

If iparm(60) = 1 - then IC PARDISO is used if the total memory of RAM (in megabytes) needed for storing the matrix factors is less than sum of two values of the environment variables: MKL_PARDISO_OOC_MAX_CORE_SIZE (its default value is 2000 MB) and MKL_PARDISO_OOC_MAX_SWAP_SIZE (its default value is 0 MB); otherwise OOC PARDISO is used. In this case amount of RAM used by OOC PARDISO can not exceed the value of MKL_PARDISO_OOC_MAX_CORE_SIZE.

If iparm(60) = 2 - then OOC PARDISO is used.

If iparm(60) is equal to 1 or 2, and the total peak memory needed for storing the 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 a 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), then all files are deleted, if it is set to 0, then 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_MAX_SWAP_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_MAX_SWAP_SIZE = K

MKL_PARDISO_OOC_KEEP_FILE = 0 (or 1)

where <path> is the directory for storing data, ooc_file is the file name without any extension, N is the maximum size of RAM in megabytes available for PARDISO (default value is 2000 MB), K is the maximum swap size in megabytes available for PARDISO (default value is 0 MB). Do not set N greater than the size of the RAM and K greater than the size of the swap.

Important

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

Alternatively the environment variables can be set via command line.

For Linux* OS:

export MKL_PARDISO_OOC_PATH = <path>/ooc_file

export MKL_PARDISO_OOC_MAX_CORE_SIZE = N

export MKL_PARDISO_OOC_MAX_CORE_SIZE = K

export MKL_PARDISO_OOC_KEEP_FILE = 0 (or 1)

For Windows* OS:

set MKL_PARDISO_OOC_PATH = <path>\ooc_file

set MKL_PARDISO_OOC_MAX_CORE_SIZE = N

set MKL_PARDISO_OOC_MAX_CORE_SIZE = K

set MKL_PARDISO_OOC_KEEP_FILE = 0 (or 1)

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.

 

Note

You can switch between IC and OOC modes after the reordering phase. There are some recommendations and limitations:

• Set iparm(60) before reordering phase to get better PARDISO performance.

• Two-level factorization algorithm is not supported in the OOC mode. If you set two-level algorithm in the OOC mode then PARDISO returns error -1.

• Switching between IC and OOC modes after reordering phase is not available in sequential mode. The program returns error -1.

 

iparm(61), iparm(62), iparm(64) - these parameters are reserved for future use. Their values must be set to 0.

  

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 vector/matrix B, which is placed in memory contiguously. The 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

(See also PARDISO parameters in tabular form.)

pt

This parameter contains internal address pointers.

iparm

On output, some iparm values report information such as the numbers of non-zero elements in the factors.

  

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 kilobytes 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 kilobytes 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.

The parameter iparm(17) provides the size in kilobytes of the total memory consumed by IC PARDISO for internal float point arrays. This parameter is computed in phase 1. See iparm(63) for the OOC mode.

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 MFLOPS (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 and 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 failure.

Description of cgs_error is given in the table below:

cgs_error

Description

1

fluctuations of the residue 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 the equation where PARDISO detects zero or negative pivot

If the solver detects a zero or negative pivot for matrix types mtype = 2 (real positive definite matrix) and mtype = 4 (complex and Hermitian positive definite matrices), 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(63) - size of the minimum OOC memory for numerical factorization and solution.

The parameter iparm(63) provides the size in kilobytes of the minimum memory required by OOC PARDISO for internal float point arrays. This parameter is computed in phase 1.

Total peak memory consumption of OOC PARDISO can be estimated as max(iparm(15), iparm(16)+iparm(63))

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). If iparm(6)=0 it contains solution vector/matrix X, which is placed contiguously in memory. The x(i+(k-1)× nrhs) element must hold the i-th component of the k-th solution vector. Note that x is only accessed in the solution phase.

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

reordering failed (matrix types 11 and 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