SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pssepsubtst.f
Go to the documentation of this file.
1 SUBROUTINE pssepsubtst( 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 REAL ABSTOL, QTQNRM, THRESH, TSTNRM, VL, VU
19* ..
20* .. Array Arguments ..
21 INTEGER DESCA( * ), ICLUSTR( * ), IFAIL( * ),
22 $ IWORK( * )
23 REAL A( * ), COPYA( * ), GAP( * ), WIN( * ),
24 $ WNEW( * ), WORK( * ), Z( * )
25* ..
26*
27* Purpose
28* =======
29*
30* PSSEPSUBTST calls PSSYEVX and then tests the output of
31* PSSYEVX
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* PSSYEVX 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 PSSEPSUBTST
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 PSSEPSUBTST
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) REAL
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) REAL
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) REAL
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) REAL
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 ("PSSTEIN"),
107* "abstol" MUST NOT be bigger than ulp*|T|.
108*
109* A (local workspace) REAL 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 PSSYEVX for a description of block cyclic layout.
114* The test matrix, which is then modified by PSSYEVX
115* A has already been padded front and back, use A(1+IPREPAD)
116*
117* COPYA (local input) REAL 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) REAL 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* PSSEPCHK and PSSEPQTQ.
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) REAL 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) REAL array, dimension (N)
145* The eigenvalues as copmuted by this call to PSSYEVX
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) REAL array,
162* dimension (NPROW*NPCOL)
163*
164* WORK (local workspace) REAL 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 PSSYEVX
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 PSSYEVX
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) REAL
189* |AQ- QL| / (ABSTOL+EPS*|A|)*N
190*
191* QTQNRM (global output) REAL
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 REAL PADVAL, FIVE, NEGONE
202 PARAMETER ( PADVAL = 13.5285e+0, five = 5.0e+0,
203 $ negone = -1.0e+0 )
204 INTEGER IPADVAL
205 PARAMETER ( IPADVAL = 927 )
206* ..
207* .. Local Scalars ..
208 INTEGER I, IAM, INDIWRK, INFO, ISIZESUBTST, ISIZESYEVX,
209 $ isizetst, j, m, maxeigs, maxil, maxiu, maxsize,
210 $ minil, mq, mycol, myil, myrow, nclusters, np,
211 $ npcol, nprow, nq, nz, oldil, oldiu, oldnz, res,
212 $ sizechk, sizemqrleft, sizemqrright, sizeqrf,
213 $ sizeqtq, sizesubtst, sizesyevx, sizetms,
214 $ sizetst, valsize, vecsize
215 REAL EPS, EPSNORMA, ERROR, MAXERROR, MAXVU,
216 $ MINERROR, MINVL, NORMWIN, OLDVL, OLDVU, ORFAC,
217 $ SAFMIN
218* ..
219* .. Local Arrays ..
220 INTEGER DESCZ( DLEN_ ), DSEED( 4 ), ITMP( 2 )
221* ..
222* .. External Functions ..
223*
224 LOGICAL LSAME
225 INTEGER NUMROC
226 REAL PSLAMCH, PSLANSY
227 EXTERNAL LSAME, NUMROC, PSLAMCH, PSLANSY
228* ..
229* .. External Subroutines ..
230 EXTERNAL blacs_gridinfo, descinit, igamn2d, igamx2d,
231 $ pichekpad, pifillpad, pschekpad, pselset,
232 $ psfillpad, pslasizesep, pslasizesyevx,
233 $ pssepchk, pssepqtq, pssyevx, sgamn2d, sgamx2d,
234 $ slacpy, slboot, sltimer
235* ..
236* .. Intrinsic Functions ..
237 INTRINSIC abs, max, min, mod
238* ..
239* .. Executable Statements ..
240* This is just to keep ftnchek happy
241 IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
242 $ rsrc_.LT.0 )RETURN
243 CALL pslasizesep( desca, iprepad, ipostpad, sizemqrleft,
244 $ sizemqrright, sizeqrf, sizetms, sizeqtq,
245 $ sizechk, sizesyevx, isizesyevx, sizesubtst,
246 $ isizesubtst, sizetst, isizetst )
247*
248 tstnrm = negone
249 qtqnrm = negone
250 eps = pslamch( desca( ctxt_ ), 'Eps' )
251 safmin = pslamch( desca( ctxt_ ), 'Safe min' )
252 normwin = safmin / eps
253 IF( n.GE.1 )
254 $ normwin = max( abs( win( 1 ) ), abs( win( n ) ), normwin )
255*
256*
257* Make sure that we aren't using information from previous calls
258*
259 nz = -13
260 oldnz = nz
261 oldil = il
262 oldiu = iu
263 oldvl = vl
264 oldvu = vu
265*
266 DO 10 i = 1, lwork1, 1
267 work( i+iprepad ) = 14.3e+0
268 10 CONTINUE
269 DO 20 i = 1, liwork, 1
270 iwork( i+iprepad ) = 14
271 20 CONTINUE
272*
273 DO 30 i = 1, n
274 wnew( i+iprepad ) = 3.14159e+0
275 30 CONTINUE
276*
277 iclustr( 1+iprepad ) = 139
278*
279 IF( lsame( jobz, 'N' ) ) THEN
280 maxeigs = 0
281 ELSE
282 IF( lsame( range, 'A' ) ) THEN
283 maxeigs = n
284 ELSE IF( lsame( range, 'I' ) ) THEN
285 maxeigs = iu - il + 1
286 ELSE
287 minvl = vl - normwin*five*eps - abstol
288 maxvu = vu + normwin*five*eps + abstol
289 minil = 1
290 maxiu = 0
291 DO 40 i = 1, n
292 IF( win( i ).LT.minvl )
293 $ minil = minil + 1
294 IF( win( i ).LE.maxvu )
295 $ maxiu = maxiu + 1
296 40 CONTINUE
297*
298 maxeigs = maxiu - minil + 1
299 END IF
300 END IF
301*
302*
303 CALL descinit( descz, desca( m_ ), desca( n_ ), desca( mb_ ),
304 $ desca( nb_ ), desca( rsrc_ ), desca( csrc_ ),
305 $ desca( ctxt_ ), desca( lld_ ), info )
306*
307 CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
308 indiwrk = 1 + iprepad + nprow*npcol + 1
309*
310 iam = 1
311 IF( myrow.EQ.0 .AND. mycol.EQ.0 )
312 $ iam = 0
313*
314* If this process is not involved in this test, bail out now
315*
316 result = -3
317 IF( myrow.GE.nprow .OR. myrow.LT.0 )
318 $ GO TO 150
319 result = 0
320*
321*
322* DSEED is not used in this call to PSLASIZESYEVX, the
323* following line just makes ftnchek happy.
324*
325 dseed( 1 ) = 1
326*
327 CALL pslasizesyevx( wknown, range, n, desca, vl, vu, il, iu,
328 $ dseed, win, maxsize, vecsize, valsize )
329*
330 np = numroc( n, desca( mb_ ), myrow, 0, nprow )
331 nq = numroc( n, desca( nb_ ), mycol, 0, npcol )
332 mq = numroc( maxeigs, desca( nb_ ), mycol, 0, npcol )
333*
334 CALL slacpy( 'A', np, nq, copya, desca( lld_ ), a( 1+iprepad ),
335 $ desca( lld_ ) )
336*
337 CALL psfillpad( desca( ctxt_ ), np, nq, a, desca( lld_ ), iprepad,
338 $ ipostpad, padval )
339*
340 CALL psfillpad( descz( ctxt_ ), np, mq, z, descz( lld_ ), iprepad,
341 $ ipostpad, padval+1.0e+0 )
342*
343 CALL psfillpad( desca( ctxt_ ), n, 1, wnew, n, iprepad, ipostpad,
344 $ padval+2.0e+0 )
345*
346 CALL psfillpad( desca( ctxt_ ), nprow*npcol, 1, gap, nprow*npcol,
347 $ iprepad, ipostpad, padval+3.0e+0 )
348*
349 CALL psfillpad( desca( ctxt_ ), lwork1, 1, work, lwork1, iprepad,
350 $ ipostpad, padval+4.0e+0 )
351*
352 CALL pifillpad( desca( ctxt_ ), liwork, 1, iwork, liwork, iprepad,
353 $ ipostpad, ipadval )
354*
355 CALL pifillpad( desca( ctxt_ ), n, 1, ifail, n, iprepad, ipostpad,
356 $ ipadval )
357*
358 CALL pifillpad( desca( ctxt_ ), 2*nprow*npcol, 1, iclustr,
359 $ 2*nprow*npcol, iprepad, ipostpad, ipadval )
360*
361* Make sure that PSSYEVX does not cheat (i.e. use answers
362* already computed.)
363*
364 DO 60 i = 1, n, 1
365 DO 50 j = 1, maxeigs, 1
366 CALL pselset( z( 1+iprepad ), i, j, desca, 13.0e+0 )
367 50 CONTINUE
368 60 CONTINUE
369*
370 orfac = -1.0e+0
371*
372 CALL slboot
373 CALL sltimer( 1 )
374 CALL sltimer( 6 )
375 CALL pssyevx( jobz, range, uplo, n, a( 1+iprepad ), ia, ja, desca,
376 $ vl, vu, il, iu, abstol, m, nz, wnew( 1+iprepad ),
377 $ orfac, z( 1+iprepad ), ia, ja, desca,
378 $ work( 1+iprepad ), lwork1, iwork( 1+iprepad ),
379 $ liwork, ifail( 1+iprepad ), iclustr( 1+iprepad ),
380 $ gap( 1+iprepad ), info )
381 CALL sltimer( 6 )
382 CALL sltimer( 1 )
383*
384 IF( thresh.LE.0 ) THEN
385 result = 0
386 ELSE
387 CALL pschekpad( desca( ctxt_ ), 'PSSYEVX-A', np, nq, a,
388 $ desca( lld_ ), iprepad, ipostpad, padval )
389*
390 CALL pschekpad( descz( ctxt_ ), 'PSSYEVX-Z', np, mq, z,
391 $ descz( lld_ ), iprepad, ipostpad,
392 $ padval+1.0e+0 )
393*
394 CALL pschekpad( desca( ctxt_ ), 'PSSYEVX-WNEW', n, 1, wnew, n,
395 $ iprepad, ipostpad, padval+2.0e+0 )
396*
397 CALL pschekpad( desca( ctxt_ ), 'PSSYEVX-GAP', nprow*npcol, 1,
398 $ gap, nprow*npcol, iprepad, ipostpad,
399 $ padval+3.0e+0 )
400*
401 CALL pschekpad( desca( ctxt_ ), 'PSSYEVX-WORK', lwork1, 1,
402 $ work, lwork1, iprepad, ipostpad,
403 $ padval+4.0e+0 )
404*
405 CALL pichekpad( desca( ctxt_ ), 'PSSYEVX-IWORK', liwork, 1,
406 $ iwork, liwork, iprepad, ipostpad, ipadval )
407*
408 CALL pichekpad( desca( ctxt_ ), 'PSSYEVX-IFAIL', n, 1, ifail,
409 $ n, iprepad, ipostpad, ipadval )
410*
411 CALL pichekpad( desca( ctxt_ ), 'PSSYEVX-ICLUSTR',
412 $ 2*nprow*npcol, 1, iclustr, 2*nprow*npcol,
413 $ iprepad, ipostpad, ipadval )
414*
415*
416* Since we now know the spectrum, we can potentially reduce MAXSIZE.
417*
418 IF( lsame( range, 'A' ) ) THEN
419 CALL pslasizesyevx( .true., range, n, desca, vl, vu, il, iu,
420 $ dseed, wnew( 1+iprepad ), maxsize,
421 $ vecsize, valsize )
422 END IF
423*
424*
425* Check INFO
426*
427*
428* Make sure that all processes return the same value of INFO
429*
430 itmp( 1 ) = info
431 itmp( 2 ) = info
432*
433 CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
434 $ -1, -1, 0 )
435 CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1, 1,
436 $ 1, -1, -1, 0 )
437*
438*
439 IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
440 IF( iam.EQ.0 )
441 $ WRITE( nout, fmt = * )
442 $ 'Different processes return different INFO'
443 result = 1
444 ELSE IF( mod( info, 2 ).EQ.1 .OR. info.GT.7 .OR. info.LT.0 )
445 $ THEN
446 IF( iam.EQ.0 )
447 $ WRITE( nout, fmt = 9999 )info
448 result = 1
449 ELSE IF( mod( info / 2, 2 ).EQ.1 .AND. lwork1.GE.maxsize ) THEN
450 IF( iam.EQ.0 )
451 $ WRITE( nout, fmt = 9996 )info
452 result = 1
453 ELSE IF( mod( info / 4, 2 ).EQ.1 .AND. lwork1.GE.vecsize ) THEN
454 IF( iam.EQ.0 )
455 $ WRITE( nout, fmt = 9996 )info
456 result = 1
457 END IF
458*
459*
460 IF( lsame( jobz, 'V' ) .AND. ( iclustr( 1+iprepad ).NE.
461 $ 0 ) .AND. ( mod( info / 2, 2 ).NE.1 ) ) THEN
462 IF( iam.EQ.0 )
463 $ WRITE( nout, fmt = 9995 )
464 result = 1
465 END IF
466*
467* Check M
468*
469 IF( ( m.LT.0 ) .OR. ( m.GT.n ) ) THEN
470 IF( iam.EQ.0 )
471 $ WRITE( nout, fmt = 9994 )
472 result = 1
473 ELSE IF( lsame( range, 'A' ) .AND. ( m.NE.n ) ) THEN
474 IF( iam.EQ.0 )
475 $ WRITE( nout, fmt = 9993 )
476 result = 1
477 ELSE IF( lsame( range, 'I' ) .AND. ( m.NE.iu-il+1 ) ) THEN
478 IF( iam.EQ.0 )
479 $ WRITE( nout, fmt = 9992 )
480 result = 1
481 ELSE IF( lsame( jobz, 'V' ) .AND.
482 $ ( .NOT.( lsame( range, 'V' ) ) ) .AND. ( m.NE.nz ) )
483 $ THEN
484 IF( iam.EQ.0 )
485 $ WRITE( nout, fmt = 9991 )
486 result = 1
487 END IF
488*
489* Check NZ
490*
491 IF( lsame( jobz, 'V' ) ) THEN
492 IF( lsame( range, 'V' ) ) THEN
493 IF( nz.GT.m ) THEN
494 IF( iam.EQ.0 )
495 $ WRITE( nout, fmt = 9990 )
496 result = 1
497 END IF
498 IF( nz.LT.m .AND. mod( info / 4, 2 ).NE.1 ) THEN
499 IF( iam.EQ.0 )
500 $ WRITE( nout, fmt = 9989 )
501 result = 1
502 END IF
503 ELSE
504 IF( nz.NE.m ) THEN
505 IF( iam.EQ.0 )
506 $ WRITE( nout, fmt = 9988 )
507 result = 1
508 END IF
509 END IF
510 END IF
511 IF( result.EQ.0 ) THEN
512*
513* Make sure that all processes return the same # of eigenvalues
514*
515 itmp( 1 ) = m
516 itmp( 2 ) = m
517*
518 CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
519 $ -1, -1, 0 )
520 CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1,
521 $ 1, 1, -1, -1, 0 )
522*
523 IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
524 IF( iam.EQ.0 )
525 $ WRITE( nout, fmt = 9987 )
526 result = 1
527 ELSE
528*
529* Make sure that different processes return the same eigenvalues
530*
531 DO 70 i = 1, m
532 work( i ) = wnew( i+iprepad )
533 work( i+m ) = wnew( i+iprepad )
534 70 CONTINUE
535*
536 CALL sgamn2d( desca( ctxt_ ), 'a', ' ', m, 1, work, m, 1,
537 $ 1, -1, -1, 0 )
538 CALL sgamx2d( desca( ctxt_ ), 'a', ' ', m, 1,
539 $ work( 1+m ), m, 1, 1, -1, -1, 0 )
540*
541 DO 80 i = 1, m
542 IF( result.EQ.0 .AND. work( i ).NE.work( m+i ) ) THEN
543 IF( iam.EQ.0 )
544 $ WRITE( nout, fmt = 9986 )
545 result = 1
546 END IF
547 80 CONTINUE
548 END IF
549 END IF
550*
551* Make sure that all processes return the same # of clusters
552*
553 IF( lsame( jobz, 'V' ) ) THEN
554 nclusters = 0
555 DO 90 i = 0, nprow*npcol - 1
556 IF( iclustr( 1+iprepad+2*i ).EQ.0 )
557 $ GO TO 100
558 nclusters = nclusters + 1
559 90 CONTINUE
560 100 CONTINUE
561 itmp( 1 ) = nclusters
562 itmp( 2 ) = nclusters
563*
564 CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
565 $ -1, -1, 0 )
566 CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1,
567 $ 1, 1, -1, -1, 0 )
568*
569 IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
570 IF( iam.EQ.0 )
571 $ WRITE( nout, fmt = 9985 )
572 result = 1
573 ELSE
574*
575* Make sure that different processes return the same clusters
576*
577 DO 110 i = 1, nclusters
578 iwork( indiwrk+i ) = iclustr( i+iprepad )
579 iwork( indiwrk+i+nclusters ) = iclustr( i+iprepad )
580 110 CONTINUE
581 CALL igamn2d( desca( ctxt_ ), 'a', ' ', nclusters*2+1, 1,
582 $ iwork( indiwrk+1 ), nclusters*2+1, 1, 1,
583 $ -1, -1, 0 )
584 CALL igamx2d( desca( ctxt_ ), 'a', ' ', nclusters*2+1, 1,
585 $ iwork( indiwrk+1+nclusters ),
586 $ nclusters*2+1, 1, 1, -1, -1, 0 )
587*
588*
589 DO 120 i = 1, nclusters
590 IF( result.EQ.0 .AND. iwork( indiwrk+i ).NE.
591 $ iwork( indiwrk+nclusters+i ) ) THEN
592 IF( iam.EQ.0 )
593 $ WRITE( nout, fmt = 9984 )
594 result = 1
595 END IF
596 120 CONTINUE
597*
598 IF( iclustr( 1+iprepad+nclusters*2 ).NE.0 ) THEN
599 IF( iam.EQ.0 )
600 $ WRITE( nout, fmt = 9983 )
601 result = 1
602 END IF
603 END IF
604 END IF
605*
606*
607 CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1,
608 $ -1, -1, 0 )
609 IF( result.NE.0 )
610 $ GO TO 150
611*
612* Compute eps * norm(A)
613*
614 IF( n.EQ.0 ) THEN
615 epsnorma = eps
616 ELSE
617 epsnorma = pslansy( 'I', uplo, n, copya, ia, ja, desca,
618 $ work )*eps
619 END IF
620*
621* Note that a couple key variables get redefined in PSSEPCHK
622* as described by this table:
623*
624* PSSEPTST name PSSEPCHK name
625* ------------- -------------
626* COPYA A
627* Z Q
628* A C
629*
630*
631 IF( lsame( jobz, 'V' ) ) THEN
632*
633* Perform the |AQ - QE| test
634*
635 CALL psfillpad( desca( ctxt_ ), sizechk, 1, work, sizechk,
636 $ iprepad, ipostpad, 4.3e+0 )
637*
638 CALL pssepchk( n, nz, copya, ia, ja, desca,
639 $ max( abstol+epsnorma, safmin ), thresh,
640 $ z( 1+iprepad ), ia, ja, descz,
641 $ a( 1+iprepad ), ia, ja, desca,
642 $ wnew( 1+iprepad ), work( 1+iprepad ),
643 $ sizechk, tstnrm, res )
644*
645 CALL pschekpad( desca( ctxt_ ), 'PSSEPCHK-WORK', sizechk, 1,
646 $ work, sizechk, iprepad, ipostpad, 4.3e+0 )
647*
648 IF( res.NE.0 )
649 $ result = 1
650*
651* Perform the |QTQ - I| test
652*
653 CALL psfillpad( desca( ctxt_ ), sizeqtq, 1, work, sizeqtq,
654 $ iprepad, ipostpad, 4.3e+0 )
655*
656*
657 CALL pssepqtq( n, nz, thresh, z( 1+iprepad ), ia, ja, descz,
658 $ a( 1+iprepad ), ia, ja, desca,
659 $ iwork( 1+iprepad+1 ), iclustr( 1+iprepad ),
660 $ gap( 1+iprepad ), work( iprepad+1 ), sizeqtq,
661 $ qtqnrm, info, res )
662*
663 CALL pschekpad( desca( ctxt_ ), 'PSSEPQTQ-WORK', sizeqtq, 1,
664 $ work, sizeqtq, iprepad, ipostpad, 4.3e+0 )
665*
666 IF( res.NE.0 )
667 $ result = 1
668*
669 IF( info.NE.0 ) THEN
670 IF( iam.EQ.0 )
671 $ WRITE( nout, fmt = 9998 )info
672 result = 1
673 END IF
674 END IF
675*
676* Check to make sure that we have the right eigenvalues
677*
678 IF( wknown ) THEN
679*
680* Set up MYIL if necessary
681*
682 myil = il
683*
684 IF( lsame( range, 'V' ) ) THEN
685 myil = 1
686 minil = 1
687 maxil = n - m + 1
688 ELSE
689 IF( lsame( range, 'A' ) ) THEN
690 myil = 1
691 END IF
692 minil = myil
693 maxil = myil
694 END IF
695*
696* Find the largest difference between the computed
697* and expected eigenvalues
698*
699 minerror = normwin
700*
701 DO 140 myil = minil, maxil
702 maxerror = 0
703*
704* Make sure that we aren't skipping any important eigenvalues
705*
706 IF( .NOT.lsame( range, 'V' ) .OR.
707 $ ( myil.EQ.1 .OR. ( win( myil-1 ).LT.vl+normwin*five*
708 $ thresh*eps ) ) ) THEN
709 IF( .NOT.lsame( range, 'V' ) .OR.
710 $ ( myil.EQ.n-m+1 .OR. ( win( myil+m ).GT.vu-
711 $ normwin*five*thresh*eps ) ) ) THEN
712*
713* Make sure that the eigenvalues that we report are OK
714*
715 DO 130 i = 1, m
716 error = abs( win( i+myil-1 )-wnew( i+iprepad ) )
717 maxerror = max( maxerror, error )
718 130 CONTINUE
719*
720 minerror = min( maxerror, minerror )
721 END IF
722 END IF
723 140 CONTINUE
724*
725*
726* If JOBZ = 'V' and RANGE='A', we might be comparing
727* against our estimate of what the eigenvalues ought to
728* be, rather than comparing against what PxSYEVX computed
729* last time around, so we have to be more generous.
730*
731 IF( lsame( jobz, 'V' ) .AND. lsame( range, 'A' ) ) THEN
732 IF( minerror.GT.normwin*five*five*thresh*eps ) THEN
733 IF( iam.EQ.0 )
734 $ WRITE( nout, fmt = 9997 )minerror, normwin
735 result = 1
736 END IF
737 ELSE
738 IF( minerror.GT.normwin*five*thresh*eps ) THEN
739 IF( iam.EQ.0 )
740 $ WRITE( nout, fmt = 9997 )minerror, normwin
741 result = 1
742 END IF
743 END IF
744 END IF
745*
746*
747* Make sure that the IL, IU, VL and VU were not altered
748*
749 IF( il.NE.oldil .OR. iu.NE.oldiu .OR. vl.NE.oldvl .OR. vu.NE.
750 $ oldvu ) THEN
751 IF( iam.EQ.0 )
752 $ WRITE( nout, fmt = 9982 )
753 result = 1
754 END IF
755*
756 IF( lsame( jobz, 'N' ) .AND. ( nz.NE.oldnz ) ) THEN
757 IF( iam.EQ.0 )
758 $ WRITE( nout, fmt = 9981 )
759 result = 1
760 END IF
761*
762 END IF
763*
764* All processes should report the same result
765*
766 CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1, -1,
767 $ -1, 0 )
768*
769 150 CONTINUE
770*
771*
772 RETURN
773*
774 9999 FORMAT( 'PSSYEVX returned INFO=', i7 )
775 9998 FORMAT( 'PSSEPQTQ returned INFO=', i7 )
776 9997 FORMAT( 'PSSEPSUBTST minerror =', d11.2, ' normwin=', d11.2 )
777 9996 FORMAT( 'PSSYEVX returned INFO=', i7,
778 $ ' despite adequate workspace' )
779 9995 FORMAT( .NE..NE.'ICLUSTR(1)0 but mod(INFO/2,2)1' )
780 9994 FORMAT( 'M not in the range 0 to N' )
781 9993 FORMAT( 'M not equal to N' )
782 9992 FORMAT( 'M not equal to IU-IL+1' )
783 9991 FORMAT( 'M not equal to NZ' )
784 9990 FORMAT( 'NZ > M' )
785 9989 FORMAT( 'NZ < M' )
786 9988 FORMAT( 'NZ not equal to M' )
787 9987 FORMAT( 'Different processes return different values for M' )
788 9986 FORMAT( 'Different processes return different eigenvalues' )
789 9985 FORMAT( 'Different processes return ',
790 $ 'different numbers of clusters' )
791 9984 FORMAT( 'Different processes return different clusters' )
792 9983 FORMAT( 'ICLUSTR not zero terminated' )
793 9982 FORMAT( 'IL, IU, VL or VU altered by PSSYEVX' )
794 9981 FORMAT( 'NZ altered by PSSYEVX with JOBZ=N' )
795*
796* End of PSSEPSUBTST
797*
798 END
descinit
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
Definition descinit.f:3
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
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
pschekpad
subroutine pschekpad(ictxt, mess, m, n, a, lda, ipre, ipost, chkval)
Definition pschekpad.f:3
pselset
subroutine pselset(a, ia, ja, desca, alpha)
Definition pselset.f:2
psfillpad
subroutine psfillpad(ictxt, m, n, a, lda, ipre, ipost, chkval)
Definition psfillpad.f:2
pslasizesep
subroutine pslasizesep(desca, iprepad, ipostpad, sizemqrleft, sizemqrright, sizeqrf, sizetms, sizeqtq, sizechk, sizesyevx, isizesyevx, sizesubtst, isizesubtst, sizetst, isizetst)
Definition pslasizesep.f:8
pslasizesyevx
subroutine pslasizesyevx(wknown, range, n, desca, vl, vu, il, iu, iseed, win, maxsize, vecsize, valsize)
Definition pslasizesyevx.f:5
pssepchk
subroutine pssepchk(ms, nv, a, ia, ja, desca, epsnorma, thresh, q, iq, jq, descq, c, ic, jc, descc, w, work, lwork, tstnrm, result)
Definition pssepchk.f:6
pssepqtq
subroutine pssepqtq(ms, nv, thresh, q, iq, jq, descq, c, ic, jc, descc, procdist, iclustr, gap, work, lwork, qtqnrm, info, res)
Definition pssepqtq.f:6
pssepsubtst
subroutine pssepsubtst(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 pssepsubtst.f:7
pssyevx
subroutine pssyevx(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 pssyevx.f:5
slboot
subroutine slboot()
Definition sltimer.f:2
sltimer
subroutine sltimer(i)
Definition sltimer.f:47