ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
pssqpsubtst.f
Go to the documentation of this file.
1 *
2 *
3  SUBROUTINE pssqpsubtst( 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  REAL ABSTOL, QTQNRM, THRESH, TSTNRM
19 * ..
20 * .. Array Arguments ..
21  INTEGER DESCA( * )
22  REAL A( * ), COPYA( * ), WIN( * ), WNEW( * ),
23  \$ WORK( * ), Z( * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * PSSQPSUBTST calls PSSYEV and then tests the output of
30 * PSSYEV
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 * PSSYEV 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 PSSQPSUBTST
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) REAL
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) REAL
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) REAL 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 PSSYEV for a description of block cyclic layout.
86 * The test matrix, which is then modified by PSSYEV
88 *
89 * COPYA (local input) REAL 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) REAL 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 * PSSEPCHK and PSSEPQTQ.
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) REAL 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) REAL array, dimension (N)
117 * The eigenvalues as computed by this call to PSSYEV.
120 *
121 * WORK (local workspace) REAL 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 PSSYEV
131 *
132 * RESULT (global output) INTEGER
133 * The result of this call to PSSYEV
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) REAL
139 * |AQ- QL| / (ABSTOL+EPS*|A|)*N
140 *
141 * QTQNRM (global output) REAL
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  REAL FIVE, NEGONE, PADVAL, ZERO
152  PARAMETER ( PADVAL = 13.5285e+0, five = 5.0e+0,
153  \$ negone = -1.0e+0, zero = 0.0e+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  REAL 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  REAL PSLAMCH, PSLANSY
174  EXTERNAL lsame, numroc, pslamch, pslansy
175 * ..
176 * .. External Subroutines ..
177  EXTERNAL blacs_gridinfo, descinit, igamn2d, igamx2d,
179  \$ pssepchk, pssepqtq, pssyev, sgamn2d, sgamx2d,
180  \$ slacpy, 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 = pslamch( desca( ctxt_ ), 'Eps' )
198  safmin = pslamch( 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.3e+0
209  10 CONTINUE
210 *
211  DO 30 i = 1, n
212  wnew( i+iprepad ) = 3.14159e+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 pslasizesyev( jobz, n, desca, minsize )
250 *
251  CALL slacpy( 'A', np, nq, copya, desca( lld_ ), a( 1+iprepad ),
252  \$ desca( lld_ ) )
253 *
254  CALL psfillpad( desca( ctxt_ ), np, nq, a, desca( lld_ ), iprepad,
256 *
257  CALL psfillpad( descz( ctxt_ ), np, mq, z, descz( lld_ ), iprepad,
259 *
262 *
265 *
266 * Make sure that PSSYEV does not cheat (i.e. use answers
268 *
269  DO 60 i = 1, n, 1
270  DO 50 j = 1, eigs, 1
271  CALL pselset( z( 1+iprepad ), i, j, desca, 13.0e+0 )
272  50 CONTINUE
273  60 CONTINUE
274 *
275  CALL slboot
276  CALL sltimer( 1 )
277  CALL sltimer( 6 )
278  CALL pssyev( 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 pschekpad( desca( ctxt_ ), 'PSSYEV-A', np, nq, a,
289 *
290  CALL pschekpad( descz( ctxt_ ), 'PSSYEV-Z', np, mq, z,
293 *
294  CALL pschekpad( desca( ctxt_ ), 'PSSYEV-WNEW', n, 1, wnew, n,
296 *
297  CALL pschekpad( desca( ctxt_ ), 'PSSYEV-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 PSSYEV.
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 sgamn2d( desca( ctxt_ ), 'a', ' ', n, 1, work, n, 1,
344  \$ 1, -1, -1, 0 )
345  CALL sgamx2d( 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 = pslansy( 'I', uplo, n, copya, ia, ja, desca,
369  \$ work )*eps
370  END IF
371 *
372 * Note that a couple key variables get redefined in PSSEPCHK
373 * as described by this table:
374 *
375 * PSSEPTST name PSSEPCHK 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 psfillpad( desca( ctxt_ ), sizechk, 1, work, sizechk,
388 *
389  resaq = 0
390 *
391  CALL pssepchk( 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 pschekpad( desca( ctxt_ ), 'PSSEPCHK-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 psfillpad( desca( ctxt_ ), sizeqtq, 1, work, sizeqtq,
410 *
411  resqtq = 0
412 *
413  CALL pssepqtq( 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 pschekpad( desca( ctxt_ ), 'PSSEPQTQ-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( 'PSSYEV returned INFO=', i7 )
469  9998 FORMAT( 'PSSEPQTQ in PSSQPSUBTST returned INFO=', i7 )
470  9997 FORMAT( 'PSSQPSUBTST minerror =', d11.2, ' normwin=', d11.2 )
471  9996 FORMAT( 'PSSYEV returned INFO=', i7,
472  \$ ' despite adequate workspace' )
473  9995 FORMAT( 'Different processes return different eigenvalues' )
474  9994 FORMAT( 'Heterogeneity detected by PSSYEV' )
475  9993 FORMAT( 'PSSYEV failed the |AQ -QE| test' )
476  9992 FORMAT( 'PSSYEV failed the |QTQ -I| test' )
477 *
478 * End of PSSQPSUBTST
479 *
480  END
max
#define max(A, B)
Definition: pcgemr.c:180
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
pslasizesyev
subroutine pslasizesyev(JOBZ, N, DESCA, MINSIZE)
Definition: pslasizesyev.f:4
pslasizesqp
subroutine pslasizesqp(DESCA, IPREPAD, IPOSTPAD, SIZEMQRLEFT, SIZEMQRRIGHT, SIZEQRF, SIZETMS, SIZEQTQ, SIZECHK, SIZESYEVX, ISIZESYEVX, SIZESYEV, SIZESYEVD, ISIZESYEVD, SIZESUBTST, ISIZESUBTST, SIZETST, ISIZETST)
Definition: pslasizesqp.f:7
pssepqtq
subroutine pssepqtq(MS, NV, THRESH, Q, IQ, JQ, DESCQ, C, IC, JC, DESCC, PROCDIST, ICLUSTR, GAP, WORK, LWORK, QTQNRM, INFO, RES)
Definition: pssepqtq.f:6
pselset
subroutine pselset(A, IA, JA, DESCA, ALPHA)
Definition: pselset.f:2
subroutine pschekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
pssqpsubtst
subroutine pssqpsubtst(WKNOWN, JOBZ, UPLO, N, THRESH, ABSTOL, A, COPYA, Z, IA, JA, DESCA, WIN, WNEW, IPREPAD, IPOSTPAD, WORK, LWORK, LWORK1, RESULT, TSTNRM, QTQNRM, NOUT)
Definition: pssqpsubtst.f:7
pssyev
subroutine pssyev(JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, DESCZ, WORK, LWORK, INFO)
Definition: pssyev.f:3
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
slboot
subroutine slboot()
Definition: sltimer.f:2
pssepchk
subroutine pssepchk(MS, NV, A, IA, JA, DESCA, EPSNORMA, THRESH, Q, IQ, JQ, DESCQ, C, IC, JC, DESCC, W, WORK, LWORK, TSTNRM, RESULT)
Definition: pssepchk.f:6