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,
3  $ IPREPAD, IPOSTPAD, WORK, LWORK, LWORK1,
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
15  INTEGER IA, IPOSTPAD, IPREPAD, JA, LWORK, LWORK1, N,
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
80 * A has already been padded front and back, use A(1+IPREPAD)
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.
91 * Z has already been padded front and back, use Z(1+IPREPAD)
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.
111 * WNEW has already been padded front and back,
112 * use WNEW(1+IPREPAD)
113 *
114 * WORK (local workspace) DOUBLE PRECISION array, dimension (LWORK)
115 * WORK has already been padded front and back,
116 * use WORK(1+IPREPAD)
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)
126 * IWORK has already been padded front and back,
127 * use IWORK(1+IPREPAD)
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,
178  $ pdchekpad, pdelset, pdfillpad, pdlasizesqp,
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
189  CALL pdlasizesqp( desca, iprepad, ipostpad, sizemqrleft,
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,
244  $ ipostpad, padval )
245 *
246  CALL pdfillpad( descz( ctxt_ ), np, mq, z, descz( lld_ ), iprepad,
247  $ ipostpad, padval+1.0d+0 )
248 *
249  CALL pdfillpad( desca( ctxt_ ), n, 1, wnew, n, iprepad, ipostpad,
250  $ padval+2.0d+0 )
251 *
252  CALL pdfillpad( desca( ctxt_ ), lwork1, 1, work, lwork1, iprepad,
253  $ ipostpad, padval+4.0d+0 )
254 *
255 * Make sure that PDSYEVD does not cheat (i.e. use answers
256 * already computed.)
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,
268  $ wnew( 1+iprepad ), z( 1+iprepad ), ia, ja, desca,
269  $ work( 1+iprepad ), lwork1, iwork( 1+iprepad ),
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,
278  $ desca( lld_ ), iprepad, ipostpad, padval )
279 *
280  CALL pdchekpad( descz( ctxt_ ), 'PDSYEVD-Z', np, mq, z,
281  $ descz( lld_ ), iprepad, ipostpad,
282  $ padval+1.0d+0 )
283 *
284  CALL pdchekpad( desca( ctxt_ ), 'PDSYEVD-WNEW', n, 1, wnew, n,
285  $ iprepad, ipostpad, padval+2.0d+0 )
286 *
287  CALL pdchekpad( desca( ctxt_ ), 'PDSYEVD-WORK', lwork1, 1,
288  $ work, lwork1, iprepad, ipostpad,
289  $ padval+4.0d+0 )
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,
376  $ iprepad, ipostpad, 4.3d+0 )
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,
384  $ wnew( 1+iprepad ), work( 1+iprepad ),
385  $ sizechk, tstnrm, resaq )
386 *
387  CALL pdchekpad( desca( ctxt_ ), 'PDSEPCHK-WORK', sizechk, 1,
388  $ work, sizechk, iprepad, ipostpad, 4.3d+0 )
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,
398  $ iprepad, ipostpad, 4.3d+0 )
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,
413  $ work, sizeqtq, iprepad, ipostpad, 4.3d+0 )
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
438  error = abs( win( i+iprepad )-wnew( i+iprepad ) )
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
pdchekpad
subroutine pdchekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pdchekpad.f:3
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
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
pdelset
subroutine pdelset(A, IA, JA, DESCA, ALPHA)
Definition: pdelset.f:2
min
#define min(A, B)
Definition: pcgemr.c:181