psstedc.f

Go to the documentation of this file.

      SUBROUTINE psstedc( COMPZ, N, D, E, Q, IQ, JQ, DESCQ, WORK, LWORK,
     $                    IWORK, LIWORK, INFO )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     March 13, 2000
*
*     .. Scalar Arguments ..
      CHARACTER          COMPZ
      INTEGER            INFO, IQ, JQ, LIWORK, LWORK, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCQ( * ), IWORK( * )
      REAL               D( * ), E( * ), Q( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*  PSSTEDC computes all eigenvalues and eigenvectors of a
*  symmetric tridiagonal matrix in parallel, using the divide and
*  conquer algorithm.
*
*  This code makes very mild assumptions about floating point
*  arithmetic. It will work on machines with a guard digit in
*  add/subtract, or on those binary machines without guard digits
*  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
*  It could conceivably fail on hexadecimal or decimal machines
*  without guard digits, but we know of none.  See SLAED3 for details.
*
*  Arguments
*  =========
*
*  COMPZ   (input) CHARACTER*1
*          = 'N':  Compute eigenvalues only.    (NOT IMPLEMENTED YET)
*          = 'I':  Compute eigenvectors of tridiagonal matrix also.
*          = 'V':  Compute eigenvectors of original dense symmetric
*                  matrix also.  On entry, Z contains the orthogonal
*                  matrix used to reduce the original matrix to
*                  tridiagonal form.            (NOT IMPLEMENTED YET)
*
*  N       (global input) INTEGER
*          The order of the tridiagonal matrix T.  N >= 0.
*
*  D       (global input/output) REAL array, dimension (N)
*          On entry, the diagonal elements of the tridiagonal matrix.
*          On exit, if INFO = 0, the eigenvalues in descending order.
*
*  E       (global input/output) REAL array, dimension (N-1)
*          On entry, the subdiagonal elements of the tridiagonal matrix.
*          On exit, E has been destroyed.
*
*  Q       (local output) REAL array,
*          local dimension ( LLD_Q, LOCc(JQ+N-1))
*          Q  contains the orthonormal eigenvectors of the symmetric
*          tridiagonal matrix.
*          On output, Q is distributed across the P processes in block
*          cyclic format.
*
*  IQ      (global input) INTEGER
*          Q's global row index, which points to the beginning of the
*          submatrix which is to be operated on.
*
*  JQ      (global input) INTEGER
*          Q's global column index, which points to the beginning of
*          the submatrix which is to be operated on.
*
*  DESCQ   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix Z.
*
*
*  WORK    (local workspace/output) REAL array,
*          dimension (LWORK)
*          On output, WORK(1) returns the workspace needed.
*
*  LWORK   (local input/output) INTEGER,
*          the dimension of the array WORK.
*          LWORK = 6*N + 2*NP*NQ
*          NP = NUMROC( N, NB, MYROW, DESCQ( RSRC_ ), NPROW )
*          NQ = NUMROC( N, NB, MYCOL, DESCQ( CSRC_ ), 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 = 2 + 7*N + 8*NPCOL
*
*  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).
*
*  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 )
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY
      INTEGER            ICOFFQ, IIQ, IPQ, IQCOL, IQROW, IROFFQ, JJQ,
     $                   ldq, liwmin, lwmin, mycol, myrow, nb, np,
     $                   npcol, nprow, nq
      REAL               ORGNRM
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            INDXG2P, NUMROC
      REAL               SLANST
      EXTERNAL           indxg2p, lsame, numroc, slanst
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, chk1mat, infog2l, pslaed0,
     $                   pslasrt, pxerbla, slascl, sstedc
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          mod, real
*     ..
*     .. 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
*
*     Test the input parameters.
*
      CALL blacs_gridinfo( descq( ctxt_ ), nprow, npcol, myrow, mycol )
      ldq = descq( lld_ )
      nb = descq( nb_ )
      np = numroc( n, nb, myrow, descq( rsrc_ ), nprow )
      nq = numroc( n, nb, mycol, descq( csrc_ ), npcol )
      info = 0
      IF( nprow.EQ.-1 ) THEN
         info = -( 600+ctxt_ )
      ELSE
         CALL chk1mat( n, 2, n, 2, iq, jq, descq, 8, info )
         IF( info.EQ.0 ) THEN
            nb = descq( nb_ )
            iroffq = mod( iq-1, descq( mb_ ) )
            icoffq = mod( jq-1, descq( nb_ ) )
            iqrow = indxg2p( iq, nb, myrow, descq( rsrc_ ), nprow )
            iqcol = indxg2p( jq, nb, mycol, descq( csrc_ ), npcol )
            lwmin = 6*n + 2*np*nq
            liwmin = 2 + 7*n + 8*npcol
            work( 1 ) = real( lwmin )
            iwork( 1 ) = liwmin
            lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
            IF( .NOT.lsame( compz, 'I' ) ) THEN
               info = -1
            ELSE IF( n.LT.0 ) THEN
               info = -2
            ELSE IF( iroffq.NE.icoffq .OR. icoffq.NE.0 ) THEN
               info = -5
            ELSE IF( descq( mb_ ).NE.descq( nb_ ) ) THEN
               info = -( 700+nb_ )
            ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
               info = -10
            ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
               info = -12
            END IF
         END IF
      END IF
      IF( info.NE.0 ) THEN
         CALL pxerbla( descq( ctxt_ ), 'PSSTEDC', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
*     Quick return
*
      IF( n.EQ.0 )
     $   GO TO 10
      CALL infog2l( iq, jq, descq, nprow, npcol, myrow, mycol, iiq, jjq,
     $              iqrow, iqcol )
      IF( n.EQ.1 ) THEN
         IF( myrow.EQ.iqrow .AND. mycol.EQ.iqcol )
     $      q( 1 ) = one
         GO TO 10
      END IF
*
*     If N is smaller than the minimum divide size NB, then
*     solve the problem with the serial divide and conquer
*     code locally.
*
      IF( n.LE.nb ) THEN
         IF( ( myrow.EQ.iqrow ) .AND. ( mycol.EQ.iqcol ) ) THEN
            ipq = iiq + ( jjq-1 )*ldq
            CALL sstedc( 'I', n, d, e, q( ipq ), ldq, work, lwork,
     $                   iwork, liwork, info )
            IF( info.NE.0 ) THEN
               info = ( n+1 ) + n
               GO TO 10
            END IF
         END IF
         GO TO 10
      END IF
*
*     If P=NPROW*NPCOL=1, solve the problem with SSTEDC.
*
      IF( npcol*nprow.EQ.1 ) THEN
         ipq = iiq + ( jjq-1 )*ldq
         CALL sstedc( 'I', n, d, e, q( ipq ), ldq, work, lwork, iwork,
     $                liwork, info )
         GO TO 10
      END IF
*
*     Scale matrix to allowable range, if necessary.
*
      orgnrm = slanst( 'M', n, d, e )
      IF( orgnrm.NE.zero ) THEN
         CALL slascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
         CALL slascl( 'G', 0, 0, orgnrm, one, n-1, 1, e, n-1, info )
      END IF
*
      CALL pslaed0( n, d, e, q, iq, jq, descq, work, iwork, info )
*
*     Sort eigenvalues and corresponding eigenvectors
*
      CALL pslasrt( 'I', n, d, q, iq, jq, descq, work, lwork, iwork,
     $              liwork, info )
*
*           Scale back.
*
      IF( orgnrm.NE.zero )
     $   CALL slascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
*
   10 CONTINUE
*
      IF( lwork.GT.0 )
     $   work( 1 ) = real( lwmin )
      IF( liwork.GT.0 )
     $   iwork( 1 ) = liwmin
      RETURN
*
*     End of PSSTEDC
*
      END