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
15  INTEGER IA, IPOSTPAD, IPREPAD, JA, LIWORK, LRWORK,
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
82 * A has already been padded front and back, use A(1+IPREPAD)
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.
93 * Z has already been padded front and back, use Z(1+IPREPAD)
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().
115 * WNEW has already been padded front and back,
116 * use WNEW(1+IPREPAD)
117 *
118 * WORK (local workspace) COMPLEX*16 array, dimension (LWORK)
119 * WORK has already been padded front and back,
120 * use WORK(1+IPREPAD)
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)
126 * RWORK has already been padded front and back,
127 * use RWORK(1+IPREPAD)
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)
136 * IWORK has already been padded front and back,
137 * use IWORK(1+IPREPAD)
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 )
164  COMPLEX*16 CPADVAL
165  parameter( cpadval = ( 13.989d+0, 1.93d+0 ) )
166  INTEGER IPADVAL
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,
194  $ pzchekpad, pzfillpad, pzgemm, pzheevd, pzlaset,
195  $ pzlasizesep, pzsepchk, pichekpad, pifillpad,
196  $ pdchekpad, pdfillpad, slboot, sltimer
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 *
206  CALL pzlasizesep( desca, iprepad, ipostpad, sizemqrleft,
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,
257  $ ipostpad, cpadval )
258 *
259  CALL pzfillpad( desca( ctxt_ ), np, nq, z, desca( lld_ ), iprepad,
260  $ ipostpad, cpadval+1.0d+0 )
261 *
262  CALL pdfillpad( desca( ctxt_ ), n, 1, wnew, n, iprepad, ipostpad,
263  $ padval+2.0d+0 )
264 *
265  CALL pdfillpad( desca( ctxt_ ), lwork1, 1, rwork, lwork1, iprepad,
266  $ ipostpad, padval+4.0d+0 )
267 *
268  CALL pifillpad( desca( ctxt_ ), liwork, 1, iwork, liwork, iprepad,
269  $ ipostpad, ipadval )
270 *
271  CALL pzfillpad( desca( ctxt_ ), lwork, 1, work, lwork, iprepad,
272  $ ipostpad, cpadval+4.1d+0 )
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,
279  $ wnew( 1+iprepad ), z( 1+iprepad ), ia, ja, desca,
280  $ work( 1+iprepad ), sizeheevd, rwork( 1+iprepad ),
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,
289  $ desca( lld_ ), iprepad, ipostpad, cpadval )
290 *
291  CALL pzchekpad( desca( ctxt_ ), 'PZHEEVD-Z', np, nq, z,
292  $ desca( lld_ ), iprepad, ipostpad,
293  $ cpadval+1.0d+0 )
294 *
295  CALL pdchekpad( desca( ctxt_ ), 'PZHEEVD-WNEW', n, 1, wnew, n,
296  $ iprepad, ipostpad, padval+2.0d+0 )
297 *
298  CALL pdchekpad( desca( ctxt_ ), 'PZHEEVD-rWORK', lwork1, 1,
299  $ rwork, lwork1, iprepad, ipostpad,
300  $ padval+4.0d+0 )
301 *
302  CALL pzchekpad( desca( ctxt_ ), 'PZHEEVD-WORK', lwork, 1, work,
303  $ lwork, iprepad, ipostpad, cpadval+4.1d+0 )
304 *
305  CALL pichekpad( desca( ctxt_ ), 'PZHEEVD-IWORK', liwork, 1,
306  $ iwork, liwork, iprepad, ipostpad, ipadval )
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,
353  $ iprepad, ipostpad, 4.3d+0 )
354 *
355  CALL pzsepchk( n, n, copya, ia, ja, desca,
356  $ max( abstol+epsnorma, safmin ), thresh,
357  $ z( 1+iprepad ), ia, ja, desca, a( 1+iprepad ),
358  $ ia, ja, desca, wnew( 1+iprepad ),
359  $ rwork( 1+iprepad ), rsizechk, tstnrm, res )
360 *
361  CALL pdchekpad( desca( ctxt_ ), 'PZSDPCHK-rWORK', rsizechk, 1,
362  $ rwork, rsizechk, iprepad, ipostpad, 4.3d+0 )
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,
372  $ iprepad, ipostpad, 4.3d+0 )
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,
389  $ rwork, rsizeqtq, iprepad, ipostpad, 4.3d+0 )
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
420  error = abs( win( i+iprepad )-wnew( i+iprepad ) )
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
pichekpad
subroutine pichekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pichekpad.f:3
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
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
pdfillpad
subroutine pdfillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pdfillpad.f:2
pzchekpad
subroutine pzchekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pzchekpad.f:3
pifillpad
subroutine pifillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pifillpad.f:2
min
#define min(A, B)
Definition: pcgemr.c:181
pzfillpad
subroutine pzfillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pzfillpad.f:2