Table Of Contents

Previous topic

2. Getting started

Next topic

4. Advanced features

3. General use of z-Pares

3.1. Naming convention of z-Pares routines

In the use of z-Pares, the user only uses one z-Pares main subroutine other than zpares_init and zpares_finalize. z-Pares main subroutines are named in the form of

zpares_TXXXYYZZ

T={s | c | d | z} (this means that T is replaced by s, c, d, or z) specifies the data type (real or complex) of the matrix data and the precision. XXX={rci | dns | mps} specifies the interface which is implemented by the routine. YY={ge | sy | he} specifies the symmetry and/or positive definiteness of the input matrices. ZZ={gv | ev} specifies the kind of eigenvalue problem which is solved (standard or generalized).

The correponding precisions and the matrix data types of T={s | c | d | z} are shown in Table 3.1. Symbols _REAL_, _COMPLEX__, and _MAT_TYPE_ represent the data types and are used in the following explanation about the data types of the variables. The user should replace them to the actual types listed in Table 3.1.

Table 3.1: Precisions and data types. double precision may be replaced by real(8) or real*8 and complex(kind(0d0)) may be replaced by complex(8), complex*16, or double complex

T precision \(A,B\) _REAL_ type _COMPLEX_ type _MAT_TYPE_ type
s single real real complex real
c single complex real complex complex
d double real double precision complex(kind(0d0)) double precision
z double complex double precision complex(kind(0d0)) complex(kind(0d0))

The meanings of XXX={rci | dns | mps} are described in Table 3.2.

Table 3.2: Interfaces of z-Pares.

XXX Description
rci Reverse communication interface
dns Dense interface using LAPACK
mps Sparse CSR interface using MUMPS

The conditions that YY={ge | sy | he} assumes are shown in Table 3.3 .

Table 3.3: Matrix symmetry and positive definiteness.

YY Property of matrix \(A\) Property of matrix \(B\)
sy Real symmetric Real symmetric and positive definite
he Complex Hermitian Complex Hermitian and positive definite
ge Other than above  

Note here that, as LAPACK users expect, YY=sy is valid with T=s or T=d, and YY=he is valid with T=c or T=z.

The explanation of ZZ={gv | ev} is shown in Table 3.4.

Table 3.4: Problem types.

ZZ Description
ev Standard eigenvalue problem : \(Ax = \lambda x\)
gv Generalized eigenvalue problem : \(Ax = \lambda Bx\)

3.2. Input and output arguments in z-Pares routines

The form of the sequence of arguments of the z-Pares routines differ depending on XXX, YY, and ZZ.

If YY is ge, the sequence of input arguments is in the form of

zpares_TXXXgeZZ(prm, nrow_local, {XXX dependent sequence}, left, right, num_ev, eigval, X, res, info)

The descriptions of the arguments are shown in Table 3.5. {XXX dependent sequence} is an argument sequence which depends on XXX. The actual sequences are described in the subsections of this section.

Table 3.5: Common arguments for zpares_TXXXgeZZ.

Name Description Type Input or output
prm Derived type zpares_prm zpares_prm IN/OUT
nrow_local Number of rows of X and/or the working spaces. See Section 2.12.3 integer IN
left Position of left edge of the search circle. See Section 2.9. _COMPLEX_ IN
right Real part of position of right edge of the search circle. See Section 2.9. _REAL_ IN
num_ev Number of calculated eigenvalues integer (IN)/OUT
eigval Array of eigenvalues _COMPLEX_, dimension(*) OUT
X Array of eigenvectors _MAT_TYPE_, dimension(nrow_local,*) (IN)/OUT
res Array of relative residual norm _REAL_, dimension(*) OUT
info Diagnostic information integer OUT

If YY is sy or he, the sequence of input arguments is in the form of

zpares_TXXX{sy|he}ZZ(prm, nrow_local, {XXX dependent sequence}, emin, emax, num_ev, eigval, X, res, info)

The descriptions of the arguments are shown in Table 3.6.

Table 3.6: Common arguments for zpares_TXXXsyZZ and zpares_TXXXheZZ.

Name Description Type Input or output
prm Derived type zpares_prm zpares_prm IN/OUT
nrow_local Number of rows of X and/or the working spaces. See Section 2.12.3 integer IN
emin Minimum value of the search interval. See Section 2.9. _REAL_ IN
emax Maximum value of the search interval. See Section 2.9. _REAL_ IN
num_ev Number of calculated eigenvalues integer (IN)/OUT
eigval Array of eigenvalues _REAL_, dimension(*) OUT
X Array of eigenvectors _MAT_TYPE_, dimension(nrow_local,*) (IN)/OUT
res Array of relative residual norms _REAL_, dimension(*) OUT
info Diagnostic information integer OUT

Note that the argument num_ev is basically OUTPUT. It is the number of approximate eigenvalues that have been computed, but is NOT the number of eigenvalues that are required by the user. The number of eigenvalues that are required by the user should be taken care by the setting of prm%L, prm%Lmax and prm%M. prm%L should be greater than or equal to the maximum multiplicity of eigenvalue inside the contour, and prm%L * prm%M should be greater than or equal to the number of eigenvalues inside the contour including the multiplicity. prm%L <= prm%Lmax should be maintained.

num_ev and X also become input arguments when the user already has approximate eigenvectors and gives them to the z-Pares routine (prm%user_source == .true.). In this case:

  • num_ev should be the number of approximate eigenpairs
  • X should contain approximate eigenvectors in leading num_ev columns

Note that num_ev >= prm%L should be maintained and in this case the estimation of the eigenvalue count in prm%num_ev_est becomes meaningless.

Real non-symmetric case

Only in case of T= {s | d} and YY=ge (real non-symmetric matrices), eigenvalues and eigenvectors will be returned with conjugate pairs if they are complex value. The computed eigenvalues are returned in complex typed 1D array eigval. Complex conjugate pairs of eigenvalues appear consecutively and the eigenvalue with the positive imaginary part appears first. The eigenvectors are returned in real typed 2D array X. X is returned, and the computed eigenvectors \(\boldsymbol{x}_j\) of eigval(j) are represented in the following manner:

  • If eigval(j) has zero imaginary part, then \(\boldsymbol{x}_j=\) V(:,j)
  • If eigval(j) and eigval(j+1) form a complex conjugate pair, then \(\boldsymbol{x}_j=\) V(:,j) + i*V(:,j+1) and \(\boldsymbol{x}_{j+1}=\) V(:,j) - i*V(:,j+1), where i is the imaginary unit.

3.2.1. Arguments in reverse communication interface

If XXX is rci, the routine implements a reverse communication interface (RCI), For RCI routines, {XXX dependent sequence} is {z, mwork, cwork}. Table 3.7 describes the arguments.

Table 3.7: {XXX dependent sequence} for zpares_TrciYYZZ.

Name Description Type Input or output
z Quadrature point _COMPLEX_ OUT
mwork Working space for performing matrix-vector multiplecation _MAT_TYPE_ -
cwork Working space for solving linear system _COMPLEX_ -

3.2.2. Arguments in dense interface

If XXX is dns, the routine implements a dense interface. The dense routines assume the matrices are given in a dense format (2D array) and solves the inner linear systems by LAPACK routines. For the dense routines, parallel implementations of the higher-level parallism is only supported. Thus prm%low_comm is always MPI_COMM_SELF.

If YY=sy or YY=he, and ZZ=gv, {XXX dependent sequence} is {UPLO, A, LDA, B, LDB}. Table 3.8 describes the arguments.

Table 3.8: {XXX dependent sequence} for zpares_Tdnssygv and zpares_Tdnshegv.

Name Description Type Input or output
UPLO UPLO={'U'|'L'}: {Upper | Lower} triangle parts of matrices are stored in A and B character IN
A Array represents matrix \(A\) _MAT_TYPE_, dimension(LDA,*) IN
LDA Leading dimension of A integer IN
B Array represents matrix \(B\) _MAT_TYPE_, dimension(LDB,*) IN
LDB Leading dimension of B integer IN

If YY=ge, {XXX dependent sequence} is {A, LDA, B, LDB}. Table 3.9 describes the arguments.

Table 3.9: {XXX dependent sequence} for zpares_Tdnsgegv.

Name Description Type Input or output
A Array represents matrix \(A\) _MAT_TYPE_, dimension(LDA,*) IN
LDA Leading dimension of A integer IN
B Array represents matrix \(B\) _MAT_TYPE_, dimension(LDB,*) IN
LDB Leading dimension of B integer IN

For the standard eigenproblem routines (i.e. ZZ=ev), B and LDB should be omitted.

In the dense interfaces, argument nrow_local should be the matrix size \(n\), since the lower-level parallelism is not supported.

3.2.3. Arguments in sparse CSR MUMPS interface

If XXX is mps, the routine implements a sparse CSR MUMPS interface. The sparse CSR MUMPS routines assumes matrices are given in Compressed Sparse Row (CSR) format and solves inner linear systems by MUMPS. For the sparse CSR MUMPS routines, parallel implementations of the higher-level parallism is only supported. Thus prm%low_comm is always MPI_COMM_SELF.

For sparse CSR MUMPS routines, if ZZ=gv, {XXX dependent sequence} is {rowptr_A, colind_A, val_A, rowptr_B, colind_B, val_B}. Table 3.10 describes the arguments.

Table 3.10: {XXX specific sequence} for zpares_TmpsYYgv.

Name Description Type Input or output
rowptr_A Array of row pointers of \(A\) integer, dimension(*) IN
colind_A Array of column indices of \(A\) integer, dimension(*) IN
val_A Array of value of \(A\) _MAT_TYPE_, dimension(*) IN
rowptr_B Array of row pointers of \(B\) integer, dimension(*) IN
colind_B Array of column indices of \(B\) integer, dimension(*) IN
val_B Array of value of \(B\) _MAT_TYPE_, dimension(*) IN

For YY=sy and YY=he, only the entries of either upper or lower triangular parts of the matrices should be stored.

For the standard eigenproblem routines (i.e. YY=ev), rowptr_B, colind_B and val_B should be omitted.

In the sparse CSR routines, argument nrow_local should be the matrix size \(n\), since the lower-level parallelism is not supported.

3.3. Parameters in zpares_prm

The derived type zpares_prm plays a central role in the use of z-Pares. It is passed to the first argument in a z-Pares main routine. Components of zpares_prm indicate input parameters for fine-tuning of the eigensolver or detailed information about the result. In the following subsections, we desribe the components of zpares_prm.

3.3.1. Integer input parameters

Table 3.11 shows integer input parameters.

Table 3.11: Integer input parameters in the components of zpares_prm.

Component Description Condition Default value
N Number of quadrature points N > 0 and even value 32
L Number of source vectors L > 0 16
Lmax Maxmum value of L Lmax \(\le\) L 64
M Maximum moment degree 0 < M \(\le\) N 16
imax Maximum refinement iteration imax \(\ge\) 0 0
n_orth Number of iteration for orthonormalization step n_orth \(\ge\) 0 3
extract Extract method. See 3.4. ZPARES_EXTRACT_RR or ZPARES_EXTRACT_EM ZPARES_EXTRACT_RR
Lstep Step size for increasing L Lstep > 0 8
high_comm Higher-level communicator - MPI_COMM_SELF
low_comm Lower-level communicator - MPI_COMM_SELF
write_unit File unit for verbosing A valid number 6 (Standard output)
verbose Verbose level {0|1|2} 0..
quad_type See Section 4.1 for description ZPARES_QUAD_ELL_TRAP or ZPARES_QUAD_USER ZPARES_QUAD_ELL_TRAP

The symbols in Table 3.11 that start with ZPARES_ are module variables typed integer,parameter of the zpares module.

3.3.2. Double precision real input parameters

Table 3.12 shows double precision input real parameters. These double precision parameters are used even if the user use a single precsion routine.

Table 3.12: Double precision real input parameters in the components of zpares_prm.

Component Description Condition Default value
delta Threshold for determining numerical rank 0 \(\le\) delta \(\le\) 1 1d-12
spu_thres Threshold for trancating eigenvalue with respect to spurious indicator 0 \(\le\) spu_thres \(\le\) 1 1d-4
asp_ratio Aspect ratio of ellipse asp_ratio > 0 1.0

3.3.3. Logical input parameters

Table 3.13 shows logical input parameters.

Table 3.13: Logical input parameters in the componets of zpares_prm.

Component If .true. Default value
calc_res Calculate residual norm .true.
trim_out Trim eigenvalues that are outside of the circle .true.
trim_spu Trim eigenvalues whose corresponding spurious indicators are smaller than spu_thres .false.
trim_res Trim eigenvalues whose corresponding residual norms are larger than tol .false.
Hermitian Assume both \(A\) and \(B\) are Hermitian .false.
B_pos_def Assume \(B\) is positive definite .false.
sym_contour See Section 4.1 for description .false.

Note that Hermitian and B_pos_def treated as .true., if YY = {sy | he}.

3.3.4. Output parameters

Table 3.14 and Table 3.15 show integer output parameters and double precsion output parameters. For the single precision routines, the single precsion results are upcasted to double precision.

Table 3.14: Integer output parameters in the components of zpares_prm.

Component Description
num_ev_est Estimation of the number of eigenvalues in the interval or the circle
iter Number of refiment iterations
num_basis Number of basis vectors used for projection.
L Auto-tuned value of the source size

Table 3.15: Double precision, pointer typed output parameters in the components of zpares_prm.

Component Description
indi_spu(:) The i-th component is the indicator of spuriousness corresponds the i-th eigenvalue. The length is num_ev. See Section 3.7.
sigval(:) Singular values used for a low-rank approximation. The i-th component is the i-th largest singular value. The length is prm%num_basis

3.4. Extracting method for eigenpairs

z-Pares extracts eigenpairs from the approximate eigensubspace by two different ways:

  • Rayleigh-Ritz procedure
  • Hankel method

In the Rayleigh-Ritz procedure, linear combinations of the solutions of the linear systems used as the basis for the procedure. The Hankel method computes eigenpairs by solving a reduced generalized eigenvalue problem of (block) Hankel matrices. In our experience, the Rayleigh-Ritz procedure gives better accuracy than the Hankel method. But the Hankel method is still attractive because no matrix-vector multiplecation of \(A\) and \(B\) is needed if the residual is not needed ( prm%calc_res == .false. ). The Rayleigh-Ritz procedure is enabled if prm%extract == ZPARES_EXTRACT_RR (default) and the Hankel method is enabled if prm%extract == ZPARES_EXTRACT_EM.

3.5. Reverse communicaton interface in complex Hermitian case

Unless YY=he, the reverse communication loop is written in the form described in Section 2.11. If YY=he, the reverse communication loop should be written as following:

  do while ( prm%itask /= ZPARES_TASK_FINISH )
     call zpares_zrcihegv &
          (prm, nrow_local, z, mwork, cwork, emin, emax, num_ev, eigval, X, res, info)

     select case (prm%itask)
     case(ZPARES_TASK_FACTO)

        ! Here the user factorizes (z*B - A)
        ! At the next return from zpares_zrcihegv,
        ! prm%itask==ZPARES_TASK_SOLVE with the same z
        
     case(ZPARES_TASK_SOLVE)

        ! i = prm%ws; j = prm%ws+prm%nc-1
        ! Here the user solves (z*B - A) X = cwork(:,i:j) 
        ! The solution X should be stored in cwork(:,i:j)
        ! At the next return from zpares_zrcihegv,
        ! prm%itask==ZPARES_TASK_FACTO_H with the same z is returned

     case(ZPARES_TASK_FACTO_H)

        ! Here the user can factorize (z*B - A)^{H}
        ! In most cases, this task is unnecessary
        ! At the next return from zpares_zrcihegv,
        ! prm%itask==ZPARES_TASK_SOLVE_H with the same z is returned
        
     case(ZPARES_TASK_SOLVE_H)

        ! i = prm%ws; j = prm%ws+prm%nc-1
        ! Here the user solves (z*B - A) X = cwork(:,i:j)
        ! The solution X should be stored in cwork(:,i:j)

     case(ZPARES_TASK_MULT_A)

        ! iw = prm%ws; jw = prm%ws+prm%nc-1
        ! ix = prm%xs; jx = prm%xs+prm%nc-1
        ! Here the user performs matrix-vector multiplications:
        ! mwork(:,iw:jw) = A*X(:,ix:jx)

     case(ZPARES_TASK_MULT_B)
        
        ! iw = prm%ws; jw = prm%ws+prm%nc-1
        ! ix = prm%xs; jx = prm%xs+prm%nc-1
        ! Here the user performs matrix-vector multiplications:
        ! mwork(:,iw:jw) = B*X(:,ix:jx)

     end select
  end do

The symbols start with ZPARES_TASK are module variables typed integer,parameter of the zpares module.

3.6. Estimations of eigenvalue count

Along with eigenpairs z-Pares computes a stochastic estimation of the eigenvalue count. This value is used for auto-tuning for the parameter prm%L. The estimated eigenvalue count is returned in prm%num_ev_est. Note that the accuracy of the estimation is problem dependent. The error may become arbitrary large.

3.7. Indicator for spurious eigenvalues

z-Pares provides an indicator of the spuriousness of an eigenvalue. The indicators are returned in prm%indi_spu(:). prm%indi_spu(i) indicates spuriousness of eigval(i). A value of the indicator takes \([0,1]\). The smaller indicator, the more spurious the eigenvalue is. The spurious eigenvalues can be discarded using threshold value prm%spu_thres.

3.8. Error diagnostics in info

The output parameter info of a z-Pares main routine tells diagnostics information on exit. It takes values listed in Table 3.16.

Table 3.16: Diagnostics information

Value Description
ZPARES_INFO_SUCCESS Successful exit
ZPARES_INFO_INVALID_PARAM Invalid parameter
ZPARES_INFO_LAPACK_ERROR LAPACK error
ZPARES_INFO_MPI_ERROR MPI error