SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ pzseprsubtst()

subroutine pzseprsubtst ( logical  wknown,
character  jobz,
character  range,
character  uplo,
integer  n,
double precision  vl,
double precision  vu,
integer  il,
integer  iu,
double precision  thresh,
double precision  abstol,
complex*16, dimension( * )  a,
complex*16, dimension( * )  copya,
complex*16, dimension( * )  z,
integer  ia,
integer  ja,
integer, dimension( * )  desca,
double precision, dimension( * )  win,
double precision, dimension( * )  wnew,
integer, dimension( * )  ifail,
integer, dimension( * )  iclustr,
double precision, dimension( * )  gap,
integer  iprepad,
integer  ipostpad,
complex*16, dimension( * )  work,
integer  lwork,
double precision, dimension( * )  rwork,
integer  lrwork,
integer  lwork1,
integer, dimension( * )  iwork,
integer  liwork,
integer  result,
double precision  tstnrm,
double precision  qtqnrm,
integer  nout 
)

Definition at line 1 of file pzseprsubtst.f.

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