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