pdlaed0.f

Go to the documentation of this file.

      SUBROUTINE pdlaed0( N, D, E, Q, IQ, JQ, DESCQ, WORK, IWORK, INFO )
*
*  -- ScaLAPACK auxiliary routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     December 31, 1998
*
*     .. Scalar Arguments ..
      INTEGER            INFO, IQ, JQ, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCQ( * ), IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), Q( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  PDLAED0 computes all eigenvalues and corresponding eigenvectors of a
*  symmetric tridiagonal matrix using the divide and conquer method.
*
*
*  Arguments
*  =========
*
*  N       (global input) INTEGER
*          The order of the tridiagonal matrix T.  N >= 0.
*
*  D       (global input/output) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (N-1)
*          On entry, the subdiagonal elements of the tridiagonal matrix.
*          On exit, E has been destroyed.
*
*  Q       (local output) DOUBLE PRECISION array,
*          global dimension (N, N),
*          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 ) DOUBLE PRECISION array, dimension (LWORK)
*          LWORK = 6*N + 2*NP*NQ, with
*          NP = NUMROC( N, MB_Q, MYROW, IQROW, NPROW )
*          NQ = NUMROC( N, NB_Q, MYCOL, IQCOL, NPCOL )
*          IQROW = INDXG2P( IQ, NB_Q, MYROW, RSRC_Q, NPROW )
*          IQCOL = INDXG2P( JQ, MB_Q, MYCOL, CSRC_Q, NPCOL )
*
*  IWORK   (local workspace/output) INTEGER array, dimension (LIWORK)
*          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).
*
*  =====================================================================
*
*     .. 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 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, ID, IDCOL, IDROW, IID, IINFO, IIQ, IM1, IM2,
     $                   IPQ, IQCOL, IQROW, J, JJD, JJQ, LDQ, MATSIZ,
     $                   MYCOL, MYROW, N1, NB, NBL, NBL1, NPCOL, NPROW,
     $                   SUBPBS, TSUBPBS
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, dgebr2d, dgebs2d, dgerv2d,
     $                   dgesd2d, dsteqr, infog2l, pdlaed1, pxerbla
*     ..
*     .. External Functions ..
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, min
*     ..
*     .. 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 )
      info = 0
      IF( descq( nb_ ).GT.n .OR. n.LT.2 )
     $   info = -1
      IF( info.NE.0 ) THEN
         CALL pxerbla( descq( ctxt_ ), 'PDLAED0', -info )
         RETURN
      END IF
*
      nb = descq( nb_ )
      ldq = descq( lld_ )
      CALL infog2l( iq, jq, descq, nprow, npcol, myrow, mycol, iiq, jjq,
     $              iqrow, iqcol )
*
*     Determine the size and placement of the submatrices, and save in
*     the leading elements of IWORK.
*
      tsubpbs = ( n-1 ) / nb + 1
      iwork( 1 ) = tsubpbs
      subpbs = 1
   10 CONTINUE
      IF( iwork( subpbs ).GT.1 ) THEN
         DO 20 j = subpbs, 1, -1
            iwork( 2*j ) = ( iwork( j )+1 ) / 2
            iwork( 2*j-1 ) = iwork( j ) / 2
   20    CONTINUE
         subpbs = 2*subpbs
         GO TO 10
      END IF
      DO 30 j = 2, subpbs
         iwork( j ) = iwork( j ) + iwork( j-1 )
   30 CONTINUE
*
*     Divide the matrix into TSUBPBS submatrices of size at most NB
*     using rank-1 modifications (cuts).
*
      DO 40 i = nb + 1, n, nb
         im1 = i - 1
         d( im1 ) = d( im1 ) - abs( e( im1 ) )
         d( i ) = d( i ) - abs( e( im1 ) )
   40 CONTINUE
*
*     Solve each submatrix eigenproblem at the bottom of the divide and
*     conquer tree. D is the same on each process.
*
      DO 50 id = 1, n, nb
         CALL infog2l( iq-1+id, jq-1+id, descq, nprow, npcol, myrow,
     $                 mycol, iid, jjd, idrow, idcol )
         matsiz = min( nb, n-id+1 )
         IF( myrow.EQ.idrow .AND. mycol.EQ.idcol ) THEN
            ipq = iid + ( jjd-1 )*ldq
            CALL dsteqr( 'I', matsiz, d( id ), e( id ), q( ipq ), ldq,
     $                   work, info )
            IF( info.NE.0 ) THEN
               CALL pxerbla( descq( ctxt_ ), 'DSTEQR', -info )
               RETURN
            END IF
            IF( myrow.NE.iqrow .OR. mycol.NE.iqcol ) THEN
               CALL dgesd2d( descq( ctxt_ ), matsiz, 1, d( id ), matsiz,
     $                       iqrow, iqcol )
            END IF
         ELSE IF( myrow.EQ.iqrow .AND. mycol.EQ.iqcol ) THEN
            CALL dgerv2d( descq( ctxt_ ), matsiz, 1, d( id ), matsiz,
     $                    idrow, idcol )
         END IF
   50 CONTINUE
*
      IF( myrow.EQ.iqrow .AND. mycol.EQ.iqcol ) THEN
         CALL dgebs2d( descq( ctxt_ ), 'A', ' ', n, 1, d, n )
      ELSE
         CALL dgebr2d( descq( ctxt_ ), 'A', ' ', n, 1, d, n, iqrow,
     $                 iqcol )
      END IF
*
*     Successively merge eigensystems of adjacent submatrices
*     into eigensystem for the corresponding larger matrix.
*
*     while ( SUBPBS > 1 )
*
   60 CONTINUE
      IF( subpbs.GT.1 ) THEN
         im2 = subpbs - 2
         DO 80 i = 0, im2, 2
            IF( i.EQ.0 ) THEN
               nbl = iwork( 2 )
               nbl1 = iwork( 1 )
               IF( nbl1.EQ.0 )
     $            GO TO 70
               id = 1
               matsiz = min( n, nbl*nb )
               n1 = nbl1*nb
            ELSE
               nbl = iwork( i+2 ) - iwork( i )
               nbl1 = nbl / 2
               IF( nbl1.EQ.0 )
     $            GO TO 70
               id = iwork( i )*nb + 1
               matsiz = min( nb*nbl, n-id+1 )
               n1 = nbl1*nb
            END IF
*
*     Merge lower order eigensystems (of size N1 and MATSIZ - N1)
*     into an eigensystem of size MATSIZ.
*
            CALL pdlaed1( matsiz, n1, d( id ), id, q, iq, jq, descq,
     $                    e( id+n1-1 ), work, iwork( subpbs+1 ), iinfo )
            IF( iinfo.NE.0 ) THEN
               info = iinfo*( n+1 ) + id
            END IF
*
   70       CONTINUE
            iwork( i / 2+1 ) = iwork( i+2 )
   80    CONTINUE
         subpbs = subpbs / 2
*
         GO TO 60
      END IF
*
*     end while
*
   90 CONTINUE
      RETURN
*
*     End of PDLAED0
*
      END