pdsvdtst.f

Go to the documentation of this file.

      SUBROUTINE pdsvdtst( 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
      DOUBLE PRECISION   THRESH
*     ..
*     .. Array Arguments ..
      INTEGER            ISEED( 4 ), RESULT( 9 )
      DOUBLE PRECISION   WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  PDSVDTST checks the singular value decomposition (SVD) routine
*  PDGESVD. PDGESVD 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 PDSVDTST 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) DOUBLE PRECISION
*         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) DOUBLE PRECISION 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(WPDLAGGE, LDU*SIZEQ + LDVT*NQ + MAX(LDU*SIZEQ, LDVT*NQ)
*         + WPDGESVD + MAX( WPDSVDCHK, WPDSVDCMP)),
*         where WPDLAGGE, WPDGESVD, WPDSVDCHK,  WPDSVDCMP are amounts
*         of workspace required respectively by PDLAGGE, PDGESVD,
*         PDSVDCHK, PDSVDCMP.
*         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 WPDLAGGE, WPDGESVD, WPDSVDCHK,
*         WPDSVDCMP 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 )
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d0, one = 1.0d0 )
*     ..
*     .. 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, ptrwork, ptrs, ptrsc, ptru,
     $                   ptruc, ptrvt, ptrvtc, sethet, SIZE, sizeq,
     $                   wpdgesvd, wpdlagge, wpdsvdchk, wpdsvdcmp
      DOUBLE PRECISION   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, dgamn2d, dgamx2d, dlabad, dscal,
     $                   igamn2d, igamx2d, igebr2d, igebs2d, pdelset,
     $                   pdgesvd, pdlacpy, pdlagge, pdlaset, pdsvdchk,
     $                   pdsvdcmp, pxerbla, slboot, slcombine, sltimer
*     ..
*     .. External Functions ..
      INTEGER            NUMROC
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           numroc, pdlamch
*     ..
*     .. 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 pdlagge( m, n, work( ptrd ), work( ptra ), ia, ja, desca,
     $              iseed, SIZE, work( ptrwork ), -1, dinfo )
      wpdlagge = int( work( ptrwork ) )
*
      CALL pdgesvd( '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 )
      wpdgesvd = int( work( ptrwork ) )
*
      CALL pdsvdchk( 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 )
      wpdsvdchk = int( work( ptrwork ) )
*
      CALL pdsvdcmp( 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 )
      wpdsvdcmp = int( work( ptrwork ) )
*
*     Calculation of workspace at last.
*
      lwmin = 1 + 2*lda*nq + 3*SIZE +
     $        max( wpdlagge, ldu*sizeq+ldvt*nq+max( ldu*sizeq,
     $        ldvt*nq )+wpdgesvd+max( wpdsvdchk, wpdsvdcmp ) )
      work( 1 ) = lwmin
*
*     If this is a "dummy" call, return.
*
      IF( lwork.EQ.-1 )
     $   GO TO 120
      IF( info.EQ.0 ) THEN
         IF( lwork.LT.lwmin ) THEN
            info = -10
         END IF
      END IF
*
      IF( info.NE.0 ) THEN
         CALL pxerbla( desca( ctxt_ ), 'PDSVDTST', -info )
         RETURN
      END IF
*
      ulp = pdlamch( context, 'P' )
      unfl = pdlamch( context, 'Safe min' )
      ovfl = one / unfl
      CALL dlabad( unfl, ovfl )
      rtunfl = sqrt( unfl )
      rtovfl = sqrt( ovfl )
*
*     This ensures that everyone starts out with the same seed.
*
      IF( myrow.EQ.0 .AND. mycol.EQ.0 ) 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.
*
         IF( itype.EQ.1 ) THEN
*
*           Zero Matrix.
*
            DO 10 i = 1, SIZE
               work( ptrd+i-1 ) = zero
   10       CONTINUE
*
            CALL pdlaset( 'All', m, n, zero, zero, work( ptra ), 
     $                    ia, ja, desca )
*
         ELSE IF( itype.EQ.2 ) THEN
*
*           Identity Matrix.
*
            DO 20 i = 1, SIZE
               work( ptrd+i-1 ) = one
   20       CONTINUE
*
            CALL pdlaset( 'All', m, n, zero, one, work( ptra ), 
     $                    ia, ja, desca )
*
         ELSE IF( itype.GT.2 ) THEN
*
*           Preset Singular Values.
*
            IF( size.NE.1 ) 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
*
            IF( itype.EQ.3 ) THEN
*
*              Diagonal Matrix with specified singular values.
*
               CALL pdlaset( 'All', m, n, zero, zero, work( ptra ), 
     $                       ia, ja, desca )
*
               DO 40 i = 1, SIZE
                  CALL pdelset( work( ptra ), i, i, desca,
     $                          work( ptrd+i-1 ) )
   40          CONTINUE
*
            ELSE IF( itype.EQ.4 ) THEN
*
*              General matrix with specified singular values.
*
               CALL pdlagge( m, n, work( ptrd ), work( ptra ), ia, ja,
     $                       desca, iseed, SIZE, work( ptrwork ),
     $                       llwork, info )
*
            ELSE IF( itype.EQ.5 ) THEN
*
*              Singular values scaled by overflow.
*
               CALL dscal( SIZE, rtovfl, work( ptrd ), 1 )
*
               CALL pdlagge( m, n, work( ptrd ), work( ptra ), ia, ja,
     $                       desca, iseed, SIZE, work( ptrwork ),
     $                       llwork, info )
*
            ELSE IF( itype.EQ.6 ) THEN
*
*              Singular values scaled by underflow.
*
               CALL dscal( SIZE, rtunfl, work( ptrd ), 1 )
               CALL pdlagge( 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
*        PDGESVD, reset pointers to WORK array to save space.
*
         DO 80 jobtype = 1, 4
*
            IF( jobtype.EQ.1 ) THEN
               jobu = 'V'
               jobvt = 'V'
               ptrvt = ptru + ldu*sizeq
               ptruc = ptrvt + ldvt*nq
               ptrwork = ptruc + ldu*sizeq
               llwork = lwork - ptrwork + 1
            ELSE IF( jobtype.EQ.2 ) THEN
               jobu = 'V'
               jobvt = 'N'
            ELSE IF( jobtype.EQ.3 ) THEN
               jobu = 'N'
               jobvt = 'V'
               ptrvtc = ptruc
               ptrwork = ptrvtc + ldvt*nq
               llwork = lwork - ptrwork + 1
            ELSE IF( jobtype.EQ.4 ) THEN
               jobu = 'N'
               jobvt = 'N'
               ptrwork = ptruc
               llwork = lwork - ptrwork + 1
            END IF
*
*           Duplicate matrix A.
*
            CALL pdlacpy( 'A', m, n, work( ptra ), ia, ja, desca,
     $                    work( ptrac ), ia, ja, desca )
*
*           Test SVD calculation with minimum amount of workspace
*           calculated earlier.
*
            IF( jobtype.EQ.1 ) THEN
*
*              Run SVD.
*
               CALL slboot
               CALL blacs_barrier( context, 'All' )
               CALL sltimer( 1 )
*
               CALL pdgesvd( jobu, jobvt, m, n, work( ptrac ), ia, ja,
     $                       desca, work( ptrs ), work( ptru ), iu, ju,
     $                       descu, work( ptrvt ), ivt, jvt, descvt,
     $                       work( ptrwork ), wpdgesvd, 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 )
*
               IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
                  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) 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.
*
               IF( info.EQ.( size+1 ) ) THEN
                  hetero = 'P'
                  sethet = 1
               END IF
*
*              If INFO was fine do more exhaustive check.
*
               IF( info.EQ.zero ) THEN
*
                  DO 50 i = 1, SIZE
                     work( i+ptrwork ) = work( i+ptrs-1 )
                     work( i+size+ptrwork ) = work( i+ptrs-1 )
   50             CONTINUE
*
                  CALL dgamn2d( desca( ctxt_ ), 'a', ' ', SIZE, 1,
     $                          work( 1+ptrwork ), SIZE, 1, 1, -1, -1,
     $                          0 )
                  CALL dgamx2d( 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+
     $                   ptrwork ) ).GT.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
*
               IF( sethet.NE.1 )
     $            hetero = 'N'
*
*              Need to copy A again since AC was overwritten by PDGESVD.
*
               CALL pdlacpy( 'A', m, n, work( ptra ), ia, ja, desca,
     $                       work( ptrac ), ia, ja, desca )
*
*              PDSVDCHK overwrites U. So before the call to PDSVDCHK
*              U is copied to UC and a pointer to UC is passed to
*              PDSVDCHK.
*
               CALL pdlacpy( 'A', m, SIZE, work( ptru ), iu, ju, descu,
     $                       work( ptruc ), iu, ju, descu )
*
*              Run tests 1 - 4.
*
               CALL pdsvdchk( 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 PDGESVD with min workspace.
*
               CALL pdgesvd( jobu, jobvt, m, n, work( ptrac ), ia, ja,
     $                       desca, work( ptrsc ), work( ptruc ), iu,
     $                       ju, descu, work( ptrvtc ), ivt, jvt, 
     $                       descvt, work( ptrwork ), wpdgesvd, info )
*
               CALL pdsvdcmp( 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
*
         IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
            DO 90 i = 1, 9
               IF( result( i ).EQ.1 ) THEN
                  pass = 1
                  WRITE( nout, fmt = * )'Test I = ', i, 'has failed'
                  WRITE( nout, fmt = * )' '
               END IF
   90       CONTINUE
            IF( pass.EQ.0 ) 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 PDSVDTST
*
      RETURN
      END