pssvdtst.f

Go to the documentation of this file.

      SUBROUTINE pssvdtst( M, N, NPROW, NPCOL, NB, ISEED, THRESH, WORK,
     $                     RESULT, LWORK, NOUT )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997          
*
*     .. Scalar Arguments ..
      INTEGER            LWORK, M, N, NB, NOUT, NPCOL, NPROW
      REAL               THRESH
*     ..
*     .. Array Arguments ..
      INTEGER            ISEED( 4 ), RESULT( 9 )
      REAL               WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  PSSVDTST checks the singular value decomposition (SVD) routine
*  PSGESVD. PSGESVD factors A = U diag(S) VT, where U and VT are
*  orthogonal and diag(S) is diagonal with the entries of the array
*  S on its diagonal. The entries of S are the singular values, stored
*  in decreasing order. U and VT can be optionally not computed,
*  computed and overwritten on A, or computed partially.
*
*  A is M by N. Let SIZE = min( M, N ). S has dimension SIZE by SIZE.
*  U is M by SIZE and VT is SIZE by  N. PDGESVD optionally calculates
*  U and VT, depending on the values of its parameters JOBU and JOBVT.
*  There are four possible  combinations of "job" parameters for a call
*  to PDGESVD, that correspond to four values of internal index JOBTYPE.
*  The table below shows the mapping between "job" parameters of 
*  PDGESVD and respective values of the index JOBTYPE together 
*  with matrices computed for each type of the job.
*
*
*              |   JOBU = 'V'     |       JOBU = 'N'
*   ---------- -------------------------------------------
*   JOBVT = 'V'|  JOBTYPE = 1     |   JOBTYPE = 3
*              |  U1, S1, VT1     |   S3, VT3
*   ---------- ------------------------------------------
*   JOBVT = 'N'|  JOBTYPE = 2     |   JOBTYPE = 4
*              |  U2, S2          |   S4
*
*
*  When PSSVDTST is called, a number of matrix "types" are specified.
*  For each type of matrix, and for the minimal workspace as well as
*  for larger than minimal workspace an M x N matrix "A" with known
*  singular values is generated and used to test the SVD routines.
*  For each matrix, A will be factored as A = U diag(S) VT and the
*  following 9 tests computed:
*
*  (1)   | A - U1 diag(S1) VT1 | / ( |A| max(M,N) ulp )
*
*  (2)   | I - U1'U1 | / ( M ulp )
*
*  (3)   | I - VT1 VT1' | / ( N ulp ),
*
*  (4)   S1 contains SIZE  nonnegative values in decreasing order.
*        (Return 0 if true, 1/ULP if false.)
*
*  (5)   | S1 - S2 | / ( SIZE ulp |S| )
*
*  (6)   | U1 - U2 | / ( M ulp )
*
*  (7)   | S1 - S3 | / ( SIZE ulp |S| )
*
*  (8)   | VT1 - VT3 | / ( N ulp )
*
*  (9)   | S1 - S4 | / (  SIZE ulp |S| )
*
*  Currently, the list of possible  matrix types is:
*
*  (1)  The zero matrix.
*
*  (2)  The identity matrix.
*
*  (3)  A diagonal matrix with evenly spaced entries
*       1, ..., ULP.
*       (ULP = (first number larger than 1) - 1 )
*
*  (4)  A matrix of the form  U D VT, where U, VT are orthogonal and
*       D has evenly spaced entries 1, ..., ULP.
*
*  (5) Same as (4), but multiplied by SQRT( overflow threshold )
*
*  (6) Same as (4), but multiplied by SQRT( underflow threshold )
*
*
*  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 p x q.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the p 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 q 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
*  ==========
*
*  M      (global input) INTEGER dimension
*         The value of the matrix row dimension.
*
*  N      (global input) INTEGER dimension
*         The value of the matrix column dimension.
*
*  NPROW  (global input) INTEGER
*         Number of process rows
*
*  NPCOL  (global input) INTEGER
*         Number of process columns
*
*  NB     (global input) INTEGER
*         The block size of the matrix A. NB >=1.
*
*  ISEED  (global input/local output) INTEGER array, dimension (4)
*         On entry, the seed of the random number generator.  The array
*         elements should be between 0 and 4095; if not they will be
*         reduced mod 4096.  Also, ISEED(4) must be odd.
*         On exit, ISEED is changed and can be used in the next call to
*         SDRVBD to continue the same random number sequence.
*
*  THRESH (global input) REAL
*         The threshold value for the test ratios.  A result is
*         included in the output file if RESULT >= THRESH.  The test
*         ratios are scaled to be O(1), so THRESH should be a small
*         multiple of 1, e.g., 10 or 100.  To have every test ratio
*         printed, use THRESH = 0.
*
*  RESULT (global input/output) INTEGER array of dimension 9. Initially
*         RESULT( I ) = 0. On the output, RESULT ( I ) = 1 if test I
*         ( see above ) wasn't passed.
*
*  WORK   (local workspace) REAL array, dimension (LWORK)
*
*  LWORK  (local input) INTEGER
*         Dimension of the array WORK. It is defined as follows
*         LWORK = 1 + 2*LDA*NQ + 3*SIZE +
*         MAX(WPSLAGGE, LDU*SIZEQ + LDVT*NQ + MAX(LDU*SIZEQ, LDVT*NQ)
*         + WPSGESVD + MAX( WPSSVDCHK, WPSSVDCMP)),
*         where WPSLAGGE, WPSGESVD, WPSSVDCHK,  WPSSVDCMP are amounts
*         of workspace required respectively by PSLAGGE, PSGESVD,
*         PSSVDCHK, PSSVDCMP.
*         Here
*         LDA = NUMROC( M, NB, MYROW, 0, NPROW ), LDU = LDA,
*         LDVT = NUMROC( SIZE, NB, MYROW, 0, NPROW ),
*         NQ = NUMROC( N, NB, MYCOL, 0, NPCOL ),
*         SIZEQ = NUMROC( SIZE, NB, MYCOL, 0, NPCOL ).
*         Values of the variables WPSLAGGE, WPSGESVD, WPSSVDCHK,
*         WPSSVDCMP are found by "dummy" calls to
*         the respective routines. In every "dummy" call, variable
*         LWORK is set to -1, thus causing respective routine
*         immediately return required workspace in WORK(1) without
*         executing any calculations
*
*  NOUT   (local input) INTEGER
*         The unit number for output file.  Only used on node 0.
*         NOUT = 6, output to screen,
*         NOUT = 0, output to stderr.
*  =====================================================================
*
*     .. Parameters ..
      INTEGER             BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
     $                    mb_, nb_, rsrc_, csrc_, lld_, ntypes
      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, ntypes = 6 )
      REAL               ZERO, ONE
      parameter( zero = 0.0e0, one = 1.0e0 )
*     ..
*     .. Local Scalars ..
      CHARACTER          HETERO, JOBU, JOBVT
      INTEGER            CONTEXT, DINFO, I, IA, IAM, INFO, ITYPE, IU,
     $                   ivt, ja, jobtype, ju, jvt, lda, ldu, ldvt,
     $                   llwork, lwmin, mycol, myrow, nnodes, nq, pass,
     $                   ptra, ptrac, ptrd, ptrs, ptrsc, ptru, ptruc, 
     $                   ptrvt, ptrvtc, ptrwork, sethet, SIZE, sizeq,
     $                   wpsgesvd, wpslagge, wpssvdchk, wpssvdcmp
      REAL               CHK, DELTA, H, MTM, OVFL, RTOVFL, RTUNFL, ULP,
     $                   unfl
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_barrier, blacs_get, blacs_gridexit,
     $                   blacs_gridinfo, blacs_gridinit, blacs_pinfo,
     $                   blacs_set,
     $                   descinit, sgamn2d, sgamx2d, slabad, sscal,
     $                   igamn2d, igamx2d, igebr2d, igebs2d, pselset,
     $                   psgesvd, pslacpy, pslagge, pslaset, pssvdchk,
     $                   pssvdcmp, pxerbla, slboot, slcombine, sltimer
*     ..
*     .. External Functions ..
      INTEGER            NUMROC
      REAL               PSLAMCH
      EXTERNAL           numroc, pslamch
*     ..
*     .. Local Arrays ..
      INTEGER            DESCA( DLEN_ ), DESCU( DLEN_ ),
     $                   descvt( dlen_ ), itmp( 2 )
      DOUBLE PRECISION   CTIME( 1 ), WTIME( 1 )
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, int, max, min, sqrt
*     ..
*     .. Executable Statements ..
*     This is just to keep ftnchek happy
      IF( block_cyclic_2d*csrc_*dtype_*lld_*mb_*m_*nb_*n_*rsrc_.LT.0 )
     $   RETURN
*
      CALL blacs_pinfo( iam, nnodes )
      CALL blacs_get( -1, 0, context )
      CALL blacs_gridinit( context, 'R', nprow, npcol )
      CALL blacs_gridinfo( context, nprow, npcol, myrow, mycol )
*
*     If this process is not a part of the contex, bail out now.
*
      IF( ( myrow.GE.nprow ) .OR. ( myrow.LT.0 ) .OR.
     $    ( mycol.GE.npcol ) .OR. ( mycol.LT.0 ) )GO TO 110
      CALL blacs_set( context, 15, 1 )
      info = 0
*
*     Check input parameters.
*
      IF( m.LE.0 ) THEN
         info = -1
      ELSE IF( n.LE.0 ) THEN
         info = -2
      ELSE IF( nprow.LE.0 ) THEN
         info = -3
      ELSE IF( npcol.LE.0 ) THEN
         info = -4
      ELSE IF( nb.LE.0 ) THEN
         info = -5
      ELSE IF( thresh.LE.0 ) THEN
         info = -7
      END IF
*
      SIZE = min( m, n )
*
*     Initialize matrix descriptors.
*
      ia  = 1
      ja  = 1
      iu  = 1
      ju  = 1
      ivt = 1
      jvt = 1
*
      lda = numroc( m, nb, myrow, 0, nprow )
      lda = max( 1, lda )
      nq = numroc( n, nb, mycol, 0, npcol )
      ldu = lda
      sizeq = numroc( SIZE, nb, mycol, 0, npcol )
      ldvt = numroc( SIZE, nb, myrow, 0, nprow )
      ldvt = max( 1, ldvt )
      CALL descinit( desca, m, n, nb, nb, 0, 0, context, lda, dinfo )
      CALL descinit( descu, m, SIZE, nb, nb, 0, 0, context, ldu, dinfo )
      CALL descinit( descvt, SIZE, n, nb, nb, 0, 0, context, ldvt,
     $               dinfo )
*
*     Set some pointers to work array in order to do "dummy" calls.
*
      ptra = 2
      ptrac = ptra + lda*nq
      ptrd = ptrac + lda*nq
      ptrs = ptrd + SIZE
      ptrsc = ptrs + SIZE
      ptrwork = ptrsc + SIZE
*
      ptru = ptrwork
      ptrvt = ptrwork
      ptruc = ptrwork
      ptrvtc = ptrwork
*
*     "Dummy" calls -- return required workspace in work(1) without
*     any calculation. 
*
      CALL pslagge( m, n, work( ptrd ), work( ptra ), ia, ja, desca,
     $              iseed, SIZE, work( ptrwork ), -1, dinfo )
      wpslagge = int( work( ptrwork ) )
*
      CALL psgesvd( 'v', 'v', M, N, WORK( PTRA ), IA, JA, DESCA, 
     $              WORK( PTRS ), WORK( PTRU ), IU, JU, DESCU, 
     $              WORK( PTRVT ), IVT, JVT, DESCVT,
     $              WORK( PTRWORK ), -1, DINFO )
      WPSGESVD = INT( WORK( PTRWORK ) )
*
      CALL PSSVDCHK( M, N, WORK( PTRAC ), IA, JA, DESCA, WORK( PTRUC ),
     $               IU, JU, DESCU, WORK( PTRVT ), IVT, JVT, DESCVT, 
     $               WORK( PTRS ), THRESH, WORK( PTRWORK ), -1, 
     $               RESULT, CHK, MTM )
      WPSSVDCHK = INT( WORK( PTRWORK ) )
*
      CALL PSSVDCMP( M, N, 1, WORK( PTRS ), WORK( PTRSC ), WORK( PTRU ),
     $               WORK( PTRUC ), IU, JU, DESCU, WORK( PTRVT ),
     $               WORK( PTRVTC ), IVT, JVT, DESCVT, THRESH, 
     $               RESULT, DELTA, WORK( PTRWORK ), -1 )
      WPSSVDCMP = INT( WORK( PTRWORK ) )
*
*     Calculation of workspace at last.
*
      LWMIN = 1 + 2*LDA*NQ + 3*SIZE +
     $        MAX( WPSLAGGE, LDU*SIZEQ+LDVT*NQ+MAX( LDU*SIZEQ,
     $        LDVT*NQ )+WPSGESVD+MAX( WPSSVDCHK, WPSSVDCMP ) )
      WORK( 1 ) = LWMIN
*
*     If this is a "dummy" call, return.
*
.EQ.      IF( LWORK-1 )
     $   GO TO 120
.EQ.      IF( INFO0 ) THEN
.LT.         IF( LWORKLWMIN ) THEN
            INFO = -10
         END IF
      END IF
*
.NE.      IF( INFO0 ) THEN
         CALL PXERBLA( DESCA( CTXT_ ), 'pssvdtst', -INFO )
         RETURN
      END IF
*
      ULP = PSLAMCH( CONTEXT, 'p' )
      UNFL = PSLAMCH( CONTEXT, 'safe min' )
      OVFL = ONE / UNFL
      CALL SLABAD( UNFL, OVFL )
      RTUNFL = SQRT( UNFL )
      RTOVFL = SQRT( OVFL )
*
*     This ensures that everyone starts out with the same seed.
*
.EQ..AND..EQ.      IF( MYROW0  MYCOL0 ) THEN
         CALL IGEBS2D( CONTEXT, 'a', ' ', 4, 1, ISEED, 4 )
      ELSE
         CALL IGEBR2D( CONTEXT, 'a', ' ', 4, 1, ISEED, 4, 0, 0 )
      END IF
*
*     Loop over matrix types.
*
      DO 100 ITYPE = 1, NTYPES
*
         PASS = 0
         SETHET = 0
         PTRWORK = PTRSC + SIZE
         LLWORK = LWORK - PTRWORK + 1
*
*        Compute A.
*
.EQ.         IF( ITYPE1 ) THEN
*
*           Zero Matrix.
*
            DO 10 I = 1, SIZE
               WORK( PTRD+I-1 ) = ZERO
   10       CONTINUE
*
            CALL PSLASET( 'all', M, N, ZERO, ZERO, WORK( PTRA ), 
     $                    IA, JA, DESCA )
*
.EQ.         ELSE IF( ITYPE2 ) THEN
*
*           Identity Matrix.
*
            DO 20 I = 1, SIZE
               WORK( PTRD+I-1 ) = ONE
   20       CONTINUE
*
            CALL PSLASET( 'all', M, N, ZERO, ONE, WORK( PTRA ),
     $                    IA, JA, DESCA )
*
.GT.         ELSE IF( ITYPE2 ) THEN
*
*           Preset Singular Values.
*
.NE.            IF( SIZE1 ) THEN
               H = ( ULP-1 ) / ( SIZE-1 )
               DO 30 I = 1, SIZE
                  WORK( PTRD+I-1 ) = 1 + H*( I-1 )
   30          CONTINUE
            ELSE
               WORK( PTRD ) = 1
            END IF
*
.EQ.            IF( ITYPE3 ) THEN
*
*              Diagonal Matrix with specified singular values.
*
               CALL PSLASET( 'all', M, N, ZERO, ZERO, WORK( PTRA ),
     $                       IA, JA, DESCA )
*
               DO 40 I = 1, SIZE
                  CALL PSELSET( WORK( PTRA ), I, I, DESCA,
     $                          WORK( PTRD+I-1 ) )
   40          CONTINUE
*
.EQ.            ELSE IF( ITYPE4 ) THEN
*
*              General matrix with specified singular values.
*
               CALL PSLAGGE( M, N, WORK( PTRD ), WORK( PTRA ), IA, JA,
     $                       DESCA, ISEED, SIZE, WORK( PTRWORK ),
     $                       LLWORK, INFO )
*
.EQ.            ELSE IF( ITYPE5 ) THEN
*
*              Singular values scaled by overflow.
*
               CALL SSCAL( SIZE, RTOVFL, WORK( PTRD ), 1 )
*
               CALL PSLAGGE( M, N, WORK( PTRD ), WORK( PTRA ), IA, JA,
     $                       DESCA, ISEED, SIZE, WORK( PTRWORK ),
     $                       LLWORK, INFO )
*
.EQ.            ELSE IF( ITYPE6 ) THEN
*
*              Singular values scaled by underflow.
*
               CALL SSCAL( SIZE, RTUNFL, WORK( PTRD ), 1 )
               CALL PSLAGGE( M, N, WORK( PTRD ), WORK( PTRA ), IA, JA,
     $                       DESCA, ISEED, SIZE, WORK( PTRWORK ),
     $                       LLWORK, INFO )
*
            END IF
*
         END IF
*
*        Set mapping between JOBTYPE and calling parameters of
*        PSGESVD, reset pointers to WORK array to save space.
*
         DO 80 JOBTYPE = 1, 4
*
.EQ.            IF( JOBTYPE1 ) THEN
               JOBU = 'v'
               JOBVT = 'v'
               PTRVT = PTRU + LDU*SIZEQ
               PTRUC = PTRVT + LDVT*NQ
               PTRWORK = PTRUC + LDU*SIZEQ
               LLWORK = LWORK - PTRWORK + 1
.EQ.            ELSE IF( JOBTYPE2 ) THEN
               JOBU = 'v'
               JOBVT = 'n'
.EQ.            ELSE IF( JOBTYPE3 ) THEN
               JOBU = 'n'
               JOBVT = 'v'
               PTRVTC = PTRUC
               PTRWORK = PTRVTC + LDVT*NQ
               LLWORK = LWORK - PTRWORK + 1
.EQ.            ELSE IF( JOBTYPE4 ) THEN
               JOBU = 'n'
               JOBVT = 'n'
               PTRWORK = PTRUC
               LLWORK = LWORK - PTRWORK + 1
            END IF
*
*           Duplicate matrix A.
*
            CALL PSLACPY( 'a', M, N, WORK( PTRA ), IA, JA, DESCA,
     $                    WORK( PTRAC ), IA, JA, DESCA )
*
*           Test SVD calculation with minimum amount of workspace
*           calculated earlier.
*
.EQ.            IF( JOBTYPE1 ) THEN
*
*              Run SVD.
               CALL SLBOOT
               CALL BLACS_BARRIER( CONTEXT, 'all' )
               CALL SLTIMER( 1 )
*
               CALL PSGESVD( JOBU, JOBVT, M, N, WORK( PTRAC ), IA, JA,
     $                       DESCA, WORK( PTRS ), WORK( PTRU ), IU, JU,
     $                       DESCU, WORK( PTRVT ), IVT, JVT, DESCVT,
     $                       WORK( PTRWORK ), WPSGESVD, INFO )
*
               CALL SLTIMER( 1 )
               CALL SLCOMBINE( CONTEXT, 'all', '>', 'w', 1, 1, WTIME )
               CALL SLCOMBINE( CONTEXT, 'all', '>', 'c', 1, 1, CTIME )
*
*              Check INFO. Different INFO for  different processes mean
*              something went wrong.
*
               ITMP( 1 ) = INFO
               ITMP( 2 ) = INFO
*
               CALL IGAMN2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, ITMP, 1, 1,
     $                       1, -1, -1, 0 )
               CALL IGAMX2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, ITMP( 2 ),
     $                       1, 1, 1, -1, -1, 0 )
*
.NE.               IF( ITMP( 1 )ITMP( 2 ) ) THEN
.EQ..AND..EQ.                  IF( MYROW0  MYCOL0 ) THEN
                     WRITE( NOUT, FMT = * )
     $                  'different processes return different info'
                     GO TO 120
                  END IF
               END IF
*
*              If INFO is  negative PXERBLA tells you. So the only thing
*              is to check for positive INFO -- detected heterogeneous
*              system.
*
.EQ.               IF( INFO( SIZE+1 ) ) THEN
                  HETERO = 'p'
                  SETHET = 1
               END IF
*
*              If INFO was fine do more exhaustive check.
*
.EQ.               IF( INFOZERO ) THEN
*
                  DO 50 I = 1, SIZE
                     WORK( I+PTRWORK ) = WORK( I+PTRS-1 )
                     WORK( I+SIZE+PTRWORK ) = WORK( I+PTRS-1 )
   50             CONTINUE
*
                  CALL SGAMN2D( DESCA( CTXT_ ), 'a', ' ', SIZE, 1,
     $                          WORK( 1+PTRWORK ), SIZE, 1, 1, -1, -1,
     $                          0 )
                  CALL SGAMX2D( DESCA( CTXT_ ), 'a', ' ', SIZE, 1,
     $                          WORK( 1+SIZE+PTRWORK ), SIZE, 1, 1, -1,
     $                          -1, 0 )
*
                  DO 60 I = 1, SIZE
                     IF( ABS( WORK( I+PTRWORK )-WORK( SIZE+I+
.GT.     $                   PTRWORK ) )ZERO ) THEN
                        WRITE( NOUT, FMT = * )'i= ', I, ' min=',
     $                     WORK( I+PTRWORK ), ' max=',
     $                     WORK( SIZE+I+PTRWORK )
                        HETERO = 't'
                        SETHET = 1
                        GO TO 70
                     END IF
*
   60             CONTINUE
   70             CONTINUE
*
               END IF
*
.NE.               IF( SETHET1 )
     $            HETERO = 'n'
*
*              After PSGESVD AC got screwed up -- need to copy again.
*
               CALL PSLACPY( 'a', M, N, WORK( PTRA ), IA, JA, DESCA,
     $                       WORK( PTRAC ), IA, JA, DESCA )
*
*              PSSVDCHK screws up U. So before the call to PSSVDCHK
*              U is copied to UC and a pointer to UC is passed to
*              PSSVDCHK.
*
               CALL PSLACPY( 'a', M, SIZE, WORK( PTRU ), IU, JU, DESCU,
     $                       WORK( PTRUC ), IU, JU, DESCU )
*
*              Run tests 1 - 4.
*
               CALL PSSVDCHK( M, N, WORK( PTRAC ), IA, JA, DESCA,
     $                        WORK( PTRUC ), IU, JU, DESCU, 
     $                        WORK( PTRVT ), IVT, JVT, DESCVT, 
     $                        WORK( PTRS ), THRESH, WORK( PTRWORK ), 
     $                        LLWORK, RESULT, CHK, MTM )
*
            ELSE
*
*              Once again test PSGESVD with min workspace.
*
               CALL PSGESVD( JOBU, JOBVT, M, N, WORK( PTRAC ), IA, JA,
     $                       DESCA, WORK( PTRSC ), WORK( PTRUC ), IU,
     $                       JU, DESCU, WORK( PTRVTC ), IVT, JVT, 
     $                       DESCVT, WORK( PTRWORK ), WPSGESVD, INFO )
*
               CALL PSSVDCMP( M, N, JOBTYPE, WORK( PTRS ),
     $                        WORK( PTRSC ), WORK( PTRU ),
     $                        WORK( PTRUC ), IU, JU, DESCU,
     $                        WORK( PTRVT ), WORK( PTRVTC ), IVT, JVT,
     $                        DESCVT, THRESH, RESULT, DELTA,
     $                        WORK( PTRWORK ), LLWORK )
*
            END IF
*
   80    CONTINUE
*
.EQ..AND..EQ.         IF( MYROW0  MYCOL0 ) THEN
            DO 90 I = 1, 9
.EQ.               IF( RESULT( I )1 ) THEN
                  PASS = 1
                  WRITE( NOUT, FMT = * )'test i = ', I, 'has failed'
                  WRITE( NOUT, FMT = * )' '
               END IF
   90       CONTINUE
.EQ.            IF( PASS0 ) THEN
               WRITE( NOUT, FMT = 9999 )'passed', WTIME( 1 ),
     $            CTIME( 1 ), M, N, NPROW, NPCOL, NB, ITYPE, CHK, MTM,
     $            DELTA, HETERO
            END IF
         END IF
  100 CONTINUE
      CALL BLACS_GRIDEXIT( CONTEXT )
  110 CONTINUE
*
 9999 FORMAT( A6, 2E10.3, 2I6, 2I4, I5, I6, 3F6.2, 4X, A1 )
  120 CONTINUE
*
*     End of PSSVDTST
*
      RETURN
      END