pslaed3()

subroutine pslaed3	(	integer	ictxt,
		integer	k,
		integer	n,
		integer	nb,
		real, dimension( * )	d,
		integer	drow,
		integer	dcol,
		real	rho,
		real, dimension( * )	dlamda,
		real, dimension( * )	w,
		real, dimension( * )	z,
		real, dimension( ldu, * )	u,
		integer	ldu,
		real, dimension( * )	buf,
		integer, dimension( * )	indx,
		integer, dimension( * )	indcol,
		integer, dimension( * )	indrow,
		integer, dimension( * )	indxr,
		integer, dimension( * )	indxc,
		integer, dimension( 0: npcol-1, 4 )	ctot,
		integer	npcol,
		integer	info
	)

Definition at line 1 of file pslaed3.f.

*
*  -- 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            DCOL, DROW, ICTXT, INFO, K, LDU, N, NB, NPCOL
      REAL               RHO
*     ..
*     .. Array Arguments ..
      INTEGER            CTOT( 0: NPCOL-1, 4 ), INDCOL( * ),
     $                   INDROW( * ), INDX( * ), INDXC( * ), INDXR( * )
      REAL               BUF( * ), D( * ), DLAMDA( * ), U( LDU, * ),
     $                   W( * ), Z( * )
*     ..
*
*  Purpose
*  =======
*
*  PSLAED3 finds the roots of the secular equation, as defined by the
*  values in D, W, and RHO, between 1 and K.  It makes the
*  appropriate calls to SLAED4
*
*  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.
*
*  Arguments
*  =========
*
*  ICTXT  (global input) INTEGER
*         The BLACS context handle, indicating the global context of
*         the operation on the matrix. The context itself is global.
*
*  K      (output) INTEGER
*         The number of non-deflated eigenvalues, and the order of the
*         related secular equation. 0 <= K <=N.
*
*  N      (input) INTEGER
*         The dimension of the symmetric tridiagonal matrix.  N >= 0.
*
*  NB      (global input) INTEGER
*          The blocking factor used to distribute the columns of the
*          matrix. NB >= 1.
*
*  D      (input/output) REAL array, dimension (N)
*         On entry, D contains the eigenvalues of the two submatrices to
*         be combined.
*         On exit, D contains the trailing (N-K) updated eigenvalues
*         (those which were deflated) sorted into increasing order.
*
*  DROW   (global input) INTEGER
*          The process row over which the first row of the matrix D is
*          distributed. 0 <= DROW < NPROW.
*
*  DCOL   (global input) INTEGER
*          The process column over which the first column of the
*          matrix D is distributed. 0 <= DCOL < NPCOL.
*
*  RHO    (global input/output) REAL
*         On entry, the off-diagonal element associated with the rank-1
*         cut which originally split the two submatrices which are now
*         being recombined.
*         On exit, RHO has been modified to the value required by
*         PSLAED3.
*
*  DLAMDA (global output) REAL array, dimension (N)
*         A copy of the first K eigenvalues which will be used by
*         SLAED4 to form the secular equation.
*
*  W      (global output) REAL array, dimension (N)
*         The first k values of the final deflation-altered z-vector
*         which will be passed to SLAED4.
*
*  Z      (global input) REAL array, dimension (N)
*         On entry, Z contains the updating vector (the last
*         row of the first sub-eigenvector matrix and the first row of
*         the second sub-eigenvector matrix).
*         On exit, the contents of Z have been destroyed by the updating
*         process.
*
*  U     (global output) REAL array
*         global dimension (N, N), local dimension (LDU, NQ).
*         (See PSLAED0 for definition of NQ.)
*         Q  contains the orthonormal eigenvectors of the symmetric
*         tridiagonal matrix.
*
*  LDU    (input) INTEGER
*         The leading dimension of the array U.
*
*  BUF    (workspace) REAL array, dimension 3*N
*
*
*  INDX   (workspace) INTEGER array, dimension (N)
*         The permutation used to sort the contents of DLAMDA into
*         ascending order.
*
*  INDCOL (workspace) INTEGER array, dimension (N)
*
*
*  INDROW (workspace) INTEGER array, dimension (N)
*
*
*  INDXR (workspace) INTEGER array, dimension (N)
*
*
*  INDXC (workspace) INTEGER array, dimension (N)
*
*  CTOT   (workspace) INTEGER array, dimension( NPCOL, 4)
*
*  NPCOL   (global input) INTEGER
*          The total number of columns over which the distributed
*           submatrix is distributed.
*
*  INFO   (output) INTEGER
*          = 0:  successful exit.
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
*          > 0:  The algorithm failed to compute the ith eigenvalue.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE
      parameter( one = 1.0e0 )
*     ..
*     .. Local Scalars ..
      INTEGER            COL, GI, I, IINFO, IIU, IPD, IU, J, JJU, JU,
     $                   KK, KL, KLC, KLR, MYCOL, MYKL, MYKLR, MYROW,
     $                   NPROW, PDC, PDR, ROW
      REAL               AUX, TEMP
*     ..
*     .. External Functions ..
      INTEGER            INDXG2L
      REAL               SLAMC3, SNRM2
      EXTERNAL           indxg2l, slamc3, snrm2
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, scopy, sgebr2d, sgebs2d,
     $                   sgerv2d, sgesd2d, slaed4
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          mod, sign, sqrt
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      iinfo = 0
*
*     Quick return if possible
*
      IF( k.EQ.0 )
     $   RETURN
*
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
      row = drow
      col = dcol
      DO 20 i = 1, n, nb
         DO 10 j = 0, nb - 1
            indrow( i+j ) = row
            indcol( i+j ) = col
   10    CONTINUE
         row = mod( row+1, nprow )
         col = mod( col+1, npcol )
   20 CONTINUE
*
      mykl = ctot( mycol, 1 ) + ctot( mycol, 2 ) + ctot( mycol, 3 )
      klr = mykl / nprow
      IF( myrow.EQ.drow ) THEN
         myklr = klr + mod( mykl, nprow )
      ELSE
         myklr = klr
      END IF
      pdc = 1
      col = dcol
   30 CONTINUE
      IF( mycol.NE.col ) THEN
         pdc = pdc + ctot( col, 1 ) + ctot( col, 2 ) + ctot( col, 3 )
         col = mod( col+1, npcol )
         GO TO 30
      END IF
      pdr = pdc
      kl = klr + mod( mykl, nprow )
      row = drow
   40 CONTINUE
      IF( myrow.NE.row ) THEN
         pdr = pdr + kl
         kl = klr
         row = mod( row+1, nprow )
         GO TO 40
      END IF
*
      DO 50 i = 1, k
         dlamda( i ) = slamc3( dlamda( i ), dlamda( i ) ) - dlamda( i )
         z( i ) = one
   50 CONTINUE
      IF( myklr.GT.0 ) THEN
         kk = pdr
         DO 80 i = 1, myklr
            CALL slaed4( k, kk, dlamda, w, buf, rho, buf( k+i ), iinfo )
            IF( iinfo.NE.0 ) THEN
               info = kk
            END IF
*
*     ..Compute part of z
*
            DO 60 j = 1, kk - 1
               z( j ) = z( j )*( buf( j ) /
     $                  ( dlamda( j )-dlamda( kk ) ) )
   60       CONTINUE
            z( kk ) = z( kk )*buf( kk )
            DO 70 j = kk + 1, k
               z( j ) = z( j )*( buf( j ) /
     $                  ( dlamda( j )-dlamda( kk ) ) )
   70       CONTINUE
            kk = kk + 1
   80    CONTINUE
*
         IF( myrow.NE.drow ) THEN
            CALL scopy( k, z, 1, buf, 1 )
            CALL sgesd2d( ictxt, k+myklr, 1, buf, k+myklr, drow, mycol )
         ELSE
            ipd = 2*k + 1
            CALL scopy( myklr, buf( k+1 ), 1, buf( ipd ), 1 )
            IF( klr.GT.0 ) THEN
               ipd = myklr + ipd
               row = mod( drow+1, nprow )
               DO 100 i = 1, nprow - 1
                  CALL sgerv2d( ictxt, k+klr, 1, buf, k+klr, row,
     $                          mycol )
                  CALL scopy( klr, buf( k+1 ), 1, buf( ipd ), 1 )
                  DO 90 j = 1, k
                     z( j ) = z( j )*buf( j )
   90             CONTINUE
                  ipd = ipd + klr
                  row = mod( row+1, nprow )
  100          CONTINUE
            END IF
         END IF
      END IF
*
      IF( myrow.EQ.drow ) THEN
         IF( mycol.NE.dcol .AND. mykl.NE.0 ) THEN
            CALL scopy( k, z, 1, buf, 1 )
            CALL scopy( mykl, buf( 2*k+1 ), 1, buf( k+1 ), 1 )
            CALL sgesd2d( ictxt, k+mykl, 1, buf, k+mykl, myrow, dcol )
         ELSE IF( mycol.EQ.dcol ) THEN
            ipd = 2*k + 1
            col = dcol
            kl = mykl
            DO 120 i = 1, npcol - 1
               ipd = ipd + kl
               col = mod( col+1, npcol )
               kl = ctot( col, 1 ) + ctot( col, 2 ) + ctot( col, 3 )
               IF( kl.NE.0 ) THEN
                  CALL sgerv2d( ictxt, k+kl, 1, buf, k+kl, myrow, col )
                  CALL scopy( kl, buf( k+1 ), 1, buf( ipd ), 1 )
                  DO 110 j = 1, k
                     z( j ) = z( j )*buf( j )
  110             CONTINUE
               END IF
  120       CONTINUE
            DO 130 i = 1, k
               z( i ) = sign( sqrt( -z( i ) ), w( i ) )
  130       CONTINUE
*
         END IF
      END IF
*
*     Diffusion
*
      IF( myrow.EQ.drow .AND. mycol.EQ.dcol ) THEN
         CALL scopy( k, z, 1, buf, 1 )
         CALL scopy( k, buf( 2*k+1 ), 1, buf( k+1 ), 1 )
         CALL sgebs2d( ictxt, 'All', ' ', 2*k, 1, buf, 2*k )
      ELSE
         CALL sgebr2d( ictxt, 'All', ' ', 2*k, 1, buf, 2*k, drow, dcol )
         CALL scopy( k, buf, 1, z, 1 )
      END IF
*
*     Copy of D at the good place
*
      klc = 0
      klr = 0
      DO 140 i = 1, k
         gi = indx( i )
         d( gi ) = buf( k+i )
         col = indcol( gi )
         row = indrow( gi )
         IF( col.EQ.mycol ) THEN
            klc = klc + 1
            indxc( klc ) = i
         END IF
         IF( row.EQ.myrow ) THEN
            klr = klr + 1
            indxr( klr ) = i
         END IF
  140 CONTINUE
*
*     Compute eigenvectors of the modified rank-1 modification.
*
      IF( mykl.NE.0 ) THEN
         DO 180 j = 1, mykl
            kk = indxc( j )
            ju = indx( kk )
            jju = indxg2l( ju, nb, j, j, npcol )
            CALL slaed4( k, kk, dlamda, w, buf, rho, aux, iinfo )
            IF( iinfo.NE.0 ) THEN
               info = kk
            END IF
            IF( k.EQ.1 .OR. k.EQ.2 ) THEN
               DO 150 i = 1, klr
                  kk = indxr( i )
                  iu = indx( kk )
                  iiu = indxg2l( iu, nb, j, j, nprow )
                  u( iiu, jju ) = buf( kk )
  150          CONTINUE
               GO TO 180
            END IF
*
            DO 160 i = 1, k
               buf( i ) = z( i ) / buf( i )
  160       CONTINUE
            temp = snrm2( k, buf, 1 )
            DO 170 i = 1, klr
               kk = indxr( i )
               iu = indx( kk )
               iiu = indxg2l( iu, nb, j, j, nprow )
               u( iiu, jju ) = buf( kk ) / temp
  170       CONTINUE
  180    CONTINUE
      END IF
*
  190 CONTINUE
      RETURN
*
*     End of PSLAED3
*

Here is the caller graph for this function: