OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
zget38.f
Go to the documentation of this file.
1*> \brief \b ZGET38
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE ZGET38( RMAX, LMAX, NINFO, KNT, NIN )
12*
13* .. Scalar Arguments ..
14* INTEGER KNT, NIN
15* ..
16* .. Array Arguments ..
17* INTEGER LMAX( 3 ), NINFO( 3 )
18* DOUBLE PRECISION RMAX( 3 )
19* ..
20*
21*
22*> \par Purpose:
23* =============
24*>
25*> \verbatim
26*>
27*> ZGET38 tests ZTRSEN, a routine for estimating condition numbers of a
28*> cluster of eigenvalues and/or its associated right invariant subspace
29*>
30*> The test matrices are read from a file with logical unit number NIN.
31*> \endverbatim
32*
33* Arguments:
34* ==========
35*
36*> \param[out] RMAX
37*> \verbatim
38*> RMAX is DOUBLE PRECISION array, dimension (3)
39*> Values of the largest test ratios.
40*> RMAX(1) = largest residuals from ZHST01 or comparing
41*> different calls to ZTRSEN
42*> RMAX(2) = largest error in reciprocal condition
43*> numbers taking their conditioning into account
44*> RMAX(3) = largest error in reciprocal condition
45*> numbers not taking their conditioning into
46*> account (may be larger than RMAX(2))
47*> \endverbatim
48*>
49*> \param[out] LMAX
50*> \verbatim
51*> LMAX is INTEGER array, dimension (3)
52*> LMAX(i) is example number where largest test ratio
53*> RMAX(i) is achieved. Also:
54*> If ZGEHRD returns INFO nonzero on example i, LMAX(1)=i
55*> If ZHSEQR returns INFO nonzero on example i, LMAX(2)=i
56*> If ZTRSEN returns INFO nonzero on example i, LMAX(3)=i
57*> \endverbatim
58*>
59*> \param[out] NINFO
60*> \verbatim
61*> NINFO is INTEGER array, dimension (3)
62*> NINFO(1) = No. of times ZGEHRD returned INFO nonzero
63*> NINFO(2) = No. of times ZHSEQR returned INFO nonzero
64*> NINFO(3) = No. of times ZTRSEN returned INFO nonzero
65*> \endverbatim
66*>
67*> \param[out] KNT
68*> \verbatim
69*> KNT is INTEGER
70*> Total number of examples tested.
71*> \endverbatim
72*>
73*> \param[in] NIN
74*> \verbatim
75*> NIN is INTEGER
76*> Input logical unit number.
77*> \endverbatim
78*
79* Authors:
80* ========
81*
82*> \author Univ. of Tennessee
83*> \author Univ. of California Berkeley
84*> \author Univ. of Colorado Denver
85*> \author NAG Ltd.
86*
87*> \ingroup complex16_eig
88*
89* =====================================================================
90 SUBROUTINE zget38( RMAX, LMAX, NINFO, KNT, NIN )
91*
92* -- LAPACK test routine --
93* -- LAPACK is a software package provided by Univ. of Tennessee, --
94* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
95*
96* .. Scalar Arguments ..
97 INTEGER KNT, NIN
98* ..
99* .. Array Arguments ..
100 INTEGER LMAX( 3 ), NINFO( 3 )
101 DOUBLE PRECISION RMAX( 3 )
102* ..
103*
104* =====================================================================
105*
106* .. Parameters ..
107 INTEGER LDT, LWORK
108 parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
109 DOUBLE PRECISION ZERO, ONE, TWO
110 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
111 DOUBLE PRECISION EPSIN
112 parameter( epsin = 5.9605d-8 )
113 COMPLEX*16 CZERO
114 parameter( czero = ( 0.0d+0, 0.0d+0 ) )
115* ..
116* .. Local Scalars ..
117 INTEGER I, INFO, ISCL, ISRT, ITMP, J, KMIN, M, N, NDIM
118 DOUBLE PRECISION BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN,
119 $ SMLNUM, STMP, TNRM, TOL, TOLIN, V, VMAX, VMIN,
120 $ VMUL
121* ..
122* .. Local Arrays ..
123 LOGICAL SELECT( LDT )
124 INTEGER IPNT( LDT ), ISELEC( LDT )
125 DOUBLE PRECISION RESULT( 2 ), RWORK( LDT ), VAL( 3 ),
126 $ WSRT( LDT )
127 COMPLEX*16 Q( LDT, LDT ), QSAV( LDT, LDT ),
128 $ QTMP( LDT, LDT ), T( LDT, LDT ),
129 $ TMP( LDT, LDT ), TSAV( LDT, LDT ),
130 $ TSAV1( LDT, LDT ), TTMP( LDT, LDT ), W( LDT ),
131 $ WORK( LWORK ), WTMP( LDT )
132* ..
133* .. External Functions ..
134 DOUBLE PRECISION DLAMCH, ZLANGE
135 EXTERNAL dlamch, zlange
136* ..
137* .. External Subroutines ..
138 EXTERNAL dlabad, zdscal, zgehrd, zhseqr, zhst01, zlacpy,
139 $ ztrsen, zunghr
140* ..
141* .. Intrinsic Functions ..
142 INTRINSIC dble, dimag, max, sqrt
143* ..
144* .. Executable Statements ..
145*
146 eps = dlamch( 'P' )
147 smlnum = dlamch( 'S' ) / eps
148 bignum = one / smlnum
149 CALL dlabad( smlnum, bignum )
150*
151* EPSIN = 2**(-24) = precision to which input data computed
152*
153 eps = max( eps, epsin )
154 rmax( 1 ) = zero
155 rmax( 2 ) = zero
156 rmax( 3 ) = zero
157 lmax( 1 ) = 0
158 lmax( 2 ) = 0
159 lmax( 3 ) = 0
160 knt = 0
161 ninfo( 1 ) = 0
162 ninfo( 2 ) = 0
163 ninfo( 3 ) = 0
164 val( 1 ) = sqrt( smlnum )
165 val( 2 ) = one
166 val( 3 ) = sqrt( sqrt( bignum ) )
167*
168* Read input data until N=0. Assume input eigenvalues are sorted
169* lexicographically (increasing by real part, then decreasing by
170* imaginary part)
171*
172 10 CONTINUE
173 READ( nin, fmt = * )n, ndim, isrt
174 IF( n.EQ.0 )
175 $ RETURN
176 READ( nin, fmt = * )( iselec( i ), i = 1, ndim )
177 DO 20 i = 1, n
178 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
179 20 CONTINUE
180 READ( nin, fmt = * )sin, sepin
181*
182 tnrm = zlange( 'M', n, n, tmp, ldt, rwork )
183 DO 200 iscl = 1, 3
184*
185* Scale input matrix
186*
187 knt = knt + 1
188 CALL zlacpy( 'F', n, n, tmp, ldt, t, ldt )
189 vmul = val( iscl )
190 DO 30 i = 1, n
191 CALL zdscal( n, vmul, t( 1, i ), 1 )
192 30 CONTINUE
193 IF( tnrm.EQ.zero )
194 $ vmul = one
195 CALL zlacpy( 'F', n, n, t, ldt, tsav, ldt )
196*
197* Compute Schur form
198*
199 CALL zgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
200 $ info )
201 IF( info.NE.0 ) THEN
202 lmax( 1 ) = knt
203 ninfo( 1 ) = ninfo( 1 ) + 1
204 GO TO 200
205 END IF
206*
207* Generate unitary matrix
208*
209 CALL zlacpy( 'L', n, n, t, ldt, q, ldt )
210 CALL zunghr( n, 1, n, q, ldt, work( 1 ), work( n+1 ), lwork-n,
211 $ info )
212*
213* Compute Schur form
214*
215 DO 50 j = 1, n - 2
216 DO 40 i = j + 2, n
217 t( i, j ) = czero
218 40 CONTINUE
219 50 CONTINUE
220 CALL zhseqr( 'S', 'V', n, 1, n, t, ldt, w, q, ldt, work, lwork,
221 $ info )
222 IF( info.NE.0 ) THEN
223 lmax( 2 ) = knt
224 ninfo( 2 ) = ninfo( 2 ) + 1
225 GO TO 200
226 END IF
227*
228* Sort, select eigenvalues
229*
230 DO 60 i = 1, n
231 ipnt( i ) = i
232 SELECT( i ) = .false.
233 60 CONTINUE
234 IF( isrt.EQ.0 ) THEN
235 DO 70 i = 1, n
236 wsrt( i ) = dble( w( i ) )
237 70 CONTINUE
238 ELSE
239 DO 80 i = 1, n
240 wsrt( i ) = dimag( w( i ) )
241 80 CONTINUE
242 END IF
243 DO 100 i = 1, n - 1
244 kmin = i
245 vmin = wsrt( i )
246 DO 90 j = i + 1, n
247 IF( wsrt( j ).LT.vmin ) THEN
248 kmin = j
249 vmin = wsrt( j )
250 END IF
251 90 CONTINUE
252 wsrt( kmin ) = wsrt( i )
253 wsrt( i ) = vmin
254 itmp = ipnt( i )
255 ipnt( i ) = ipnt( kmin )
256 ipnt( kmin ) = itmp
257 100 CONTINUE
258 DO 110 i = 1, ndim
259 SELECT( ipnt( iselec( i ) ) ) = .true.
260 110 CONTINUE
261*
262* Compute condition numbers
263*
264 CALL zlacpy( 'F', n, n, q, ldt, qsav, ldt )
265 CALL zlacpy( 'F', n, n, t, ldt, tsav1, ldt )
266 CALL ztrsen( 'B', 'V', SELECT, n, t, ldt, q, ldt, wtmp, m, s,
267 $ sep, work, lwork, info )
268 IF( info.NE.0 ) THEN
269 lmax( 3 ) = knt
270 ninfo( 3 ) = ninfo( 3 ) + 1
271 GO TO 200
272 END IF
273 septmp = sep / vmul
274 stmp = s
275*
276* Compute residuals
277*
278 CALL zhst01( n, 1, n, tsav, ldt, t, ldt, q, ldt, work, lwork,
279 $ rwork, result )
280 vmax = max( result( 1 ), result( 2 ) )
281 IF( vmax.GT.rmax( 1 ) ) THEN
282 rmax( 1 ) = vmax
283 IF( ninfo( 1 ).EQ.0 )
284 $ lmax( 1 ) = knt
285 END IF
286*
287* Compare condition number for eigenvalue cluster
288* taking its condition number into account
289*
290 v = max( two*dble( n )*eps*tnrm, smlnum )
291 IF( tnrm.EQ.zero )
292 $ v = one
293 IF( v.GT.septmp ) THEN
294 tol = one
295 ELSE
296 tol = v / septmp
297 END IF
298 IF( v.GT.sepin ) THEN
299 tolin = one
300 ELSE
301 tolin = v / sepin
302 END IF
303 tol = max( tol, smlnum / eps )
304 tolin = max( tolin, smlnum / eps )
305 IF( eps*( sin-tolin ).GT.stmp+tol ) THEN
306 vmax = one / eps
307 ELSE IF( sin-tolin.GT.stmp+tol ) THEN
308 vmax = ( sin-tolin ) / ( stmp+tol )
309 ELSE IF( sin+tolin.LT.eps*( stmp-tol ) ) THEN
310 vmax = one / eps
311 ELSE IF( sin+tolin.LT.stmp-tol ) THEN
312 vmax = ( stmp-tol ) / ( sin+tolin )
313 ELSE
314 vmax = one
315 END IF
316 IF( vmax.GT.rmax( 2 ) ) THEN
317 rmax( 2 ) = vmax
318 IF( ninfo( 2 ).EQ.0 )
319 $ lmax( 2 ) = knt
320 END IF
321*
322* Compare condition numbers for invariant subspace
323* taking its condition number into account
324*
325 IF( v.GT.septmp*stmp ) THEN
326 tol = septmp
327 ELSE
328 tol = v / stmp
329 END IF
330 IF( v.GT.sepin*sin ) THEN
331 tolin = sepin
332 ELSE
333 tolin = v / sin
334 END IF
335 tol = max( tol, smlnum / eps )
336 tolin = max( tolin, smlnum / eps )
337 IF( eps*( sepin-tolin ).GT.septmp+tol ) THEN
338 vmax = one / eps
339 ELSE IF( sepin-tolin.GT.septmp+tol ) THEN
340 vmax = ( sepin-tolin ) / ( septmp+tol )
341 ELSE IF( sepin+tolin.LT.eps*( septmp-tol ) ) THEN
342 vmax = one / eps
343 ELSE IF( sepin+tolin.LT.septmp-tol ) THEN
344 vmax = ( septmp-tol ) / ( sepin+tolin )
345 ELSE
346 vmax = one
347 END IF
348 IF( vmax.GT.rmax( 2 ) ) THEN
349 rmax( 2 ) = vmax
350 IF( ninfo( 2 ).EQ.0 )
351 $ lmax( 2 ) = knt
352 END IF
353*
354* Compare condition number for eigenvalue cluster
355* without taking its condition number into account
356*
357 IF( sin.LE.dble( 2*n )*eps .AND. stmp.LE.dble( 2*n )*eps ) THEN
358 vmax = one
359 ELSE IF( eps*sin.GT.stmp ) THEN
360 vmax = one / eps
361 ELSE IF( sin.GT.stmp ) THEN
362 vmax = sin / stmp
363 ELSE IF( sin.LT.eps*stmp ) THEN
364 vmax = one / eps
365 ELSE IF( sin.LT.stmp ) THEN
366 vmax = stmp / sin
367 ELSE
368 vmax = one
369 END IF
370 IF( vmax.GT.rmax( 3 ) ) THEN
371 rmax( 3 ) = vmax
372 IF( ninfo( 3 ).EQ.0 )
373 $ lmax( 3 ) = knt
374 END IF
375*
376* Compare condition numbers for invariant subspace
377* without taking its condition number into account
378*
379 IF( sepin.LE.v .AND. septmp.LE.v ) THEN
380 vmax = one
381 ELSE IF( eps*sepin.GT.septmp ) THEN
382 vmax = one / eps
383 ELSE IF( sepin.GT.septmp ) THEN
384 vmax = sepin / septmp
385 ELSE IF( sepin.LT.eps*septmp ) THEN
386 vmax = one / eps
387 ELSE IF( sepin.LT.septmp ) THEN
388 vmax = septmp / sepin
389 ELSE
390 vmax = one
391 END IF
392 IF( vmax.GT.rmax( 3 ) ) THEN
393 rmax( 3 ) = vmax
394 IF( ninfo( 3 ).EQ.0 )
395 $ lmax( 3 ) = knt
396 END IF
397*
398* Compute eigenvalue condition number only and compare
399* Update Q
400*
401 vmax = zero
402 CALL zlacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
403 CALL zlacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
404 septmp = -one
405 stmp = -one
406 CALL ztrsen( 'E', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
407 $ m, stmp, septmp, work, lwork, info )
408 IF( info.NE.0 ) THEN
409 lmax( 3 ) = knt
410 ninfo( 3 ) = ninfo( 3 ) + 1
411 GO TO 200
412 END IF
413 IF( s.NE.stmp )
414 $ vmax = one / eps
415 IF( -one.NE.septmp )
416 $ vmax = one / eps
417 DO 130 i = 1, n
418 DO 120 j = 1, n
419 IF( ttmp( i, j ).NE.t( i, j ) )
420 $ vmax = one / eps
421 IF( qtmp( i, j ).NE.q( i, j ) )
422 $ vmax = one / eps
423 120 CONTINUE
424 130 CONTINUE
425*
426* Compute invariant subspace condition number only and compare
427* Update Q
428*
429 CALL zlacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
430 CALL zlacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
431 septmp = -one
432 stmp = -one
433 CALL ztrsen( 'V', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
434 $ m, stmp, septmp, work, lwork, info )
435 IF( info.NE.0 ) THEN
436 lmax( 3 ) = knt
437 ninfo( 3 ) = ninfo( 3 ) + 1
438 GO TO 200
439 END IF
440 IF( -one.NE.stmp )
441 $ vmax = one / eps
442 IF( sep.NE.septmp )
443 $ vmax = one / eps
444 DO 150 i = 1, n
445 DO 140 j = 1, n
446 IF( ttmp( i, j ).NE.t( i, j ) )
447 $ vmax = one / eps
448 IF( qtmp( i, j ).NE.q( i, j ) )
449 $ vmax = one / eps
450 140 CONTINUE
451 150 CONTINUE
452*
453* Compute eigenvalue condition number only and compare
454* Do not update Q
455*
456 CALL zlacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
457 CALL zlacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
458 septmp = -one
459 stmp = -one
460 CALL ztrsen( 'E', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
461 $ m, stmp, septmp, work, lwork, info )
462 IF( info.NE.0 ) THEN
463 lmax( 3 ) = knt
464 ninfo( 3 ) = ninfo( 3 ) + 1
465 GO TO 200
466 END IF
467 IF( s.NE.stmp )
468 $ vmax = one / eps
469 IF( -one.NE.septmp )
470 $ vmax = one / eps
471 DO 170 i = 1, n
472 DO 160 j = 1, n
473 IF( ttmp( i, j ).NE.t( i, j ) )
474 $ vmax = one / eps
475 IF( qtmp( i, j ).NE.qsav( i, j ) )
476 $ vmax = one / eps
477 160 CONTINUE
478 170 CONTINUE
479*
480* Compute invariant subspace condition number only and compare
481* Do not update Q
482*
483 CALL zlacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
484 CALL zlacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
485 septmp = -one
486 stmp = -one
487 CALL ztrsen( 'V', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
488 $ m, stmp, septmp, work, lwork, info )
489 IF( info.NE.0 ) THEN
490 lmax( 3 ) = knt
491 ninfo( 3 ) = ninfo( 3 ) + 1
492 GO TO 200
493 END IF
494 IF( -one.NE.stmp )
495 $ vmax = one / eps
496 IF( sep.NE.septmp )
497 $ vmax = one / eps
498 DO 190 i = 1, n
499 DO 180 j = 1, n
500 IF( ttmp( i, j ).NE.t( i, j ) )
501 $ vmax = one / eps
502 IF( qtmp( i, j ).NE.qsav( i, j ) )
503 $ vmax = one / eps
504 180 CONTINUE
505 190 CONTINUE
506 IF( vmax.GT.rmax( 1 ) ) THEN
507 rmax( 1 ) = vmax
508 IF( ninfo( 1 ).EQ.0 )
509 $ lmax( 1 ) = knt
510 END IF
511 200 CONTINUE
512 GO TO 10
513*
514* End of ZGET38
515*
516 END
dlabad
subroutine dlabad(small, large)
DLABAD
Definition dlabad.f:74
zgehrd
subroutine zgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
ZGEHRD
Definition zgehrd.f:167
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
zhseqr
subroutine zhseqr(job, compz, n, ilo, ihi, h, ldh, w, z, ldz, work, lwork, info)
ZHSEQR
Definition zhseqr.f:299
ztrsen
subroutine ztrsen(job, compq, select, n, t, ldt, q, ldq, w, m, s, sep, work, lwork, info)
ZTRSEN
Definition ztrsen.f:264
zunghr
subroutine zunghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
ZUNGHR
Definition zunghr.f:126
zdscal
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
zhst01
subroutine zhst01(n, ilo, ihi, a, lda, h, ldh, q, ldq, work, lwork, rwork, result)
ZHST01
Definition zhst01.f:140
zget38
subroutine zget38(rmax, lmax, ninfo, knt, nin)
ZGET38
Definition zget38.f:91
max
#define max(a, b)
Definition macros.h:21