      SUBROUTINE PDGESVD( JOBU, JOBVT, M, N, A, IA, JA, DESCA, S, U, IU,
     $                    JU, DESCU, VT, IVT, JVT, DESCVT, WORK, LWORK,
     $                    INFO )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 15, 1997
*
*     .. Scalar Arguments ..
      CHARACTER          JOBU, JOBVT
      INTEGER            IA, INFO, IU, IVT, JA, JU, JVT, LWORK, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCU( * ), DESCVT( * )
      DOUBLE PRECISION   A( * ), S( * ), U( * ), VT( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  PDGESVD computes the singular value decomposition (SVD) of an
*  M-by-N matrix A, optionally computing the left and/or right
*  singular vectors. The SVD is written as
*
*       A = U * SIGMA * transpose(V)
*
*  where SIGMA is an M-by-N matrix which is zero except for its
*  min(M,N) diagonal elements, U is an M-by-M orthogonal matrix, and
*  V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
*  are the singular values of A and the columns of U and V are the
*  corresponding right and left singular vectors, respectively. The
*  singular values are returned in array S in decreasing order and
*  only the first min(M,N) columns of U and rows of VT = V**T are
*  computed.
*
*  Notes
*  =====
*  Each global data object is described by an associated description
*  vector. This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector DESCA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
*                                 DTYPE_A = 1.
*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) DESCA( M_ )    The number of rows in the global
*                                 array A.
*  N_A    (global) DESCA( N_ )    The number of columns in the global
*                                 array A.
*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
*                                 the rows of the array.
*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
*                                 the columns of the array.
*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix, and
*  assume that its process grid has dimension r x c. LOCr( K ) denotes
*  the number of elements of K that a process would receive if K were
*  distributed over the r processes of its process column. Similarly,
*  LOCc( K ) denotes the number of elements of K that a process would
*  receive if K were distributed over the c processes of its process
*  row. The values of LOCr() and LOCc() may be determined via a call
*  to the ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*  Arguments
*  =========
*
*          MP = number of local rows in A and U
*          NQ = number of local columns in A and VT
*          SIZE = min( M, N )
*          SIZEQ = number of local columns in U
*          SIZEP = number of local rows in VT
*
*  JOBU    (global input) CHARACTER*1
*          Specifies options for computing U:
*          = 'V':  the first SIZE columns of U (the left singular
*                  vectors) are returned in the array U;
*          = 'N':  no columns of U (no left singular vectors) are
*                  computed.
*
*  JOBVT   (global input) CHARACTER*1
*          Specifies options for computing V**T:
*          = 'V':  the first SIZE rows of V**T (the right singular
*                  vectors) are returned in the array VT;
*          = 'N':  no rows of V**T (no right singular vectors) are
*                  computed.
*
*  M       (global input) INTEGER
*          The number of rows of the input matrix A.  M >= 0.
*
*  N       (global input) INTEGER
*          The number of columns of the input matrix A.  N >= 0.
*
*  A       (local input/workspace) block cyclic DOUBLE PRECISION array,
*          global dimension (M, N), local dimension (MP, NQ)
*          On exit, the contents of A are destroyed.
*
*  IA      (global input) INTEGER
*          The row index in the global array A indicating the first
*          row of sub( A ).
*
*  JA      (global input) INTEGER
*          The column index in the global array A indicating the
*          first column of sub( A ).
*
*  DESCA   (global input) INTEGER array of dimension DLEN_
*          The array descriptor for the distributed matrix A.
*
*  S       (global output) DOUBLE PRECISION array, dimension SIZE
*          The singular values of A, sorted so that S(i) >= S(i+1).
*
*  U       (local output) DOUBLE PRECISION array, local dimension
*          (MP, SIZEQ), global dimension (M, SIZE)
*          if JOBU = 'V', U contains the first min(m,n) columns of U
*          if JOBU = 'N', U is not referenced.
*
*  IU      (global input) INTEGER
*          The row index in the global array U indicating the first
*          row of sub( U ).
*
*  JU      (global input) INTEGER
*          The column index in the global array U indicating the
*          first column of sub( U ).
*
*  DESCU   (global input) INTEGER array of dimension DLEN_
*          The array descriptor for the distributed matrix U.
*
*  VT      (local output) DOUBLE PRECISION array, local dimension
*          (SIZEP, NQ), global dimension (SIZE, N).
*          If JOBVT = 'V', VT contains the first SIZE rows of
*          V**T. If JOBVT = 'N', VT is not referenced.
*
*  IVT     (global input) INTEGER
*          The row index in the global array VT indicating the first
*          row of sub( VT ).
*
*  JVT     (global input) INTEGER
*          The column index in the global array VT indicating the
*          first column of sub( VT ).
*
*  DESCVT   (global input) INTEGER array of dimension DLEN_
*          The array descriptor for the distributed matrix VT.
*
*  WORK    (local workspace/output) DOUBLE PRECISION array, dimension
*          (LWORK)
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
*  LWORK   (local input) INTEGER
*          The dimension of the array WORK.
*
*          LWORK > 2 + 6*SIZEB + MAX(WATOBD, WBDTOSVD),
*
*          where SIZEB = MAX(M,N), and WATOBD and WBDTOSVD refer,
*          respectively, to the workspace required to bidiagonalize
*          the matrix A and to go from the bidiagonal matrix to the
*          singular value decomposition U*S*VT.
*
*          For WATOBD, the following holds:
*
*          WATOBD = MAX(MAX(WPDLANGE,WPDGEBRD),
*                       MAX(WPDLARED2D,WPDLARED1D)),
*
*          where WPDLANGE, WPDLARED1D, WPDLARED2D, WPDGEBRD are the
*          workspaces required respectively for the subprograms
*          PDLANGE, PDLARED1D, PDLARED2D, PDGEBRD. Using the
*          standard notation
*
*          MP = NUMROC( M, MB, MYROW, DESCA( CTXT_ ), NPROW),
*          NQ = NUMROC( N, NB, MYCOL, DESCA( LLD_ ), NPCOL),
*
*          the workspaces required for the above subprograms are
*
*          WPDLANGE = MP,
*          WPDLARED1D = NQ0,
*          WPDLARED2D = MP0,
*          WPDGEBRD = NB*(MP + NQ + 1) + NQ,
*
*          where NQ0 and MP0 refer, respectively, to the values obtained
*          at MYCOL = 0 and MYROW = 0. In general, the upper limit for
*          the workspace is given by a workspace required on
*          processor (0,0):
*
*          WATOBD <= NB*(MP0 + NQ0 + 1) + NQ0.
*
*          In case of a homogeneous process grid this upper limit can
*          be used as an estimate of the minimum workspace for every
*          processor.
*
*          For WBDTOSVD, the following holds:
*
*          WBDTOSVD = SIZE*(WANTU*NRU + WANTVT*NCVT) +
*                     MAX(WDBDSQR,
*                         MAX(WANTU*WPDORMBRQLN, WANTVT*WPDORMBRPRT)),
*
*          where
*
*                          1, if left(right) singular vectors are wanted
*          WANTU(WANTVT) =
*                          0, otherwise
*
*          and WDBDSQR, WPDORMBRQLN and WPDORMBRPRT refer respectively
*          to the workspace required for the subprograms DBDSQR,
*          PDORMBR(QLN), and PDORMBR(PRT), where QLN and PRT are the
*          values of the arguments VECT, SIDE, and TRANS in the call
*          to PDORMBR. NRU is equal to the local number of rows of
*          the matrix U when distributed 1-dimensional "column" of
*          processes. Analogously, NCVT is equal to the local number
*          of columns of the matrix VT when distributed across
*          1-dimensional "row" of processes. Calling the LAPACK
*          procedure DBDSQR requires
*
*          WDBDSQR = MAX(1, 2*SIZE + (2*SIZE - 4)*MAX(WANTU, WANTVT))
*
*          on every processor. Finally,
*
*          WPDORMBRQLN = MAX( (NB*(NB-1))/2, (SIZEQ+MP)*NB)+NB*NB,
*          WPDORMBRPRT = MAX( (MB*(MB-1))/2, (SIZEP+NQ)*MB )+MB*MB,
*
*          If LWORK = -1, then 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.
*
*  INFO    (global output) INTEGER
*          = 0:  successful exit.
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
*          > 0:  if SBDSQR did not converge
*                If INFO = MIN(M,N) + 1, then PDGESVD has detected
*                heterogeneity by finding that eigenvalues were not
*                identical across the process grid. In this case, the
*                accuracy of the results from PDGESVD cannot be
*                guaranteed.
*
*  =====================================================================
*
*  The results of PDGEBRD, and therefore PDGESVD, may vary slightly
*  from run to run with the same input data. If repeatability is an
*  issue, call BLACS_SET with the appropriate option after defining
*  the process grid.
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
     $                   MB_, NB_, RSRC_, CSRC_, LLD_, ITHVAL
      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, ITHVAL = 10 )
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
*     ..
*     .. Local Scalars ..
      CHARACTER          UPLO
      INTEGER            CONTEXTC, CONTEXTR, I, INDD, INDD2, INDE,
     $                   INDE2, INDTAUP, INDTAUQ, INDU, INDV, INDWORK,
     $                   IOFFD, IOFFE, ISCALE, J, K, LDU, LDVT, LLWORK,
     $                   LWMIN, MAXIM, MB, MP, MYPCOL, MYPCOLC, MYPCOLR,
     $                   MYPROW, MYPROWC, MYPROWR, NB, NCVT, NPCOL,
     $                   NPCOLC, NPCOLR, NPROCS, NPROW, NPROWC, NPROWR,
     $                   NQ, NRU, SIZE, SIZEB, SIZEP, SIZEPOS, SIZEQ,
     $                   WANTU, WANTVT, WATOBD, WBDTOSVD, WDBDSQR,
     $                   WPDGEBRD, WPDLANGE, WPDORMBRPRT, WPDORMBRQLN
      DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
     $                   SMLNUM
*     ..
*     .. Local Arrays ..
      INTEGER            DESCTU( DLEN_ ), DESCTVT( DLEN_ ), IDUM1( 3 ),
     $                   IDUM2( 3 )
      DOUBLE PRECISION   C( 1, 1 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            NUMROC
      DOUBLE PRECISION   PDLAMCH, PDLANGE
      EXTERNAL           LSAME, NUMROC, PDLAMCH, PDLANGE
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GET, BLACS_GRIDEXIT, BLACS_GRIDINFO,
     $                   BLACS_GRIDINIT, CHK1MAT, DBDSQR, DESCINIT,
     $                   DGAMN2D, DGAMX2D, DSCAL, IGAMX2D, IGEBR2D,
     $                   IGEBS2D, PCHK1MAT, PDGEBRD, PDGEMR2D,
     $                   PDLARED1D, PDLARED2D, PDLASCL, PDLASET,
     $                   PDORMBR, PXERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          DBLE, MAX, MIN, SQRT
*     ..
*     .. Executable Statements ..
*     This is just to keep ftnchek happy
      IF( BLOCK_CYCLIC_2D*DTYPE_*LLD_*MB_*M_*NB_*N_.LT.0 )
     $   RETURN
*
      CALL BLACS_GRIDINFO( DESCA( CTXT_ ), NPROW, NPCOL, MYPROW,
     $                     MYPCOL )
      ISCALE = 0
      INFO = 0
*
      IF( NPROW.EQ.-1 ) THEN
         INFO = -( 800+CTXT_ )
      ELSE
*
         SIZE = MIN( M, N )
         SIZEB = MAX( M, N )
         NPROCS = NPROW*NPCOL
         IF( M.GE.N ) THEN
            IOFFD = JA - 1
            IOFFE = IA - 1
            SIZEPOS = 1
         ELSE
            IOFFD = IA - 1
            IOFFE = JA - 1
            SIZEPOS = 3
         END IF
*
         IF( LSAME( JOBU, 'V' ) ) THEN
            WANTU = 1
         ELSE
            WANTU = 0
         END IF
         IF( LSAME( JOBVT, 'V' ) ) THEN
            WANTVT = 1
         ELSE
            WANTVT = 0
         END IF
*
         CALL CHK1MAT( M, 3, N, 4, IA, JA, DESCA, 8, INFO )
         IF( WANTU.EQ.1 ) THEN
            CALL CHK1MAT( M, 3, SIZE, SIZEPOS, IU, JU, DESCU, 13, INFO )
         END IF
         IF( WANTVT.EQ.1 ) THEN
            CALL CHK1MAT( SIZE, SIZEPOS, N, 4, IVT, JVT, DESCVT, 17,
     $                    INFO )
         END IF
         CALL IGAMX2D( DESCA( CTXT_ ), 'A', ' ', 1, 1, INFO, 1, 1, 1,
     $                 -1, -1, 0 )
*
         IF( INFO.EQ.0 ) THEN
*
*           Set up pointers into the WORK array.
*
            INDD = 2
            INDE = INDD + SIZEB + IOFFD
            INDD2 = INDE + SIZEB + IOFFE
            INDE2 = INDD2 + SIZEB + IOFFD
            INDTAUQ = INDE2 + SIZEB + IOFFE
            INDTAUP = INDTAUQ + SIZEB + JA - 1
            INDWORK = INDTAUP + SIZEB + IA - 1
            LLWORK = LWORK - INDWORK + 1
*
*           Initialize contexts for "column" and "row" process matrices.
*
            CALL BLACS_GET( DESCA( CTXT_ ), 10, CONTEXTC )
            CALL BLACS_GRIDINIT( CONTEXTC, 'R', NPROCS, 1 )
            CALL BLACS_GRIDINFO( CONTEXTC, NPROWC, NPCOLC, MYPROWC,
     $                           MYPCOLC )
            CALL BLACS_GET( DESCA( CTXT_ ), 10, CONTEXTR )
            CALL BLACS_GRIDINIT( CONTEXTR, 'R', 1, NPROCS )
            CALL BLACS_GRIDINFO( CONTEXTR, NPROWR, NPCOLR, MYPROWR,
     $                           MYPCOLR )
*
*           Set local dimensions of matrices (this is for MB=NB=1).
*
            NRU = NUMROC( SIZE, 1, MYPROWC, 0, NPROCS )
            NCVT = NUMROC( SIZE, 1, MYPCOLR, 0, NPROCS )
            NB = DESCA( NB_ )
            MB = DESCA( MB_ )
            MP = NUMROC( M, MB, MYPROW, DESCA( RSRC_ ), NPROW )
            NQ = NUMROC( N, NB, MYPCOL, DESCA( CSRC_ ), NPCOL )
            IF( WANTVT.EQ.1 ) THEN
               SIZEP = NUMROC( SIZE, DESCVT( MB_ ), MYPROW,
     $                 DESCVT( RSRC_ ), NPROW )
            ELSE
               SIZEP = 0
            END IF
            IF( WANTU.EQ.1 ) THEN
               SIZEQ = NUMROC( SIZE, DESCU( NB_ ), MYPCOL,
     $                 DESCU( CSRC_ ), NPCOL )
            ELSE
               SIZEQ = 0
            END IF
*
*           Transmit MAX(NQ0, MP0).
*
            IF( MYPROW.EQ.0 .AND. MYPCOL.EQ.0 ) THEN
               MAXIM = MAX( NQ, MP )
               CALL IGEBS2D( DESCA( CTXT_ ), 'All', ' ', 1, 1, MAXIM,
     $                       1 )
            ELSE
               CALL IGEBR2D( DESCA( CTXT_ ), 'All', ' ', 1, 1, MAXIM, 1,
     $                       0, 0 )
            END IF
*
            WPDLANGE = MP
            WPDGEBRD = NB*( MP+NQ+1 ) + NQ
            WATOBD = MAX( MAX( WPDLANGE, WPDGEBRD ), MAXIM )
*
            WDBDSQR = MAX( 1, 2*SIZE+( 2*SIZE-4 )*MAX( WANTU, WANTVT ) )
            WPDORMBRQLN = MAX( ( NB*( NB-1 ) ) / 2, ( SIZEQ+MP )*NB ) +
     $                    NB*NB
            WPDORMBRPRT = MAX( ( MB*( MB-1 ) ) / 2, ( SIZEP+NQ )*MB ) +
     $                    MB*MB
            WBDTOSVD = SIZE*( WANTU*NRU+WANTVT*NCVT ) +
     $                 MAX( WDBDSQR, MAX( WANTU*WPDORMBRQLN,
     $                 WANTVT*WPDORMBRPRT ) )
*
*           Finally, calculate required workspace.
*
            LWMIN = 2 + 6*SIZEB + MAX( WATOBD, WBDTOSVD )
            WORK( 1 ) = DBLE( LWMIN )
*
            IF( WANTU.NE.1 .AND. .NOT.( LSAME( JOBU, 'N' ) ) ) THEN
               INFO = -1
            ELSE IF( WANTVT.NE.1 .AND. .NOT.( LSAME( JOBVT, 'N' ) ) )
     $                THEN
               INFO = -2
            ELSE IF( LWORK.LT.LWMIN .AND. LWORK.NE.-1 ) THEN
               INFO = -19
            END IF
*
         END IF
*
         IDUM1( 1 ) = WANTU
         IDUM1( 2 ) = WANTVT
         IF( LWORK.EQ.-1 ) THEN
            IDUM1( 3 ) = -1
         ELSE
            IDUM1( 3 ) = 1
         END IF
         IDUM2( 1 ) = 1
         IDUM2( 2 ) = 2
         IDUM2( 3 ) = 19
         CALL PCHK1MAT( M, 3, N, 4, IA, JA, DESCA, 8, 3, IDUM1, IDUM2,
     $                  INFO )
         IF( INFO.EQ.0 ) THEN
            IF( WANTU.EQ.1 ) THEN
               CALL PCHK1MAT( M, 3, N, 4, IU, JU, DESCU, 13, 0, IDUM1,
     $                        IDUM2, INFO )
            END IF
            IF( WANTVT.EQ.1 ) THEN
               CALL PCHK1MAT( M, 3, N, 4, IVT, JVT, DESCVT, 17, 0,
     $                        IDUM1, IDUM2, INFO )
            END IF
         END IF
*
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL PXERBLA( DESCA( CTXT_ ), 'PDGESVD', -INFO )
         RETURN
      ELSE IF( LWORK.EQ.-1 ) THEN
         GO TO 40
      END IF
*
*     Quick return if possible.
*
      IF( M.LE.0 .OR. N.LE.0 )
     $   GO TO 40
*
*     Get machine constants.
*
      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 ) ) )
*
*     Scale matrix to allowable range, if necessary.
*
      ANRM = PDLANGE( '1', M, 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( 'G', ONE, SIGMA, M, N, A, IA, JA, DESCA, INFO )
      END IF
*
      CALL PDGEBRD( M, N, A, IA, JA, DESCA, WORK( INDD ), WORK( INDE ),
     $              WORK( INDTAUQ ), WORK( INDTAUP ), WORK( INDWORK ),
     $              LLWORK, INFO )
*
*     Copy D and E to all processes.
*     Array D is in local array of dimension:
*     LOCc(JA+MIN(M,N)-1) if M >= N; LOCr(IA+MIN(M,N)-1) otherwise.
*     Array E is in local array of dimension
*     LOCr(IA+MIN(M,N)-1) if M >= N; LOCc(JA+MIN(M,N)-2) otherwise.
*
      IF( M.GE.N ) THEN
*        Distribute D
         CALL PDLARED1D( N+IOFFD, IA, JA, DESCA, WORK( INDD ),
     $                   WORK( INDD2 ), WORK( INDWORK ), LLWORK )
*        Distribute E
         CALL PDLARED2D( M+IOFFE, IA, JA, DESCA, WORK( INDE ),
     $                   WORK( INDE2 ), WORK( INDWORK ), LLWORK )
      ELSE
*        Distribute D
         CALL PDLARED2D( M+IOFFD, IA, JA, DESCA, WORK( INDD ),
     $                   WORK( INDD2 ), WORK( INDWORK ), LLWORK )
*        Distribute E
         CALL PDLARED1D( N+IOFFE, IA, JA, DESCA, WORK( INDE ),
     $                   WORK( INDE2 ), WORK( INDWORK ), LLWORK )
      END IF
*
*     Prepare for calling PDBDSQR.
*
      IF( M.GE.N ) THEN
         UPLO = 'U'
      ELSE
         UPLO = 'L'
      END IF
*
      INDU = INDWORK
      INDV = INDU + SIZE*NRU
      INDWORK = INDV + SIZE*NCVT
*
      LDU = MAX( 1, NRU )
      LDVT = MAX( 1, SIZE )
*
      CALL DESCINIT( DESCTU, SIZE, SIZE, 1, 1, 0, 0, CONTEXTC, LDU,
     $               INFO )
      CALL DESCINIT( DESCTVT, SIZE, SIZE, 1, 1, 0, 0, CONTEXTR, LDVT,
     $               INFO )
*
      IF( WANTU.EQ.1 ) THEN
         CALL PDLASET( 'Full', M, SIZE, ZERO, ONE, WORK( INDU ), 1, 1,
     $                 DESCTU )
      ELSE
         NRU = 0
      END IF
*
      IF( WANTVT.EQ.1 ) THEN
         CALL PDLASET( 'Full', SIZE, N, ZERO, ONE, WORK( INDV ), 1, 1,
     $                 DESCTVT )
      ELSE
         NCVT = 0
      END IF
*
      CALL DBDSQR( UPLO, SIZE, NCVT, NRU, 0, WORK( INDD2+IOFFD ),
     $             WORK( INDE2+IOFFE ), WORK( INDV ), SIZE,
     $             WORK( INDU ), LDU, C, 1, WORK( INDWORK ), INFO )
*
*     Redistribute elements of U and VT in the block-cyclic fashion.
*
      IF( WANTU.EQ.1 )
     $   CALL PDGEMR2D( SIZE, SIZE, WORK( INDU ), 1, 1, DESCTU, U, IU,
     $                  JU, DESCU, DESCU( CTXT_ ) )
*
      IF( WANTVT.EQ.1 )
     $   CALL PDGEMR2D( SIZE, SIZE, WORK( INDV ), 1, 1, DESCTVT, VT,
     $                  IVT, JVT, DESCVT, DESCVT( CTXT_ ) )
*
*     Set to ZERO "non-square" elements of the larger matrices U, VT.
*
      IF( M.GT.N .AND. WANTU.EQ.1 ) THEN
         CALL PDLASET( 'Full', M-SIZE, SIZE, ZERO, ZERO, U, IA+SIZE, JU,
     $                 DESCU )
      ELSE IF( N.GT.M .AND. WANTVT.EQ.1 ) THEN
         CALL PDLASET( 'Full', SIZE, N-SIZE, ZERO, ZERO, VT, IVT,
     $                 JVT+SIZE, DESCVT )
      END IF
*
*     Multiply Householder rotations from bidiagonalized matrix.
*
      IF( WANTU.EQ.1 )
     $   CALL PDORMBR( 'Q', 'L', 'N', M, SIZE, N, A, IA, JA, DESCA,
     $                 WORK( INDTAUQ ), U, IU, JU, DESCU,
     $                 WORK( INDWORK ), LLWORK, INFO )
*
      IF( WANTVT.EQ.1 )
     $   CALL PDORMBR( 'P', 'R', 'T', SIZE, N, M, A, IA, JA, DESCA,
     $                 WORK( INDTAUP ), VT, IVT, JVT, DESCVT,
     $                 WORK( INDWORK ), LLWORK, INFO )
*
*     Copy singular values into output array S.
*
      DO 10 I = 1, SIZE
         S( I ) = WORK( INDD2+IOFFD+I-1 )
   10 CONTINUE
*
*     If matrix was scaled, then rescale singular values appropriately.
*
      IF( ISCALE.EQ.1 ) THEN
         CALL DSCAL( SIZE, ONE / SIGMA, S, 1 )
      END IF
*
*     Compare every ith eigenvalue, or all if there are only a few,
*     across the process grid to check for heterogeneity.
*
      IF( SIZE.LE.ITHVAL ) THEN
         J = SIZE
         K = 1
      ELSE
         J = SIZE / ITHVAL
         K = ITHVAL
      END IF
*
      DO 20 I = 1, J
         WORK( I+INDE ) = S( ( I-1 )*K+1 )
         WORK( I+INDD2 ) = S( ( I-1 )*K+1 )
   20 CONTINUE
*
      CALL DGAMN2D( DESCA( CTXT_ ), 'a', ' ', J, 1, WORK( 1+INDE ), J,
     $              1, 1, -1, -1, 0 )
      CALL DGAMX2D( DESCA( CTXT_ ), 'a', ' ', J, 1, WORK( 1+INDD2 ), J,
     $              1, 1, -1, -1, 0 )
*
      DO 30 I = 1, J
         IF( ( WORK( I+INDE )-WORK( I+INDD2 ) ).NE.ZERO ) THEN
            INFO = SIZE + 1
         END IF
   30 CONTINUE
*
   40 CONTINUE
*
      CALL BLACS_GRIDEXIT( CONTEXTC )
      CALL BLACS_GRIDEXIT( CONTEXTR )
*
*     End of PDGESVD
*
      RETURN
      END
