psgesvd.f

Go to the documentation of this file.

 
      SUBROUTINE psgesvd(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) --
*     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(*)
      REAL A(*),U(*),VT(*),WORK(*)
      REAL S(*)
*     ..
*
*  Purpose
*  =======
*
*  PSGESVD 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 REAL
*          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) REAL               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) REAL               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) REAL               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 + 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(WPSLANGE,WPSGEBRD),
*                       MAX(WPSLARED2D,WP(pre)LARED1D)),
*
*          where WPSLANGE, WPSLARED1D, WPSLARED2D, WPSGEBRD are the
*          workspaces required respectively for the subprograms
*          PSLANGE, PSLARED1D, PSLARED2D, PSGEBRD. 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
*
*          WPSLANGE = MP,
*          WPSLARED1D = NQ0,
*          WPSLARED2D = MP0,
*          WPSGEBRD = 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(WSBDSQR,
*                         MAX(WANTU*WPSORMBRQLN, WANTVT*WPSORMBRPRT)),
*
*          where
*
*                          1, if left(right) singular vectors are wanted
*          WANTU(WANTVT) =
*                          0, otherwise
*
*          and WSBDSQR, WPSORMBRQLN and WPSORMBRPRT refer respectively
*          to the workspace required for the subprograms SBDSQR,
*          PSORMBR(QLN), and PSORMBR(PRT), where QLN and PRT are the
*          values of the arguments VECT, SIDE, and TRANS in the call
*          to PSORMBR. 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 SBDSQR requires
*
*          WSBDSQR = MAX(1, 4*SIZE )
*
*          on every processor. Finally,
*
*          WPSORMBRQLN = MAX( (NB*(NB-1))/2, (SIZEQ+MP)*NB)+NB*NB,
*          WPSORMBRPRT = 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    (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 PSGESVD has detected
*                heterogeneity by finding that eigenvalues were not
*                identical across the process grid. In this case, the
*                accuracy of the results from PSGESVD cannot be
*                guaranteed.
*
*  =====================================================================
*
*  The results of PSGEBRD, and therefore PSGESVD, 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 PSGESVD inherits the same alignement requirement as
*  the routine PSGEBRD, 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)
      REAL ZERO,ONE
      parameter(zero= (0.0e+0),one= (1.0e+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,wsbdsqr,wpsgebrd,wpslange,wpsormbrprt,
     +        wpsormbrqln
      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,PSLANGE
      EXTERNAL lsame,numroc,pdlamch,pzlange
*     ..
*     .. External Subroutines ..
      EXTERNAL blacs_get,blacs_gridexit,blacs_gridinfo,blacs_gridinit,
     +         chk1mat,sbdsqr,descinit,sgamn2d,sgamx2d,sscal,igamx2d,
     +         igebr2d,igebs2d,pchk1mat,psgebrd,psgemr2d,pslared1d,
     +         pslared2d,pslascl,pslaset,psormbr,pxerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC max,min,sqrt,real
*     ..
*     .. 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(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)
              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
*
              wpslange = mp
              wpsgebrd = nb* (mp+nq+1) + nq
              watobd = max(max(wpslange,wpsgebrd),maxim)
*
              wsbdsqr = max(1,4*size)
              wpsormbrqln = max((nb* (nb-1))/2, (sizeq+mp)*nb) + nb*nb
              wpsormbrprt = max((mb* (mb-1))/2, (sizep+nq)*mb) + mb*mb
              wbdtosvd = size* (wantu*nru+wantvt*ncvt) +
     +                   max(wsbdsqr,max(wantu*wpsormbrqln,
     +                   wantvt*wpsormbrprt))
*
*           Finally, calculate required workspace.
*
              lwmin = 1 + 6*sizeb + max(watobd,wbdtosvd)
              work(1) = real(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,SIZE,4,iu,ju,descu,13,0,idum1,idum2,
     +                          info)
              END IF
              IF (wantvt.EQ.1) THEN
                  CALL pchk1mat(SIZE,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_),'PSGESVD',-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 = pslamch(desca(ctxt_),'Safe minimum')
      eps = pslamch(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 = pslange('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 pslascl('G',one,sigma,m,n,a,ia,ja,desca,info)
      END IF
*
      CALL psgebrd(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 pslared1d(n+ioffd,ia,ja,desca,work(indd),work(indd2),
     +                   work(indwork),llwork)
*        Distribute E
          CALL pslared2d(m+ioffe,ia,ja,desca,work(inde),work(inde2),
     +                   work(indwork),llwork)
      ELSE
*        Distribute D
          CALL pslared2d(m+ioffd,ia,ja,desca,work(indd),work(indd2),
     +                   work(indwork),llwork)
*        Distribute E
          CALL pslared1d(n+ioffe,ia,ja,desca,work(inde),work(inde2),
     +                   work(indwork),llwork)
      END IF
*
*     Prepare for calling PSBDSQR.
*
      IF (m.GE.n) 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)
*
      IF (wantu.EQ.1) THEN
          CALL pslaset('Full',m,SIZE,zero,one,work(indu),1,1,desctu)
      ELSE
          nru = 0
      END IF
*
      IF (wantvt.EQ.1) THEN
          CALL pslaset('Full',SIZE,n,zero,one,work(indv),1,1,desctvt)
      ELSE
          ncvt = 0
      END IF
*
      CALL sbdsqr(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 psgemr2d(m,SIZE,work(indu),1,1,desctu,u,iu,
     +                              ju,descu,descu(ctxt_))
*
      IF (wantvt.EQ.1) CALL psgemr2d(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.
*
      IF (m.GT.n .AND. wantu.EQ.1) THEN
          CALL pslaset('Full',m-SIZE,SIZE,zero,zero,u,ia+SIZE,ju,descu)
      ELSE IF (n.GT.m .AND. wantvt.EQ.1) THEN
          CALL pslaset('Full',SIZE,n-SIZE,zero,zero,vt,ivt,jvt+SIZE,
     +                 descvt)
      END IF
*
*     Multiply Householder rotations from bidiagonalized matrix.
*
      IF (wantu.EQ.1) CALL psormbr('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 psormbr('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 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.
*
      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 sgamn2d(desca(ctxt_),'a',' ',j,1,work(1+inde),j,1,1,-1,-1,0)
      CALL sgamx2d(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 PSGESVD
*
      RETURN
      END