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