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