ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pssvdtst.f
Go to the documentation of this file.
1  SUBROUTINE pssvdtst( 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  REAL THRESH
12 * ..
13 * .. Array Arguments ..
14  INTEGER ISEED( 4 ), RESULT( 9 )
15  REAL WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * PSSVDTST checks the singular value decomposition (SVD) routine
22 * PSGESVD. PSGESVD 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 PSSVDTST 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) REAL
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) REAL 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(WPSLAGGE, LDU*SIZEQ + LDVT*NQ + MAX(LDU*SIZEQ, LDVT*NQ)
187 * + WPSGESVD + MAX( WPSSVDCHK, WPSSVDCMP)),
188 * where WPSLAGGE, WPSGESVD, WPSSVDCHK, WPSSVDCMP are amounts
189 * of workspace required respectively by PSLAGGE, PSGESVD,
190 * PSSVDCHK, PSSVDCMP.
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 WPSLAGGE, WPSGESVD, WPSSVDCHK,
197 * WPSSVDCMP 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  REAL ZERO, ONE
216  parameter( zero = 0.0e0, one = 1.0e0 )
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, ptrs, ptrsc, ptru, ptruc,
224  $ ptrvt, ptrvtc, ptrwork, sethet, SIZE, sizeq,
225  $ wpsgesvd, wpslagge, wpssvdchk, wpssvdcmp
226  REAL 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, sgamn2d, sgamx2d, slabad, sscal,
234  $ igamn2d, igamx2d, igebr2d, igebs2d, pselset,
235  $ psgesvd, pslacpy, pslagge, pslaset, pssvdchk,
236  $ pssvdcmp, pxerbla, slboot, slcombine, sltimer
237 * ..
238 * .. External Functions ..
239  INTEGER NUMROC
240  REAL PSLAMCH
241  EXTERNAL numroc, pslamch
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 pslagge( m, n, work( ptrd ), work( ptra ), ia, ja, desca,
325  $ iseed, SIZE, work( ptrwork ), -1, dinfo )
326  wpslagge = int( work( ptrwork ) )
327 *
328  CALL psgesvd( '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  wpsgesvd = int( work( ptrwork ) )
333 *
334  CALL pssvdchk( 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  wpssvdchk = int( work( ptrwork ) )
339 *
340  CALL pssvdcmp( 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  wpssvdcmp = int( work( ptrwork ) )
345 *
346 * Calculation of workspace at last.
347 *
348  lwmin = 1 + 2*lda*nq + 3*SIZE +
349  $ max( wpslagge, ldu*sizeq+ldvt*nq+max( ldu*sizeq,
350  $ ldvt*nq )+wpsgesvd+max( wpssvdchk, wpssvdcmp ) )
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_ ), 'PSSVDTST', -info )
365  RETURN
366  END IF
367 *
368  ulp = pslamch( context, 'P' )
369  unfl = pslamch( context, 'Safe min' )
370  ovfl = one / unfl
371  CALL slabad( 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 pslaset( '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 pslaset( '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 pslaset( 'All', m, n, zero, zero, work( ptra ),
434  $ ia, ja, desca )
435 *
436  DO 40 i = 1, SIZE
437  CALL pselset( 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 pslagge( 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 sscal( SIZE, rtovfl, work( ptrd ), 1 )
454 *
455  CALL pslagge( 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 sscal( SIZE, rtunfl, work( ptrd ), 1 )
464  CALL pslagge( 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 * PSGESVD, 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 pslacpy( '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  CALL slboot
512  CALL blacs_barrier( context, 'All' )
513  CALL sltimer( 1 )
514 *
515  CALL psgesvd( jobu, jobvt, m, n, work( ptrac ), ia, ja,
516  $ desca, work( ptrs ), work( ptru ), iu, ju,
517  $ descu, work( ptrvt ), ivt, jvt, descvt,
518  $ work( ptrwork ), wpsgesvd, info )
519 *
520  CALL sltimer( 1 )
521  CALL slcombine( context, 'All', '>', 'W', 1, 1, wtime )
522  CALL slcombine( context, 'All', '>', 'C', 1, 1, ctime )
523 *
524 * Check INFO. Different INFO for different processes mean
525 * something went wrong.
526 *
527  itmp( 1 ) = info
528  itmp( 2 ) = info
529 *
530  CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1,
531  $ 1, -1, -1, 0 )
532  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ),
533  $ 1, 1, 1, -1, -1, 0 )
534 *
535  IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
536  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
537  WRITE( nout, fmt = * )
538  $ 'Different processes return different INFO'
539  GO TO 120
540  END IF
541  END IF
542 *
543 * If INFO is negative PXERBLA tells you. So the only thing
544 * is to check for positive INFO -- detected heterogeneous
545 * system.
546 *
547  IF( info.EQ.( size+1 ) ) THEN
548  hetero = 'P'
549  sethet = 1
550  END IF
551 *
552 * If INFO was fine do more exhaustive check.
553 *
554  IF( info.EQ.zero ) THEN
555 *
556  DO 50 i = 1, SIZE
557  work( i+ptrwork ) = work( i+ptrs-1 )
558  work( i+size+ptrwork ) = work( i+ptrs-1 )
559  50 CONTINUE
560 *
561  CALL sgamn2d( desca( ctxt_ ), 'a', ' ', SIZE, 1,
562  $ work( 1+ptrwork ), SIZE, 1, 1, -1, -1,
563  $ 0 )
564  CALL sgamx2d( desca( ctxt_ ), 'a', ' ', SIZE, 1,
565  $ work( 1+size+ptrwork ), SIZE, 1, 1, -1,
566  $ -1, 0 )
567 *
568  DO 60 i = 1, SIZE
569  IF( abs( work( i+ptrwork )-work( size+i+
570  $ ptrwork ) ).GT.zero ) THEN
571  WRITE( nout, fmt = * )'I= ', i, ' MIN=',
572  $ work( i+ptrwork ), ' MAX=',
573  $ work( size+i+ptrwork )
574  hetero = 'T'
575  sethet = 1
576  GO TO 70
577  END IF
578 *
579  60 CONTINUE
580  70 CONTINUE
581 *
582  END IF
583 *
584  IF( sethet.NE.1 )
585  $ hetero = 'N'
586 *
587 * After PSGESVD AC got screwed up -- need to copy again.
588 *
589  CALL pslacpy( 'A', m, n, work( ptra ), ia, ja, desca,
590  $ work( ptrac ), ia, ja, desca )
591 *
592 * PSSVDCHK screws up U. So before the call to PSSVDCHK
593 * U is copied to UC and a pointer to UC is passed to
594 * PSSVDCHK.
595 *
596  CALL pslacpy( 'A', m, SIZE, work( ptru ), iu, ju, descu,
597  $ work( ptruc ), iu, ju, descu )
598 *
599 * Run tests 1 - 4.
600 *
601  CALL pssvdchk( m, n, work( ptrac ), ia, ja, desca,
602  $ work( ptruc ), iu, ju, descu,
603  $ work( ptrvt ), ivt, jvt, descvt,
604  $ work( ptrs ), thresh, work( ptrwork ),
605  $ llwork, result, chk, mtm )
606 *
607  ELSE
608 *
609 * Once again test PSGESVD with min workspace.
610 *
611  CALL psgesvd( jobu, jobvt, m, n, work( ptrac ), ia, ja,
612  $ desca, work( ptrsc ), work( ptruc ), iu,
613  $ ju, descu, work( ptrvtc ), ivt, jvt,
614  $ descvt, work( ptrwork ), wpsgesvd, info )
615 *
616  CALL pssvdcmp( m, n, jobtype, work( ptrs ),
617  $ work( ptrsc ), work( ptru ),
618  $ work( ptruc ), iu, ju, descu,
619  $ work( ptrvt ), work( ptrvtc ), ivt, jvt,
620  $ descvt, thresh, result, delta,
621  $ work( ptrwork ), llwork )
622 *
623  END IF
624 *
625  80 CONTINUE
626 *
627  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
628  DO 90 i = 1, 9
629  IF( result( i ).EQ.1 ) THEN
630  pass = 1
631  WRITE( nout, fmt = * )'Test I = ', i, 'has failed'
632  WRITE( nout, fmt = * )' '
633  END IF
634  90 CONTINUE
635  IF( pass.EQ.0 ) THEN
636  WRITE( nout, fmt = 9999 )'Passed', wtime( 1 ),
637  $ ctime( 1 ), m, n, nprow, npcol, nb, itype, chk, mtm,
638  $ delta, hetero
639  END IF
640  END IF
641  100 CONTINUE
642  CALL blacs_gridexit( context )
643  110 CONTINUE
644 *
645  9999 FORMAT( a6, 2e10.3, 2i6, 2i4, i5, i6, 3f6.2, 4x, a1 )
646  120 CONTINUE
647 *
648 * End of PSSVDTST
649 *
650  RETURN
651  END
max
#define max(A, B)
Definition: pcgemr.c:180
pssvdtst
subroutine pssvdtst(M, N, NPROW, NPCOL, NB, ISEED, THRESH, WORK, RESULT, LWORK, NOUT)
Definition: pssvdtst.f:3
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
pssvdcmp
subroutine pssvdcmp(M, N, JOBTYPE, S, SC, U, UC, IU, JU, DESCU, VT, VTC, IVT, JVT, DESCVT, THRESH, RESULT, DELTA, WORK, LWORK)
Definition: pssvdcmp.f:4
pselset
subroutine pselset(A, IA, JA, DESCA, ALPHA)
Definition: pselset.f:2
psgesvd
subroutine psgesvd(JOBU, JOBVT, M, N, A, IA, JA, DESCA, S, U, IU, JU, DESCU, VT, IVT, JVT, DESCVT, WORK, LWORK, INFO)
Definition: psgesvd.f:4
pslacpy
subroutine pslacpy(UPLO, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB)
Definition: pslacpy.f:3
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
slboot
subroutine slboot()
Definition: sltimer.f:2
pssvdchk
subroutine pssvdchk(M, N, A, IA, JA, DESCA, U, IU, JU, DESCU, VT, IVT, JVT, DESCVT, S, THRESH, WORK, LWORK, RESULT, CHK, MTM)
Definition: pssvdchk.f:4
pslagge
subroutine pslagge(M, N, D, A, IA, JA, DESCA, ISEED, ORDER, WORK, LWORK, INFO)
Definition: pslagge.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
pslaset
subroutine pslaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: psblastst.f:6863
slcombine
subroutine slcombine(ICTXT, SCOPE, OP, TIMETYPE, N, IBEG, TIMES)
Definition: sltimer.f:267
min
#define min(A, B)
Definition: pcgemr.c:181