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