ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdsepsubtst.f
Go to the documentation of this file.
1  SUBROUTINE pdsepsubtst( WKNOWN, JOBZ, RANGE, UPLO, N, VL, VU, IL,
2  $ IU, THRESH, ABSTOL, A, COPYA, Z, IA, JA,
3  $ DESCA, WIN, WNEW, IFAIL, ICLUSTR, GAP,
4  $ IPREPAD, IPOSTPAD, WORK, LWORK, LWORK1,
5  $ IWORK, LIWORK, RESULT, TSTNRM, QTQNRM,
6  $ NOUT )
7 *
8 * -- ScaLAPACK routine (version 1.7) --
9 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
10 * and University of California, Berkeley.
11 * March 15, 2000
12 *
13 * .. Scalar Arguments ..
14  LOGICAL WKNOWN
15  CHARACTER JOBZ, RANGE, UPLO
16  INTEGER IA, IL, IPOSTPAD, IPREPAD, IU, JA, LIWORK,
17  $ LWORK, LWORK1, N, NOUT, RESULT
18  DOUBLE PRECISION ABSTOL, QTQNRM, THRESH, TSTNRM, VL, VU
19 * ..
20 * .. Array Arguments ..
21  INTEGER DESCA( * ), ICLUSTR( * ), IFAIL( * ),
22  $ IWORK( * )
23  DOUBLE PRECISION A( * ), COPYA( * ), GAP( * ), WIN( * ),
24  $ WNEW( * ), WORK( * ), Z( * )
25 * ..
26 *
27 * Purpose
28 * =======
29 *
30 * PDSEPSUBTST calls PDSYEVX and then tests the output of
31 * PDSYEVX
32 * If JOBZ = 'V' then the following two tests are performed:
33 * |AQ -QL| / (abstol + eps * norm(A) ) < N*THRESH
34 * |QT * Q - I| / eps * norm(A) < N*THRESH
35 * If WKNOWN then
36 * we check to make sure that the eigenvalues match expectations
37 * i.e. |WIN - WNEW(1+IPREPAD)| / (eps * |WIN|) < THRESH
38 * where WIN is the array of eigenvalues as computed by
39 * PDSYEVX when eigenvectors are requested
40 *
41 * Arguments
42 * =========
43 *
44 * NP = the number of rows local to a given process.
45 * NQ = the number of columns local to a given process.
46 *
47 * WKNOWN (global input) INTEGER
48 * .FALSE.: WIN does not contain the eigenvalues
49 * .TRUE.: WIN does contain the eigenvalues
50 *
51 * JOBZ (global input) CHARACTER*1
52 * Specifies whether or not to compute the eigenvectors:
53 * = 'N': Compute eigenvalues only.
54 * = 'V': Compute eigenvalues and eigenvectors.
55 * Must be 'V' on first call to PDSEPSUBTST
56 *
57 * RANGE (global input) CHARACTER*1
58 * = 'A': all eigenvalues will be found.
59 * = 'V': all eigenvalues in the interval [VL,VU]
60 * will be found.
61 * = 'I': the IL-th through IU-th eigenvalues will be found.
62 * Must be 'A' on first call to PDSEPSUBTST
63 *
64 * UPLO (global input) CHARACTER*1
65 * Specifies whether the upper or lower triangular part of the
66 * symmetric matrix A is stored:
67 * = 'U': Upper triangular
68 * = 'L': Lower triangular
69 *
70 * N (global input) INTEGER
71 * Size of the matrix to be tested. (global size)
72 *
73 * VL (global input) DOUBLE PRECISION
74 * If RANGE='V', the lower bound of the interval to be searched
75 * for eigenvalues. Not referenced if RANGE = 'A' or 'I'.
76 *
77 * VU (global input) DOUBLE PRECISION
78 * If RANGE='V', the upper bound of the interval to be searched
79 * for eigenvalues. Not referenced if RANGE = 'A' or 'I'.
80 *
81 * IL (global input) INTEGER
82 * If RANGE='I', the index (from smallest to largest) of the
83 * smallest eigenvalue to be returned. IL >= 1.
84 * Not referenced if RANGE = 'A' or 'V'.
85 *
86 * IU (global input) INTEGER
87 * If RANGE='I', the index (from smallest to largest) of the
88 * largest eigenvalue to be returned. min(IL,N) <= IU <= N.
89 * Not referenced if RANGE = 'A' or 'V'.
90 *
91 * THRESH (global input) DOUBLE PRECISION
92 * A test will count as "failed" if the "error", computed as
93 * described below, exceeds THRESH. Note that the error
94 * is scaled to be O(1), so THRESH should be a reasonably
95 * small multiple of 1, e.g., 10 or 100. In particular,
96 * it should not depend on the precision (single vs. double)
97 * or the size of the matrix. It must be at least zero.
98 *
99 * ABSTOL (global input) DOUBLE PRECISION
100 * The absolute tolerance for the eigenvalues. An
101 * eigenvalue is considered to be located if it has
102 * been determined to lie in an interval whose width
103 * is "abstol" or less. If "abstol" is less than or equal
104 * to zero, then ulp*|T| will be used, where |T| is
105 * the 1-norm of the matrix. If eigenvectors are
106 * desired later by inverse iteration ("PDSTEIN"),
107 * "abstol" MUST NOT be bigger than ulp*|T|.
108 *
109 * A (local workspace) DOUBLE PRECISION array
110 * global dimension (N, N), local dimension (DESCA(DLEN_), NQ)
111 * A is distributed in a block cyclic manner over both rows
112 * and columns.
113 * See PDSYEVX for a description of block cyclic layout.
114 * The test matrix, which is then modified by PDSYEVX
115 * A has already been padded front and back, use A(1+IPREPAD)
116 *
117 * COPYA (local input) DOUBLE PRECISION array, dimension(N*N)
118 * COPYA holds a copy of the original matrix A
119 * identical in both form and content to A
120 *
121 * Z (local workspace) DOUBLE PRECISION array, dim (N*N)
122 * Z is distributed in the same manner as A
123 * Z contains the eigenvector matrix
124 * Z is used as workspace by the test routines
125 * PDSEPCHK and PDSEPQTQ.
126 * Z has already been padded front and back, use Z(1+IPREPAD)
127 *
128 * IA (global input) INTEGER
129 * On entry, IA specifies the global row index of the submatrix
130 * of the global matrix A, COPYA and Z to operate on.
131 *
132 * JA (global input) INTEGER
133 * On entry, IA specifies the global column index of the submat
134 * of the global matrix A, COPYA and Z to operate on.
135 *
136 * DESCA (global/local input) INTEGER array of dimension 8
137 * The array descriptor for the matrix A, COPYA and Z.
138 *
139 * WIN (global input) DOUBLE PRECISION array, dimension (N)
140 * If .not. WKNOWN, WIN is ignored on input
141 * Otherwise, WIN() is taken as the standard by which the
142 * eigenvalues are to be compared against.
143 *
144 * WNEW (global workspace) DOUBLE PRECISION array, dimension (N)
145 * The eigenvalues as copmuted by this call to PDSYEVX
146 * If JOBZ <> 'V' or RANGE <> 'A' these eigenvalues are
147 * compared against those in WIN().
148 * WNEW has already been padded front and back,
149 * use WNEW(1+IPREPAD)
150 *
151 * IFAIL (global output) INTEGER array, dimension (N)
152 * If JOBZ = 'V', then on normal exit, the first M elements of
153 * IFAIL are zero. If INFO > 0 on exit, then IFAIL contains the
154 * indices of the eigenvectors that failed to converge.
155 * If JOBZ = 'N', then IFAIL is not referenced.
156 * IFAIL has already been padded front and back,
157 * use IFAIL(1+IPREPAD)
158 *
159 * ICLUSTR (global workspace) integer array, dimension (2*NPROW*NPCOL)
160 *
161 * GAP (global workspace) DOUBLE PRECISION array,
162 * dimension (NPROW*NPCOL)
163 *
164 * WORK (local workspace) DOUBLE PRECISION array, dimension (LWORK)
165 * WORK has already been padded front and back,
166 * use WORK(1+IPREPAD)
167 *
168 * LWORK (local input) INTEGER
169 * The actual length of the array WORK after padding.
170 *
171 *
172 * LWORK1 (local input) INTEGER
173 * The amount of real workspace to pass to PDSYEVX
174 *
175 * IWORK (local workspace) INTEGER array, dimension (LIWORK)
176 * IWORK has already been padded front and back,
177 * use IWORK(1+IPREPAD)
178 *
179 * LIWORK (local input) INTEGER
180 * The length of the array IWORK after padding.
181 *
182 * RESULT (global output) INTEGER
183 * The result of this call to PDSYEVX
184 * RESULT = -3 => This process did not participate
185 * RESULT = 0 => All tests passed
186 * RESULT = 1 => ONe or more tests failed
187 *
188 * TSTNRM (global output) DOUBLE PRECISION
189 * |AQ- QL| / (ABSTOL+EPS*|A|)*N
190 *
191 * QTQNRM (global output) DOUBLE PRECISION
192 * |QTQ -I| / N*EPS
193 *
194 * .. Parameters ..
195 *
196  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
197  $ MB_, NB_, RSRC_, CSRC_, LLD_
198  PARAMETER ( BLOCK_CYCLIC_2D = 1, dlen_ = 9, dtype_ = 1,
199  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
200  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
201  DOUBLE PRECISION PADVAL, FIVE, NEGONE
202  PARAMETER ( PADVAL = 13.5285d+0, five = 5.0d+0,
203  $ negone = -1.0d+0 )
204  INTEGER IPADVAL
205  PARAMETER ( IPADVAL = 927 )
206 * ..
207 * .. Local Scalars ..
208  LOGICAL MISSLARGEST, MISSSMALLEST
209  INTEGER I, IAM, INDIWRK, INFO, ISIZESUBTST, ISIZESYEVX,
210  $ isizetst, j, m, maxeigs, maxil, maxiu, maxsize,
211  $ minil, mq, mycol, myil, myrow, nclusters, np,
212  $ npcol, nprow, nq, nz, oldil, oldiu, oldnz, res,
213  $ sizechk, sizemqrleft, sizemqrright, sizeqrf,
214  $ sizeqtq, sizesubtst, sizesyevx, sizetms,
215  $ sizetst, valsize, vecsize
216  DOUBLE PRECISION EPS, EPSNORMA, ERROR, MAXERROR, MAXVU,
217  $ MINERROR, MINVL, NORMWIN, OLDVL, OLDVU, ORFAC,
218  $ SAFMIN
219 * ..
220 * .. Local Arrays ..
221  INTEGER DESCZ( DLEN_ ), DSEED( 4 ), ITMP( 2 )
222 * ..
223 * .. External Functions ..
224 *
225  LOGICAL LSAME
226  INTEGER NUMROC
227  DOUBLE PRECISION PDLAMCH, PDLANSY
228  EXTERNAL LSAME, NUMROC, PDLAMCH, PDLANSY
229 * ..
230 * .. External Subroutines ..
231  EXTERNAL blacs_gridinfo, descinit, dgamn2d, dgamx2d,
232  $ dlacpy, igamn2d, igamx2d, pdchekpad, pdelset,
233  $ pdfillpad, pdlasizesep, pdlasizesyevx,
234  $ pdsepchk, pdsepqtq, pdsyevx, pichekpad,
235  $ pifillpad, slboot, sltimer
236 * ..
237 * .. Intrinsic Functions ..
238  INTRINSIC abs, max, min, mod
239 * ..
240 * .. Executable Statements ..
241 * This is just to keep ftnchek happy
242  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
243  $ rsrc_.LT.0 )RETURN
244  CALL pdlasizesep( desca, iprepad, ipostpad, sizemqrleft,
245  $ sizemqrright, sizeqrf, sizetms, sizeqtq,
246  $ sizechk, sizesyevx, isizesyevx, sizesubtst,
247  $ isizesubtst, sizetst, isizetst )
248 *
249  tstnrm = negone
250  qtqnrm = negone
251  eps = pdlamch( desca( ctxt_ ), 'Eps' )
252  safmin = pdlamch( desca( ctxt_ ), 'Safe min' )
253 *
254  normwin = safmin / eps
255  IF( n.GE.1 )
256  $ normwin = max( abs( win( 1 ) ), abs( win( n ) ), normwin )
257 *
258 * Make sure that we aren't using information from previous calls
259 *
260  nz = -13
261  oldnz = nz
262  oldil = il
263  oldiu = iu
264  oldvl = vl
265  oldvu = vu
266 *
267  DO 10 i = 1, lwork1, 1
268  work( i+iprepad ) = 14.3d+0
269  10 CONTINUE
270  DO 20 i = 1, liwork, 1
271  iwork( i+iprepad ) = 14
272  20 CONTINUE
273 *
274  DO 30 i = 1, n
275  wnew( i+iprepad ) = 3.14159d+0
276  30 CONTINUE
277 *
278  iclustr( 1+iprepad ) = 139
279 *
280  IF( lsame( jobz, 'N' ) ) THEN
281  maxeigs = 0
282  ELSE
283  IF( lsame( range, 'A' ) ) THEN
284  maxeigs = n
285  ELSE IF( lsame( range, 'I' ) ) THEN
286  maxeigs = iu - il + 1
287  ELSE
288  minvl = vl - normwin*five*eps - abstol
289  maxvu = vu + normwin*five*eps + abstol
290  minil = 1
291  maxiu = 0
292  DO 40 i = 1, n
293  IF( win( i ).LT.minvl )
294  $ minil = minil + 1
295  IF( win( i ).LE.maxvu )
296  $ maxiu = maxiu + 1
297  40 CONTINUE
298 *
299  maxeigs = maxiu - minil + 1
300  END IF
301  END IF
302 *
303 *
304  CALL descinit( descz, desca( m_ ), desca( n_ ), desca( mb_ ),
305  $ desca( nb_ ), desca( rsrc_ ), desca( csrc_ ),
306  $ desca( ctxt_ ), desca( lld_ ), info )
307 *
308  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
309  indiwrk = 1 + iprepad + nprow*npcol + 1
310 *
311  iam = 1
312  IF( myrow.EQ.0 .AND. mycol.EQ.0 )
313  $ iam = 0
314 *
315 * If this process is not involved in this test, bail out now
316 *
317  result = -3
318  IF( myrow.GE.nprow .OR. myrow.LT.0 )
319  $ GO TO 150
320  result = 0
321 *
322 *
323 * DSEED is not used in this call to PDLASIZESYEVX, the
324 * following line just makes ftnchek happy.
325 *
326  dseed( 1 ) = 1
327 *
328  CALL pdlasizesyevx( wknown, range, n, desca, vl, vu, il, iu,
329  $ dseed, win, maxsize, vecsize, valsize )
330 *
331  np = numroc( n, desca( mb_ ), myrow, 0, nprow )
332  nq = numroc( n, desca( nb_ ), mycol, 0, npcol )
333  mq = numroc( maxeigs, desca( nb_ ), mycol, 0, npcol )
334 *
335  CALL dlacpy( 'A', np, nq, copya, desca( lld_ ), a( 1+iprepad ),
336  $ desca( lld_ ) )
337 *
338  CALL pdfillpad( desca( ctxt_ ), np, nq, a, desca( lld_ ), iprepad,
339  $ ipostpad, padval )
340 *
341  CALL pdfillpad( descz( ctxt_ ), np, mq, z, descz( lld_ ), iprepad,
342  $ ipostpad, padval+1.0d+0 )
343 *
344  CALL pdfillpad( desca( ctxt_ ), n, 1, wnew, n, iprepad, ipostpad,
345  $ padval+2.0d+0 )
346 *
347  CALL pdfillpad( desca( ctxt_ ), nprow*npcol, 1, gap, nprow*npcol,
348  $ iprepad, ipostpad, padval+3.0d+0 )
349 *
350  CALL pdfillpad( desca( ctxt_ ), lwork1, 1, work, lwork1, iprepad,
351  $ ipostpad, padval+4.0d+0 )
352 *
353  CALL pifillpad( desca( ctxt_ ), liwork, 1, iwork, liwork, iprepad,
354  $ ipostpad, ipadval )
355 *
356  CALL pifillpad( desca( ctxt_ ), n, 1, ifail, n, iprepad, ipostpad,
357  $ ipadval )
358 *
359  CALL pifillpad( desca( ctxt_ ), 2*nprow*npcol, 1, iclustr,
360  $ 2*nprow*npcol, iprepad, ipostpad, ipadval )
361 *
362 * Make sure that PDSYEVX does not cheat (i.e. use answers
363 * already computed.)
364 *
365  DO 60 i = 1, n, 1
366  DO 50 j = 1, maxeigs, 1
367  CALL pdelset( z( 1+iprepad ), i, j, desca, 13.0d+0 )
368  50 CONTINUE
369  60 CONTINUE
370 *
371  orfac = -1.0d+0
372 *
373  CALL slboot
374  CALL sltimer( 1 )
375  CALL sltimer( 6 )
376  CALL pdsyevx( jobz, range, uplo, n, a( 1+iprepad ), ia, ja, desca,
377  $ vl, vu, il, iu, abstol, m, nz, wnew( 1+iprepad ),
378  $ orfac, z( 1+iprepad ), ia, ja, desca,
379  $ work( 1+iprepad ), lwork1, iwork( 1+iprepad ),
380  $ liwork, ifail( 1+iprepad ), iclustr( 1+iprepad ),
381  $ gap( 1+iprepad ), info )
382  CALL sltimer( 6 )
383  CALL sltimer( 1 )
384 *
385  IF( thresh.LE.0 ) THEN
386  result = 0
387  ELSE
388  CALL pdchekpad( desca( ctxt_ ), 'PDSYEVX-A', np, nq, a,
389  $ desca( lld_ ), iprepad, ipostpad, padval )
390 *
391  CALL pdchekpad( descz( ctxt_ ), 'PDSYEVX-Z', np, mq, z,
392  $ descz( lld_ ), iprepad, ipostpad,
393  $ padval+1.0d+0 )
394 *
395  CALL pdchekpad( desca( ctxt_ ), 'PDSYEVX-WNEW', n, 1, wnew, n,
396  $ iprepad, ipostpad, padval+2.0d+0 )
397 *
398  CALL pdchekpad( desca( ctxt_ ), 'PDSYEVX-GAP', nprow*npcol, 1,
399  $ gap, nprow*npcol, iprepad, ipostpad,
400  $ padval+3.0d+0 )
401 *
402  CALL pdchekpad( desca( ctxt_ ), 'PDSYEVX-WORK', lwork1, 1,
403  $ work, lwork1, iprepad, ipostpad,
404  $ padval+4.0d+0 )
405 *
406  CALL pichekpad( desca( ctxt_ ), 'PDSYEVX-IWORK', liwork, 1,
407  $ iwork, liwork, iprepad, ipostpad, ipadval )
408 *
409  CALL pichekpad( desca( ctxt_ ), 'PDSYEVX-IFAIL', n, 1, ifail,
410  $ n, iprepad, ipostpad, ipadval )
411 *
412  CALL pichekpad( desca( ctxt_ ), 'PDSYEVX-ICLUSTR',
413  $ 2*nprow*npcol, 1, iclustr, 2*nprow*npcol,
414  $ iprepad, ipostpad, ipadval )
415 *
416 *
417 * Since we now know the spectrum, we can potentially reduce MAXSIZE.
418 *
419  IF( lsame( range, 'A' ) ) THEN
420  CALL pdlasizesyevx( .true., range, n, desca, vl, vu, il, iu,
421  $ dseed, wnew( 1+iprepad ), maxsize,
422  $ vecsize, valsize )
423  END IF
424 *
425 *
426 * Check INFO
427 *
428 *
429 * Make sure that all processes return the same value of INFO
430 *
431  itmp( 1 ) = info
432  itmp( 2 ) = info
433 *
434  CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
435  $ -1, -1, 0 )
436  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1, 1,
437  $ 1, -1, -1, 0 )
438 *
439 *
440  IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
441  IF( iam.EQ.0 )
442  $ WRITE( nout, fmt = * )
443  $ 'Different processes return different INFO'
444  result = 1
445  ELSE IF( mod( info, 2 ).EQ.1 .OR. info.GT.7 .OR. info.LT.0 )
446  $ THEN
447  IF( iam.EQ.0 )
448  $ WRITE( nout, fmt = 9999 )info
449  result = 1
450  ELSE IF( mod( info / 2, 2 ).EQ.1 .AND. lwork1.GE.maxsize ) THEN
451  IF( iam.EQ.0 )
452  $ WRITE( nout, fmt = 9996 )info
453  result = 1
454  ELSE IF( mod( info / 4, 2 ).EQ.1 .AND. lwork1.GE.vecsize ) THEN
455  IF( iam.EQ.0 )
456  $ WRITE( nout, fmt = 9996 )info
457  result = 1
458  END IF
459 *
460 *
461  IF( lsame( jobz, 'V' ) .AND. ( iclustr( 1+iprepad ).NE.
462  $ 0 ) .AND. ( mod( info / 2, 2 ).NE.1 ) ) THEN
463  IF( iam.EQ.0 )
464  $ WRITE( nout, fmt = 9995 )
465  result = 1
466  END IF
467 *
468 * Check M
469 *
470  IF( ( m.LT.0 ) .OR. ( m.GT.n ) ) THEN
471  IF( iam.EQ.0 )
472  $ WRITE( nout, fmt = 9994 )
473  result = 1
474  ELSE IF( lsame( range, 'A' ) .AND. ( m.NE.n ) ) THEN
475  IF( iam.EQ.0 )
476  $ WRITE( nout, fmt = 9993 )
477  result = 1
478  ELSE IF( lsame( range, 'I' ) .AND. ( m.NE.iu-il+1 ) ) THEN
479  IF( iam.EQ.0 )
480  $ WRITE( nout, fmt = 9992 )
481  result = 1
482  ELSE IF( lsame( jobz, 'V' ) .AND.
483  $ ( .NOT.( lsame( range, 'V' ) ) ) .AND. ( m.NE.nz ) )
484  $ THEN
485  IF( iam.EQ.0 )
486  $ WRITE( nout, fmt = 9991 )
487  result = 1
488  END IF
489 *
490 * Check NZ
491 *
492  IF( lsame( jobz, 'V' ) ) THEN
493  IF( lsame( range, 'V' ) ) THEN
494  IF( nz.GT.m ) THEN
495  IF( iam.EQ.0 )
496  $ WRITE( nout, fmt = 9990 )
497  result = 1
498  END IF
499  IF( nz.LT.m .AND. mod( info / 4, 2 ).NE.1 ) THEN
500  IF( iam.EQ.0 )
501  $ WRITE( nout, fmt = 9989 )
502  result = 1
503  END IF
504  ELSE
505  IF( nz.NE.m ) THEN
506  IF( iam.EQ.0 )
507  $ WRITE( nout, fmt = 9988 )
508  result = 1
509  END IF
510  END IF
511  END IF
512  IF( result.EQ.0 ) THEN
513 *
514 * Make sure that all processes return the same # of eigenvalues
515 *
516  itmp( 1 ) = m
517  itmp( 2 ) = m
518 *
519  CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
520  $ -1, -1, 0 )
521  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1,
522  $ 1, 1, -1, -1, 0 )
523 *
524  IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
525  IF( iam.EQ.0 )
526  $ WRITE( nout, fmt = 9987 )
527  result = 1
528  ELSE
529 *
530 * Make sure that different processes return the same eigenvalues
531 *
532  DO 70 i = 1, m
533  work( i ) = wnew( i+iprepad )
534  work( i+m ) = wnew( i+iprepad )
535  70 CONTINUE
536 *
537  CALL dgamn2d( desca( ctxt_ ), 'a', ' ', m, 1, work, m, 1,
538  $ 1, -1, -1, 0 )
539  CALL dgamx2d( desca( ctxt_ ), 'a', ' ', m, 1,
540  $ work( 1+m ), m, 1, 1, -1, -1, 0 )
541 *
542  DO 80 i = 1, m
543 *
544  IF( result.EQ.0 .AND. ( abs( work( i )-work( m+
545  $ i ) ).GT.five*eps*abs( work( i ) ) ) ) THEN
546  IF( iam.EQ.0 )
547  $ WRITE( nout, fmt = 9986 )
548  result = 1
549  END IF
550  80 CONTINUE
551  END IF
552  END IF
553 *
554 * Make sure that all processes return the same # of clusters
555 *
556  IF( lsame( jobz, 'V' ) ) THEN
557  nclusters = 0
558  DO 90 i = 0, nprow*npcol - 1
559  IF( iclustr( 1+iprepad+2*i ).EQ.0 )
560  $ GO TO 100
561  nclusters = nclusters + 1
562  90 CONTINUE
563  100 CONTINUE
564  itmp( 1 ) = nclusters
565  itmp( 2 ) = nclusters
566 *
567  CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
568  $ -1, -1, 0 )
569  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1,
570  $ 1, 1, -1, -1, 0 )
571 *
572  IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
573  IF( iam.EQ.0 )
574  $ WRITE( nout, fmt = 9985 )
575  result = 1
576  ELSE
577 *
578 * Make sure that different processes return the same clusters
579 *
580  DO 110 i = 1, nclusters
581  iwork( indiwrk+i ) = iclustr( i+iprepad )
582  iwork( indiwrk+i+nclusters ) = iclustr( i+iprepad )
583  110 CONTINUE
584  CALL igamn2d( desca( ctxt_ ), 'a', ' ', nclusters*2+1, 1,
585  $ iwork( indiwrk+1 ), nclusters*2+1, 1, 1,
586  $ -1, -1, 0 )
587  CALL igamx2d( desca( ctxt_ ), 'a', ' ', nclusters*2+1, 1,
588  $ iwork( indiwrk+1+nclusters ),
589  $ nclusters*2+1, 1, 1, -1, -1, 0 )
590 *
591 *
592  DO 120 i = 1, nclusters
593  IF( result.EQ.0 .AND. iwork( indiwrk+i ).NE.
594  $ iwork( indiwrk+nclusters+i ) ) THEN
595  IF( iam.EQ.0 )
596  $ WRITE( nout, fmt = 9984 )
597  result = 1
598  END IF
599  120 CONTINUE
600 *
601  IF( iclustr( 1+iprepad+nclusters*2 ).NE.0 ) THEN
602  IF( iam.EQ.0 )
603  $ WRITE( nout, fmt = 9983 )
604  result = 1
605  END IF
606  END IF
607  END IF
608 *
609 *
610  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1,
611  $ -1, -1, 0 )
612  IF( result.NE.0 )
613  $ GO TO 150
614 *
615 * Compute eps * norm(A)
616 *
617  IF( n.EQ.0 ) THEN
618  epsnorma = eps
619  ELSE
620  epsnorma = pdlansy( 'I', uplo, n, copya, ia, ja, desca,
621  $ work )*eps
622  END IF
623 *
624 * Note that a couple key variables get redefined in PDSEPCHK
625 * as described by this table:
626 *
627 * PDSEPTST name PDSEPCHK name
628 * ------------- -------------
629 * COPYA A
630 * Z Q
631 * A C
632 *
633 *
634  IF( lsame( jobz, 'V' ) ) THEN
635 *
636 * Perform the |AQ - QE| test
637 *
638  CALL pdfillpad( desca( ctxt_ ), sizechk, 1, work, sizechk,
639  $ iprepad, ipostpad, 4.3d+0 )
640 *
641  CALL pdsepchk( n, nz, copya, ia, ja, desca,
642  $ max( abstol+epsnorma, safmin ), thresh,
643  $ z( 1+iprepad ), ia, ja, descz,
644  $ a( 1+iprepad ), ia, ja, desca,
645  $ wnew( 1+iprepad ), work( 1+iprepad ),
646  $ sizechk, tstnrm, res )
647 *
648  CALL pdchekpad( desca( ctxt_ ), 'PDSEPCHK-WORK', sizechk, 1,
649  $ work, sizechk, iprepad, ipostpad, 4.3d+0 )
650 *
651  IF( res.NE.0 )
652  $ result = 1
653 *
654 * Perform the |QTQ - I| test
655 *
656  CALL pdfillpad( desca( ctxt_ ), sizeqtq, 1, work, sizeqtq,
657  $ iprepad, ipostpad, 4.3d+0 )
658 *
659 *
660  CALL pdsepqtq( n, nz, thresh, z( 1+iprepad ), ia, ja, descz,
661  $ a( 1+iprepad ), ia, ja, desca,
662  $ iwork( 1+iprepad+1 ), iclustr( 1+iprepad ),
663  $ gap( 1+iprepad ), work( iprepad+1 ), sizeqtq,
664  $ qtqnrm, info, res )
665 *
666  CALL pdchekpad( desca( ctxt_ ), 'PDSEPQTQ-WORK', sizeqtq, 1,
667  $ work, sizeqtq, iprepad, ipostpad, 4.3d+0 )
668 *
669  IF( res.NE.0 )
670  $ result = 1
671 *
672  IF( info.NE.0 ) THEN
673  IF( iam.EQ.0 )
674  $ WRITE( nout, fmt = 9998 )info
675  result = 1
676  END IF
677  END IF
678 *
679 * Check to make sure that we have the right eigenvalues
680 *
681  IF( wknown ) THEN
682 *
683 * Set up MYIL if necessary
684 *
685  myil = il
686 *
687  IF( lsame( range, 'V' ) ) THEN
688  myil = 1
689  minil = 1
690  maxil = n - m + 1
691  ELSE
692  IF( lsame( range, 'A' ) ) THEN
693  myil = 1
694  END IF
695  minil = myil
696  maxil = myil
697  END IF
698 *
699 * Find the largest difference between the computed
700 * and expected eigenvalues
701 *
702  minerror = normwin
703 *
704  DO 140 myil = minil, maxil
705  maxerror = 0
706 *
707 * Make sure that we aren't skipping any important eigenvalues
708 *
709  misssmallest = .true.
710  IF( .NOT.lsame( range, 'V' ) .OR. ( myil.EQ.1 ) )
711  $ misssmallest = .false.
712  IF( misssmallest .AND. ( win( myil-1 ).LT.vl+normwin*
713  $ five*thresh*eps ) )misssmallest = .false.
714  misslargest = .true.
715  IF( .NOT.lsame( range, 'V' ) .OR. ( myil.EQ.maxil ) )
716  $ misslargest = .false.
717  IF( misslargest .AND. ( win( myil+m ).GT.vu-normwin*five*
718  $ thresh*eps ) )misslargest = .false.
719  IF( .NOT.misssmallest ) THEN
720  IF( .NOT.misslargest ) THEN
721 *
722 * Make sure that the eigenvalues that we report are OK
723 *
724  DO 130 i = 1, m
725  error = abs( win( i+myil-1 )-wnew( i+iprepad ) )
726  maxerror = max( maxerror, error )
727  130 CONTINUE
728 *
729  minerror = min( maxerror, minerror )
730  END IF
731  END IF
732  140 CONTINUE
733 *
734 *
735 * If JOBZ = 'V' and RANGE='A', we might be comparing
736 * against our estimate of what the eigenvalues ought to
737 * be, rather than comparing against what PxSYEVX computed
738 * last time around, so we have to be more generous.
739 *
740  IF( lsame( jobz, 'V' ) .AND. lsame( range, 'A' ) ) THEN
741  IF( minerror.GT.normwin*five*five*thresh*eps ) THEN
742  IF( iam.EQ.0 )
743  $ WRITE( nout, fmt = 9997 )minerror, normwin
744  result = 1
745  END IF
746  ELSE
747  IF( minerror.GT.normwin*five*thresh*eps ) THEN
748  IF( iam.EQ.0 )
749  $ WRITE( nout, fmt = 9997 )minerror, normwin
750  result = 1
751  END IF
752  END IF
753  END IF
754 *
755 *
756 * Make sure that the IL, IU, VL and VU were not altered
757 *
758  IF( il.NE.oldil .OR. iu.NE.oldiu .OR. vl.NE.oldvl .OR. vu.NE.
759  $ oldvu ) THEN
760  IF( iam.EQ.0 )
761  $ WRITE( nout, fmt = 9982 )
762  result = 1
763  END IF
764 *
765  IF( lsame( jobz, 'N' ) .AND. ( nz.NE.oldnz ) ) THEN
766  IF( iam.EQ.0 )
767  $ WRITE( nout, fmt = 9981 )
768  result = 1
769  END IF
770 *
771  END IF
772 *
773 * All processes should report the same result
774 *
775  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1, -1,
776  $ -1, 0 )
777 *
778  150 CONTINUE
779 *
780 *
781  RETURN
782 *
783  9999 FORMAT( 'PDSYEVX returned INFO=', i7 )
784  9998 FORMAT( 'PDSEPQTQ returned INFO=', i7 )
785  9997 FORMAT( 'PDSEPSUBTST minerror =', d11.2, ' normwin=', d11.2 )
786  9996 FORMAT( 'PDSYEVX returned INFO=', i7,
787  $ ' despite adequate workspace' )
788  9995 FORMAT( .NE..NE.'ICLUSTR(1)0 but mod(INFO/2,2)1' )
789  9994 FORMAT( 'M not in the range 0 to N' )
790  9993 FORMAT( 'M not equal to N' )
791  9992 FORMAT( 'M not equal to IU-IL+1' )
792  9991 FORMAT( 'M not equal to NZ' )
793  9990 FORMAT( 'NZ > M' )
794  9989 FORMAT( 'NZ < M' )
795  9988 FORMAT( 'NZ not equal to M' )
796  9987 FORMAT( 'Different processes return different values for M' )
797  9986 FORMAT( 'Different processes return different eigenvalues' )
798  9985 FORMAT( 'Different processes return ',
799  $ 'different numbers of clusters' )
800  9984 FORMAT( 'Different processes return different clusters' )
801  9983 FORMAT( 'ICLUSTR not zero terminated' )
802  9982 FORMAT( 'IL, IU, VL or VU altered by PDSYEVX' )
803  9981 FORMAT( 'NZ altered by PDSYEVX with JOBZ=N' )
804 *
805 * End of PDSEPSUBTST
806 *
807  END
pichekpad
subroutine pichekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pichekpad.f:3
pdsepqtq
subroutine pdsepqtq(MS, NV, THRESH, Q, IQ, JQ, DESCQ, C, IC, JC, DESCC, PROCDIST, ICLUSTR, GAP, WORK, LWORK, QTQNRM, INFO, RES)
Definition: pdsepqtq.f:6
max
#define max(A, B)
Definition: pcgemr.c:180
pdlasizesyevx
subroutine pdlasizesyevx(WKNOWN, RANGE, N, DESCA, VL, VU, IL, IU, ISEED, WIN, MAXSIZE, VECSIZE, VALSIZE)
Definition: pdlasizesyevx.f:5
pdchekpad
subroutine pdchekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pdchekpad.f:3
pdsepsubtst
subroutine pdsepsubtst(WKNOWN, JOBZ, RANGE, UPLO, N, VL, VU, IL, IU, THRESH, ABSTOL, A, COPYA, Z, IA, JA, DESCA, WIN, WNEW, IFAIL, ICLUSTR, GAP, IPREPAD, IPOSTPAD, WORK, LWORK, LWORK1, IWORK, LIWORK, RESULT, TSTNRM, QTQNRM, NOUT)
Definition: pdsepsubtst.f:7
pdlasizesep
subroutine pdlasizesep(DESCA, IPREPAD, IPOSTPAD, SIZEMQRLEFT, SIZEMQRRIGHT, SIZEQRF, SIZETMS, SIZEQTQ, SIZECHK, SIZESYEVX, ISIZESYEVX, SIZESUBTST, ISIZESUBTST, SIZETST, ISIZETST)
Definition: pdlasizesep.f:8
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
slboot
subroutine slboot()
Definition: sltimer.f:2
pdfillpad
subroutine pdfillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pdfillpad.f:2
pdsyevx
subroutine pdsyevx(JOBZ, RANGE, UPLO, N, A, IA, JA, DESCA, VL, VU, IL, IU, ABSTOL, M, NZ, W, ORFAC, Z, IZ, JZ, DESCZ, WORK, LWORK, IWORK, LIWORK, IFAIL, ICLUSTR, GAP, INFO)
Definition: pdsyevx.f:5
pdsepchk
subroutine pdsepchk(MS, NV, A, IA, JA, DESCA, EPSNORMA, THRESH, Q, IQ, JQ, DESCQ, C, IC, JC, DESCC, W, WORK, LWORK, TSTNRM, RESULT)
Definition: pdsepchk.f:6
pdelset
subroutine pdelset(A, IA, JA, DESCA, ALPHA)
Definition: pdelset.f:2
pifillpad
subroutine pifillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pifillpad.f:2
min
#define min(A, B)
Definition: pcgemr.c:181