ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
pdsdpsubtst.f
Go to the documentation of this file.
1  SUBROUTINE pdsdpsubtst( WKNOWN, UPLO, N, THRESH, ABSTOL, A,
2  \$ COPYA, Z, IA, JA, DESCA, WIN, WNEW,
4  \$ IWORK, LIWORK,
5  \$ RESULT, TSTNRM, QTQNRM, NOUT )
6 *
7 * -- ScaLAPACK testing routine (version 1.7) --
8 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
9 * and University of California, Berkeley.
10 * March 16, 2000
11 *
12 * .. Scalar Arguments ..
13  LOGICAL WKNOWN
14  CHARACTER UPLO
16  \$ NOUT, RESULT, LIWORK
17  DOUBLE PRECISION ABSTOL, QTQNRM, THRESH, TSTNRM
18 * ..
19 * .. Array Arguments ..
20  INTEGER DESCA( * ), IWORK( * )
21  DOUBLE PRECISION A( * ), COPYA( * ), WIN( * ), WNEW( * ),
22  \$ WORK( * ), Z( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * PDSDPSUBTST calls PDSYEVD and then tests the output of
29 * PDSYEVD
30 * The following two tests are performed:
31 * |AQ -QL| / (abstol + eps * norm(A) ) < N*THRESH
32 * |QT * Q - I| / eps * norm(A) < N*THRESH
33 * If WKNOWN then
34 * we check to make sure that the eigenvalues match expectations
35 * i.e. |WIN - WNEW(1+IPREPAD)| / (eps * |WIN|) < THRESH
36 * where WIN is the array of eigenvalues as computed by
37 * PDSYEVD when eigenvectors are requested
38 *
39 * Arguments
40 * =========
41 *
42 * NP = the number of rows local to a given process.
43 * NQ = the number of columns local to a given process.
44 *
45 * WKNOWN (global input) INTEGER
46 * .FALSE.: WIN does not contain the eigenvalues
47 * .TRUE.: WIN does contain the eigenvalues
48 *
49 * UPLO (global input) CHARACTER*1
50 * Specifies whether the upper or lower triangular part of the
51 * symmetric matrix A is stored:
52 * = 'U': Upper triangular
53 * = 'L': Lower triangular
54 *
55 * N (global input) INTEGER
56 * Size of the matrix to be tested. (global size)
57 *
58 * THRESH (global input) DOUBLE PRECISION
59 * A test will count as "failed" if the "error", computed as
60 * described below, exceeds THRESH. Note that the error
61 * is scaled to be O(1), so THRESH should be a reasonably
62 * small multiple of 1, e.g., 10 or 100. In particular,
63 * it should not depend on the precision (single vs. double)
64 * or the size of the matrix. It must be at least zero.
65 *
66 * ABSTOL (global input) DOUBLE PRECISION
67 * The absolute tolerance for the eigenvalues. An
68 * eigenvalue is considered to be located if it has
69 * been determined to lie in an interval whose width
70 * is "abstol" or less. If "abstol" is less than or equal
71 * to zero, then ulp*|T| will be used, where |T| is
72 * the 1-norm of the matrix.
73 *
74 * A (local workspace) DOUBLE PRECISION array
75 * global dimension (N, N), local dimension (DESCA(DLEN_), NQ)
76 * A is distributed in a block cyclic manner over both rows
77 * and columns.
78 * See PDSYEVD for a description of block cyclic layout.
79 * The test matrix, which is then modified by PDSYEVD
81 *
82 * COPYA (local input) DOUBLE PRECISION array, dimension(N*N)
83 * COPYA holds a copy of the original matrix A
84 * identical in both form and content to A
85 *
86 * Z (local workspace) DOUBLE PRECISION array, dim (N*N)
87 * Z is distributed in the same manner as A
88 * Z contains the eigenvector matrix
89 * Z is used as workspace by the test routines
90 * PDSEPCHK and PDSEPQTQ.
92 *
93 * IA (global input) INTEGER
94 * On entry, IA specifies the global row index of the submatrix
95 * of the global matrix A, COPYA and Z to operate on.
96 *
97 * JA (global input) INTEGER
98 * On entry, IA specifies the global column index of the submat
99 * of the global matrix A, COPYA and Z to operate on.
100 *
101 * DESCA (global/local input) INTEGER array of dimension 8
102 * The array descriptor for the matrix A, COPYA and Z.
103 *
104 * WIN (global input) DOUBLE PRECISION array, dimension (N)
105 * If .not. WKNOWN, WIN is ignored on input
106 * Otherwise, WIN() is taken as the standard by which the
107 * eigenvalues are to be compared against.
108 *
109 * WNEW (global workspace) DOUBLE PRECISION array, dimension (N)
110 * The eigenvalues as computed by this call to PDSYEVD.
113 *
114 * WORK (local workspace) DOUBLE PRECISION array, dimension (LWORK)
117 *
118 * LWORK (local input) INTEGER
119 * The actual length of the array WORK after padding.
120 *
121 *
122 * LWORK1 (local input) INTEGER
123 * The amount of real workspace to pass to PDSYEVD
124 *
125 * IWORK (local workspace) INTEGER array, dimension (LIWORK)
128 *
129 * LIWORK (local input) INTEGER
130 * The length of the array IWORK after padding.
131 *
132 * RESULT (global output) INTEGER
133 * The result of this call to PDSYEVD
134 * RESULT = -3 => This process did not participate
135 * RESULT = 0 => All tests passed
136 * RESULT = 1 => ONe or more tests failed
137 *
138 * TSTNRM (global output) DOUBLE PRECISION
139 * |AQ- QL| / (ABSTOL+EPS*|A|)*N
140 *
141 * QTQNRM (global output) DOUBLE PRECISION
142 * |QTQ -I| / N*EPS
143 *
144 * .. Parameters ..
145 *
146  INTEGER BLOCK_CYCLIC_2D, DLEN_, DT_, CTXT_, M_, N_,
147  \$ MB_, NB_, RSRC_, CSRC_, LLD_
148  PARAMETER ( BLOCK_CYCLIC_2D = 1, dlen_ = 9, dt_ = 1,
149  \$ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
150  \$ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
151  DOUBLE PRECISION FIVE, NEGONE, PADVAL, ZERO
152  parameter( padval = 13.5285d+0, five = 5.0d+0,
153  \$ negone = -1.0d+0, zero = 0.0d+0 )
154 * ..
155 * .. Local Scalars ..
156  INTEGER I, IAM, INFO, ISIZESUBTST, ISIZESYEVX,
157  \$ ISIZETST, J, MINSIZE, MQ, MYCOL, MYROW,
158  \$ NP, NPCOL, NPROW, NQ, RESAQ, RESQTQ,
159  \$ SIZECHK, SIZEMQRLEFT, SIZEMQRRIGHT, SIZEQRF,
160  \$ sizeqtq, sizesubtst, sizesyev, sizesyevx,
161  \$ sizetms, sizetst, sizesyevd, isizesyevd,
162  \$ trilwmin
163  DOUBLE PRECISION EPS, EPSNORMA, ERROR, MAXERROR, MINERROR,
164  \$ NORMWIN, SAFMIN
165 * ..
166 * .. Local Arrays ..
167  INTEGER DESCZ( DLEN_ ), ITMP( 2 )
168 * ..
169 * .. External Functions ..
170 *
171  LOGICAL LSAME
172  INTEGER NUMROC
173  DOUBLE PRECISION PDLAMCH, PDLANSY
174  EXTERNAL LSAME, NUMROC, PDLAMCH, PDLANSY
175 * ..
176 * .. External Subroutines ..
177  EXTERNAL blacs_gridinfo, descinit, igamn2d, igamx2d,
179  \$ pdsepchk, pdsepqtq, pdsyevd, dgamn2d,
180  \$ dgamx2d, dlacpy, slboot, sltimer
181 * ..
182 * .. Intrinsic Functions ..
183  INTRINSIC abs, max, min, mod
184 * ..
185 * .. Executable Statements ..
186 * This is just to keep ftnchek happy
187  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dt_*lld_*mb_*m_*nb_*n_*
188  \$ rsrc_.LT.0 )RETURN
190  \$ sizemqrright, sizeqrf, sizetms, sizeqtq,
191  \$ sizechk, sizesyevx, isizesyevx, sizesyev,
192  \$ sizesyevd, isizesyevd, sizesubtst, isizesubtst,
193  \$ sizetst, isizetst )
194 *
195  tstnrm = negone
196  qtqnrm = negone
197  eps = pdlamch( desca( ctxt_ ), 'Eps' )
198  safmin = pdlamch( desca( ctxt_ ), 'Safe min' )
199 *
200  normwin = safmin / eps
201  IF( n.GE.1 )
202  \$ normwin = max( abs( win( 1+iprepad ) ),
203  \$ abs( win( n+iprepad ) ), normwin )
204 *
205 * Make sure that we aren't using information from previous calls
206 *
207  DO 10 i = 1, lwork1, 1
208  work( i+iprepad ) = 14.3d+0
209  10 CONTINUE
210 *
211  DO 30 i = 1, n
212  wnew( i+iprepad ) = 3.14159d+0
213  30 CONTINUE
214 *
215  CALL descinit( descz, desca( m_ ), desca( n_ ), desca( mb_ ),
216  \$ desca( nb_ ), desca( rsrc_ ), desca( csrc_ ),
217  \$ desca( ctxt_ ), desca( lld_ ), info )
218 *
219  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
220 *
221  iam = 1
222  IF( myrow.EQ.0 .AND. mycol.EQ.0 )
223  \$ iam = 0
224 *
225 * If this process is not involved in this test, bail out now
226 *
227  IF( myrow.GE.nprow .OR. myrow.LT.0 )
228  \$ GO TO 150
229  result = 0
230 *
231  np = numroc( n, desca( mb_ ), myrow, 0, nprow )
232  nq = numroc( n, desca( nb_ ), mycol, 0, npcol )
233  mq = numroc( n, desca( nb_ ), mycol, 0, npcol )
234 *
235 * Find the amount of workspace needed with or without eigenvectors.
236 *
237  trilwmin = 3*n + max( desca( nb_ )*( np+1 ), 3*desca( nb_ ) )
238  minsize = max( 1 + 6*n + 2*np*nq, trilwmin ) + 2*n
239 *
240  CALL dlacpy( 'A', np, nq, copya, desca( lld_ ), a( 1+iprepad ),
241  \$ desca( lld_ ) )
242 *
243  CALL pdfillpad( desca( ctxt_ ), np, nq, a, desca( lld_ ), iprepad,
245 *
246  CALL pdfillpad( descz( ctxt_ ), np, mq, z, descz( lld_ ), iprepad,
248 *
251 *
254 *
255 * Make sure that PDSYEVD does not cheat (i.e. use answers
257 *
258  DO 60 i = 1, n, 1
259  DO 50 j = 1, n, 1
260  CALL pdelset( z( 1+iprepad ), i, j, desca, 13.0d+0 )
261  50 CONTINUE
262  60 CONTINUE
263 *
264  CALL slboot
265  CALL sltimer( 1 )
266  CALL sltimer( 6 )
267  CALL pdsyevd( 'V', uplo, n, a( 1+iprepad ), ia, ja, desca,
270  \$ liwork, info )
271  CALL sltimer( 6 )
272  CALL sltimer( 1 )
273 *
274  IF( thresh.LE.0 ) THEN
275  result = 0
276  ELSE
277  CALL pdchekpad( desca( ctxt_ ), 'PDSYEVD-A', np, nq, a,
279 *
280  CALL pdchekpad( descz( ctxt_ ), 'PDSYEVD-Z', np, mq, z,
283 *
284  CALL pdchekpad( desca( ctxt_ ), 'PDSYEVD-WNEW', n, 1, wnew, n,
286 *
287  CALL pdchekpad( desca( ctxt_ ), 'PDSYEVD-WORK', lwork1, 1,
290 *
291 * Check INFO
292 *
293 *
294 * Make sure that all processes return the same value of INFO
295 *
296  itmp( 1 ) = info
297  itmp( 2 ) = info
298 *
299  CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
300  \$ -1, -1, 0 )
301  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1, 1,
302  \$ 1, -1, -1, 0 )
303 *
304 *
305  IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
306  IF( iam.EQ.0 )
307  \$ WRITE( nout, fmt = * )
308  \$ 'Different processes return different INFO'
309  result = 1
310  ELSE IF( info.NE.0 ) THEN
311  IF( iam.EQ.0 ) THEN
312  WRITE( nout, fmt = 9999 )info
313  IF( info.EQ.(n+1) )
314  \$ WRITE( nout, fmt = 9994 )
315  result = 1
316  END IF
317  ELSE IF( info.EQ.14 .AND. lwork1.GE.minsize ) THEN
318  IF( iam.EQ.0 )
319  \$ WRITE( nout, fmt = 9996 )info
320  result = 1
321  END IF
322 *
323  IF( result.EQ.0 .OR. info.GT.n ) THEN
324 *
325 * Make sure that different processes return the same eigenvalues.
326 * This is a more exhaustive check that provided by PDSYEVD.
327 *
328  DO 70 i = 1, n
329  work( i ) = wnew( i+iprepad )
330  work( i+n ) = wnew( i+iprepad )
331  70 CONTINUE
332 *
333  CALL dgamn2d( desca( ctxt_ ), 'a', ' ', n, 1, work, n, 1,
334  \$ 1, -1, -1, 0 )
335  CALL dgamx2d( desca( ctxt_ ), 'a', ' ', n, 1,
336  \$ work( 1+n ), n, 1, 1, -1, -1, 0 )
337 *
338  DO 80 i = 1, n
339 *
340  IF( abs( work( i )-work( n+i ) ).GT.zero ) THEN
341  IF( iam.EQ.0 )
342  \$ WRITE( nout, fmt = 9995 )
343  result = 1
344  GO TO 90
345  END IF
346  80 CONTINUE
347  90 CONTINUE
348  END IF
349 *
350  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1,
351  \$ -1, -1, 0 )
352 *
353 * Compute eps * norm(A)
354 *
355  IF( n.EQ.0 ) THEN
356  epsnorma = eps
357  ELSE
358  epsnorma = pdlansy( 'I', uplo, n, copya, ia, ja, desca,
359  \$ work )*eps
360  END IF
361 *
362 * Note that a couple key variables get redefined in PDSEPCHK
363 * as described by this table:
364 *
365 * PDSEPTST name PDSEPCHK name
366 * ------------- -------------
367 * COPYA A
368 * Z Q
369 * A C
370 *
371 *
372 *
373 * Perform the |AQ - QE| test
374 *
375  CALL pdfillpad( desca( ctxt_ ), sizechk, 1, work, sizechk,
377 *
378  resaq = 0
379 *
380  CALL pdsepchk( n, n, copya, ia, ja, desca,
381  \$ max( abstol+epsnorma, safmin ), thresh,
382  \$ z( 1+iprepad ), ia, ja, descz,
383  \$ a( 1+iprepad ), ia, ja, desca,
385  \$ sizechk, tstnrm, resaq )
386 *
387  CALL pdchekpad( desca( ctxt_ ), 'PDSEPCHK-WORK', sizechk, 1,
389 *
390  IF( resaq.NE.0 ) THEN
391  result = 1
392  WRITE( nout, fmt = 9993 )
393  END IF
394 *
395 * Perform the |QTQ - I| test
396 *
397  CALL pdfillpad( desca( ctxt_ ), sizeqtq, 1, work, sizeqtq,
399 *
400  resqtq = 0
401 *
402 *
403  DO 40 i = 1, 2
404  iwork( iprepad + i ) = 0
405  40 CONTINUE
406  CALL pdsepqtq( n, n, thresh, z( 1+iprepad ), ia, ja, descz,
407  \$ a( 1+iprepad ), ia, ja, desca,
408  \$ iwork( 1 ), iwork( 1 ), work( 1 ),
409  \$ work( iprepad+1 ), sizeqtq, qtqnrm, info,
410  \$ resqtq )
411 *
412  CALL pdchekpad( desca( ctxt_ ), 'PDSEPQTQ-WORK', sizeqtq, 1,
414 *
415  IF( resqtq.NE.0 ) THEN
416  result = 1
417  WRITE( nout, fmt = 9992 )
418  END IF
419 *
420  IF( info.NE.0 ) THEN
421  IF( iam.EQ.0 )
422  \$ WRITE( nout, fmt = 9998 )info
423  result = 1
424  END IF
425  ENDIF
426 *
427 * Check to make sure that we have the right eigenvalues
428 *
429  IF( wknown .AND. n.GT.0 ) THEN
430 *
431 * Find the largest difference between the computed
432 * and expected eigenvalues
433 *
434  minerror = normwin
435  maxerror = 0
436 *
437  DO 140 i = 1, n
439  maxerror = max( maxerror, error )
440  140 CONTINUE
441  minerror = min( maxerror, minerror )
442 *
443  IF( minerror.GT.normwin*five*thresh*eps ) THEN
444  IF( iam.EQ.0 )
445  \$ WRITE( nout, fmt = 9997 )minerror, normwin
446  result = 1
447  END IF
448  END IF
449 *
450 * All processes should report the same result
451 *
452  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1, -1,
453  \$ -1, 0 )
454 *
455  150 CONTINUE
456 *
457 *
458  RETURN
459 *
460  9999 FORMAT( 'PDSYEVD returned INFO=', i7 )
461  9998 FORMAT( 'PDSEPQTQ in PDSDPSUBTST returned INFO=', i7 )
462  9997 FORMAT( 'PDSDPSUBTST minerror =', d11.2, ' normwin=', d11.2 )
463  9996 FORMAT( 'PDSYEVD returned INFO=', i7,
464  \$ ' despite adequate workspace' )
465  9995 FORMAT( 'Different processes return different eigenvalues' )
466  9994 FORMAT( 'Heterogeneity detected by PDSYEVD' )
467  9993 FORMAT( 'PDSYEVD failed the |AQ -QE| test' )
468  9992 FORMAT( 'PDSYEVD failed the |QTQ -I| test' )
469 *
470 * End of PDSDPSUBTST
471 *
472  END
473
474
475
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
pdsdpsubtst
subroutine pdsdpsubtst(WKNOWN, UPLO, N, THRESH, ABSTOL, A, COPYA, Z, IA, JA, DESCA, WIN, WNEW, IPREPAD, IPOSTPAD, WORK, LWORK, LWORK1, IWORK, LIWORK, RESULT, TSTNRM, QTQNRM, NOUT)
Definition: pdsdpsubtst.f:6
pdsyevd
subroutine pdsyevd(JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, DESCZ, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdsyevd.f:3
subroutine pdchekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
pdlasizesqp
subroutine pdlasizesqp(DESCA, IPREPAD, IPOSTPAD, SIZEMQRLEFT, SIZEMQRRIGHT, SIZEQRF, SIZETMS, SIZEQTQ, SIZECHK, SIZESYEVX, ISIZESYEVX, SIZESYEV, SIZESYEVD, ISIZESYEVD, SIZESUBTST, ISIZESUBTST, SIZETST, ISIZETST)
Definition: pdlasizesqp.f:7
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
slboot
subroutine slboot()
Definition: sltimer.f:2