ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pcseprsubtst.f
Go to the documentation of this file.
1  SUBROUTINE pcseprsubtst( 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  REAL ABSTOL, QTQNRM, THRESH, TSTNRM, VL, VU
21  INTEGER LRWORK
22 * ..
23 * .. Array Arguments ..
24  INTEGER DESCA( * ), ICLUSTR( * ), IFAIL( * ),
25  $ iwork( * )
26  COMPLEX A( * ), COPYA( * ), WORK( * ), Z( * )
27  REAL GAP( * ), RWORK( * ), WIN( * ), WNEW( * )
28 * ..
29 *
30 * Purpose
31 * =======
32 *
33 * PCSEPRSUBTST calls PCSYEVR 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) REAL
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) REAL
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) REAL
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) REAL
101 * The absolute tolerance for the residual test.
102 *
103 * A (local workspace) COMPLEX 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 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 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 * PCSEPCHK and PCSEPQTQ.
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) REAL 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) REAL 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) REAL array,
155 * dimension (NPROW*NPCOL)
156 *
157 * WORK (local workspace) COMPLEX 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) REAL
188 * |AQ- QL| / (ABSTOL+EPS*|A|)*N
189 *
190 * QTQNRM (global output) REAL
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  REAL PADVAL, FIVE, NEGONE
201  PARAMETER ( PADVAL = 13.5285e0, five = 5.0e0,
202  $ negone = -1.0e0 )
203  COMPLEX ZPADVAL
204  PARAMETER ( ZPADVAL = ( 13.989e0, 1.93e0 ) )
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  REAL 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  REAL PSLAMCH, PCLANHE
230  EXTERNAL LSAME, NUMROC, PSLAMCH, PCLANHE
231 * ..
232 * .. External Subroutines ..
233  EXTERNAL blacs_gridinfo, clacpy, descinit, igamn2d,
234  $ igamx2d, pcchekpad, pcelset, pcfillpad,
235  $ pcheevr, pclasizeheevr, pclasizesepr, pcsepchk,
236  $ pcsepqtq, pichekpad, pifillpad, pschekpad,
237  $ psfillpad, sgamn2d, sgamx2d, slboot, sltimer
238 * ..
239 * .. Intrinsic Functions ..
240  INTRINSIC abs, max, min, mod
241 * ..
242 * .. Executable Statements ..
243 *
244  CALL pclasizesepr( desca, iprepad, ipostpad, sizemqrleft,
245  $ sizemqrright, sizeqrf, sizetms, sizeqtq,
246  $ sizechk, sizeevr, rsizeevr, isizeevr,
247  $ sizesubtst, rsizesubtst, isizesubtst,
248  $ sizetst, rsizetst, isizetst )
249 *
250  tstnrm = negone
251  qtqnrm = negone
252  eps = pslamch( desca( ctxt_ ), 'Eps' )
253  safmin = pslamch( desca( ctxt_ ), 'Safe min' )
254 *
255  normwin = safmin / eps
256  IF( n.GE.1 )
257  $ normwin = max( abs( win( 1 ) ), abs( win( n ) ), normwin )
258 *
259 * Make sure that no information from previous calls is used
260 *
261  nz = -13
262  oldnz = nz
263  oldil = il
264  oldiu = iu
265  oldvl = vl
266  oldvu = vu
267 *
268  DO 10 i = 1, lwork1, 1
269  rwork( i+iprepad ) = 14.3e0
270  10 CONTINUE
271 *
272  DO 15 i = 1, lwork, 1
273  work( i+iprepad ) = ( 15.63e0, 1.1e0 )
274  15 CONTINUE
275 *
276  DO 20 i = 1, liwork, 1
277  iwork( i+iprepad ) = 14
278  20 CONTINUE
279 *
280  DO 30 i = 1, n
281  wnew( i+iprepad ) = 3.14159e0
282  30 CONTINUE
283 *
284  iclustr( 1+iprepad ) = 139
285 *
286  IF (lsame( range, 'V' ) ) THEN
287 * WRITE(*,*) 'VL VU = ', VL, ' ', VU
288  END IF
289 
290  IF( lsame( jobz, 'N' ) ) THEN
291  maxeigs = 0
292  ELSE
293  IF( lsame( range, 'A' ) ) THEN
294  maxeigs = n
295  ELSE IF( lsame( range, 'I' ) ) THEN
296  maxeigs = iu - il + 1
297  ELSE
298  minvl = vl - normwin*five*eps - abstol
299  maxvu = vu + normwin*five*eps + abstol
300 * WRITE(*,*) 'MINVL = ', MINVL, ' MAXVU = ', MAXVU
301 * WRITE(*,*) 'WIN = ', WIN( 1 )
302  minil = 1
303  maxiu = 0
304  DO 40 i = 1, n
305  IF( win( i ).LT.minvl )
306  $ minil = minil + 1
307  IF( win( i ).LE.maxvu )
308  $ maxiu = maxiu + 1
309  40 CONTINUE
310 *
311  maxeigs = maxiu - minil + 1
312  END IF
313  END IF
314 *
315 *
316  CALL descinit( descz, desca( m_ ), desca( n_ ), desca( mb_ ),
317  $ desca( nb_ ), desca( rsrc_ ), desca( csrc_ ),
318  $ desca( ctxt_ ), desca( lld_ ), info )
319 *
320  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
321  indiwrk = 1 + iprepad + nprow*npcol + 1
322 *
323  iam = 1
324  IF( myrow.EQ.0 .AND. mycol.EQ.0 )
325  $ iam = 0
326 *
327 * If this process is not involved in this test, bail out now
328 *
329  result = -3
330  IF( myrow.GE.nprow .OR. myrow.LT.0 )
331  $ GO TO 150
332  result = 0
333 *
334  iseed( 1 ) = 1
335 *
336  CALL pclasizeheevr( wknown, range, n, desca, vl, vu, il, iu,
337  $ iseed, win, maxsize, vecsize, valsize )
338 *
339  np = numroc( n, desca( mb_ ), myrow, 0, nprow )
340  nq = numroc( n, desca( nb_ ), mycol, 0, npcol )
341  mq = numroc( maxeigs, desca( nb_ ), mycol, 0, npcol )
342 *
343  CALL clacpy( 'A', np, nq, copya, desca( lld_ ), a( 1+iprepad ),
344  $ desca( lld_ ) )
345 *
346  CALL pcfillpad( desca( ctxt_ ), np, nq, a, desca( lld_ ), iprepad,
347  $ ipostpad, zpadval )
348 *
349  CALL pcfillpad( descz( ctxt_ ), np, mq, z, descz( lld_ ), iprepad,
350  $ ipostpad, zpadval+1.0e0 )
351 *
352 * WRITE(*,*) ' NP = ', NP, ' MQ = ', MQ, ' LDZ = ', DESCZ( LLD_ ),
353 * $ ' IPREPAD = ', IPREPAD, ' IPOSTPAD = ', IPOSTPAD,
354 * $ ' MAXEIGS = ', MAXEIGS
355 * WRITE(*,*) ' PADZ( 1 ) = ', Z( 1 ), ' PADZ( 2 ) = ', Z( 2 ),
356 * $ ' PADZ( 3 ) = ', Z( 3 ), ' PADZ( 4 ) = ', Z( 4 )
357 *
358  CALL psfillpad( desca( ctxt_ ), n, 1, wnew, n, iprepad, ipostpad,
359  $ padval+2.0e0 )
360 *
361  CALL psfillpad( desca( ctxt_ ), nprow*npcol, 1, gap, nprow*npcol,
362  $ iprepad, ipostpad, padval+3.0e0 )
363 *
364  CALL psfillpad( desca( ctxt_ ), lwork1, 1, rwork,lwork1, iprepad,
365  $ ipostpad, padval+4.0e0 )
366 *
367  CALL pifillpad( desca( ctxt_ ), liwork, 1, iwork, liwork, iprepad,
368  $ ipostpad, ipadval )
369 *
370  CALL pifillpad( desca( ctxt_ ), n, 1, ifail, n, iprepad, ipostpad,
371  $ ipadval )
372 *
373  CALL pifillpad( desca( ctxt_ ), 2*nprow*npcol, 1, iclustr,
374  $ 2*nprow*npcol, iprepad, ipostpad, ipadval )
375 *
376  CALL pcfillpad( desca( ctxt_ ), lwork, 1, work, lwork, iprepad,
377  $ ipostpad, zpadval+4.1e0 )
378 *
379 * Make sure that PCHEEVR does not cheat (i.e. use answers
380 * already computed.)
381 *
382  DO 60 i = 1, n, 1
383  DO 50 j = 1, maxeigs, 1
384  CALL pcelset( z( 1+iprepad ), i, j, desca,
385  $ ( 13.0e0, 1.34e0 ) )
386  50 CONTINUE
387  60 CONTINUE
388 *
389 * Reset and start the timer
390 *
391  CALL slboot
392  CALL sltimer( 1 )
393  CALL sltimer( 6 )
394 
395 *********************************
396 *
397 * Main call to PCHEEVR
398 *
399  CALL pcheevr( jobz, range, uplo, n, a( 1+iprepad ), ia, ja, desca,
400  $ vl, vu, il, iu, m, nz, wnew( 1+iprepad ),
401  $ z( 1+iprepad ), ia, ja, desca,
402  $ work( 1+iprepad ), sizeevr,
403  $ rwork( 1+iprepad ), lwork1,
404  $ iwork( 1+iprepad ), liwork, info )
405 *
406 *********************************
407 *
408 * Stop timer
409 *
410  CALL sltimer( 6 )
411  CALL sltimer( 1 )
412 *
413 * Indicate that there are no unresolved clusters.
414 * This is necessary so that the tester
415 * (adapted from the one originally made for PSSYEVX)
416 * works correctly.
417  iclustr( 1+iprepad ) = 0
418 *
419  IF( thresh.LE.0 ) THEN
420  result = 0
421  ELSE
422  CALL pcchekpad( desca( ctxt_ ), 'PCHEEVR-A', np, nq, a,
423  $ desca( lld_ ), iprepad, ipostpad, zpadval )
424 *
425  CALL pcchekpad( descz( ctxt_ ), 'PCHEEVR-Z', np, mq, z,
426  $ descz( lld_ ), iprepad, ipostpad,
427  $ zpadval+1.0e0 )
428 *
429  CALL pschekpad( desca( ctxt_ ), 'PCHEEVR-WNEW', n, 1, wnew, n,
430  $ iprepad, ipostpad, padval+2.0e0 )
431 *
432  CALL pschekpad( desca( ctxt_ ), 'PCHEEVR-GAP', nprow*npcol, 1,
433  $ gap, nprow*npcol, iprepad, ipostpad,
434  $ padval+3.0e0 )
435 *
436  CALL pschekpad( desca( ctxt_ ), 'PCHEEVR-RWORK',lwork1, 1,
437  $ rwork, lwork1, iprepad, ipostpad,
438  $ padval+4.0e0 )
439 *
440  CALL pcchekpad( desca( ctxt_ ), 'PCHEEVR-WORK',lwork, 1,
441  $ work, lwork, iprepad, ipostpad,
442  $ zpadval+4.1e0 )
443 *
444  CALL pichekpad( desca( ctxt_ ), 'PCHEEVR-IWORK', liwork, 1,
445  $ iwork, liwork, iprepad, ipostpad, ipadval )
446 *
447  CALL pichekpad( desca( ctxt_ ), 'PCHEEVR-IFAIL', n, 1, ifail,
448  $ n, iprepad, ipostpad, ipadval )
449 *
450  CALL pichekpad( desca( ctxt_ ), 'PCHEEVR-ICLUSTR',
451  $ 2*nprow*npcol, 1, iclustr, 2*nprow*npcol,
452  $ iprepad, ipostpad, ipadval )
453 *
454 * If we now know the spectrum, we can potentially reduce MAXSIZE.
455 *
456  IF( lsame( range, 'A' ) ) THEN
457  CALL pclasizeheevr( .true., range, n, desca, vl, vu, il, iu,
458  $ iseed, wnew( 1+iprepad ), maxsize,
459  $ vecsize, valsize )
460  END IF
461 *
462 * Check INFO
463 * Make sure that all processes return the same value of INFO
464 *
465  itmp( 1 ) = info
466  itmp( 2 ) = info
467 *
468  CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
469  $ -1, -1, 0 )
470  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1, 1,
471  $ 1, -1, -1, 0 )
472 *
473 *
474  IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
475  IF( iam.EQ.0 )
476  $ WRITE( nout, fmt = * )
477  $ 'Different processes return different INFO'
478  result = 1
479  ELSE IF( mod( info, 2 ).EQ.1 .OR. info.GT.7 .OR. info.LT.0 )
480  $ THEN
481  IF( iam.EQ.0 )
482  $ WRITE( nout, fmt = 9999 )info
483  result = 1
484  ELSE IF( mod( info / 2, 2 ).EQ.1 .AND. lwork1.GE.maxsize ) THEN
485  IF( iam.EQ.0 )
486  $ WRITE( nout, fmt = 9996 )info
487  result = 1
488  ELSE IF( mod( info / 4, 2 ).EQ.1 .AND. lwork1.GE.vecsize ) THEN
489  IF( iam.EQ.0 )
490  $ WRITE( nout, fmt = 9996 )info
491  result = 1
492  END IF
493 *
494  IF( lsame( jobz, 'V' ) .AND. ( iclustr( 1+iprepad ).NE.
495  $ 0 ) .AND. ( mod( info / 2, 2 ).NE.1 ) ) THEN
496  IF( iam.EQ.0 )
497  $ WRITE( nout, fmt = 9995 )
498  result = 1
499  END IF
500 *
501 * Check M
502 *
503  IF( ( m.LT.0 ) .OR. ( m.GT.n ) ) THEN
504  IF( iam.EQ.0 )
505  $ WRITE( nout, fmt = 9994 )
506  WRITE( nout,*) 'M = ', m, '\n', 'N = ', n
507  result = 1
508  ELSE IF( lsame( range, 'A' ) .AND. ( m.NE.n ) ) THEN
509  IF( iam.EQ.0 )
510  $ WRITE( nout, fmt = 9993 )
511  result = 1
512  ELSE IF( lsame( range, 'I' ) .AND. ( m.NE.iu-il+1 ) ) THEN
513  IF( iam.EQ.0 ) THEN
514  WRITE( nout, fmt = 9992 )
515  WRITE( nout,*) 'IL = ', il, ' IU = ', iu, ' M = ', m
516  END IF
517  result = 1
518  ELSE IF( lsame( jobz, 'V' ) .AND.
519  $ ( .NOT.( lsame( range, 'V' ) ) ) .AND. ( m.NE.nz ) )
520  $ THEN
521  IF( iam.EQ.0 )
522  $ WRITE( nout, fmt = 9991 )
523  result = 1
524  END IF
525 *
526 * Check NZ
527 *
528  IF( lsame( jobz, 'V' ) ) THEN
529  IF( lsame( range, 'V' ) ) THEN
530  IF( nz.GT.m ) THEN
531  IF( iam.EQ.0 )
532  $ WRITE( nout, fmt = 9990 )
533  result = 1
534  END IF
535  IF( nz.LT.m .AND. mod( info / 4, 2 ).NE.1 ) THEN
536  IF( iam.EQ.0 )
537  $ WRITE( nout, fmt = 9989 )
538  result = 1
539  END IF
540  ELSE
541  IF( nz.NE.m ) THEN
542  IF( iam.EQ.0 )
543  $ WRITE( nout, fmt = 9988 )
544  result = 1
545  END IF
546  END IF
547  END IF
548  IF( result.EQ.0 ) THEN
549 *
550 * Make sure that all processes return the same # of eigenvalues
551 *
552  itmp( 1 ) = m
553  itmp( 2 ) = m
554 *
555  CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
556  $ -1, -1, 0 )
557  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1,
558  $ 1, 1, -1, -1, 0 )
559 *
560  IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
561  IF( iam.EQ.0 )
562  $ WRITE( nout, fmt = 9987 )
563  result = 1
564  ELSE
565 *
566 * Ensure that different processes return the same eigenvalues
567 *
568  DO 70 i = 1, m
569  rwork( i ) = wnew( i+iprepad )
570  rwork( i+m ) = wnew( i+iprepad )
571  70 CONTINUE
572 *
573  CALL sgamn2d( desca( ctxt_ ), 'a', ' ', m, 1, rwork, m,
574  $ 1, 1, -1, -1, 0 )
575  CALL sgamx2d( desca( ctxt_ ), 'a', ' ', m, 1,
576  $ rwork( 1+m ), m, 1, 1, -1, -1, 0 )
577 *
578  DO 80 i = 1, m
579  IF( result.EQ.0 .AND. ( abs( rwork( i )-rwork( m+
580  $ i ) ).GT.five*eps*abs( rwork( i ) ) ) ) THEN
581  IF( iam.EQ.0 )
582  $ WRITE( nout, fmt = 9986 )
583  result = 1
584  END IF
585  80 CONTINUE
586  END IF
587  END IF
588 *
589 * Make sure that all processes return the same # of clusters
590 *
591  IF( lsame( jobz, 'V' ) ) THEN
592  nclusters = 0
593  DO 90 i = 0, nprow*npcol - 1
594  IF( iclustr( 1+iprepad+2*i ).EQ.0 )
595  $ GO TO 100
596  nclusters = nclusters + 1
597  90 CONTINUE
598  100 CONTINUE
599  itmp( 1 ) = nclusters
600  itmp( 2 ) = nclusters
601 *
602  CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
603  $ -1, -1, 0 )
604  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1,
605  $ 1, 1, -1, -1, 0 )
606 *
607  IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
608  IF( iam.EQ.0 )
609  $ WRITE( nout, fmt = 9985 )
610  result = 1
611  ELSE
612 *
613 * Make sure that different processes return the same clusters
614 *
615  DO 110 i = 1, nclusters
616  iwork( indiwrk+i ) = iclustr( i+iprepad )
617  iwork( indiwrk+i+nclusters ) = iclustr( i+iprepad )
618  110 CONTINUE
619  CALL igamn2d( desca( ctxt_ ), 'a', ' ', nclusters*2+1, 1,
620  $ iwork( indiwrk+1 ), nclusters*2+1, 1, 1,
621  $ -1, -1, 0 )
622  CALL igamx2d( desca( ctxt_ ), 'a', ' ', nclusters*2+1, 1,
623  $ iwork( indiwrk+1+nclusters ),
624  $ nclusters*2+1, 1, 1, -1, -1, 0 )
625 *
626  DO 120 i = 1, nclusters
627  IF( result.EQ.0 .AND. iwork( indiwrk+i ).NE.
628  $ iwork( indiwrk+nclusters+i ) ) THEN
629  IF( iam.EQ.0 )
630  $ WRITE( nout, fmt = 9984 )
631  result = 1
632  END IF
633  120 CONTINUE
634 *
635  IF( iclustr( 1+iprepad+nclusters*2 ).NE.0 ) THEN
636  IF( iam.EQ.0 )
637  $ WRITE( nout, fmt = 9983 )
638  result = 1
639  END IF
640  END IF
641  END IF
642 *
643  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1,
644  $ -1, -1, 0 )
645  IF( result.NE.0 )
646  $ GO TO 150
647 *
648 * Compute eps * norm(A)
649 *
650  IF( n.EQ.0 ) THEN
651  epsnorma = eps
652  ELSE
653  epsnorma = pclanhe( 'I', uplo, n, copya, ia, ja, desca,
654  $ rwork )*eps
655  END IF
656 *
657  IF( lsame( jobz, 'V' ) ) THEN
658 *
659 * Perform the |A Z - Z W| test
660 *
661  CALL psfillpad( desca( ctxt_ ), sizechk, 1, rwork,sizechk,
662  $ iprepad, ipostpad, 4.3e0 )
663 *
664  CALL pcsepchk( n, nz, copya, ia, ja, desca,
665  $ max( abstol+epsnorma, safmin ), thresh,
666  $ z( 1+iprepad ), ia, ja, descz,
667  $ a( 1+iprepad ), ia, ja, desca,
668  $ wnew( 1+iprepad ), rwork( 1+iprepad ),
669  $ sizechk, tstnrm, res )
670 *
671  CALL pschekpad( desca( ctxt_ ), 'PCSEPCHK-RWORK',sizechk, 1,
672  $ rwork,sizechk, iprepad, ipostpad, 4.3e0 )
673 *
674  IF( res.NE.0 )
675  $ result = 1
676 *
677 * Perform the |QTQ - I| test
678 *
679  CALL psfillpad( desca( ctxt_ ), sizeqtq, 1,rwork, sizeqtq,
680  $ iprepad, ipostpad, 4.3e0 )
681 *
682 *
683  CALL pcsepqtq( n, nz, thresh, z( 1+iprepad ), ia, ja, descz,
684  $ a( 1+iprepad ), ia, ja, desca,
685  $ iwork( 1+iprepad+1 ), iclustr( 1+iprepad ),
686  $ gap( 1+iprepad ),rwork( iprepad+1 ), sizeqtq,
687  $ qtqnrm, info, res )
688 *
689  CALL pschekpad( desca( ctxt_ ), 'PSSEPQTQ-RWORK',sizeqtq, 1,
690  $ rwork,sizeqtq, iprepad, ipostpad, 4.3e0 )
691 *
692  IF( res.NE.0 )
693  $ result = 1
694 *
695  IF( info.NE.0 ) THEN
696  IF( iam.EQ.0 )
697  $ WRITE( nout, fmt = 9998 )info
698  result = 1
699  END IF
700  END IF
701 *
702 * Check to make sure that the right eigenvalues have been obtained
703 *
704  IF( wknown ) THEN
705 * Set up MYIL if necessary
706  myil = il
707 *
708  IF( lsame( range, 'V' ) ) THEN
709  myil = 1
710  minil = 1
711  maxil = n - m + 1
712  ELSE
713  IF( lsame( range, 'A' ) ) THEN
714  myil = 1
715  END IF
716  minil = myil
717  maxil = myil
718  END IF
719 *
720 * Find the largest difference between the computed
721 * and expected eigenvalues
722 *
723  minerror = normwin
724 *
725  DO 140 myil = minil, maxil
726  maxerror = 0
727 *
728 * Make sure that we aren't skipping any important eigenvalues
729 *
730  misssmallest = .true.
731  IF( .NOT.lsame( range, 'V' ) .OR. ( myil.EQ.1 ) )
732  $ misssmallest = .false.
733  IF( misssmallest .AND. ( win( myil-1 ).LT.vl+normwin*
734  $ five*thresh*eps ) )misssmallest = .false.
735  misslargest = .true.
736  IF( .NOT.lsame( range, 'V' ) .OR. ( myil.EQ.maxil ) )
737  $ misslargest = .false.
738  IF( misslargest .AND. ( win( myil+m ).GT.vu-normwin*five*
739  $ thresh*eps ) )misslargest = .false.
740  IF( .NOT.misssmallest ) THEN
741  IF( .NOT.misslargest ) THEN
742 *
743 * Make sure that the eigenvalues that we report are OK
744 *
745  DO 130 i = 1, m
746 * WRITE(*,*) 'WIN WNEW = ', WIN( I+MYIL-1 ),
747 * $ WNEW( I+IPREPAD )
748  error = abs( win( i+myil-1 )-wnew( i+iprepad ) )
749  maxerror = max( maxerror, error )
750  130 CONTINUE
751 *
752  minerror = min( maxerror, minerror )
753  END IF
754  END IF
755  140 CONTINUE
756 *
757 * If JOBZ = 'V' and RANGE='A', we might be comparing
758 * against our estimate of what the eigenvalues ought to
759 * be, rather than comparing against what was computed
760 * last time around, so we have to be more generous.
761 *
762  IF( lsame( jobz, 'V' ) .AND. lsame( range, 'A' ) ) THEN
763  IF( minerror.GT.normwin*five*five*thresh*eps ) THEN
764  IF( iam.EQ.0 )
765  $ WRITE( nout, fmt = 9997 )minerror, normwin
766  result = 1
767  END IF
768  ELSE
769  IF( minerror.GT.normwin*five*thresh*eps ) THEN
770  IF( iam.EQ.0 )
771  $ WRITE( nout, fmt = 9997 )minerror, normwin
772  result = 1
773  END IF
774  END IF
775  END IF
776 *
777 * Make sure that the IL, IU, VL and VU were not altered
778 *
779  IF( il.NE.oldil .OR. iu.NE.oldiu .OR. vl.NE.oldvl .OR. vu.NE.
780  $ oldvu ) THEN
781  IF( iam.EQ.0 )
782  $ WRITE( nout, fmt = 9982 )
783  result = 1
784  END IF
785 *
786  IF( lsame( jobz, 'N' ) .AND. ( nz.NE.oldnz ) ) THEN
787  IF( iam.EQ.0 )
788  $ WRITE( nout, fmt = 9981 )
789  result = 1
790  END IF
791 *
792  END IF
793 *
794 * All processes should report the same result
795 *
796  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1, -1,
797  $ -1, 0 )
798 *
799  150 CONTINUE
800 *
801  RETURN
802 *
803  9999 FORMAT( 'PCHEEVR returned INFO=', i7 )
804  9998 FORMAT( 'PCSEPQTQ returned INFO=', i7 )
805  9997 FORMAT( 'PCSEPRSUBTST minerror =', d11.2, ' normwin=', d11.2 )
806  9996 FORMAT( 'PCHEEVR returned INFO=', i7,
807  $ ' despite adequate workspace' )
808  9995 FORMAT( .NE..NE.'ICLUSTR(1)0 but mod(INFO/2,2)1' )
809  9994 FORMAT( 'M not in the range 0 to N' )
810  9993 FORMAT( 'M not equal to N' )
811  9992 FORMAT( 'M not equal to IU-IL+1' )
812  9991 FORMAT( 'M not equal to NZ' )
813  9990 FORMAT( 'NZ > M' )
814  9989 FORMAT( 'NZ < M' )
815  9988 FORMAT( 'NZ not equal to M' )
816  9987 FORMAT( 'Different processes return different values for M' )
817  9986 FORMAT( 'Different processes return different eigenvalues' )
818  9985 FORMAT( 'Different processes return ',
819  $ 'different numbers of clusters' )
820  9984 FORMAT( 'Different processes return different clusters' )
821  9983 FORMAT( 'ICLUSTR not zero terminated' )
822  9982 FORMAT( 'IL, IU, VL or VU altered by PCHEEVR' )
823  9981 FORMAT( 'NZ altered by PCHEEVR with JOBZ=N' )
824 *
825 * End of PCSEPRSUBTST
826 *
827  END
pichekpad
subroutine pichekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pichekpad.f:3
max
#define max(A, B)
Definition: pcgemr.c:180
pcsepqtq
subroutine pcsepqtq(MS, NV, THRESH, Q, IQ, JQ, DESCQ, C, IC, JC, DESCC, PROCDIST, ICLUSTR, GAP, WORK, LWORK, QTQNRM, INFO, RES)
Definition: pcsepqtq.f:6
pcheevr
subroutine pcheevr(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: pcheevr.f:6
pclasizeheevr
subroutine pclasizeheevr(WKNOWN, RANGE, N, DESCA, VL, VU, IL, IU, ISEED, WIN, MAXSIZE, VECSIZE, VALSIZE)
Definition: pclasizeheevr.f:3
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
pcsepchk
subroutine pcsepchk(MS, NV, A, IA, JA, DESCA, EPSNORMA, THRESH, Q, IQ, JQ, DESCQ, C, IC, JC, DESCC, W, WORK, LWORK, TSTNRM, RESULT)
Definition: pcsepchk.f:6
pschekpad
subroutine pschekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pschekpad.f:3
pclasizesepr
subroutine pclasizesepr(DESCA, IPREPAD, IPOSTPAD, SIZEMQRLEFT, SIZEMQRRIGHT, SIZEQRF, SIZETMS, SIZEQTQ, SIZECHK, SIZEHEEVR, RSIZEHEEVR, ISIZEHEEVR, SIZESUBTST, RSIZESUBTST, ISIZESUBTST, SIZETST, RSIZETST, ISIZETST)
Definition: pclasizesepr.f:7
pcchekpad
subroutine pcchekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pcchekpad.f:3
pcelset
subroutine pcelset(A, IA, JA, DESCA, ALPHA)
Definition: pcelset.f:2
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
slboot
subroutine slboot()
Definition: sltimer.f:2
pcseprsubtst
subroutine pcseprsubtst(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: pcseprsubtst.f:7
pcfillpad
subroutine pcfillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pcfillpad.f:2
psfillpad
subroutine psfillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: psfillpad.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