pdsyevd.f

Go to the documentation of this file.

      SUBROUTINE pdsyevd( JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ,
     $                    DESCZ, WORK, LWORK, IWORK, LIWORK, INFO )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     March 14, 2000
*
*     .. Scalar Arguments ..
      CHARACTER          JOBZ, UPLO
      INTEGER            IA, INFO, IZ, JA, JZ, LIWORK, LWORK, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCZ( * ), IWORK( * )
      DOUBLE PRECISION   A( * ), W( * ), WORK( * ), Z( * )
*     ..
*
*  Purpose
*  =======
*
*  PDSYEVD computes  all the eigenvalues and eigenvectors
*  of a real symmetric matrix A by calling the recommended sequence
*  of ScaLAPACK routines.
*
*  In its present form, PDSYEVD assumes a homogeneous system and makes
*  no checks for consistency of the eigenvalues or eigenvectors across
*  the different processes.  Because of this, it is possible that a
*  heterogeneous system may return incorrect results without any error
*  messages.
*
*  Arguments
*  =========
*
*     NP = the number of rows local to a given process.
*     NQ = the number of columns local to a given process.
*
*  JOBZ    (input) CHARACTER*1
*          = 'N':  Compute eigenvalues only;     (NOT IMPLEMENTED YET)
*          = 'V':  Compute eigenvalues and eigenvectors.
*
*  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 to be operated on, i.e. the
*          order of the distributed submatrix sub( A ). N >= 0.
*
*  A       (local input/workspace) block cyclic DOUBLE PRECISION 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.
*
*  W       (global output) DOUBLE PRECISION array, dimension (N)
*          If INFO=0, the eigenvalues in ascending order.
*
*  Z       (local output) DOUBLE PRECISION array,
*          global dimension (N, N),
*          local dimension ( LLD_Z, LOCc(JZ+N-1) )
*          Z contains the orthonormal eigenvectors
*          of the symmetric matrix A.
*
*  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) DOUBLE PRECISION array,
*          dimension (LWORK)
*          On output, WORK(1) returns the workspace required.
*
*  LWORK   (local input) INTEGER
*          LWORK >= MAX( 1+6*N+2*NP*NQ, TRILWMIN ) + 2*N
*          TRILWMIN = 3*N + MAX( NB*( NP+1 ), 3*NB )
*          NP = NUMROC( N, NB, MYROW, IAROW, NPROW )
*          NQ = NUMROC( N, NB, MYCOL, IACOL, NPCOL )
*
*          If LWORK = -1, the LWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          size for the WORK array.  The required workspace is returned
*          as the first element of WORK and no error message is issued
*          by PXERBLA.
*
*  IWORK   (local workspace/output) INTEGER array, dimension (LIWORK)
*          On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK.
*
*  LIWORK  (input) INTEGER
*          The dimension of the array IWORK.
*          LIWORK = 7*N + 8*NPCOL + 2
*
*  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:  The algorithm failed to compute the INFO/(N+1) th
*                eigenvalue while working on the submatrix lying in
*                global rows and columns mod(INFO,N+1).
*
*  Alignment requirements
*  ======================
*
*  The distributed submatrices sub( A ), sub( Z ) must verify
*  some alignment properties, namely the following expression
*  should be true:
*  ( MB_A.EQ.NB_A.EQ.MB_Z.EQ.NB_Z .AND. IROFFA.EQ.ICOFFA .AND.
*    IROFFA.EQ.0 .AND.IROFFA.EQ.IROFFZ. AND. IAROW.EQ.IZROW)
*    with IROFFA = MOD( IA-1, MB_A )
*     and ICOFFA = MOD( JA-1, NB_A ).
*
*  Further Details
*  ======= =======
*
*  Contributed by Francoise Tisseur, University of Manchester.
*
*  Reference:  F. Tisseur and J. Dongarra, "A Parallel Divide and
*              Conquer Algorithm for the Symmetric Eigenvalue Problem
*              on Distributed Memory Architectures",
*              SIAM J. Sci. Comput., 6:20 (1999), pp. 2223--2236.
*              (see also LAPACK Working Note 132)
*                http://www.netlib.org/lapack/lawns/lawn132.ps
*
*  =====================================================================
*
*     .. Parameters ..
*
      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 )
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY, UPPER
      INTEGER            IACOL, IAROW, ICOFFA, ICOFFZ, ICTXT, IINFO,
     $                   indd, inde, inde2, indtau, indwork, indwork2,
     $                   iroffa, iroffz, iscale, liwmin, llwork,
     $                   llwork2, lwmin, mycol, myrow, nb, np, npcol,
     $                   nprow, nq, offset, trilwmin
      DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
     $                   smlnum
*     ..
*     .. Local Arrays ..
*     ..
      INTEGER            IDUM1( 2 ), IDUM2( 2 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            INDXG2P, NUMROC
      DOUBLE PRECISION   PDLAMCH, PDLANSY
      EXTERNAL           lsame, indxg2p, numroc, pdlamch, pdlansy
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, chk1mat, dscal, pchk1mat,
     $                   pdlared1d, pdlascl, pdlaset, pdormtr, pdstedc,
     $                   pdsytrd, pxerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dble, ichar, max, min, mod, sqrt
*     ..
*     .. Executable Statements ..
*       This is just to keep ftnchek and toolpack/1 happy
      IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
     $    rsrc_.LT.0 )RETURN
*
*     Quick return
*
      IF( n.EQ.0 )
     $   RETURN
*
*     Test the input arguments.
*
      ictxt = descz( ctxt_ )
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      info = 0
      IF( nprow.EQ.-1 ) THEN
         info = -( 600+ctxt_ )
      ELSE
         CALL chk1mat( n, 3, n, 3, ia, ja, desca, 7, info )
         CALL chk1mat( n, 3, n, 3, iz, jz, descz, 12, info )
         IF( info.EQ.0 ) THEN
            upper = lsame( uplo, 'U' )
            nb = desca( nb_ )
            iroffa = mod( ia-1, desca( mb_ ) )
            icoffa = mod( ja-1, desca( nb_ ) )
            iroffz = mod( iz-1, descz( mb_ ) )
            icoffz = mod( jz-1, descz( nb_ ) )
            iarow = indxg2p( ia, nb, myrow, desca( rsrc_ ), nprow )
            iacol = indxg2p( ja, nb, mycol, desca( csrc_ ), npcol )
            np = numroc( n, nb, myrow, iarow, nprow )
            nq = numroc( n, nb, mycol, iacol, npcol )
*
            lquery = ( lwork.EQ.-1 )
            trilwmin = 3*n + max( nb*( np+1 ), 3*nb )
            lwmin = max( 1+6*n+2*np*nq, trilwmin ) + 2*n
            liwmin = 7*n + 8*npcol + 2
            work( 1 ) = dble( lwmin )
            iwork( 1 ) = liwmin
            IF( .NOT.lsame( jobz, 'V' ) ) THEN
               info = -1
            ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
               info = -2
            ELSE IF( iroffa.NE.icoffa .OR. icoffa.NE.0 ) THEN
               info = -6
            ELSE IF( iroffa.NE.iroffz .OR. icoffa.NE.icoffz ) THEN
               info = -10
            ELSE IF( desca( m_ ).NE.descz( m_ ) ) THEN
               info = -( 1200+m_ )
            ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
               info = -( 700+nb_ )
            ELSE IF( descz( mb_ ).NE.descz( nb_ ) ) THEN
               info = -( 1200+nb_ )
            ELSE IF( desca( mb_ ).NE.descz( mb_ ) ) THEN
               info = -( 1200+mb_ )
            ELSE IF( desca( ctxt_ ).NE.descz( ctxt_ ) ) THEN
               info = -( 1200+ctxt_ )
            ELSE IF( desca( rsrc_ ).NE.descz( rsrc_ ) ) THEN
               info = -( 1200+rsrc_ )
            ELSE IF( desca( csrc_ ).NE.descz( csrc_ ) ) THEN
               info = -( 1200+csrc_ )
            ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
               info = -14
            ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
               info = -16
            END IF
         END IF
         IF( upper ) THEN
            idum1( 1 ) = ichar( 'U' )
         ELSE
            idum1( 1 ) = ichar( 'L' )
         END IF
         idum2( 1 ) = 2
         IF( lwork.EQ.-1 ) THEN
            idum1( 2 ) = -1
         ELSE
            idum1( 2 ) = 1
         END IF
         idum2( 2 ) = 14
         CALL pchk1mat( n, 3, n, 3, ia, ja, desca, 7, 2, idum1, idum2,
     $                  info )
      END IF
      IF( info.NE.0 ) THEN
         CALL pxerbla( ictxt, 'PDSYEVD', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
*     Set up pointers into the WORK array
*
      indtau = 1
      inde = indtau + n
      indd = inde + n
      inde2 = indd + n
      indwork = inde2 + n
      llwork = lwork - indwork + 1
      indwork2 = indd
      llwork2 = lwork - indwork2 + 1
*
*     Scale matrix to allowable range, if necessary.
*
      iscale = 0
      safmin = pdlamch( desca( ctxt_ ), 'Safe minimum' )
      eps = pdlamch( desca( ctxt_ ), 'Precision' )
      smlnum = safmin / eps
      bignum = one / smlnum
      rmin = sqrt( smlnum )
      rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
      anrm = pdlansy( 'M', uplo, n, a, ia, ja, desca, work( indwork ) )
*
      IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
         iscale = 1
         sigma = rmin / anrm
      ELSE IF( anrm.GT.rmax ) THEN
         iscale = 1
         sigma = rmax / anrm
      END IF
*
      IF( iscale.EQ.1 ) THEN
         CALL pdlascl( uplo, one, sigma, n, n, a, ia, ja, desca, iinfo )
      END IF
*
*     Reduce symmetric matrix to tridiagonal form.
*
*
      CALL pdsytrd( uplo, n, a, ia, ja, desca, work( indd ),
     $              work( inde2 ), work( indtau ), work( indwork ),
     $              llwork, iinfo )
*
*     Copy the values of D, E to all processes.
*
      CALL pdlared1d( n, ia, ja, desca, work( indd ), w,
     $                work( indwork ), llwork )
*
      CALL pdlared1d( n, ia, ja, desca, work( inde2 ), work( inde ),
     $                work( indwork ), llwork )
*
      CALL pdlaset( 'Full', n, n, zero, one, z, 1, 1, descz )
*
      IF( upper ) THEN
         offset = 1
      ELSE
         offset = 0
      END IF
      CALL pdstedc( 'I', n, w, work( inde+offset ), z, iz, jz, descz,
     $              work( indwork2 ), llwork2, iwork, liwork, info )
*
      CALL pdormtr( 'L', uplo, 'N', n, n, a, ia, ja, desca,
     $              work( indtau ), z, iz, jz, descz, work( indwork2 ),
     $              llwork2, iinfo )
*
*     If matrix was scaled, then rescale eigenvalues appropriately.
*
      IF( iscale.EQ.1 ) THEN
         CALL dscal( n, one / sigma, w, 1 )
      END IF
*
      RETURN
*
*     End of PDSYEVD
*
      END