pcgesvd.f

Go to the documentation of this file.

 
      SUBROUTINE pcgesvd(JOBU,JOBVT,M,N,A,IA,JA,DESCA,S,U,IU,JU,DESCU,
     +                   VT,IVT,JVT,DESCVT,WORK,LWORK,RWORK,INFO)
*
*  -- ScaLAPACK routine (version 1.7) --
*     Univ. of Tennessee, Oak Ridge National Laboratory
*     and Univ. of California Berkeley.
*     Jan 2006
 
*
*     .. Scalar Arguments ..
      CHARACTER JOBU,JOBVT
      INTEGER IA,INFO,IU,IVT,JA,JU,JVT,LWORK,M,N
*     ..
*     .. Array Arguments ..
      INTEGER DESCA(*),DESCU(*),DESCVT(*)
      COMPLEX A(*),U(*),VT(*),WORK(*)
      REAL S(*)
      REAL RWORK(*)
*     ..
*
*  Purpose
*  =======
*
*  PCGESVD 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 COMPLEX
*          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) REAL               array, dimension SIZE
*          The singular values of A, sorted so that S(i) >= S(i+1).
*
*  U       (local output) COMPLEX            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) COMPLEX            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) COMPLEX            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 >= 1 + 2*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(WPCLANGE,WPCGEBRD),
*                       MAX(WPCLARED2D,WP(pre)LARED1D)),
*
*          where WPCLANGE, WPCLARED1D, WPCLARED2D, WPCGEBRD are the
*          workspaces required respectively for the subprograms
*          PCLANGE, PSLARED1D, PSLARED2D, PCGEBRD. 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
*
*          WPCLANGE = MP,
*          WPSLARED1D = NQ0,
*          WPSLARED2D = MP0,
*          WPCGEBRD = 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(WCBDSQR,
*                         MAX(WANTU*WPCORMBRQLN, WANTVT*WPCORMBRPRT)),
*
*          where
*
*                          1, if left(right) singular vectors are wanted
*          WANTU(WANTVT) =
*                          0, otherwise
*
*          and WCBDSQR, WPCORMBRQLN and WPCORMBRPRT refer respectively
*          to the workspace required for the subprograms CBDSQR,
*          PCUNMBR(QLN), and PCUNMBR(PRT), where QLN and PRT are the
*          values of the arguments VECT, SIDE, and TRANS in the call
*          to PCUNMBR. 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 CBDSQR requires
*
*          WCBDSQR = MAX(1, 4*SIZE )
*
*          on every processor. Finally,
*
*          WPCORMBRQLN = MAX( (NB*(NB-1))/2, (SIZEQ+MP)*NB)+NB*NB,
*          WPCORMBRPRT = 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.
*
*  RWORK   (workspace) REAL             array, dimension (1+4*SIZEB)
*          On exit, if INFO = 0, RWORK(1) returns the necessary size
*          for RWORK.
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
 
*          > 0:  if CBDSQR did not converge
*                If INFO = MIN(M,N) + 1, then PCGESVD has detected
*                heterogeneity by finding that eigenvalues were not
*                identical across the process grid. In this case, the
*                accuracy of the results from PCGESVD cannot be
*                guaranteed.
*
*  =====================================================================
*
*  The results of PCGEBRD, and therefore PCGESVD, 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.
*
*  Alignment requirements
*  ======================
*
*  The routine PCGESVD inherits the same alignement requirement as
*  the routine PCGEBRD, namely:
*
*  The distributed submatrix sub( A ) must verify some alignment proper-
*  ties, namely the following expressions should be true:
*  ( MB_A.EQ.NB_A .AND. IROFFA.EQ.ICOFFA )
*          where NB = MB_A = NB_A,
*          IROFFA = MOD( IA-1, NB ), ICOFFA = MOD( JA-1, NB ),
*
*  =====================================================================
*
*
*     .. 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)
      COMPLEX ZERO,ONE
      parameter(zero= ((0.0e+0,0.0e+0)),one= ((1.0e+0,0.0e+0)))
      REAL DZERO,DONE
      parameter(dzero=0.0d+0,done=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,wcbdsqr,wpcgebrd,wpclange,wpcormbrprt,
     +        wpcormbrqln
      REAL ANRM,BIGNUM,EPS,RMAX,RMIN,SAFMIN,SIGMA,SMLNUM
*     ..
*     .. Local Arrays ..
      INTEGER DESCTU(DLEN_),DESCTVT(DLEN_),IDUM1(3),IDUM2(3)
      REAL C(1,1)
*     ..
*     .. External Functions ..
      LOGICAL LSAME
      INTEGER NUMROC
      REAL PSLAMCH,PCLANGE
      EXTERNAL lsame,numroc,pdlamch,pzlange
*     ..
*     .. External Subroutines ..
      EXTERNAL blacs_get,blacs_gridexit,blacs_gridinfo,blacs_gridinit,
     +         chk1mat,cbdsqr,descinit,sgamn2d,sgamx2d,sscal,igamx2d,
     +         igebr2d,igebs2d,pchk1mat,pcgebrd,pcgemr2d,pslared1d,
     +         pslared2d,pclascl,pclaset,pcunmbr,pxerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC max,min,sqrt,real
      INTRINSIC cmplx
*     ..
*     .. 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 = 2
              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(M,1,MYPROWC,0,NPROCS)
              NCVT = NUMROC(N,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)
.EQ.              IF (WANTVT1) THEN
                  SIZEP = NUMROC(SIZE,DESCVT(MB_),MYPROW,DESCVT(RSRC_),
     +                    NPROW)
              ELSE
                  SIZEP = 0
              END IF
.EQ.              IF (WANTU1) THEN
                  SIZEQ = NUMROC(SIZE,DESCU(NB_),MYPCOL,DESCU(CSRC_),
     +                    NPCOL)
              ELSE
                  SIZEQ = 0
              END IF
*
*           Transmit MAX(NQ0, MP0).
*
.EQ..AND..EQ.              IF (MYPROW0  MYPCOL0) 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
*
              WPCLANGE = MP
              WPCGEBRD = NB* (MP+NQ+1) + NQ
              WATOBD = MAX(MAX(WPCLANGE,WPCGEBRD),MAXIM)
*
              WCBDSQR = MAX(1,4*SIZE)
              WPCORMBRQLN = MAX((NB* (NB-1))/2, (SIZEQ+MP)*NB) + NB*NB
              WPCORMBRPRT = MAX((MB* (MB-1))/2, (SIZEP+NQ)*MB) + MB*MB
              WBDTOSVD = SIZE* (WANTU*NRU+WANTVT*NCVT) +
     +                   MAX(WCBDSQR,MAX(WANTU*WPCORMBRQLN,
     +                   WANTVT*WPCORMBRPRT))
*
*           Finally, calculate required workspace.
*
              LWMIN = 1 + 2*SIZEB + MAX(WATOBD,WBDTOSVD)
              WORK(1) = CMPLX(LWMIN,0D+00)
              RWORK(1) = REAL(1+4*SIZEB)
*
.NE..AND..NOT.              IF (WANTU1   (LSAME(JOBU,'n'))) THEN
                  INFO = -1
.NE..AND..NOT.              ELSE IF (WANTVT1   (LSAME(JOBVT,'n'))) THEN
                  INFO = -2
.LT..AND..NE.              ELSE IF (LWORKLWMIN  LWORK-1) THEN
                  INFO = -19
              END IF
*
          END IF
*
          IDUM1(1) = WANTU
          IDUM1(2) = WANTVT
.EQ.          IF (LWORK-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)
.EQ.          IF (INFO0) THEN
.EQ.              IF (WANTU1) THEN
                  CALL PCHK1MAT(M,3,SIZE,4,IU,JU,DESCU,13,0,IDUM1,IDUM2,
     +                          INFO)
              END IF
.EQ.              IF (WANTVT1) THEN
                  CALL PCHK1MAT(SIZE,3,N,4,IVT,JVT,DESCVT,17,0,IDUM1,
     +                          IDUM2,INFO)
              END IF
          END IF
*
      END IF
*
.NE.      IF (INFO0) THEN
          CALL PXERBLA(DESCA(CTXT_),'pcgesvd',-INFO)
          RETURN
.EQ.      ELSE IF (LWORK-1) THEN
          GO TO 40
      END IF
*
*     Quick return if possible.
*
.LE..OR..LE.      IF (M0  N0) GO TO 40
*
*     Get machine constants.
*
      SAFMIN = PSLAMCH(DESCA(CTXT_),'safe minimum')
      EPS = PSLAMCH(DESCA(CTXT_),'precision')
      SMLNUM = SAFMIN/EPS
      BIGNUM = DONE/SMLNUM
      RMIN = SQRT(SMLNUM)
      RMAX = MIN(SQRT(BIGNUM),DONE/SQRT(SQRT(SAFMIN)))
*
*     Scale matrix to allowable range, if necessary.
*
      ANRM = PCLANGE('1',M,N,A,IA,JA,DESCA,WORK(INDWORK))
.GT..AND..LT.      IF (ANRMDZERO  ANRMRMIN) THEN
          ISCALE = 1
          SIGMA = RMIN/ANRM
.GT.      ELSE IF (ANRMRMAX) THEN
          ISCALE = 1
          SIGMA = RMAX/ANRM
      END IF
*
.EQ.      IF (ISCALE1) THEN
          CALL PCLASCL('g',DONE,SIGMA,M,N,A,IA,JA,DESCA,INFO)
      END IF
*
      CALL PCGEBRD(M,N,A,IA,JA,DESCA,RWORK(INDD),RWORK(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.
*
.GE.      IF (MN) THEN
*        Distribute D
          CALL PSLARED1D(N+IOFFD,IA,JA,DESCA,RWORK(INDD),RWORK(INDD2),
     +                   WORK(INDWORK),LLWORK)
*        Distribute E
          CALL PSLARED2D(M+IOFFE,IA,JA,DESCA,RWORK(INDE),RWORK(INDE2),
     +                   WORK(INDWORK),LLWORK)
      ELSE
*        Distribute D
          CALL PSLARED2D(M+IOFFD,IA,JA,DESCA,RWORK(INDD),RWORK(INDD2),
     +                   WORK(INDWORK),LLWORK)
*        Distribute E
          CALL PSLARED1D(N+IOFFE,IA,JA,DESCA,RWORK(INDE),RWORK(INDE2),
     +                   WORK(INDWORK),LLWORK)
      END IF
*
*     Prepare for calling PCBDSQR.
*
.GE.      IF (MN) THEN
          UPLO = 'u'
      ELSE
          UPLO = 'l'
      END IF
*
      INDU = INDWORK
      INDV = INDU + SIZE*NRU*WANTU
      INDWORK = INDV + SIZE*NCVT*WANTVT
*
      LDU = MAX(1,NRU)
      LDVT = MAX(1,SIZE)
*
      CALL DESCINIT(DESCTU,M,SIZE,1,1,0,0,CONTEXTC,LDU,INFO)
      CALL DESCINIT(DESCTVT,SIZE,N,1,1,0,0,CONTEXTR,LDVT,INFO)
*
.EQ.      IF (WANTU1) THEN
          CALL PCLASET('full',M,SIZE,ZERO,ONE,WORK(INDU),1,1,DESCTU)
      ELSE
          NRU = 0
      END IF
*
.EQ.      IF (WANTVT1) THEN
          CALL PCLASET('full',SIZE,N,ZERO,ONE,WORK(INDV),1,1,DESCTVT)
      ELSE
          NCVT = 0
      END IF
*
      CALL CBDSQR(UPLO,SIZE,NCVT,NRU,0,RWORK(INDD2+IOFFD),
     +            RWORK(INDE2+IOFFE),WORK(INDV),SIZE,WORK(INDU),LDU,C,1,
     +            WORK(INDWORK),INFO)
*
*     Redistribute elements of U and VT in the block-cyclic fashion.
*
.EQ.      IF (WANTU1) CALL PCGEMR2D(M,SIZE,WORK(INDU),1,1,DESCTU,U,IU,
     +                              JU,DESCU,DESCU(CTXT_))
*
.EQ.      IF (WANTVT1) CALL PCGEMR2D(SIZE,N,WORK(INDV),1,1,DESCTVT,VT,
     +                               IVT,JVT,DESCVT,DESCVT(CTXT_))
*
*     Set to ZERO "non-square" elements of the larger matrices U, VT.
*
.GT..AND..EQ.      IF (MN  WANTU1) THEN
          CALL PCLASET('full',M-SIZE,SIZE,ZERO,ZERO,U,IA+SIZE,JU,DESCU)
.GT..AND..EQ.      ELSE IF (NM  WANTVT1) THEN
          CALL PCLASET('full',SIZE,N-SIZE,ZERO,ZERO,VT,IVT,JVT+SIZE,
     +                 DESCVT)
      END IF
*
*     Multiply Householder rotations from bidiagonalized matrix.
*
.EQ.      IF (WANTU1) CALL PCUNMBR('q','l','n',M,SIZE,N,A,IA,JA,DESCA,
     +                             WORK(INDTAUQ),U,IU,JU,DESCU,
     +                             WORK(INDWORK),LLWORK,INFO)
*
.EQ.      IF (WANTVT1) CALL PCUNMBR('p','r','c',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) = RWORK(INDD2+IOFFD+I-1)
   10 CONTINUE
*
*     If matrix was scaled, then rescale singular values appropriately.
*
.EQ.      IF (ISCALE1) THEN
          CALL SSCAL(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.
*
.LE.      IF (SIZEITHVAL) THEN
          J = SIZE
          K = 1
      ELSE
          J = SIZE/ITHVAL
          K = ITHVAL
      END IF
*
      DO 20 I = 1,J
          RWORK(I+INDE) = S((I-1)*K+1)
          RWORK(I+INDD2) = S((I-1)*K+1)
   20 CONTINUE
*
      CALL SGAMN2D(DESCA(CTXT_),'a',' ',J,1,RWORK(1+INDE),J,1,1,-1,-1,0)
      CALL SGAMX2D(DESCA(CTXT_),'a',' ',J,1,RWORK(1+INDD2),J,1,1,-1,-1,
     +             0)
*
      DO 30 I = 1,J
.NE.          IF ((RWORK(I+INDE)-RWORK(I+INDD2))DZERO) THEN
              INFO = SIZE + 1
          END IF
   30 CONTINUE
*
   40 CONTINUE
*
      CALL BLACS_GRIDEXIT(CONTEXTC)
      CALL BLACS_GRIDEXIT(CONTEXTR)
*
*     End of PCGESVD
*
      RETURN
      END