pdlaed1.f

Go to the documentation of this file.

      SUBROUTINE pdlaed1( N, N1, D, ID, Q, IQ, JQ, DESCQ, RHO, 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            ID, INFO, IQ, JQ, N, N1
      DOUBLE PRECISION   RHO
*     ..
*     .. Array Arguments ..
      INTEGER            DESCQ( * ), IWORK( * )
      DOUBLE PRECISION   D( * ), Q( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  PDLAED1 computes the updated eigensystem of a diagonal
*  matrix after modification by a rank-one symmetric matrix,
*  in parallel.
*
*    T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)
*
*     where Z = Q'u, u is a vector of length N with ones in the
*     N1 and N1 + 1 th elements and zeros elsewhere.
*
*     The eigenvectors of the original matrix are stored in Q, and the
*     eigenvalues are in D.  The algorithm consists of three stages:
*
*        The first stage consists of deflating the size of the problem
*        when there are multiple eigenvalues or if there is a zero in
*        the Z vector.  For each such occurence the dimension of the
*        secular equation problem is reduced by one.  This stage is
*        performed by the routine PDLAED2.
*
*        The second stage consists of calculating the updated
*        eigenvalues. This is done by finding the roots of the secular
*        equation via the routine SLAED4 (as called by PDLAED3).
*        This routine also calculates the eigenvectors of the current
*        problem.
*
*        The final stage consists of computing the updated eigenvectors
*        directly using the updated eigenvalues.  The eigenvectors for
*        the current problem are multiplied with the eigenvectors from
*        the overall problem.
*
*  Arguments
*  =========
*
*  N       (global input) INTEGER
*          The order of the tridiagonal matrix T.  N >= 0.
*
*
*  N1      (input) INTEGER
*          The location of the last eigenvalue in the leading
*          sub-matrix.
*          min(1,N) <= N1 <= N.
*
*  D       (global input/output) DOUBLE PRECISION array, dimension (N)
*          On entry,the eigenvalues of the rank-1-perturbed matrix.
*          On exit, the eigenvalues of the repaired matrix.
*
*  ID      (global input) INTEGER
*          Q's global row/col index, which points to the beginning
*          of the submatrix which is to be operated on.
*
*  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.
*
*  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.
*
*  RHO    (input) DOUBLE PRECISION
*         The subdiagonal entry used to create the rank-1 modification.
*
*  WORK    (local workspace/output) DOUBLE PRECISION array,
*          dimension 6*N + 2*NP*NQ
*
*  IWORK   (local workspace/output) INTEGER array,
*          dimension 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 ith eigenvalue.
*
*  =====================================================================
*
*     .. 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 ..
      INTEGER            COL, COLTYP, IBUF, ICTOT, ICTXT, IDLMDA, IIQ,
     $                   indcol, indrow, indx, indxc, indxp, indxr, inq,
     $                   ipq, ipq2, ipsm, ipu, ipwork, iq1, iq2, iqcol,
     $                   iqq, iqrow, iw, iz, j, jc, jj2c, jjc, jjq, jnq,
     $                   k, ldq, ldq2, ldu, mycol, myrow, nb, nn, nn1,
     $                   nn2, np, npcol, nprow, nq
*     ..
*     .. Local Arrays ..
      INTEGER            DESCQ2( DLEN_ ), DESCU( DLEN_ )
*     ..
*     .. External Functions ..
      INTEGER            NUMROC
      EXTERNAL           numroc
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, dcopy, descinit, infog1l,
     $                   infog2l, pdgemm, pdlaed2, pdlaed3, pdlaedz,
     $                   pdlaset, pxerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max, 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( nprow.EQ.-1 ) THEN
         info = -( 600+ctxt_ )
      ELSE IF( n.LT.0 ) THEN
         info = -1
      ELSE IF( id.GT.descq( n_ ) ) THEN
         info = -4
      ELSE IF( n1.GE.n ) THEN
         info = -2
      END IF
      IF( info.NE.0 ) THEN
         CALL pxerbla( descq( ctxt_ ), 'PDLAED1', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( n.EQ.0 )
     $   RETURN
*
*     The following values are  integer pointers which indicate
*     the portion of the workspace used by a particular array
*     in PDLAED2 and PDLAED3.
*
      ictxt = descq( ctxt_ )
      nb = descq( nb_ )
      ldq = descq( lld_ )
*
      CALL infog2l( iq-1+id, jq-1+id, descq, nprow, npcol, myrow, mycol,
     $              iiq, jjq, iqrow, iqcol )
*
      np = numroc( n, descq( mb_ ), myrow, iqrow, nprow )
      nq = numroc( n, descq( nb_ ), mycol, iqcol, npcol )
*
      ldq2 = max( np, 1 )
      ldu = ldq2
*
      iz = 1
      idlmda = iz + n
      iw = idlmda + n
      ipq2 = iw + n
      ipu = ipq2 + ldq2*nq
      ibuf = ipu + ldu*nq
*     (IBUF est de taille 3*N au maximum)
*
      ictot = 1
      ipsm = ictot + npcol*4
      indx = ipsm + npcol*4
      indxc = indx + n
      indxp = indxc + n
      indcol = indxp + n
      coltyp = indcol + n
      indrow = coltyp + n
      indxr = indrow + n
*
      CALL descinit( descq2, n, n, nb, nb, iqrow, iqcol, ictxt, ldq2,
     $               info )
      CALL descinit( descu, n, n, nb, nb, iqrow, iqcol, ictxt, ldu,
     $               info )
*
*     Form the z-vector which consists of the last row of Q_1 and the
*     first row of Q_2.
*
      ipwork = idlmda
      CALL pdlaedz( n, n1, id, q, iq, jq, ldq, descq, work( iz ),
     $              work( ipwork ) )
*
*     Deflate eigenvalues.
*
      ipq = iiq + ( jjq-1 )*ldq
      CALL pdlaed2( ictxt, k, n, n1, nb, d, iqrow, iqcol, q( ipq ), ldq,
     $              rho, work( iz ), work( iw ), work( idlmda ),
     $              work( ipq2 ), ldq2, work( ibuf ), iwork( ictot ),
     $              iwork( ipsm ), npcol, iwork( indx ), iwork( indxc ),
     $              iwork( indxp ), iwork( indcol ), iwork( coltyp ),
     $              nn, nn1, nn2, iq1, iq2 )
*
*
*     Solve Secular Equation.
*
      IF( k.NE.0 ) THEN
         CALL pdlaset( 'A', n, n, zero, one, work( ipu ), 1, 1, descu )
         CALL pdlaed3( ictxt, k, n, nb, d, iqrow, iqcol, rho,
     $                 work( idlmda ), work( iw ), work( iz ),
     $                 work( ipu ), ldq2, work( ibuf ), iwork( indx ),
     $                 iwork( indcol ), iwork( indrow ), iwork( indxr ),
     $                 iwork( indxc ), iwork( ictot ), npcol, info )
*
*     Compute the updated eigenvectors.
*
         iqq = min( iq1, iq2 )
         IF( nn1.GT.0 ) THEN
            inq = iq - 1 + id
            jnq = jq - 1 + id + iqq - 1
            CALL pdgemm( 'N', 'N', n1, nn, nn1, one, work( ipq2 ), 1,
     $                   iq1, descq2, work( ipu ), iq1, iqq, descu,
     $                   zero, q, inq, jnq, descq )
         END IF
         IF( nn2.GT.0 ) THEN
            inq = iq - 1 + id + n1
            jnq = jq - 1 + id + iqq - 1
            CALL pdgemm( 'N', 'N', n-n1, nn, nn2, one, work( ipq2 ),
     $                   n1+1, iq2, descq2, work( ipu ), iq2, iqq,
     $                   descu, zero, q, inq, jnq, descq )
         END IF
*
         DO 10 j = k + 1, n
            jc = iwork( indx+j-1 )
            CALL infog1l( jq-1+jc, nb, npcol, mycol, iqcol, jjc, col )
            CALL infog1l( jc, nb, npcol, mycol, iqcol, jj2c, col )
            IF( mycol.EQ.col ) THEN
               iq2 = ipq2 + ( jj2c-1 )*ldq2
               inq = ipq + ( jjc-1 )*ldq
               CALL dcopy( np, work( iq2 ), 1, q( inq ), 1 )
            END IF
   10    CONTINUE
      END IF
*
   20 CONTINUE
      RETURN
*
*     End of PDLAED1
*
      END