ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
pcsdpsubtst.f
Go to the documentation of this file.
1  SUBROUTINE pcsdpsubtst( 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  REAL ABSTOL, QTQNRM, THRESH, TSTNRM
18 * ..
19 * .. Array Arguments ..
20  INTEGER DESCA( * ), IWORK( * )
21  REAL RWORK( * ), WIN( * ), WNEW( * )
22  COMPLEX A( * ), COPYA( * ), WORK( * ), Z( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * PCSDPSUBTST calls PCHEEVD and then tests the output of
29 * PCHEEVD
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 * PCHEEVD 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) REAL
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) REAL
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 ("PCSTEIN"),
74 * "abstol" MUST NOT be bigger than ulp*|T|.
75 *
76 * A (local workspace) COMPLEX 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 PCHEEVD for a description of block cyclic layout.
81 * The test matrix, which is then modified by PCHEEVD
83 *
84 * COPYA (local input) COMPLEX 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 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 * PCSEPCHK and PCSEPQTQ.
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) REAL 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) REAL array, dimension (N)
112 * The eigenvalues as copmuted by this call to PCHEEVD
113 * If JOBZ <> 'V' or RANGE <> 'A' these eigenvalues are
114 * compared against those in WIN().
117 *
118 * WORK (local workspace) COMPLEX array, dimension (LWORK)
121 *
122 * LWORK (local input) INTEGER
123 * The actual length of the array WORK after padding.
124 *
125 * RWORK (local workspace) REAL 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 real workspace to pass to PCHEEVD
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 PCHEEVD
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) REAL
149 * |AQ- QL| / (ABSTOL+EPS*|A|)*N
150 *
151 * QTQNRM (global output) REAL
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 )
162  parameter( padval = 13.5285e+0, five = 5.0e+0,
163  \$ negone = -1.0e+0 )
165  parameter( cpadval = ( 13.989e+0, 1.93e+0 ) )
167  PARAMETER ( IPADVAL = 927 )
168  COMPLEX CZERO, CONE, CNEGONE
169  parameter( czero = 0.0e+0, cone = 1.0e+0,
170  \$ cnegone = -1.0e+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  REAL 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  REAL PCLANGE, PCLANHE, PSLAMCH
190  EXTERNAL NUMROC, PCLANGE, PCLANHE, PSLAMCH
191 * ..
192 * .. External Subroutines ..
193  EXTERNAL blacs_gridinfo, clacpy, igamn2d, igamx2d,
197 * ..
198 * .. Intrinsic Functions ..
199  INTRINSIC abs, max, min, real
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 = pslamch( desca( ctxt_ ), 'Eps' )
216  safmin = pslamch( 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.3e+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.63e+0, 1.1e+0 )
231  30 CONTINUE
232 *
233  DO 40 i = 1, n
234  wnew( i+iprepad ) = 3.14159e+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 clacpy( 'A', np, nq, copya, desca( lld_ ), a( 1+iprepad ),
254  \$ desca( lld_ ) )
255 *
256  CALL pcfillpad( desca( ctxt_ ), np, nq, a, desca( lld_ ), iprepad,
258 *
259  CALL pcfillpad( 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 pcheevd( '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 pcchekpad( desca( ctxt_ ), 'PCHEEVD-A', np, nq, a,
290 *
291  CALL pcchekpad( desca( ctxt_ ), 'PCHEEVD-Z', np, nq, z,
294 *
295  CALL pschekpad( desca( ctxt_ ), 'PCHEEVD-WNEW', n, 1, wnew, n,
297 *
298  CALL pschekpad( desca( ctxt_ ), 'PCHEEVD-rWORK', lwork1, 1,
301 *
302  CALL pcchekpad( desca( ctxt_ ), 'PCHEEVD-WORK', lwork, 1, work,
304 *
305  CALL pichekpad( desca( ctxt_ ), 'PCHEEVD-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 = pclanhe( 'I', uplo, n, copya, ia, ja, desca,
338  \$ rwork )*eps
339  END IF
340 *
341 * Note that a couple key variables get redefined in PCSEPCHK
342 * as described by this table:
343 *
344 * PCSEPTST name PCSEPCHK name
345 * ------------- -------------
346 * COPYA A
347 * Z Q
348 * A C
349 *
350 * Perform the |AQ - QE| test
351 *
352  CALL psfillpad( desca( ctxt_ ), rsizechk, 1, rwork, rsizechk,
354 *
355  CALL pcsepchk( 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 pschekpad( desca( ctxt_ ), 'PCSDPCHK-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 psfillpad( desca( ctxt_ ), rsizeqtq, 1, rwork, rsizeqtq,
373 *
374 *
375  res = 0
376  ulp = pslamch( desca( ctxt_ ), 'P' )
377  CALL pclaset( 'A', n, n, czero, cone, a( 1+iprepad ), ia, ja,
378  \$ desca )
379  CALL pcgemm( '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 = pclange( '1', n, n, a( 1+iprepad ), ia, ja, desca,
383  \$ work( 1+iprepad ) )
384  qtqnrm = norm / ( real( max( n, 1 ) )*ulp )
385  IF( qtqnrm.GT.thresh ) THEN
386  res = 1
387  END IF
388  CALL pschekpad( desca( ctxt_ ), 'PCSEPQTQ-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.0
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( 'PCHEEVD returned INFO=', i7 )
443  9998 FORMAT( 'PCSEPQTQ returned INFO=', i7 )
444  9997 FORMAT( 'PCSDPSUBTST minerror =', d11.2, ' normwin=', d11.2 )
445  9996 FORMAT( 'PCHEEVD returned INFO=', i7,
446  \$ ' despite adequate workspace' )
447  9995 FORMAT( 'PCHEEVD failed the |AQ -QE| test' )
448  9994 FORMAT( 'PCHEEVD failed the |QTQ -I| test' )
449 *
450 * End of PCSDPSUBTST
451 *
452  END
subroutine pichekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
max
#define max(A, B)
Definition: pcgemr.c:180
pcheevd
subroutine pcheevd(JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, DESCZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
Definition: pcheevd.f:4
pcsdpsubtst
subroutine pcsdpsubtst(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: pcsdpsubtst.f:6
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
pclasizesep
subroutine pclasizesep(DESCA, IPREPAD, IPOSTPAD, SIZEMQRLEFT, SIZEMQRRIGHT, SIZEQRF, SIZETMS, RSIZEQTQ, RSIZECHK, SIZEHEEVX, RSIZEHEEVX, ISIZEHEEVX, SIZEHEEVD, RSIZEHEEVD, ISIZEHEEVD, SIZESUBTST, RSIZESUBTST, ISIZESUBTST, SIZETST, RSIZETST, ISIZETST)
Definition: pclasizesep.f:7
pcsepchk
subroutine pcsepchk(MS, NV, A, IA, JA, DESCA, EPSNORMA, THRESH, Q, IQ, JQ, DESCQ, C, IC, JC, DESCC, W, WORK, LWORK, TSTNRM, RESULT)
Definition: pcsepchk.f:6
subroutine pschekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
subroutine pcchekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
pclaset
subroutine pclaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pcblastst.f:7508
slboot
subroutine slboot()
Definition: sltimer.f:2