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,
5  $ IPREPAD, IPOSTPAD, WORK, LWORK, LWORK1,
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
16  INTEGER IA, IPOSTPAD, IPREPAD, JA, LWORK, LWORK1, N,
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
87 * A has already been padded front and back, use A(1+IPREPAD)
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.
98 * Z has already been padded front and back, use Z(1+IPREPAD)
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.
118 * WNEW has already been padded front and back,
119 * use WNEW(1+IPREPAD)
120 *
121 * WORK (local workspace) DOUBLE PRECISION array, dimension (LWORK)
122 * WORK has already been padded front and back,
123 * use WORK(1+IPREPAD)
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,
179  $ pdfillpad, pdlasizesqp, pdsepchk, pdsepqtq,
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
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  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,
255  $ ipostpad, padval )
256 *
257  CALL pdfillpad( descz( ctxt_ ), np, mq, z, descz( lld_ ), iprepad,
258  $ ipostpad, padval+1.0d+0 )
259 *
260  CALL pdfillpad( desca( ctxt_ ), n, 1, wnew, n, iprepad, ipostpad,
261  $ padval+2.0d+0 )
262 *
263  CALL pdfillpad( desca( ctxt_ ), lwork1, 1, work, lwork1, iprepad,
264  $ ipostpad, padval+4.0d+0 )
265 *
266 * Make sure that PDSYEV does not cheat (i.e. use answers
267 * already computed.)
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,
279  $ wnew( 1+iprepad ), z( 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,
288  $ desca( lld_ ), iprepad, ipostpad, padval )
289 *
290  CALL pdchekpad( descz( ctxt_ ), 'PDSYEV-Z', np, mq, z,
291  $ descz( lld_ ), iprepad, ipostpad,
292  $ padval+1.0d+0 )
293 *
294  CALL pdchekpad( desca( ctxt_ ), 'PDSYEV-WNEW', n, 1, wnew, n,
295  $ iprepad, ipostpad, padval+2.0d+0 )
296 *
297  CALL pdchekpad( desca( ctxt_ ), 'PDSYEV-WORK', lwork1, 1,
298  $ work, lwork1, iprepad, ipostpad,
299  $ padval+4.0d+0 )
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,
387  $ iprepad, ipostpad, 4.3d+0 )
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,
395  $ wnew( 1+iprepad ), work( 1+iprepad ),
396  $ sizechk, tstnrm, resaq )
397 *
398  CALL pdchekpad( desca( ctxt_ ), 'PDSEPCHK-WORK', sizechk, 1,
399  $ work, sizechk, iprepad, ipostpad, 4.3d+0 )
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,
409  $ iprepad, ipostpad, 4.3d+0 )
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,
420  $ work, sizeqtq, iprepad, ipostpad, 4.3d+0 )
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
445  error = abs( win( i+iprepad )-wnew( i+iprepad ) )
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
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
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
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