SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ pcsdpsubtst()

subroutine pcsdpsubtst ( logical  wknown,
character  uplo,
integer  n,
real  thresh,
real  abstol,
complex, dimension( * )  a,
complex, dimension( * )  copya,
complex, dimension( * )  z,
integer  ia,
integer  ja,
integer, dimension( * )  desca,
real, dimension( * )  win,
real, dimension( * )  wnew,
integer  iprepad,
integer  ipostpad,
complex, dimension( * )  work,
integer  lwork,
real, dimension( * )  rwork,
integer  lrwork,
integer  lwork1,
integer, dimension( * )  iwork,
integer  liwork,
integer  result,
real  tstnrm,
real  qtqnrm,
integer  nout 
)

Definition at line 1 of file pcsdpsubtst.f.

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 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
82* A has already been padded front and back, use A(1+IPREPAD)
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.
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) 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().
115* WNEW has already been padded front and back,
116* use WNEW(1+IPREPAD)
117*
118* WORK (local workspace) COMPLEX 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) REAL 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 real workspace to pass to PCHEEVD
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 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 )
161 REAL PADVAL, FIVE, NEGONE
162 parameter( padval = 13.5285e+0, five = 5.0e+0,
163 $ negone = -1.0e+0 )
164 COMPLEX CPADVAL
165 parameter( cpadval = ( 13.989e+0, 1.93e+0 ) )
166 INTEGER IPADVAL
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,
194 $ pcchekpad, pcfillpad, pcgemm, pcheevd, pclaset,
195 $ pclasizesep, pcsepchk, pichekpad, pifillpad,
196 $ pschekpad, psfillpad, slboot, sltimer
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*
206 CALL pclasizesep( 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 = 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,
257 $ ipostpad, cpadval )
258*
259 CALL pcfillpad( desca( ctxt_ ), np, nq, z, desca( lld_ ), iprepad,
260 $ ipostpad, cpadval+1.0e+0 )
261*
262 CALL psfillpad( desca( ctxt_ ), n, 1, wnew, n, iprepad, ipostpad,
263 $ padval+2.0e+0 )
264*
265 CALL psfillpad( desca( ctxt_ ), lwork1, 1, rwork, lwork1, iprepad,
266 $ ipostpad, padval+4.0e+0 )
267*
268 CALL pifillpad( desca( ctxt_ ), liwork, 1, iwork, liwork, iprepad,
269 $ ipostpad, ipadval )
270*
271 CALL pcfillpad( desca( ctxt_ ), lwork, 1, work, lwork, iprepad,
272 $ ipostpad, cpadval+4.1e+0 )
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,
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 pcchekpad( desca( ctxt_ ), 'PCHEEVD-A', np, nq, a,
289 $ desca( lld_ ), iprepad, ipostpad, cpadval )
290*
291 CALL pcchekpad( desca( ctxt_ ), 'PCHEEVD-Z', np, nq, z,
292 $ desca( lld_ ), iprepad, ipostpad,
293 $ cpadval+1.0e+0 )
294*
295 CALL pschekpad( desca( ctxt_ ), 'PCHEEVD-WNEW', n, 1, wnew, n,
296 $ iprepad, ipostpad, padval+2.0e+0 )
297*
298 CALL pschekpad( desca( ctxt_ ), 'PCHEEVD-rWORK', lwork1, 1,
299 $ rwork, lwork1, iprepad, ipostpad,
300 $ padval+4.0e+0 )
301*
302 CALL pcchekpad( desca( ctxt_ ), 'PCHEEVD-WORK', lwork, 1, work,
303 $ lwork, iprepad, ipostpad, cpadval+4.1e+0 )
304*
305 CALL pichekpad( desca( ctxt_ ), 'PCHEEVD-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 = 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,
353 $ iprepad, ipostpad, 4.3e+0 )
354*
355 CALL pcsepchk( 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 pschekpad( desca( ctxt_ ), 'PCSDPCHK-rWORK', rsizechk, 1,
362 $ rwork, rsizechk, iprepad, ipostpad, 4.3e+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 psfillpad( desca( ctxt_ ), rsizeqtq, 1, rwork, rsizeqtq,
372 $ iprepad, ipostpad, 4.3e+0 )
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,
389 $ rwork, rsizeqtq, iprepad, ipostpad, 4.3e+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.0
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( '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*
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition numroc.f:2
pclaset
subroutine pclaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
Definition pcblastst.f:7508
pslamch
real function pslamch(ictxt, cmach)
Definition pcblastst.f:7455
pcchekpad
subroutine pcchekpad(ictxt, mess, m, n, a, lda, ipre, ipost, chkval)
Definition pcchekpad.f:3
pcfillpad
subroutine pcfillpad(ictxt, m, n, a, lda, ipre, ipost, chkval)
Definition pcfillpad.f:2
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
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
pclange
real function pclange(norm, m, n, a, ia, ja, desca, work)
Definition pclange.f:3
pclanhe
real function pclanhe(norm, uplo, n, a, ia, ja, desca, work)
Definition pclanhe.f:3
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
pichekpad
subroutine pichekpad(ictxt, mess, m, n, a, lda, ipre, ipost, chkval)
Definition pichekpad.f:3
pifillpad
subroutine pifillpad(ictxt, m, n, a, lda, ipre, ipost, chkval)
Definition pifillpad.f:2
pschekpad
subroutine pschekpad(ictxt, mess, m, n, a, lda, ipre, ipost, chkval)
Definition pschekpad.f:3
psfillpad
subroutine psfillpad(ictxt, m, n, a, lda, ipre, ipost, chkval)
Definition psfillpad.f:2
slboot
subroutine slboot()
Definition sltimer.f:2
sltimer
subroutine sltimer(i)
Definition sltimer.f:47
Here is the call graph for this function:
Here is the caller graph for this function: