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