ScaLAPACK  2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
pdsyevd.f
Go to the documentation of this file.
00001       SUBROUTINE PDSYEVD( JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ,
00002      $                    DESCZ, WORK, LWORK, IWORK, LIWORK, INFO )
00003 *
00004 *  -- ScaLAPACK routine (version 1.7) --
00005 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
00006 *     and University of California, Berkeley.
00007 *     March 14, 2000
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          JOBZ, UPLO
00011       INTEGER            IA, INFO, IZ, JA, JZ, LIWORK, LWORK, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            DESCA( * ), DESCZ( * ), IWORK( * )
00015       DOUBLE PRECISION   A( * ), W( * ), WORK( * ), Z( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  PDSYEVD computes  all the eigenvalues and eigenvectors
00022 *  of a real symmetric matrix A by calling the recommended sequence
00023 *  of ScaLAPACK routines.
00024 *
00025 *  In its present form, PDSYEVD assumes a homogeneous system and makes
00026 *  no checks for consistency of the eigenvalues or eigenvectors across
00027 *  the different processes.  Because of this, it is possible that a
00028 *  heterogeneous system may return incorrect results without any error
00029 *  messages.
00030 *
00031 *  Arguments
00032 *  =========
00033 *
00034 *     NP = the number of rows local to a given process.
00035 *     NQ = the number of columns local to a given process.
00036 *
00037 *  JOBZ    (input) CHARACTER*1
00038 *          = 'N':  Compute eigenvalues only;     (NOT IMPLEMENTED YET)
00039 *          = 'V':  Compute eigenvalues and eigenvectors.
00040 *
00041 *  UPLO    (global input) CHARACTER*1
00042 *          Specifies whether the upper or lower triangular part of the
00043 *          symmetric matrix A is stored:
00044 *          = 'U':  Upper triangular
00045 *          = 'L':  Lower triangular
00046 *
00047 *  N       (global input) INTEGER
00048 *          The number of rows and columns to be operated on, i.e. the
00049 *          order of the distributed submatrix sub( A ). N >= 0.
00050 *
00051 *  A       (local input/workspace) block cyclic DOUBLE PRECISION array,
00052 *          global dimension (N, N), local dimension ( LLD_A,
00053 *          LOCc(JA+N-1) )
00054 *          On entry, the symmetric matrix A.  If UPLO = 'U', only the
00055 *          upper triangular part of A is used to define the elements of
00056 *          the symmetric matrix.  If UPLO = 'L', only the lower
00057 *          triangular part of A is used to define the elements of the
00058 *          symmetric matrix.
00059 *          On exit, the lower triangle (if UPLO='L') or the upper
00060 *          triangle (if UPLO='U') of A, including the diagonal, is
00061 *          destroyed.
00062 *
00063 *  IA      (global input) INTEGER
00064 *          A's global row index, which points to the beginning of the
00065 *          submatrix which is to be operated on.
00066 *
00067 *  JA      (global input) INTEGER
00068 *          A's global column index, which points to the beginning of
00069 *          the submatrix which is to be operated on.
00070 *
00071 *  DESCA   (global and local input) INTEGER array of dimension DLEN_.
00072 *          The array descriptor for the distributed matrix A.
00073 *
00074 *  W       (global output) DOUBLE PRECISION array, dimension (N)
00075 *          If INFO=0, the eigenvalues in ascending order.
00076 *
00077 *  Z       (local output) DOUBLE PRECISION array,
00078 *          global dimension (N, N),
00079 *          local dimension ( LLD_Z, LOCc(JZ+N-1) )
00080 *          Z contains the orthonormal eigenvectors
00081 *          of the symmetric matrix A.
00082 *
00083 *  IZ      (global input) INTEGER
00084 *          Z's global row index, which points to the beginning of the
00085 *          submatrix which is to be operated on.
00086 *
00087 *  JZ      (global input) INTEGER
00088 *          Z's global column index, which points to the beginning of
00089 *          the submatrix which is to be operated on.
00090 *
00091 *  DESCZ   (global and local input) INTEGER array of dimension DLEN_.
00092 *          The array descriptor for the distributed matrix Z.
00093 *          DESCZ( CTXT_ ) must equal DESCA( CTXT_ )
00094 *
00095 *  WORK    (local workspace/output) DOUBLE PRECISION array,
00096 *          dimension (LWORK)
00097 *          On output, WORK(1) returns the workspace required.
00098 *
00099 *  LWORK   (local input) INTEGER
00100 *          LWORK >= MAX( 1+6*N+2*NP*NQ, TRILWMIN ) + 2*N
00101 *          TRILWMIN = 3*N + MAX( NB*( NP+1 ), 3*NB )
00102 *          NP = NUMROC( N, NB, MYROW, IAROW, NPROW )
00103 *          NQ = NUMROC( N, NB, MYCOL, IACOL, NPCOL )
00104 *
00105 *          If LWORK = -1, the LWORK is global input and a workspace
00106 *          query is assumed; the routine only calculates the minimum
00107 *          size for the WORK array.  The required workspace is returned
00108 *          as the first element of WORK and no error message is issued
00109 *          by PXERBLA.
00110 *
00111 *  IWORK   (local workspace/output) INTEGER array, dimension (LIWORK)
00112 *          On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK.
00113 *
00114 *  LIWORK  (input) INTEGER
00115 *          The dimension of the array IWORK.
00116 *          LIWORK = 7*N + 8*NPCOL + 2
00117 *
00118 *  INFO    (global output) INTEGER
00119 *          = 0:  successful exit
00120 *          < 0:  If the i-th argument is an array and the j-entry had
00121 *                an illegal value, then INFO = -(i*100+j), if the i-th
00122 *                argument is a scalar and had an illegal value, then
00123 *                INFO = -i.
00124 *          > 0:  The algorithm failed to compute the INFO/(N+1) th
00125 *                eigenvalue while working on the submatrix lying in
00126 *                global rows and columns mod(INFO,N+1).
00127 *
00128 *  Alignment requirements
00129 *  ======================
00130 *
00131 *  The distributed submatrices sub( A ), sub( Z ) must verify
00132 *  some alignment properties, namely the following expression
00133 *  should be true:
00134 *  ( MB_A.EQ.NB_A.EQ.MB_Z.EQ.NB_Z .AND. IROFFA.EQ.ICOFFA .AND.
00135 *    IROFFA.EQ.0 .AND.IROFFA.EQ.IROFFZ. AND. IAROW.EQ.IZROW)
00136 *    with IROFFA = MOD( IA-1, MB_A )
00137 *     and ICOFFA = MOD( JA-1, NB_A ).
00138 *
00139 *  Further Details
00140 *  ======= =======
00141 *
00142 *  Contributed by Francoise Tisseur, University of Manchester.
00143 *
00144 *  Reference:  F. Tisseur and J. Dongarra, "A Parallel Divide and
00145 *              Conquer Algorithm for the Symmetric Eigenvalue Problem
00146 *              on Distributed Memory Architectures",
00147 *              SIAM J. Sci. Comput., 6:20 (1999), pp. 2223--2236.
00148 *              (see also LAPACK Working Note 132)
00149 *                http://www.netlib.org/lapack/lawns/lawn132.ps
00150 *
00151 *  =====================================================================
00152 *
00153 *     .. Parameters ..
00154 *
00155       INTEGER            BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
00156      $                   MB_, NB_, RSRC_, CSRC_, LLD_
00157       PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
00158      $                   CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
00159      $                   RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
00160       DOUBLE PRECISION   ZERO, ONE
00161       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00162 *     ..
00163 *     .. Local Scalars ..
00164       LOGICAL            LQUERY, UPPER
00165       INTEGER            IACOL, IAROW, ICOFFA, ICOFFZ, ICTXT, IINFO,
00166      $                   INDD, INDE, INDE2, INDTAU, INDWORK, INDWORK2,
00167      $                   IROFFA, IROFFZ, ISCALE, LIWMIN, LLWORK,
00168      $                   LLWORK2, LWMIN, MYCOL, MYROW, NB, NP, NPCOL,
00169      $                   NPROW, NQ, OFFSET, TRILWMIN
00170       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
00171      $                   SMLNUM
00172 *     ..
00173 *     .. Local Arrays ..
00174 *     ..
00175       INTEGER            IDUM1( 2 ), IDUM2( 2 )
00176 *     ..
00177 *     .. External Functions ..
00178       LOGICAL            LSAME
00179       INTEGER            INDXG2P, NUMROC
00180       DOUBLE PRECISION   PDLAMCH, PDLANSY
00181       EXTERNAL           LSAME, INDXG2P, NUMROC, PDLAMCH, PDLANSY
00182 *     ..
00183 *     .. External Subroutines ..
00184       EXTERNAL           BLACS_GRIDINFO, CHK1MAT, DSCAL, PCHK1MAT,
00185      $                   PDLARED1D, PDLASCL, PDLASET, PDORMTR, PDSTEDC,
00186      $                   PDSYTRD, PXERBLA
00187 *     ..
00188 *     .. Intrinsic Functions ..
00189       INTRINSIC          DBLE, ICHAR, MAX, MIN, MOD, SQRT
00190 *     ..
00191 *     .. Executable Statements ..
00192 *       This is just to keep ftnchek and toolpack/1 happy
00193       IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_*
00194      $    RSRC_.LT.0 )RETURN
00195 *
00196 *     Quick return
00197 *
00198       IF( N.EQ.0 )
00199      $   RETURN
00200 *
00201 *     Test the input arguments.
00202 *
00203       CALL BLACS_GRIDINFO( DESCZ( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL )
00204 *
00205       INFO = 0
00206       IF( NPROW.EQ.-1 ) THEN
00207          INFO = -( 600+CTXT_ )
00208       ELSE
00209          CALL CHK1MAT( N, 3, N, 3, IA, JA, DESCA, 7, INFO )
00210          CALL CHK1MAT( N, 3, N, 3, IZ, JZ, DESCZ, 12, INFO )
00211          IF( INFO.EQ.0 ) THEN
00212             UPPER = LSAME( UPLO, 'U' )
00213             NB = DESCA( NB_ )
00214             IROFFA = MOD( IA-1, DESCA( MB_ ) )
00215             ICOFFA = MOD( JA-1, DESCA( NB_ ) )
00216             IROFFZ = MOD( IZ-1, DESCZ( MB_ ) )
00217             ICOFFZ = MOD( JZ-1, DESCZ( NB_ ) )
00218             IAROW = INDXG2P( IA, NB, MYROW, DESCA( RSRC_ ), NPROW )
00219             IACOL = INDXG2P( JA, NB, MYCOL, DESCA( CSRC_ ), NPCOL )
00220             NP = NUMROC( N, NB, MYROW, IAROW, NPROW )
00221             NQ = NUMROC( N, NB, MYCOL, IACOL, NPCOL )
00222 *
00223             LQUERY = ( LWORK.EQ.-1 )
00224             TRILWMIN = 3*N + MAX( NB*( NP+1 ), 3*NB )
00225             LWMIN = MAX( 1+6*N+2*NP*NQ, TRILWMIN ) + 2*N
00226             LIWMIN = 7*N + 8*NPCOL + 2
00227             WORK( 1 ) = DBLE( LWMIN )
00228             IWORK( 1 ) = LIWMIN
00229             IF( .NOT.LSAME( JOBZ, 'V' ) ) THEN
00230                INFO = -1
00231             ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00232                INFO = -2
00233             ELSE IF( IROFFA.NE.ICOFFA .OR. ICOFFA.NE.0 ) THEN
00234                INFO = -6
00235             ELSE IF( IROFFA.NE.IROFFZ .OR. ICOFFA.NE.ICOFFZ ) THEN
00236                INFO = -10
00237             ELSE IF( DESCA( M_ ).NE.DESCZ( M_ ) ) THEN
00238                INFO = -( 1200+M_ )
00239             ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
00240                INFO = -( 700+NB_ )
00241             ELSE IF( DESCZ( MB_ ).NE.DESCZ( NB_ ) ) THEN
00242                INFO = -( 1200+NB_ )
00243             ELSE IF( DESCA( MB_ ).NE.DESCZ( MB_ ) ) THEN
00244                INFO = -( 1200+MB_ )
00245             ELSE IF( DESCA( CTXT_ ).NE.DESCZ( CTXT_ ) ) THEN
00246                INFO = -( 1200+CTXT_ )
00247             ELSE IF( DESCA( RSRC_ ).NE.DESCZ( RSRC_ ) ) THEN
00248                INFO = -( 1200+RSRC_ )
00249             ELSE IF( DESCA( CSRC_ ).NE.DESCZ( CSRC_ ) ) THEN
00250                INFO = -( 1200+CSRC_ )
00251             ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00252                INFO = -14
00253             ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
00254                INFO = -16
00255             END IF
00256          END IF
00257          IF( UPPER ) THEN
00258             IDUM1( 1 ) = ICHAR( 'U' )
00259          ELSE
00260             IDUM1( 1 ) = ICHAR( 'L' )
00261          END IF
00262          IDUM2( 1 ) = 2
00263          IF( LWORK.EQ.-1 ) THEN
00264             IDUM1( 2 ) = -1
00265          ELSE
00266             IDUM1( 2 ) = 1
00267          END IF
00268          IDUM2( 2 ) = 14
00269          CALL PCHK1MAT( N, 3, N, 3, IA, JA, DESCA, 7, 2, IDUM1, IDUM2,
00270      $                  INFO )
00271       END IF
00272       IF( INFO.NE.0 ) THEN
00273          CALL PXERBLA( ICTXT, 'PDSYEVD', -INFO )
00274          RETURN
00275       ELSE IF( LQUERY ) THEN
00276          RETURN
00277       END IF
00278 *
00279 *     Set up pointers into the WORK array
00280 *
00281       INDTAU = 1
00282       INDE = INDTAU + N
00283       INDD = INDE + N
00284       INDE2 = INDD + N
00285       INDWORK = INDE2 + N
00286       LLWORK = LWORK - INDWORK + 1
00287       INDWORK2 = INDD
00288       LLWORK2 = LWORK - INDWORK2 + 1
00289 *
00290 *     Scale matrix to allowable range, if necessary.
00291 *
00292       ISCALE = 0
00293       SAFMIN = PDLAMCH( DESCA( CTXT_ ), 'Safe minimum' )
00294       EPS = PDLAMCH( DESCA( CTXT_ ), 'Precision' )
00295       SMLNUM = SAFMIN / EPS
00296       BIGNUM = ONE / SMLNUM
00297       RMIN = SQRT( SMLNUM )
00298       RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
00299       ANRM = PDLANSY( 'M', UPLO, N, A, IA, JA, DESCA, WORK( INDWORK ) )
00300 *
00301       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
00302          ISCALE = 1
00303          SIGMA = RMIN / ANRM
00304       ELSE IF( ANRM.GT.RMAX ) THEN
00305          ISCALE = 1
00306          SIGMA = RMAX / ANRM
00307       END IF
00308 *
00309       IF( ISCALE.EQ.1 ) THEN
00310          CALL PDLASCL( UPLO, ONE, SIGMA, N, N, A, IA, JA, DESCA, IINFO )
00311       END IF
00312 *
00313 *     Reduce symmetric matrix to tridiagonal form.
00314 *
00315 *
00316       CALL PDSYTRD( UPLO, N, A, IA, JA, DESCA, WORK( INDD ),
00317      $              WORK( INDE2 ), WORK( INDTAU ), WORK( INDWORK ),
00318      $              LLWORK, IINFO )
00319 *
00320 *     Copy the values of D, E to all processes.
00321 *
00322       CALL PDLARED1D( N, IA, JA, DESCA, WORK( INDD ), W,
00323      $                WORK( INDWORK ), LLWORK )
00324 *
00325       CALL PDLARED1D( N, IA, JA, DESCA, WORK( INDE2 ), WORK( INDE ),
00326      $                WORK( INDWORK ), LLWORK )
00327 *
00328       CALL PDLASET( 'Full', N, N, ZERO, ONE, Z, 1, 1, DESCZ )
00329 *
00330       IF( UPPER ) THEN
00331          OFFSET = 1
00332       ELSE
00333          OFFSET = 0
00334       END IF
00335       CALL PDSTEDC( 'I', N, W, WORK( INDE+OFFSET ), Z, IZ, JZ, DESCZ,
00336      $              WORK( INDWORK2 ), LLWORK2, IWORK, LIWORK, INFO )
00337 *
00338       CALL PDORMTR( 'L', UPLO, 'N', N, N, A, IA, JA, DESCA,
00339      $              WORK( INDTAU ), Z, IZ, JZ, DESCZ, WORK( INDWORK2 ),
00340      $              LLWORK2, IINFO )
00341 *
00342 *     If matrix was scaled, then rescale eigenvalues appropriately.
00343 *
00344       IF( ISCALE.EQ.1 ) THEN
00345          CALL DSCAL( N, ONE / SIGMA, W, 1 )
00346       END IF
00347 *
00348       RETURN
00349 *
00350 *     End of PDSYEVD
00351 *
00352       END