Man Linux: Main Page and Category List

NAME

       PSSYEVX - compute selected eigenvalues and, optionally, eigenvectors of
       a real symmetric matrix  A  by  calling  the  recommended  sequence  of
       ScaLAPACK routines

SYNOPSIS

       SUBROUTINE PSSYEVX( JOBZ, RANGE, UPLO, N, A, IA, JA, DESCA, VL, VU, IL,
                           IU, ABSTOL, M, NZ, W,  ORFAC,  Z,  IZ,  JZ,  DESCZ,
                           WORK,  LWORK,  IWORK,  LIWORK, IFAIL, ICLUSTR, GAP,
                           INFO )

           CHARACTER       JOBZ, RANGE, UPLO

           INTEGER         IA, IL, INFO, IU, IZ, JA, JZ, LIWORK, LWORK, M,  N,
                           NZ

           REAL            ABSTOL, ORFAC, VL, VU

           INTEGER         DESCA(  *  ), DESCZ( * ), ICLUSTR( * ), IFAIL( * ),
                           IWORK( * )

           REAL            A( * ), GAP( * ), W( * ), WORK( * ), Z( * )

           INTEGER         BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_, MB_,
                           NB_, RSRC_, CSRC_, LLD_

           PARAMETER       ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, CTXT_
                           = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, RSRC_  =  7,
                           CSRC_ = 8, LLD_ = 9 )

           REAL            ZERO, ONE, TEN, FIVE

           PARAMETER       (  ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 10.0E+0, FIVE
                           = 5.0E+0 )

           INTEGER         IERREIN, IERRCLS, IERRSPC, IERREBZ

           PARAMETER       ( IERREIN = 1, IERRCLS = 2, IERRSPC = 4, IERREBZ  =
                           8 )

           LOGICAL         ALLEIG, INDEIG, LOWER, LQUERY, QUICKRETURN, VALEIG,
                           WANTZ

           CHARACTER       ORDER

           INTEGER         CSRC_A,  I,  IACOL,  IAROW,  ICOFFA,  IINFO,  INDD,
                           INDD2,   INDE,   INDE2,   INDIBL,  INDISP,  INDTAU,
                           INDWORK,  IROFFA,   IROFFZ,   ISCALE,   ISIZESTEBZ,
                           ISIZESTEIN, IZROW, LALLWORK, LIWMIN, LLWORK, LWMIN,
                           MAXEIGS, MB_A, MB_Z, MQ0, MYCOL, MYROW,  NB,  NB_A,
                           NB_Z,  NEIG,  NN,  NNP,  NP0, NPCOL, NPROCS, NPROW,
                           NSPLIT, NZZ,  OFFSET,  RSRC_A,  RSRC_Z,  SIZEORMTR,
                           SIZESTEIN, SIZESYEVX

           REAL            ABSTLL,  ANRM,  BIGNUM,  EPS,  RMAX,  RMIN, SAFMIN,
                           SIGMA, SMLNUM, VLL, VUU

           INTEGER         IDUM1( 4 ), IDUM2( 4 )

           LOGICAL         LSAME

           INTEGER         ICEIL, INDXG2P, NUMROC

           REAL            PSLAMCH, PSLANSY

           EXTERNAL        LSAME, ICEIL, INDXG2P, NUMROC, PSLAMCH, PSLANSY

           EXTERNAL        BLACS_GRIDINFO,   CHK1MAT,    IGAMN2D,    PCHK2MAT,
                           PSELGET,   PSLARED1D,  PSLASCL,  PSORMTR,  PSSTEBZ,
                           PSSTEIN,  PSSYTRD,   PXERBLA,   SGEBR2D,   SGEBS2D,
                           SLASRT, SSCAL

           INTRINSIC       ABS, ICHAR, MAX, MIN, MOD, REAL, SQRT

PURPOSE

       PSSYEVX  computes selected eigenvalues and, optionally, eigenvectors of
       a real symmetric matrix  A  by  calling  the  recommended  sequence  of
       ScaLAPACK  routines.  Eigenvalues/vectors can be selected by specifying
       a range of values or a range of indices for the desired eigenvalues.

NOTES

       Each global data object  is  described  by  an  associated  description
       vector.   This  vector stores the information required to establish the
       mapping between an object element and  its  corresponding  process  and
       memory location.

       Let  A  be  a generic term for any 2D block cyclicly distributed array.
       Such a global array has an associated description vector DESCA.  In the
       following  comments,  the  character _ should be read as "of the global
       array".

       NOTATION        STORED IN      EXPLANATION
       --------------- -------------- --------------------------------------
       DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
                                      DTYPE_A = 1.
       CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
                                      the BLACS process grid A is distribu-
                                      ted over. The context itself is glo-
                                      bal, but the handle (the integer
                                      value) may vary.
       M_A    (global) DESCA( M_ )    The number of rows in the global
                                      array A.
       N_A    (global) DESCA( N_ )    The number of columns in the global
                                      array A.
       MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
                                      the rows of the array.
       NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
                                      the columns of the array.
       RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
                                      row of the array A is distributed.
       CSRC_A (global) DESCA( CSRC_ ) The process column over which the
                                      first column of the array A is
                                      distributed.
       LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
                                      array.  LLD_A >= MAX(1,LOCr(M_A)).

       Let K be the number of rows or columns of  a  distributed  matrix,  and
       assume  that  its  process grid has dimension p x q.  LOCr( K ) denotes
       the number of elements of K that a process  would  receive  if  K  were
       distributed over the p processes of its process column.
       Similarly, LOCc( K ) denotes the number of elements of K that a process
       would receive if K were distributed over the q processes of its process
       row.
       The  values  of  LOCr()  and LOCc() may be determined via a call to the
       ScaLAPACK tool function, NUMROC:
               LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
               LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
       An upper bound for these quantities may be computed by:
               LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
               LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A

       PSSYEVX assumes IEEE 754 standard compliant arithmetic.  To port  to  a
       system  which does not have IEEE 754 arithmetic, modify the appropriate
       SLmake.inc file to include the compiler switch -DNO_IEEE.  This  switch
       only affects the compilation of pslaiect.c.

ARGUMENTS

          NP = the number of rows local to a given process.
          NQ = the number of columns local to a given process.

       JOBZ    (global input) CHARACTER*1
               Specifies whether or not to compute the eigenvectors:
               = ’N’:  Compute eigenvalues only.
               = ’V’:  Compute eigenvalues and eigenvectors.

       RANGE   (global input) CHARACTER*1
               = ’A’: all eigenvalues will be found.
               = ’V’: all eigenvalues in the interval [VL,VU] will be found.
               = ’I’: the IL-th through IU-th eigenvalues will be found.

       UPLO    (global input) CHARACTER*1
               Specifies whether the upper or lower triangular part of the
               symmetric matrix A is stored:
               = ’U’:  Upper triangular
               = ’L’:  Lower triangular

       N       (global input) INTEGER
               The number of rows and columns of the matrix A.  N >= 0.

       A       (local input/workspace) block cyclic REAL array,
               global dimension (N, N),
               local dimension ( LLD_A, LOCc(JA+N-1) )

               On entry, the symmetric matrix A.  If UPLO = ’U’, only the
               upper triangular part of A is used to define the elements of
               the symmetric matrix.  If UPLO = ’L’, only the lower
               triangular part of A is used to define the elements of the
               symmetric matrix.

               On exit, the lower triangle (if UPLO=’L’) or the upper
               triangle (if UPLO=’U’) of A, including the diagonal, is
               destroyed.

       IA      (global input) INTEGER
               A’s global row index, which points to the beginning of the
               submatrix which is to be operated on.

       JA      (global input) INTEGER
               A’s global column index, which points to the beginning of
               the submatrix which is to be operated on.

       DESCA   (global and local input) INTEGER array of dimension DLEN_.
               The array descriptor for the distributed matrix A.
               If DESCA( CTXT_ ) is incorrect, PSSYEVX cannot guarantee
               correct error reporting.

       VL      (global input) REAL
               If RANGE=’V’, the lower bound of the interval to be searched
               for eigenvalues.  Not referenced if RANGE = ’A’ or ’I’.

       VU      (global input) REAL
               If RANGE=’V’, the upper bound of the interval to be searched
               for eigenvalues.  Not referenced if RANGE = ’A’ or ’I’.

       IL      (global input) INTEGER
               If RANGE=’I’, the index (from smallest to largest) of the
               smallest eigenvalue to be returned.  IL >= 1.
               Not referenced if RANGE = ’A’ or ’V’.

       IU      (global input) INTEGER
               If RANGE=’I’, the index (from smallest to largest) of the
               largest eigenvalue to be returned.  min(IL,N) <= IU <= N.
               Not referenced if RANGE = ’A’ or ’V’.

       ABSTOL  (global input) REAL
               If JOBZ=’V’, setting ABSTOL to PSLAMCH( CONTEXT, ’U’) yields
               the most orthogonal eigenvectors.

               The absolute error tolerance for the eigenvalues.
               An approximate eigenvalue is accepted as converged
               when it is determined to lie in an interval [a,b]
               of width less than or equal to

                       ABSTOL + EPS *   max( |a|,|b| ) ,

               where EPS is the machine precision.  If ABSTOL is less than
               or equal to zero, then EPS*norm(T) will be used in its place,
               where norm(T) is the 1-norm of the tridiagonal matrix
               obtained by reducing A to tridiagonal form.

               Eigenvalues will be computed most accurately when ABSTOL is
               set to twice the underflow threshold 2*PSLAMCH(’S’) not zero.
               If this routine returns with ((MOD(INFO,2).NE.0) .OR.
               (MOD(INFO/8,2).NE.0)), indicating that some eigenvalues or
               eigenvectors did not converge, try setting ABSTOL to
               2*PSLAMCH(’S’).

               See "Computing Small Singular Values of Bidiagonal Matrices
               with Guaranteed High Relative Accuracy," by Demmel and
               Kahan, LAPACK Working Note #3.

               See "On the correctness of Parallel Bisection in Floating
               Point" by Demmel, Dhillon and Ren, LAPACK Working Note #70

       M       (global output) INTEGER
               Total number of eigenvalues found.  0 <= M <= N.

       NZ      (global output) INTEGER
               Total number of eigenvectors computed.  0 <= NZ <= M.
               The number of columns of Z that are filled.
               If JOBZ .NE. ’V’, NZ is not referenced.
               If JOBZ .EQ. ’V’, NZ = M unless the user supplies
               insufficient space and PSSYEVX is not able to detect this
               before beginning computation.  To get all the eigenvectors
               requested, the user must supply both sufficient
               space to hold the eigenvectors in Z (M .LE. DESCZ(N_))
               and sufficient workspace to compute them.  (See LWORK below.)
               PSSYEVX is always able to detect insufficient space without
               computation unless RANGE .EQ. ’V’.

       W       (global output) REAL array, dimension (N)
               On normal exit, the first M entries contain the selected
               eigenvalues in ascending order.

       ORFAC   (global input) REAL
               Specifies which eigenvectors should be reorthogonalized.
               Eigenvectors that correspond to eigenvalues which are within
               tol=ORFAC*norm(A) of each other are to be reorthogonalized.
               However, if the workspace is insufficient (see LWORK),
               tol may be decreased until all eigenvectors to be
               reorthogonalized can be stored in one process.
               No reorthogonalization will be done if ORFAC equals zero.
               A default value of 10^-3 is used if ORFAC is negative.
               ORFAC should be identical on all processes.

       Z       (local output) REAL array,
               global dimension (N, N),
               local dimension ( LLD_Z, LOCc(JZ+N-1) )
               If JOBZ = ’V’, then on normal exit the first M columns of Z
               contain the orthonormal eigenvectors of the matrix
               corresponding to the selected eigenvalues.  If an eigenvector
               fails to converge, then that column of Z contains the latest
               approximation to the eigenvector, and the index of the
               eigenvector is returned in IFAIL.
               If JOBZ = ’N’, then Z is not referenced.

       IZ      (global input) INTEGER
               Z’s global row index, which points to the beginning of the
               submatrix which is to be operated on.

       JZ      (global input) INTEGER
               Z’s global column index, which points to the beginning of
               the submatrix which is to be operated on.

       DESCZ   (global and local input) INTEGER array of dimension DLEN_.
               The array descriptor for the distributed matrix Z.
               DESCZ( CTXT_ ) must equal DESCA( CTXT_ )

       WORK    (local workspace/output) REAL array,
               dimension (LWORK)
               On return, WROK(1) contains the optimal amount of
               workspace required for efficient execution.
               if JOBZ=’N’ WORK(1) = optimal amount of workspace
                  required to compute eigenvalues efficiently
               if JOBZ=’V’ WORK(1) = optimal amount of workspace
                  required to compute eigenvalues and eigenvectors
                  efficiently with no guarantee on orthogonality.
                  If RANGE=’V’, it is assumed that all eigenvectors
                  may be required.

       LWORK   (local input) INTEGER
               Size of WORK
               See below for definitions of variables used to define LWORK.
               If no eigenvectors are requested (JOBZ = ’N’) then
                  LWORK >= 5 * N + MAX( 5 * NN, NB * ( NP0 + 1 ) )
               If eigenvectors are requested (JOBZ = ’V’ ) then
                  the amount of workspace required to guarantee that all
                  eigenvectors are computed is:
                  LWORK >= 5*N + MAX( 5*NN, NP0 * MQ0 + 2 * NB * NB ) +
                    ICEIL( NEIG, NPROW*NPCOL)*NN

                  The computed eigenvectors may not be orthogonal if the
                  minimal workspace is supplied and ORFAC is too small.
                  If you want to guarantee orthogonality (at the cost
                  of potentially poor performance) you should add
                  the following to LWORK:
                     (CLUSTERSIZE-1)*N
                  where CLUSTERSIZE is the number of eigenvalues in the
                  largest cluster, where a cluster is defined as a set of
                  close eigenvalues: { W(K),...,W(K+CLUSTERSIZE-1) |
                                       W(J+1) <= W(J) + ORFAC*2*norm(A) }
               Variable definitions:
                  NEIG = number of eigenvectors requested
                  NB = DESCA( MB_ ) = DESCA( NB_ ) =
                       DESCZ( MB_ ) = DESCZ( NB_ )
                  NN = MAX( N, NB, 2 )
                  DESCA( RSRC_ ) = DESCA( NB_ ) = DESCZ( RSRC_ ) =
                                   DESCZ( CSRC_ ) = 0
                  NP0 = NUMROC( NN, NB, 0, 0, NPROW )
                  MQ0 = NUMROC( MAX( NEIG, NB, 2 ), NB, 0, 0, NPCOL )
                  ICEIL( X, Y ) is a ScaLAPACK function returning
                  ceiling(X/Y)

               When LWORK is too small:
                  If LWORK is too small to guarantee orthogonality,
                  PSSYEVX attempts to maintain orthogonality in
                  the clusters with the smallest
                  spacing between the eigenvalues.
                  If LWORK is too small to compute all the eigenvectors
                  requested, no computation is performed and INFO=-23
                  is returned.  Note that when RANGE=’V’, PSSYEVX does
                  not know how many eigenvectors are requested until
                  the eigenvalues are computed.  Therefore, when RANGE=’V’
                  and as long as LWORK is large enough to allow PSSYEVX to
                  compute the eigenvalues, PSSYEVX will compute the
                  eigenvalues and as many eigenvectors as it can.

               Relationship between workspace, orthogonality & performance:
                  Greater performance can be achieved if adequate workspace
                  is provided.  On the other hand, in some situations,
                  performance can decrease as the workspace provided
                  increases above the workspace amount shown below:

                  For optimal performance, greater workspace may be
                  needed, i.e.
                     LWORK >=  MAX( LWORK, 5*N + NSYTRD_LWOPT )
                     Where:
                       LWORK, as defined previously, depends upon the number
                          of eigenvectors requested, and
                       NSYTRD_LWOPT = N + 2*( ANB+1 )*( 4*NPS+2 ) +
                         ( NPS + 3 ) *  NPS

                       ANB = PJLAENV( DESCA( CTXT_), 3, ’PSSYTTRD’, ’L’,
                          0, 0, 0, 0)
                       SQNPC = INT( SQRT( DBLE( NPROW * NPCOL ) ) )
                       NPS = MAX( NUMROC( N, 1, 0, 0, SQNPC ), 2*ANB )

                       NUMROC is a ScaLAPACK tool functions;
                       PJLAENV is a ScaLAPACK envionmental inquiry function
                       MYROW, MYCOL, NPROW and NPCOL can be determined by
                         calling the subroutine BLACS_GRIDINFO.

                     For large N, no extra workspace is needed, however the
                     biggest boost in performance comes for small N, so it
                     is wise to provide the extra workspace (typically less
                     than a Megabyte per process).

                  If CLUSTERSIZE >= N/SQRT(NPROW*NPCOL), then providing
                  enough space to compute all the eigenvectors
                  orthogonally will cause serious degradation in
                  performance. In the limit (i.e. CLUSTERSIZE = N-1)
                  PSSTEIN will perform no better than SSTEIN on 1
                  processor.
                  For CLUSTERSIZE = N/SQRT(NPROW*NPCOL) reorthogonalizing
                  all eigenvectors will increase the total execution time
                  by a factor of 2 or more.
                  For CLUSTERSIZE > N/SQRT(NPROW*NPCOL) execution time will
                  grow as the square of the cluster size, all other factors
                  remaining equal and assuming enough workspace.  Less
                  workspace means less reorthogonalization but faster
                  execution.

               If LWORK = -1, then LWORK is global input and a workspace
               query is assumed; the routine only calculates the size
               required for optimal performance for all work arrays. Each of
               these values is returned in the first entry of the
               corresponding work arrays, and no error message is issued by
               PXERBLA.

       IWORK   (local workspace) INTEGER array
               On return, IWORK(1) contains the amount of integer workspace
               required.

       LIWORK  (local input) INTEGER
               size of IWORK
               LIWORK >= 6 * NNP
               Where:
                 NNP = MAX( N, NPROW*NPCOL + 1, 4 )
               If LIWORK = -1, then LIWORK is global input and a workspace
               query is assumed; the routine only calculates the minimum
               and optimal size for all work arrays. Each of these
               values is returned in the first entry of the corresponding
               work array, and no error message is issued by PXERBLA.

       IFAIL   (global output) INTEGER array, dimension (N)
               If JOBZ = ’V’, then on normal exit, the first M elements of
               IFAIL are zero.  If (MOD(INFO,2).NE.0) on exit, then
               IFAIL contains the
               indices of the eigenvectors that failed to converge.
               If JOBZ = ’N’, then IFAIL is not referenced.

       ICLUSTR (global output) integer array, dimension (2*NPROW*NPCOL)
               This array contains indices of eigenvectors corresponding to
               a cluster of eigenvalues that could not be reorthogonalized
               due to insufficient workspace (see LWORK, ORFAC and INFO).
               Eigenvectors corresponding to clusters of eigenvalues indexed
               ICLUSTR(2*I-1) to ICLUSTR(2*I), could not be
               reorthogonalized due to lack of workspace. Hence the
               eigenvectors corresponding to these clusters may not be
               orthogonal.  ICLUSTR() is a zero terminated array.
               (ICLUSTR(2*K).NE.0 .AND. ICLUSTR(2*K+1).EQ.0) if and only if
               K is the number of clusters
               ICLUSTR is not referenced if JOBZ = ’N’

       GAP     (global output) REAL array,
                  dimension (NPROW*NPCOL)
               This array contains the gap between eigenvalues whose
               eigenvectors could not be reorthogonalized. The output
               values in this array correspond to the clusters indicated
               by the array ICLUSTR. As a result, the dot product between
               eigenvectors correspoding to the I^th cluster may be as high
               as ( C * n ) / GAP(I) where C is a small constant.

       INFO    (global output) INTEGER
               = 0:  successful exit
               < 0:  If the i-th argument is an array and the j-entry had
                     an illegal value, then INFO = -(i*100+j), if the i-th
                     argument is a scalar and had an illegal value, then
                     INFO = -i.
               > 0:  if (MOD(INFO,2).NE.0), then one or more eigenvectors
                       failed to converge.  Their indices are stored
                       in IFAIL.  Ensure ABSTOL=2.0*PSLAMCH( ’U’ )
                       Send e-mail to scalapack@cs.utk.edu
                     if (MOD(INFO/2,2).NE.0),then eigenvectors corresponding
                       to one or more clusters of eigenvalues could not be
                       reorthogonalized because of insufficient workspace.
                       The indices of the clusters are stored in the array
                       ICLUSTR.
                     if (MOD(INFO/4,2).NE.0), then space limit prevented
                       PSSYEVX from computing all of the eigenvectors
                       between VL and VU.  The number of eigenvectors
                       computed is returned in NZ.
                     if (MOD(INFO/8,2).NE.0), then PSSTEBZ failed to compute
                       eigenvalues.  Ensure ABSTOL=2.0*PSLAMCH( ’U’ )
                       Send e-mail to scalapack@cs.utk.edu