ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
pzsdpsubtst.f
Go to the documentation of this file.
1  SUBROUTINE pzsdpsubtst( WKNOWN, UPLO, N, THRESH, ABSTOL, A, COPYA,
2  \$ Z, IA, JA, DESCA, WIN, WNEW, IPREPAD,
3  \$ IPOSTPAD, WORK, LWORK, RWORK, LRWORK,
4  \$ LWORK1, IWORK, LIWORK, RESULT, TSTNRM,
5  \$ 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 * February 28, 2000
11 *
12 * .. Scalar Arguments ..
13  LOGICAL WKNOWN
14  CHARACTER UPLO
16  \$ LWORK, LWORK1, N, NOUT, RESULT
17  DOUBLE PRECISION ABSTOL, QTQNRM, THRESH, TSTNRM
18 * ..
19 * .. Array Arguments ..
20  INTEGER DESCA( * ), IWORK( * )
21  DOUBLE PRECISION RWORK( * ), WIN( * ), WNEW( * )
22  COMPLEX*16 A( * ), COPYA( * ), WORK( * ), Z( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * PZSDPSUBTST calls PZHEEVD and then tests the output of
29 * PZHEEVD
30 * If JOBZ = 'V' then 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 * PZHEEVD 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 * Hermitian 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. If eigenvectors are
73 * desired later by inverse iteration ("PZSTEIN"),
74 * "abstol" MUST NOT be bigger than ulp*|T|.
75 *
76 * A (local workspace) COMPLEX*16 array
77 * global dimension (N, N), local dimension (DESCA(DLEN_), NQ)
78 * A is distributed in a block cyclic manner over both rows
79 * and columns.
80 * See PZHEEVD for a description of block cyclic layout.
81 * The test matrix, which is then modified by PZHEEVD
83 *
84 * COPYA (local input) COMPLEX*16 array, dimension(N*N)
85 * COPYA holds a copy of the original matrix A
86 * identical in both form and content to A
87 *
88 * Z (local workspace) COMPLEX*16 array, dim (N*N)
89 * Z is distributed in the same manner as A
90 * Z contains the eigenvector matrix
91 * Z is used as workspace by the test routines
92 * PZSEPCHK and PZSEPQTQ.
94 *
95 * IA (global input) INTEGER
96 * On entry, IA specifies the global row index of the submatrix
97 * of the global matrix A, COPYA and Z to operate on.
98 *
99 * JA (global input) INTEGER
100 * On entry, IA specifies the global column index of the submat
101 * of the global matrix A, COPYA and Z to operate on.
102 *
103 * DESCA (global/local input) INTEGER array of dimension 8
104 * The array descriptor for the matrix A, COPYA and Z.
105 *
106 * WIN (global input) DOUBLE PRECISION array, dimension (N)
107 * If .not. WKNOWN, WIN is ignored on input
108 * Otherwise, WIN() is taken as the standard by which the
109 * eigenvalues are to be compared against.
110 *
111 * WNEW (global workspace) DOUBLE PRECISION array, dimension (N)
112 * The eigenvalues as copmuted by this call to PZHEEVD
113 * If JOBZ <> 'V' or RANGE <> 'A' these eigenvalues are
114 * compared against those in WIN().
117 *
118 * WORK (local workspace) COMPLEX*16 array, dimension (LWORK)
121 *
122 * LWORK (local input) INTEGER
123 * The actual length of the array WORK after padding.
124 *
125 * RWORK (local workspace) DOUBLE PRECISION array, dimension (LRWORK)
128 *
129 * LRWORK (local input) INTEGER
130 * The actual length of the array RWORK after padding.
131 *
132 * LWORK1 (local input) INTEGER
133 * The amount of DOUBLE PRECISION workspace to pass to PZHEEVD
134 *
135 * IWORK (local workspace) INTEGER array, dimension (LIWORK)
138 *
139 * LIWORK (local input) INTEGER
140 * The length of the array IWORK after padding.
141 *
142 * RESULT (global output) INTEGER
143 * The result of this call to PZHEEVD
144 * RESULT = -3 => This process did not participate
145 * RESULT = 0 => All tests passed
146 * RESULT = 1 => ONe or more tests failed
147 *
148 * TSTNRM (global output) DOUBLE PRECISION
149 * |AQ- QL| / (ABSTOL+EPS*|A|)*N
150 *
151 * QTQNRM (global output) DOUBLE PRECISION
152 * |QTQ -I| / N*EPS
153 *
154 * .. Parameters ..
155 *
156  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
157  \$ MB_, NB_, RSRC_, CSRC_, LLD_
158  PARAMETER ( BLOCK_CYCLIC_2D = 1, dlen_ = 9, dtype_ = 1,
159  \$ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
160  \$ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
161  DOUBLE PRECISION PADVAL, FIVE, NEGONE
162  parameter( padval = 13.5285d+0, five = 5.0d+0,
163  \$ negone = -1.0d+0 )
165  parameter( cpadval = ( 13.989d+0, 1.93d+0 ) )
167  PARAMETER ( IPADVAL = 927 )
168  COMPLEX*16 CZERO, CONE, CNEGONE
169  parameter( czero = 0.0d+0, cone = 1.0d+0,
170  \$ cnegone = -1.0d+0 )
171 * ..
172 * .. Local Scalars ..
173  INTEGER I, IAM, INFO, ISIZEHEEVD, ISIZEHEEVX,
174  \$ ISIZESUBTST, ISIZETST, MYCOL, MYROW, NP, NPCOL,
175  \$ NPROW, NQ, RES, RSIZECHK, RSIZEHEEVD,
176  \$ RSIZEHEEVX, RSIZEQTQ, RSIZESUBTST, RSIZETST,
177  \$ sizeheevd, sizeheevx, sizemqrleft,
178  \$ sizemqrright, sizeqrf, sizesubtst, sizetms,
179  \$ sizetst
180  DOUBLE PRECISION EPS, EPSNORMA, ERROR, MAXERROR, MINERROR, NORM,
181  \$ NORMWIN, SAFMIN, ULP
182 * ..
183 * .. Local Arrays ..
184  INTEGER ITMP( 2 )
185 * ..
186 * .. External Functions ..
187 *
188  INTEGER NUMROC
189  DOUBLE PRECISION PZLANGE, PZLANHE, PDLAMCH
190  EXTERNAL NUMROC, PZLANGE, PZLANHE, PDLAMCH
191 * ..
192 * .. External Subroutines ..
193  EXTERNAL blacs_gridinfo, zlacpy, igamn2d, igamx2d,
197 * ..
198 * .. Intrinsic Functions ..
199  INTRINSIC abs, max, min, dble
200 * ..
201 * .. Executable Statements ..
202 * This is just to keep ftnchek happy
203  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
204  \$ rsrc_.LT.0 )RETURN
205 *
207  \$ sizemqrright, sizeqrf, sizetms, rsizeqtq,
208  \$ rsizechk, sizeheevx, rsizeheevx, isizeheevx,
209  \$ sizeheevd, rsizeheevd, isizeheevd, sizesubtst,
210  \$ rsizesubtst, isizesubtst, sizetst, rsizetst,
211  \$ isizetst )
212 *
213  tstnrm = negone
214  qtqnrm = negone
215  eps = pdlamch( desca( ctxt_ ), 'Eps' )
216  safmin = pdlamch( desca( ctxt_ ), 'Safe min' )
217 *
218  normwin = safmin / eps
219  IF( n.GE.1 )
220  \$ normwin = max( abs( win( 1+iprepad ) ),
221  \$ abs( win( n+iprepad ) ), normwin )
222 *
223  DO 10 i = 1, lwork1, 1
224  rwork( i+iprepad ) = 14.3d+0
225  10 CONTINUE
226  DO 20 i = 1, liwork, 1
227  iwork( i ) = 14
228  20 CONTINUE
229  DO 30 i = 1, lwork, 1
230  work( i+iprepad ) = ( 15.63d+0, 1.1d+0 )
231  30 CONTINUE
232 *
233  DO 40 i = 1, n
234  wnew( i+iprepad ) = 3.14159d+0
235  40 CONTINUE
236 *
237  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
238 *
239  iam = 1
240  IF( myrow.EQ.0 .AND. mycol.EQ.0 )
241  \$ iam = 0
242 *
243 * If this process is not involved in this test, bail out now
244 *
245  result = -3
246  IF( myrow.GE.nprow .OR. myrow.LT.0 )
247  \$ GO TO 60
248  result = 0
249 *
250  np = numroc( n, desca( mb_ ), myrow, 0, nprow )
251  nq = numroc( n, desca( nb_ ), mycol, 0, npcol )
252 *
253  CALL zlacpy( 'A', np, nq, copya, desca( lld_ ), a( 1+iprepad ),
254  \$ desca( lld_ ) )
255 *
256  CALL pzfillpad( desca( ctxt_ ), np, nq, a, desca( lld_ ), iprepad,
258 *
259  CALL pzfillpad( desca( ctxt_ ), np, nq, z, desca( lld_ ), iprepad,
261 *
264 *
267 *
270 *
273 *
274  CALL slboot
275  CALL sltimer( 1 )
276  CALL sltimer( 6 )
277 *
278  CALL pzheevd( 'V', uplo, n, a( 1+iprepad ), ia, ja, desca,
281  \$ lwork1, iwork( 1+iprepad ), liwork, info )
282  CALL sltimer( 6 )
283  CALL sltimer( 1 )
284 *
285  IF( thresh.LE.0 ) THEN
286  result = 0
287  ELSE
288  CALL pzchekpad( desca( ctxt_ ), 'PZHEEVD-A', np, nq, a,
290 *
291  CALL pzchekpad( desca( ctxt_ ), 'PZHEEVD-Z', np, nq, z,
294 *
295  CALL pdchekpad( desca( ctxt_ ), 'PZHEEVD-WNEW', n, 1, wnew, n,
297 *
298  CALL pdchekpad( desca( ctxt_ ), 'PZHEEVD-rWORK', lwork1, 1,
301 *
302  CALL pzchekpad( desca( ctxt_ ), 'PZHEEVD-WORK', lwork, 1, work,
304 *
305  CALL pichekpad( desca( ctxt_ ), 'PZHEEVD-IWORK', liwork, 1,
307 *
308 * Check INFO
309 *
310 * Make sure that all processes return the same value of INFO
311 *
312  itmp( 1 ) = info
313  itmp( 2 ) = info
314 *
315  CALL igamn2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp, 1, 1, 1,
316  \$ -1, -1, 0 )
317  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, itmp( 2 ), 1, 1,
318  \$ 1, -1, -1, 0 )
319 *
320 *
321  IF( itmp( 1 ).NE.itmp( 2 ) ) THEN
322  IF( iam.EQ.0 )
323  \$ WRITE( nout, fmt = * )
324  \$ 'Different processes return different INFO'
325  result = 1
326  ELSE IF( info.NE.0 ) THEN
327  IF( iam.EQ.0 )
328  \$ WRITE( nout, fmt = 9996 )info
329  result = 1
330  END IF
331 *
332 * Compute eps * norm(A)
333 *
334  IF( n.EQ.0 ) THEN
335  epsnorma = eps
336  ELSE
337  epsnorma = pzlanhe( 'I', uplo, n, copya, ia, ja, desca,
338  \$ rwork )*eps
339  END IF
340 *
341 * Note that a couple key variables get redefined in PZSEPCHK
342 * as described by this table:
343 *
344 * PZSEPTST name PZSEPCHK name
345 * ------------- -------------
346 * COPYA A
347 * Z Q
348 * A C
349 *
350 * Perform the |AQ - QE| test
351 *
352  CALL pdfillpad( desca( ctxt_ ), rsizechk, 1, rwork, rsizechk,
354 *
355  CALL pzsepchk( n, n, copya, ia, ja, desca,
356  \$ max( abstol+epsnorma, safmin ), thresh,
358  \$ ia, ja, desca, wnew( 1+iprepad ),
359  \$ rwork( 1+iprepad ), rsizechk, tstnrm, res )
360 *
361  CALL pdchekpad( desca( ctxt_ ), 'PZSDPCHK-rWORK', rsizechk, 1,
363 *
364  IF( res.NE.0 ) THEN
365  result = 1
366  WRITE( nout, fmt = 9995 )
367  END IF
368 *
369 * Perform the |QTQ - I| test
370 *
371  CALL pdfillpad( desca( ctxt_ ), rsizeqtq, 1, rwork, rsizeqtq,
373 *
374 *
375  res = 0
376  ulp = pdlamch( desca( ctxt_ ), 'P' )
377  CALL pzlaset( 'A', n, n, czero, cone, a( 1+iprepad ), ia, ja,
378  \$ desca )
379  CALL pzgemm( 'Conjugate transpose', 'N', n, n, n, cnegone,
380  \$ z( 1+iprepad ), ia, ja, desca, z( 1+iprepad ), ia,
381  \$ ja, desca, cone, a( 1+iprepad ), ia, ja, desca )
382  norm = pzlange( '1', n, n, a( 1+iprepad ), ia, ja, desca,
383  \$ work( 1+iprepad ) )
384  qtqnrm = norm / ( dble( max( n, 1 ) )*ulp )
385  IF( qtqnrm.GT.thresh ) THEN
386  res = 1
387  END IF
388  CALL pdchekpad( desca( ctxt_ ), 'PZSEPQTQ-rWORK', rsizeqtq, 1,
390 *
391  IF( res.NE.0 ) THEN
392  result = 1
393  WRITE( nout, fmt = 9994 )
394  END IF
395 *
396  IF( info.NE.0 ) THEN
397  IF( iam.EQ.0 )
398  \$ WRITE( nout, fmt = 9998 )info
399  result = 1
400  END IF
401 *
402  IF( info.NE.0 ) THEN
403  IF( iam.EQ.0 )
404  \$ WRITE( nout, fmt = 9998 )info
405  result = 1
406  END IF
407  END IF
408 *
409 * Check to make sure that we have the right eigenvalues
410 *
411  IF( wknown .AND. n.GT.0 ) THEN
412 *
413 * Find the largest difference between the computed
414 * and expected eigenvalues
415 *
416  minerror = normwin
417  maxerror = 0.0d+00
418 *
419  DO 50 i = 1, n
421  maxerror = max( maxerror, error )
422  50 CONTINUE
423  minerror = min( maxerror, minerror )
424 *
425  IF( minerror.GT.normwin*five*thresh*eps ) THEN
426  IF( iam.EQ.0 )
427  \$ WRITE( nout, fmt = 9997 )minerror, normwin
428  result = 1
429  END IF
430  END IF
431 *
432 *
433 * All processes should report the same result
434 *
435  CALL igamx2d( desca( ctxt_ ), 'a', ' ', 1, 1, result, 1, 1, 1, -1,
436  \$ -1, 0 )
437 *
438  60 CONTINUE
439 *
440  RETURN
441 *
442  9999 FORMAT( 'PZHEEVD returned INFO=', i7 )
443  9998 FORMAT( 'PZSEPQTQ returned INFO=', i7 )
444  9997 FORMAT( 'PZSDPSUBTST minerror =', d11.2, ' normwin=', d11.2 )
445  9996 FORMAT( 'PZHEEVD returned INFO=', i7,
446  \$ ' despite adequate workspace' )
447  9995 FORMAT( 'PZHEEVD failed the |AQ -QE| test' )
448  9994 FORMAT( 'PZHEEVD failed the |QTQ -I| test' )
449 *
450 * End of PZSDPSUBTST
451 *
452  END
subroutine pichekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
max
#define max(A, B)
Definition: pcgemr.c:180
subroutine pdchekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
pzsepchk
subroutine pzsepchk(MS, NV, A, IA, JA, DESCA, EPSNORMA, THRESH, Q, IQ, JQ, DESCQ, C, IC, JC, DESCC, W, WORK, LWORK, TSTNRM, RESULT)
Definition: pzsepchk.f:6
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
pzheevd
subroutine pzheevd(JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, DESCZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
Definition: pzheevd.f:4
pzlaset
subroutine pzlaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pzblastst.f:7509
pzsdpsubtst
subroutine pzsdpsubtst(WKNOWN, UPLO, N, THRESH, ABSTOL, A, COPYA, Z, IA, JA, DESCA, WIN, WNEW, IPREPAD, IPOSTPAD, WORK, LWORK, RWORK, LRWORK, LWORK1, IWORK, LIWORK, RESULT, TSTNRM, QTQNRM, NOUT)
Definition: pzsdpsubtst.f:6
pzlasizesep
subroutine pzlasizesep(DESCA, IPREPAD, IPOSTPAD, SIZEMQRLEFT, SIZEMQRRIGHT, SIZEQRF, SIZETMS, RSIZEQTQ, RSIZECHK, SIZEHEEVX, RSIZEHEEVX, ISIZEHEEVX, SIZEHEEVD, RSIZEHEEVD, ISIZEHEEVD, SIZESUBTST, RSIZESUBTST, ISIZESUBTST, SIZETST, RSIZETST, ISIZETST)
Definition: pzlasizesep.f:7
slboot
subroutine slboot()
Definition: sltimer.f:2
subroutine pdfillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)