ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
pdsqpsubtst.f
Go to the documentation of this file.
1 *
2 *
3  SUBROUTINE pdsqpsubtst( WKNOWN, JOBZ, UPLO, N, THRESH, ABSTOL, A,
4  \$ COPYA, Z, IA, JA, DESCA, WIN, WNEW,
6  \$ RESULT, TSTNRM, QTQNRM, NOUT )
7 *
8 * -- ScaLAPACK routine (version 1.7) --
9 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
10 * and University of California, Berkeley.
11 * November 15, 1997
12 *
13 * .. Scalar Arguments ..
14  LOGICAL WKNOWN
15  CHARACTER JOBZ, UPLO
17  \$ nout, result
18  DOUBLE PRECISION ABSTOL, QTQNRM, THRESH, TSTNRM
19 * ..
20 * .. Array Arguments ..
21  INTEGER DESCA( * )
22  DOUBLE PRECISION A( * ), COPYA( * ), WIN( * ), WNEW( * ),
23  \$ WORK( * ), Z( * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * PDSQPSUBTST calls PDSYEV and then tests the output of
30 * PDSYEV
31 * If JOBZ = 'V' then the following two tests are performed:
32 * |AQ -QL| / (abstol + eps * norm(A) ) < N*THRESH
33 * |QT * Q - I| / eps * norm(A) < N*THRESH
34 * If WKNOWN then
35 * we check to make sure that the eigenvalues match expectations
36 * i.e. |WIN - WNEW(1+IPREPAD)| / (eps * |WIN|) < THRESH
37 * where WIN is the array of eigenvalues as computed by
38 * PDSYEV when eigenvectors are requested
39 *
40 * Arguments
41 * =========
42 *
43 * NP = the number of rows local to a given process.
44 * NQ = the number of columns local to a given process.
45 *
46 * WKNOWN (global input) INTEGER
47 * .FALSE.: WIN does not contain the eigenvalues
48 * .TRUE.: WIN does contain the eigenvalues
49 *
50 * JOBZ (global input) CHARACTER*1
51 * Specifies whether or not to compute the eigenvectors:
52 * = 'N': Compute eigenvalues only.
53 * = 'V': Compute eigenvalues and eigenvectors.
54 * Must be 'V' on first call to PDSQPSUBTST
55 *
56 * UPLO (global input) CHARACTER*1
57 * Specifies whether the upper or lower triangular part of the
58 * symmetric matrix A is stored:
59 * = 'U': Upper triangular
60 * = 'L': Lower triangular
61 *
62 * N (global input) INTEGER
63 * Size of the matrix to be tested. (global size)
64 *
65 * THRESH (global input) DOUBLE PRECISION
66 * A test will count as "failed" if the "error", computed as
67 * described below, exceeds THRESH. Note that the error
68 * is scaled to be O(1), so THRESH should be a reasonably
69 * small multiple of 1, e.g., 10 or 100. In particular,
70 * it should not depend on the precision (single vs. double)
71 * or the size of the matrix. It must be at least zero.
72 *
73 * ABSTOL (global input) DOUBLE PRECISION
74 * The absolute tolerance for the eigenvalues. An
75 * eigenvalue is considered to be located if it has
76 * been determined to lie in an interval whose width
77 * is "abstol" or less. If "abstol" is less than or equal
78 * to zero, then ulp*|T| will be used, where |T| is
79 * the 1-norm of the matrix.
80 *
81 * A (local workspace) DOUBLE PRECISION array
82 * global dimension (N, N), local dimension (DESCA(DLEN_), NQ)
83 * A is distributed in a block cyclic manner over both rows
84 * and columns.
85 * See PDSYEV for a description of block cyclic layout.
86 * The test matrix, which is then modified by PDSYEV
88 *
89 * COPYA (local input) DOUBLE PRECISION array, dimension(N*N)
90 * COPYA holds a copy of the original matrix A
91 * identical in both form and content to A
92 *
93 * Z (local workspace) DOUBLE PRECISION array, dim (N*N)
94 * Z is distributed in the same manner as A
95 * Z contains the eigenvector matrix
96 * Z is used as workspace by the test routines
97 * PDSEPCHK and PDSEPQTQ.
99 *
100 * IA (global input) INTEGER
101 * On entry, IA specifies the global row index of the submatrix
102 * of the global matrix A, COPYA and Z to operate on.
103 *
104 * JA (global input) INTEGER
105 * On entry, IA specifies the global column index of the submat
106 * of the global matrix A, COPYA and Z to operate on.
107 *
108 * DESCA (global/local input) INTEGER array of dimension 8
109 * The array descriptor for the matrix A, COPYA and Z.
110 *
111 * WIN (global input) DOUBLE PRECISION array, dimension (N)
112 * If .not. WKNOWN, WIN is ignored on input
113 * Otherwise, WIN() is taken as the standard by which the
114 * eigenvalues are to be compared against.
115 *
116 * WNEW (global workspace) DOUBLE PRECISION array, dimension (N)
117 * The eigenvalues as computed by this call to PDSYEV.
120 *
121 * WORK (local workspace) DOUBLE PRECISION array, dimension (LWORK)
124 *
125 * LWORK (local input) INTEGER
126 * The actual length of the array WORK after padding.
127 *
128 *
129 * LWORK1 (local input) INTEGER
130 * The amount of real workspace to pass to PDSYEV
131 *
132 * RESULT (global output) INTEGER
133 * The result of this call to PDSYEV
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, EIGS, 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  DOUBLE PRECISION EPS, EPSNORMA, ERROR, MAXERROR, MINERROR,
163  \$ NORMWIN, SAFMIN
164 * ..
165 * .. Local Arrays ..
166  INTEGER DESCZ( DLEN_ ), ITMP( 2 ),
167  \$ IWORK( 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, dgamn2d, dgamx2d,
178  \$ dlacpy, igamn2d, igamx2d, pdchekpad, pdelset,
180  \$ pdsyev, 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  DO 40 i = 1, 2
216  iwork( i ) = 0
217  40 CONTINUE
218 *
219  IF( lsame( jobz, 'N' ) ) THEN
220  eigs = 0
221  ELSE
222  eigs = n
223  END IF
224 *
225 *
226  CALL descinit( descz, desca( m_ ), desca( n_ ), desca( mb_ ),
227  \$ desca( nb_ ), desca( rsrc_ ), desca( csrc_ ),
228  \$ desca( ctxt_ ), desca( lld_ ), info )
229 *
230  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
231 *
232  iam = 1
233  IF( myrow.EQ.0 .AND. mycol.EQ.0 )
234  \$ iam = 0
235 *
236 * If this process is not involved in this test, bail out now
237 *
238  result = -3
239  IF( myrow.GE.nprow .OR. myrow.LT.0 )
240  \$ GO TO 150
241  result = 0
242 *
243  np = numroc( n, desca( mb_ ), myrow, 0, nprow )
244  nq = numroc( n, desca( nb_ ), mycol, 0, npcol )
245  mq = numroc( eigs, desca( nb_ ), mycol, 0, npcol )
246 *
247 * Find the amount of workspace needed with or without eigenvectors.
248 *
249  CALL pdlasizesyev( jobz, n, desca, minsize )
250 *
251  CALL dlacpy( 'A', np, nq, copya, desca( lld_ ), a( 1+iprepad ),
252  \$ desca( lld_ ) )
253 *
254  CALL pdfillpad( desca( ctxt_ ), np, nq, a, desca( lld_ ), iprepad,
256 *
257  CALL pdfillpad( descz( ctxt_ ), np, mq, z, descz( lld_ ), iprepad,
259 *
262 *
265 *
266 * Make sure that PDSYEV does not cheat (i.e. use answers
268 *
269  DO 60 i = 1, n, 1
270  DO 50 j = 1, eigs, 1
271  CALL pdelset( z( 1+iprepad ), i, j, desca, 13.0d+0 )
272  50 CONTINUE
273  60 CONTINUE
274 *
275  CALL slboot
276  CALL sltimer( 1 )
277  CALL sltimer( 6 )
278  CALL pdsyev( jobz, uplo, n, a( 1+iprepad ), ia, ja, desca,
280  \$ work( 1+iprepad ), lwork1, info )
281  CALL sltimer( 6 )
282  CALL sltimer( 1 )
283 *
284  IF( thresh.LE.0 ) THEN
285  result = 0
286  ELSE
287  CALL pdchekpad( desca( ctxt_ ), 'PDSYEV-A', np, nq, a,
289 *
290  CALL pdchekpad( descz( ctxt_ ), 'PDSYEV-Z', np, mq, z,
293 *
294  CALL pdchekpad( desca( ctxt_ ), 'PDSYEV-WNEW', n, 1, wnew, n,
296 *
297  CALL pdchekpad( desca( ctxt_ ), 'PDSYEV-WORK', lwork1, 1,
300 *
301 * Check INFO
302 *
303 *
304 * Make sure that all processes return the same value of INFO
305 *
306  itmp( 1 ) = info
307  itmp( 2 ) = info
308 *
309  CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
310  \$ -1, -1, 0 )
311  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1, 1,
312  \$ 1, -1, -1, 0 )
313 *
314 *
315  IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
316  IF( iam.EQ.0 )
317  \$ WRITE( nout, fmt = * )
318  \$ 'Different processes return different INFO'
319  result = 1
320  ELSE IF( info.NE.0 ) THEN
321  IF( iam.EQ.0 ) THEN
322  WRITE( nout, fmt = 9999 )info
323  IF( info.EQ.(n+1) )
324  \$ WRITE( nout, fmt = 9994 )
325  result = 1
326  END IF
327  ELSE IF( info.EQ.14 .AND. lwork1.GE.minsize ) THEN
328  IF( iam.EQ.0 )
329  \$ WRITE( nout, fmt = 9996 )info
330  result = 1
331  END IF
332 *
333  IF( result.EQ.0 .OR. info.GT.n ) THEN
334 *
335 * Make sure that different processes return the same eigenvalues.
336 * This is a more exhaustive check that provided by PDSYEV.
337 *
338  DO 70 i = 1, n
339  work( i ) = wnew( i+iprepad )
340  work( i+n ) = wnew( i+iprepad )
341  70 CONTINUE
342 *
343  CALL dgamn2d( desca( ctxt_ ), 'a', ' ', n, 1, work, n, 1,
344  \$ 1, -1, -1, 0 )
345  CALL dgamx2d( desca( ctxt_ ), 'a', ' ', n, 1,
346  \$ work( 1+n ), n, 1, 1, -1, -1, 0 )
347 *
348  DO 80 i = 1, n
349 *
350  IF( abs( work( i )-work( n+i ) ).GT.zero ) THEN
351  IF( iam.EQ.0 )
352  \$ WRITE( nout, fmt = 9995 )
353  result = 1
354  GO TO 90
355  END IF
356  80 CONTINUE
357  90 CONTINUE
358  END IF
359 *
360  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1,
361  \$ -1, -1, 0 )
362 *
363 * Compute eps * norm(A)
364 *
365  IF( n.EQ.0 ) THEN
366  epsnorma = eps
367  ELSE
368  epsnorma = pdlansy( 'I', uplo, n, copya, ia, ja, desca,
369  \$ work )*eps
370  END IF
371 *
372 * Note that a couple key variables get redefined in PDSEPCHK
373 * as described by this table:
374 *
375 * PDSEPTST name PDSEPCHK name
376 * ------------- -------------
377 * COPYA A
378 * Z Q
379 * A C
380 *
381 *
382  IF( lsame( jobz, 'V' ) ) THEN
383 *
384 * Perform the |AQ - QE| test
385 *
386  CALL pdfillpad( desca( ctxt_ ), sizechk, 1, work, sizechk,
388 *
389  resaq = 0
390 *
391  CALL pdsepchk( n, n, copya, ia, ja, desca,
392  \$ max( abstol+epsnorma, safmin ), thresh,
393  \$ z( 1+iprepad ), ia, ja, descz,
394  \$ a( 1+iprepad ), ia, ja, desca,
396  \$ sizechk, tstnrm, resaq )
397 *
398  CALL pdchekpad( desca( ctxt_ ), 'PDSEPCHK-WORK', sizechk, 1,
400 *
401  IF( resaq.NE.0 ) THEN
402  result = 1
403  WRITE( nout, fmt = 9993 )
404  END IF
405 *
406 * Perform the |QTQ - I| test
407 *
408  CALL pdfillpad( desca( ctxt_ ), sizeqtq, 1, work, sizeqtq,
410 *
411  resqtq = 0
412 *
413  CALL pdsepqtq( n, n, thresh, z( 1+iprepad ), ia, ja, descz,
414  \$ a( 1+iprepad ), ia, ja, desca,
415  \$ iwork( 1 ), iwork( 1 ), work( 1 ),
416  \$ work( iprepad+1 ), sizeqtq, qtqnrm, info,
417  \$ resqtq )
418 *
419  CALL pdchekpad( desca( ctxt_ ), 'PDSEPQTQ-WORK', sizeqtq, 1,
421 *
422  IF( resqtq.NE.0 ) THEN
423  result = 1
424  WRITE( nout, fmt = 9992 )
425  END IF
426 *
427  IF( info.NE.0 ) THEN
428  IF( iam.EQ.0 )
429  \$ WRITE( nout, fmt = 9998 )info
430  result = 1
431  END IF
432  END IF
433 *
434 * Check to make sure that we have the right eigenvalues
435 *
436  IF( wknown .AND. n.GT.0 ) THEN
437 *
438 * Find the largest difference between the computed
439 * and expected eigenvalues
440 *
441  minerror = normwin
442  maxerror = 0
443 *
444  DO 140 i = 1, n
446  maxerror = max( maxerror, error )
447  140 CONTINUE
448  minerror = min( maxerror, minerror )
449 *
450  IF( minerror.GT.normwin*five*thresh*eps ) THEN
451  IF( iam.EQ.0 )
452  \$ WRITE( nout, fmt = 9997 )minerror, normwin
453  result = 1
454  END IF
455  END IF
456  END IF
457 *
458 * All processes should report the same result
459 *
460  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1, -1,
461  \$ -1, 0 )
462 *
463  150 CONTINUE
464 *
465 *
466  RETURN
467 *
468  9999 FORMAT( 'PDSYEV returned INFO=', i7 )
469  9998 FORMAT( 'PDSEPQTQ in PDSQPSUBTST returned INFO=', i7 )
470  9997 FORMAT( 'PDSQPSUBTST minerror =', d11.2, ' normwin=', d11.2 )
471  9996 FORMAT( 'PDSYEV returned INFO=', i7,
472  \$ ' despite adequate workspace' )
473  9995 FORMAT( 'Different processes return different eigenvalues' )
474  9994 FORMAT( 'Heterogeneity detected by PDSYEV' )
475  9993 FORMAT( 'PDSYEV failed the |AQ -QE| test' )
476  9992 FORMAT( 'PDSYEV failed the |QTQ -I| test' )
477 *
478 * End of PDSQPSUBTST
479 *
480  END
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
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
pdsqpsubtst
subroutine pdsqpsubtst(WKNOWN, JOBZ, UPLO, N, THRESH, ABSTOL, A, COPYA, Z, IA, JA, DESCA, WIN, WNEW, IPREPAD, IPOSTPAD, WORK, LWORK, LWORK1, RESULT, TSTNRM, QTQNRM, NOUT)
Definition: pdsqpsubtst.f:7
pdsyev
subroutine pdsyev(JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, DESCZ, WORK, LWORK, INFO)
Definition: pdsyev.f:3
pdlasizesyev
subroutine pdlasizesyev(JOBZ, N, DESCA, MINSIZE)
Definition: pdlasizesyev.f:4
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
slboot
subroutine slboot()
Definition: sltimer.f:2