ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
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
pichekpad
subroutine pichekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pichekpad.f:3
max
#define max(A, B)
Definition: pcgemr.c:180
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
pslasizesyevx
subroutine pslasizesyevx(WKNOWN, RANGE, N, DESCA, VL, VU, IL, IU, ISEED, WIN, MAXSIZE, VECSIZE, VALSIZE)
Definition: pslasizesyevx.f:5
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
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
pselset
subroutine pselset(A, IA, JA, DESCA, ALPHA)
Definition: pselset.f:2
pschekpad
subroutine pschekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pschekpad.f:3
pslasizesep
subroutine pslasizesep(DESCA, IPREPAD, IPOSTPAD, SIZEMQRLEFT, SIZEMQRRIGHT, SIZEQRF, SIZETMS, SIZEQTQ, SIZECHK, SIZESYEVX, ISIZESYEVX, SIZESUBTST, ISIZESUBTST, SIZETST, ISIZETST)
Definition: pslasizesep.f:8
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
slboot
subroutine slboot()
Definition: sltimer.f:2
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
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
psfillpad
subroutine psfillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: psfillpad.f:2
pifillpad
subroutine pifillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pifillpad.f:2
min
#define min(A, B)
Definition: pcgemr.c:181