ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdsvdtst.f
Go to the documentation of this file.
1  SUBROUTINE pdsvdtst( M, N, NPROW, NPCOL, NB, ISEED, THRESH, WORK,
2  $ RESULT, LWORK, NOUT )
3 *
4 * -- ScaLAPACK routine (version 1.7) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * May 1, 1997
8 *
9 * .. Scalar Arguments ..
10  INTEGER LWORK, M, N, NB, NOUT, NPCOL, NPROW
11  DOUBLE PRECISION THRESH
12 * ..
13 * .. Array Arguments ..
14  INTEGER ISEED( 4 ), RESULT( 9 )
15  DOUBLE PRECISION WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * PDSVDTST checks the singular value decomposition (SVD) routine
22 * PDGESVD. PDGESVD factors A = U diag(S) VT, where U and VT are
23 * orthogonal and diag(S) is diagonal with the entries of the array
24 * S on its diagonal. The entries of S are the singular values, stored
25 * in decreasing order. U and VT can be optionally not computed,
26 * computed and overwritten on A, or computed partially.
27 *
28 * A is M by N. Let SIZE = min( M, N ). S has dimension SIZE by SIZE.
29 * U is M by SIZE and VT is SIZE by N. PDGESVD optionally calculates
30 * U and VT, depending on the values of its parameters JOBU and JOBVT.
31 * There are four possible combinations of "job" parameters for a call
32 * to PDGESVD, that correspond to four values of internal index JOBTYPE.
33 * The table below shows the mapping between "job" parameters of
34 * PDGESVD and respective values of the index JOBTYPE together
35 * with matrices computed for each type of the job.
36 *
37 *
38 * | JOBU = 'V' | JOBU = 'N'
39 * ---------- -------------------------------------------
40 * JOBVT = 'V'| JOBTYPE = 1 | JOBTYPE = 3
41 * | U1, S1, VT1 | S3, VT3
42 * ---------- ------------------------------------------
43 * JOBVT = 'N'| JOBTYPE = 2 | JOBTYPE = 4
44 * | U2, S2 | S4
45 *
46 *
47 * When PDSVDTST is called, a number of matrix "types" are specified.
48 * For each type of matrix, and for the minimal workspace as well as
49 * for larger than minimal workspace an M x N matrix "A" with known
50 * singular values is generated and used to test the SVD routines.
51 * For each matrix, A will be factored as A = U diag(S) VT and the
52 * following 9 tests computed:
53 *
54 * (1) | A - U1 diag(S1) VT1 | / ( |A| max(M,N) ulp )
55 *
56 * (2) | I - U1'U1 | / ( M ulp )
57 *
58 * (3) | I - VT1 VT1' | / ( N ulp ),
59 *
60 * (4) S1 contains SIZE nonnegative values in decreasing order.
61 * (Return 0 if true, 1/ULP if false.)
62 *
63 * (5) | S1 - S2 | / ( SIZE ulp |S| )
64 *
65 * (6) | U1 - U2 | / ( M ulp )
66 *
67 * (7) | S1 - S3 | / ( SIZE ulp |S| )
68 *
69 * (8) | VT1 - VT3 | / ( N ulp )
70 *
71 * (9) | S1 - S4 | / ( SIZE ulp |S| )
72 *
73 * Currently, the list of possible matrix types is:
74 *
75 * (1) The zero matrix.
76 *
77 * (2) The identity matrix.
78 *
79 * (3) A diagonal matrix with evenly spaced entries
80 * 1, ..., ULP.
81 * (ULP = (first number larger than 1) - 1 )
82 *
83 * (4) A matrix of the form U D VT, where U, VT are orthogonal and
84 * D has evenly spaced entries 1, ..., ULP.
85 *
86 * (5) Same as (4), but multiplied by SQRT( overflow threshold )
87 *
88 * (6) Same as (4), but multiplied by SQRT( underflow threshold )
89 *
90 *
91 * Notes
92 * =====
93 *
94 * Each global data object is described by an associated description
95 * vector. This vector stores the information required to establish
96 * the mapping between an object element and its corresponding process
97 * and memory location.
98 *
99 * Let A be a generic term for any 2D block cyclicly distributed array.
100 * Such a global array has an associated description vector DESCA.
101 * In the following comments, the character _ should be read as
102 * "of the global array".
103 *
104 * NOTATION STORED IN EXPLANATION
105 * --------------- -------------- --------------------------------------
106 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
107 * DTYPE_A = 1.
108 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
109 * the BLACS process grid A is distribu-
110 * ted over. The context itself is glo-
111 * bal, but the handle (the integer
112 * value) may vary.
113 * M_A (global) DESCA( M_ ) The number of rows in the global
114 * array A.
115 * N_A (global) DESCA( N_ ) The number of columns in the global
116 * array A.
117 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
118 * the rows of the array.
119 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
120 * the columns of the array.
121 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
122 * row of the array A is distributed.
123 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
124 * first column of the array A is
125 * distributed.
126 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
127 * array. LLD_A >= MAX(1,LOCr(M_A)).
128 *
129 * Let K be the number of rows or columns of a distributed matrix,
130 * and assume that its process grid has dimension p x q.
131 * LOCr( K ) denotes the number of elements of K that a process
132 * would receive if K were distributed over the p processes of its
133 * process column.
134 * Similarly, LOCc( K ) denotes the number of elements of K that a
135 * process would receive if K were distributed over the q processes of
136 * its process row.
137 * The values of LOCr() and LOCc() may be determined via a call to the
138 * ScaLAPACK tool function, NUMROC:
139 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
140 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
141 * An upper bound for these quantities may be computed by:
142 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
143 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
144 *
145 * Arguments
146 * ==========
147 *
148 * M (global input) INTEGER dimension
149 * The value of the matrix row dimension.
150 *
151 * N (global input) INTEGER dimension
152 * The value of the matrix column dimension.
153 *
154 * NPROW (global input) INTEGER
155 * Number of process rows
156 *
157 * NPCOL (global input) INTEGER
158 * Number of process columns
159 *
160 * NB (global input) INTEGER
161 * The block size of the matrix A. NB >=1.
162 *
163 * ISEED (global input/local output) INTEGER array, dimension (4)
164 * On entry, the seed of the random number generator. The array
165 * elements should be between 0 and 4095; if not they will be
166 * reduced mod 4096. Also, ISEED(4) must be odd.
167 * On exit, ISEED is changed and can be used in the next call to
168 * SDRVBD to continue the same random number sequence.
169 *
170 * THRESH (global input) DOUBLE PRECISION
171 * The threshold value for the test ratios. A result is
172 * included in the output file if RESULT >= THRESH. The test
173 * ratios are scaled to be O(1), so THRESH should be a small
174 * multiple of 1, e.g., 10 or 100. To have every test ratio
175 * printed, use THRESH = 0.
176 *
177 * RESULT (global input/output) INTEGER array of dimension 9. Initially
178 * RESULT( I ) = 0. On the output, RESULT ( I ) = 1 if test I
179 * ( see above ) wasn't passed.
180 *
181 * WORK (local workspace) DOUBLE PRECISION array, dimension (LWORK)
182 *
183 * LWORK (local input) INTEGER
184 * Dimension of the array WORK. It is defined as follows
185 * LWORK = 1 + 2*LDA*NQ + 3*SIZE +
186 * MAX(WPDLAGGE, LDU*SIZEQ + LDVT*NQ + MAX(LDU*SIZEQ, LDVT*NQ)
187 * + WPDGESVD + MAX( WPDSVDCHK, WPDSVDCMP)),
188 * where WPDLAGGE, WPDGESVD, WPDSVDCHK, WPDSVDCMP are amounts
189 * of workspace required respectively by PDLAGGE, PDGESVD,
190 * PDSVDCHK, PDSVDCMP.
191 * Here
192 * LDA = NUMROC( M, NB, MYROW, 0, NPROW ), LDU = LDA,
193 * LDVT = NUMROC( SIZE, NB, MYROW, 0, NPROW ),
194 * NQ = NUMROC( N, NB, MYCOL, 0, NPCOL ),
195 * SIZEQ = NUMROC( SIZE, NB, MYCOL, 0, NPCOL ).
196 * Values of the variables WPDLAGGE, WPDGESVD, WPDSVDCHK,
197 * WPDSVDCMP are found by "dummy" calls to
198 * the respective routines. In every "dummy" call, variable
199 * LWORK is set to -1, thus causing respective routine
200 * immediately return required workspace in WORK(1) without
201 * executing any calculations
202 *
203 * NOUT (local input) INTEGER
204 * The unit number for output file. Only used on node 0.
205 * NOUT = 6, output to screen,
206 * NOUT = 0, output to stderr.
207 * =====================================================================
208 *
209 * .. Parameters ..
210  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
211  $ mb_, nb_, rsrc_, csrc_, lld_, ntypes
212  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
213  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
214  $ rsrc_ = 7, csrc_ = 8, lld_ = 9, ntypes = 6 )
215  DOUBLE PRECISION ZERO, ONE
216  parameter( zero = 0.0d0, one = 1.0d0 )
217 * ..
218 * .. Local Scalars ..
219  CHARACTER HETERO, JOBU, JOBVT
220  INTEGER CONTEXT, DINFO, I, IA, IAM, INFO, ITYPE, IU,
221  $ ivt, ja, jobtype, ju, jvt, lda, ldu, ldvt,
222  $ llwork, lwmin, mycol, myrow, nnodes, nq, pass,
223  $ ptra, ptrac, ptrd, ptrwork, ptrs, ptrsc, ptru,
224  $ ptruc, ptrvt, ptrvtc, sethet, SIZE, sizeq,
225  $ wpdgesvd, wpdlagge, wpdsvdchk, wpdsvdcmp
226  DOUBLE PRECISION CHK, DELTA, H, MTM, OVFL, RTOVFL, RTUNFL, ULP,
227  $ unfl
228 * ..
229 * .. External Subroutines ..
230  EXTERNAL blacs_barrier, blacs_get, blacs_gridexit,
231  $ blacs_gridinfo, blacs_gridinit, blacs_pinfo,
232  $ blacs_set,
233  $ descinit, dgamn2d, dgamx2d, dlabad, dscal,
234  $ igamn2d, igamx2d, igebr2d, igebs2d, pdelset,
235  $ pdgesvd, pdlacpy, pdlagge, pdlaset, pdsvdchk,
236  $ pdsvdcmp, pxerbla, slboot, slcombine, sltimer
237 * ..
238 * .. External Functions ..
239  INTEGER NUMROC
240  DOUBLE PRECISION PDLAMCH
241  EXTERNAL numroc, pdlamch
242 * ..
243 * .. Local Arrays ..
244  INTEGER DESCA( DLEN_ ), DESCU( DLEN_ ),
245  $ descvt( dlen_ ), itmp( 2 )
246  DOUBLE PRECISION CTIME( 1 ), WTIME( 1 )
247 * ..
248 * .. Intrinsic Functions ..
249  INTRINSIC abs, int, max, min, sqrt
250 * ..
251 * .. Executable Statements ..
252 * This is just to keep ftnchek happy
253  IF( block_cyclic_2d*csrc_*dtype_*lld_*mb_*m_*nb_*n_*rsrc_.LT.0 )
254  $ RETURN
255 *
256  CALL blacs_pinfo( iam, nnodes )
257  CALL blacs_get( -1, 0, context )
258  CALL blacs_gridinit( context, 'R', nprow, npcol )
259  CALL blacs_gridinfo( context, nprow, npcol, myrow, mycol )
260 *
261 * If this process is not a part of the contex, bail out now.
262 *
263  IF( ( myrow.GE.nprow ) .OR. ( myrow.LT.0 ) .OR.
264  $ ( mycol.GE.npcol ) .OR. ( mycol.LT.0 ) )GO TO 110
265  CALL blacs_set( context, 15, 1 )
266  info = 0
267 *
268 * Check input parameters.
269 *
270  IF( m.LE.0 ) THEN
271  info = -1
272  ELSE IF( n.LE.0 ) THEN
273  info = -2
274  ELSE IF( nprow.LE.0 ) THEN
275  info = -3
276  ELSE IF( npcol.LE.0 ) THEN
277  info = -4
278  ELSE IF( nb.LE.0 ) THEN
279  info = -5
280  ELSE IF( thresh.LE.0 ) THEN
281  info = -7
282  END IF
283 *
284  SIZE = min( m, n )
285 *
286 * Initialize matrix descriptors.
287 *
288  ia = 1
289  ja = 1
290  iu = 1
291  ju = 1
292  ivt = 1
293  jvt = 1
294 *
295  lda = numroc( m, nb, myrow, 0, nprow )
296  lda = max( 1, lda )
297  nq = numroc( n, nb, mycol, 0, npcol )
298  ldu = lda
299  sizeq = numroc( SIZE, nb, mycol, 0, npcol )
300  ldvt = numroc( SIZE, nb, myrow, 0, nprow )
301  ldvt = max( 1, ldvt )
302  CALL descinit( desca, m, n, nb, nb, 0, 0, context, lda, dinfo )
303  CALL descinit( descu, m, SIZE, nb, nb, 0, 0, context, ldu, dinfo )
304  CALL descinit( descvt, SIZE, n, nb, nb, 0, 0, context, ldvt,
305  $ dinfo )
306 *
307 * Set some pointers to work array in order to do "dummy" calls.
308 *
309  ptra = 2
310  ptrac = ptra + lda*nq
311  ptrd = ptrac + lda*nq
312  ptrs = ptrd + SIZE
313  ptrsc = ptrs + SIZE
314  ptrwork = ptrsc + SIZE
315 *
316  ptru = ptrwork
317  ptrvt = ptrwork
318  ptruc = ptrwork
319  ptrvtc = ptrwork
320 *
321 * "Dummy" calls -- return required workspace in work(1) without
322 * any calculation.
323 *
324  CALL pdlagge( m, n, work( ptrd ), work( ptra ), ia, ja, desca,
325  $ iseed, SIZE, work( ptrwork ), -1, dinfo )
326  wpdlagge = int( work( ptrwork ) )
327 *
328  CALL pdgesvd( 'V', 'V', m, n, work( ptra ), ia, ja, desca,
329  $ work( ptrs ), work( ptru ), iu, ju, descu,
330  $ work( ptrvt ), ivt, jvt, descvt,
331  $ work( ptrwork ), -1, dinfo )
332  wpdgesvd = int( work( ptrwork ) )
333 *
334  CALL pdsvdchk( m, n, work( ptrac ), ia, ja, desca, work( ptruc ),
335  $ iu, ju, descu, work( ptrvt ), ivt, jvt, descvt,
336  $ work( ptrs ), thresh, work( ptrwork ), -1,
337  $ result, chk, mtm )
338  wpdsvdchk = int( work( ptrwork ) )
339 *
340  CALL pdsvdcmp( m, n, 1, work( ptrs ), work( ptrsc ), work( ptru ),
341  $ work( ptruc ), iu, ju, descu, work( ptrvt ),
342  $ work( ptrvtc ), ivt, jvt, descvt, thresh,
343  $ result, delta, work( ptrwork ), -1 )
344  wpdsvdcmp = int( work( ptrwork ) )
345 *
346 * Calculation of workspace at last.
347 *
348  lwmin = 1 + 2*lda*nq + 3*SIZE +
349  $ max( wpdlagge, ldu*sizeq+ldvt*nq+max( ldu*sizeq,
350  $ ldvt*nq )+wpdgesvd+max( wpdsvdchk, wpdsvdcmp ) )
351  work( 1 ) = lwmin
352 *
353 * If this is a "dummy" call, return.
354 *
355  IF( lwork.EQ.-1 )
356  $ GO TO 120
357  IF( info.EQ.0 ) THEN
358  IF( lwork.LT.lwmin ) THEN
359  info = -10
360  END IF
361  END IF
362 *
363  IF( info.NE.0 ) THEN
364  CALL pxerbla( desca( ctxt_ ), 'PDSVDTST', -info )
365  RETURN
366  END IF
367 *
368  ulp = pdlamch( context, 'P' )
369  unfl = pdlamch( context, 'Safe min' )
370  ovfl = one / unfl
371  CALL dlabad( unfl, ovfl )
372  rtunfl = sqrt( unfl )
373  rtovfl = sqrt( ovfl )
374 *
375 * This ensures that everyone starts out with the same seed.
376 *
377  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
378  CALL igebs2d( context, 'a', ' ', 4, 1, iseed, 4 )
379  ELSE
380  CALL igebr2d( context, 'a', ' ', 4, 1, iseed, 4, 0, 0 )
381  END IF
382 *
383 * Loop over matrix types.
384 *
385  DO 100 itype = 1, ntypes
386 *
387  pass = 0
388  sethet = 0
389  ptrwork = ptrsc + SIZE
390  llwork = lwork - ptrwork + 1
391 *
392 * Compute A.
393 *
394  IF( itype.EQ.1 ) THEN
395 *
396 * Zero Matrix.
397 *
398  DO 10 i = 1, SIZE
399  work( ptrd+i-1 ) = zero
400  10 CONTINUE
401 *
402  CALL pdlaset( 'All', m, n, zero, zero, work( ptra ),
403  $ ia, ja, desca )
404 *
405  ELSE IF( itype.EQ.2 ) THEN
406 *
407 * Identity Matrix.
408 *
409  DO 20 i = 1, SIZE
410  work( ptrd+i-1 ) = one
411  20 CONTINUE
412 *
413  CALL pdlaset( 'All', m, n, zero, one, work( ptra ),
414  $ ia, ja, desca )
415 *
416  ELSE IF( itype.GT.2 ) THEN
417 *
418 * Preset Singular Values.
419 *
420  IF( size.NE.1 ) THEN
421  h = ( ulp-1 ) / ( size-1 )
422  DO 30 i = 1, SIZE
423  work( ptrd+i-1 ) = 1 + h*( i-1 )
424  30 CONTINUE
425  ELSE
426  work( ptrd ) = 1
427  END IF
428 *
429  IF( itype.EQ.3 ) THEN
430 *
431 * Diagonal Matrix with specified singular values.
432 *
433  CALL pdlaset( 'All', m, n, zero, zero, work( ptra ),
434  $ ia, ja, desca )
435 *
436  DO 40 i = 1, SIZE
437  CALL pdelset( work( ptra ), i, i, desca,
438  $ work( ptrd+i-1 ) )
439  40 CONTINUE
440 *
441  ELSE IF( itype.EQ.4 ) THEN
442 *
443 * General matrix with specified singular values.
444 *
445  CALL pdlagge( m, n, work( ptrd ), work( ptra ), ia, ja,
446  $ desca, iseed, SIZE, work( ptrwork ),
447  $ llwork, info )
448 *
449  ELSE IF( itype.EQ.5 ) THEN
450 *
451 * Singular values scaled by overflow.
452 *
453  CALL dscal( SIZE, rtovfl, work( ptrd ), 1 )
454 *
455  CALL pdlagge( m, n, work( ptrd ), work( ptra ), ia, ja,
456  $ desca, iseed, SIZE, work( ptrwork ),
457  $ llwork, info )
458 *
459  ELSE IF( itype.EQ.6 ) THEN
460 *
461 * Singular values scaled by underflow.
462 *
463  CALL dscal( SIZE, rtunfl, work( ptrd ), 1 )
464  CALL pdlagge( m, n, work( ptrd ), work( ptra ), ia, ja,
465  $ desca, iseed, SIZE, work( ptrwork ),
466  $ llwork, info )
467 *
468  END IF
469 *
470  END IF
471 *
472 * Set mapping between JOBTYPE and calling parameters of
473 * PDGESVD, reset pointers to WORK array to save space.
474 *
475  DO 80 jobtype = 1, 4
476 *
477  IF( jobtype.EQ.1 ) THEN
478  jobu = 'V'
479  jobvt = 'V'
480  ptrvt = ptru + ldu*sizeq
481  ptruc = ptrvt + ldvt*nq
482  ptrwork = ptruc + ldu*sizeq
483  llwork = lwork - ptrwork + 1
484  ELSE IF( jobtype.EQ.2 ) THEN
485  jobu = 'V'
486  jobvt = 'N'
487  ELSE IF( jobtype.EQ.3 ) THEN
488  jobu = 'N'
489  jobvt = 'V'
490  ptrvtc = ptruc
491  ptrwork = ptrvtc + ldvt*nq
492  llwork = lwork - ptrwork + 1
493  ELSE IF( jobtype.EQ.4 ) THEN
494  jobu = 'N'
495  jobvt = 'N'
496  ptrwork = ptruc
497  llwork = lwork - ptrwork + 1
498  END IF
499 *
500 * Duplicate matrix A.
501 *
502  CALL pdlacpy( 'A', m, n, work( ptra ), ia, ja, desca,
503  $ work( ptrac ), ia, ja, desca )
504 *
505 * Test SVD calculation with minimum amount of workspace
506 * calculated earlier.
507 *
508  IF( jobtype.EQ.1 ) THEN
509 *
510 * Run SVD.
511 *
512  CALL slboot
513  CALL blacs_barrier( context, 'All' )
514  CALL sltimer( 1 )
515 *
516  CALL pdgesvd( jobu, jobvt, m, n, work( ptrac ), ia, ja,
517  $ desca, work( ptrs ), work( ptru ), iu, ju,
518  $ descu, work( ptrvt ), ivt, jvt, descvt,
519  $ work( ptrwork ), wpdgesvd, info )
520 *
521  CALL sltimer( 1 )
522  CALL slcombine( context, 'All', '>', 'W', 1, 1, wtime )
523  CALL slcombine( context, 'All', '>', 'C', 1, 1, ctime )
524 *
525 * Check INFO. Different INFO for different processes mean
526 * something went wrong.
527 *
528  itmp( 1 ) = info
529  itmp( 2 ) = info
530 *
531  CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1,
532  $ 1, -1, -1, 0 )
533  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ),
534  $ 1, 1, 1, -1, -1, 0 )
535 *
536  IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
537  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
538  WRITE( nout, fmt = * )
539  $ 'Different processes return different INFO'
540  GO TO 120
541  END IF
542  END IF
543 *
544 * If INFO is negative PXERBLA tells you. So the only thing
545 * is to check for positive INFO -- detected heterogeneous
546 * system.
547 *
548  IF( info.EQ.( size+1 ) ) THEN
549  hetero = 'P'
550  sethet = 1
551  END IF
552 *
553 * If INFO was fine do more exhaustive check.
554 *
555  IF( info.EQ.zero ) THEN
556 *
557  DO 50 i = 1, SIZE
558  work( i+ptrwork ) = work( i+ptrs-1 )
559  work( i+size+ptrwork ) = work( i+ptrs-1 )
560  50 CONTINUE
561 *
562  CALL dgamn2d( desca( ctxt_ ), 'a', ' ', SIZE, 1,
563  $ work( 1+ptrwork ), SIZE, 1, 1, -1, -1,
564  $ 0 )
565  CALL dgamx2d( desca( ctxt_ ), 'a', ' ', SIZE, 1,
566  $ work( 1+size+ptrwork ), SIZE, 1, 1, -1,
567  $ -1, 0 )
568 *
569  DO 60 i = 1, SIZE
570  IF( abs( work( i+ptrwork )-work( size+i+
571  $ ptrwork ) ).GT.zero ) THEN
572  WRITE( nout, fmt = * )'I= ', i, ' MIN=',
573  $ work( i+ptrwork ), ' MAX=',
574  $ work( size+i+ptrwork )
575  hetero = 'T'
576  sethet = 1
577  GO TO 70
578  END IF
579 *
580  60 CONTINUE
581  70 CONTINUE
582 *
583  END IF
584 *
585  IF( sethet.NE.1 )
586  $ hetero = 'N'
587 *
588 * Need to copy A again since AC was overwritten by PDGESVD.
589 *
590  CALL pdlacpy( 'A', m, n, work( ptra ), ia, ja, desca,
591  $ work( ptrac ), ia, ja, desca )
592 *
593 * PDSVDCHK overwrites U. So before the call to PDSVDCHK
594 * U is copied to UC and a pointer to UC is passed to
595 * PDSVDCHK.
596 *
597  CALL pdlacpy( 'A', m, SIZE, work( ptru ), iu, ju, descu,
598  $ work( ptruc ), iu, ju, descu )
599 *
600 * Run tests 1 - 4.
601 *
602  CALL pdsvdchk( m, n, work( ptrac ), ia, ja, desca,
603  $ work( ptruc ), iu, ju, descu,
604  $ work( ptrvt ), ivt, jvt, descvt,
605  $ work( ptrs ), thresh, work( ptrwork ),
606  $ llwork, result, chk, mtm )
607 *
608  ELSE
609 *
610 * Once again test PDGESVD with min workspace.
611 *
612  CALL pdgesvd( jobu, jobvt, m, n, work( ptrac ), ia, ja,
613  $ desca, work( ptrsc ), work( ptruc ), iu,
614  $ ju, descu, work( ptrvtc ), ivt, jvt,
615  $ descvt, work( ptrwork ), wpdgesvd, info )
616 *
617  CALL pdsvdcmp( m, n, jobtype, work( ptrs ),
618  $ work( ptrsc ), work( ptru ),
619  $ work( ptruc ), iu, ju, descu,
620  $ work( ptrvt ), work( ptrvtc ), ivt, jvt,
621  $ descvt, thresh, result, delta,
622  $ work( ptrwork ), llwork )
623 *
624  END IF
625 *
626  80 CONTINUE
627 *
628  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
629  DO 90 i = 1, 9
630  IF( result( i ).EQ.1 ) THEN
631  pass = 1
632  WRITE( nout, fmt = * )'Test I = ', i, 'has failed'
633  WRITE( nout, fmt = * )' '
634  END IF
635  90 CONTINUE
636  IF( pass.EQ.0 ) THEN
637  WRITE( nout, fmt = 9999 )'Passed', wtime( 1 ),
638  $ ctime( 1 ), m, n, nprow, npcol, nb, itype, chk, mtm,
639  $ delta, hetero
640  END IF
641  END IF
642  100 CONTINUE
643  CALL blacs_gridexit( context )
644  110 CONTINUE
645 *
646  9999 FORMAT( a6, 2e10.3, 2i6, 2i4, i5, i6, 3f6.2, 4x, a1 )
647  120 CONTINUE
648 *
649 * End of PDSVDTST
650 *
651  RETURN
652  END
pdlagge
subroutine pdlagge(M, N, D, A, IA, JA, DESCA, ISEED, ORDER, WORK, LWORK, INFO)
Definition: pdlagge.f:3
pdsvdchk
subroutine pdsvdchk(M, N, A, IA, JA, DESCA, U, IU, JU, DESCU, VT, IVT, JVT, DESCVT, S, THRESH, WORK, LWORK, RESULT, CHK, MTM)
Definition: pdsvdchk.f:4
max
#define max(A, B)
Definition: pcgemr.c:180
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
pdgesvd
subroutine pdgesvd(JOBU, JOBVT, M, N, A, IA, JA, DESCA, S, U, IU, JU, DESCU, VT, IVT, JVT, DESCVT, WORK, LWORK, INFO)
Definition: pdgesvd.f:4
pdsvdcmp
subroutine pdsvdcmp(M, N, JOBTYPE, S, SC, U, UC, IU, JU, DESCU, VT, VTC, IVT, JVT, DESCVT, THRESH, RESULT, DELTA, WORK, LWORK)
Definition: pdsvdcmp.f:4
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
slboot
subroutine slboot()
Definition: sltimer.f:2
pdlaset
subroutine pdlaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pdblastst.f:6862
pdsvdtst
subroutine pdsvdtst(M, N, NPROW, NPCOL, NB, ISEED, THRESH, WORK, RESULT, LWORK, NOUT)
Definition: pdsvdtst.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
pdlacpy
subroutine pdlacpy(UPLO, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB)
Definition: pdlacpy.f:3
slcombine
subroutine slcombine(ICTXT, SCOPE, OP, TIMETYPE, N, IBEG, TIMES)
Definition: sltimer.f:267
pdelset
subroutine pdelset(A, IA, JA, DESCA, ALPHA)
Definition: pdelset.f:2
min
#define min(A, B)
Definition: pcgemr.c:181