OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dget38.f
Go to the documentation of this file.
1*> \brief \b DGET38
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 DGET38( 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*> DGET38 tests DTRSEN, 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 DHST01 or comparing
41*> different calls to DTRSEN
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 DGEHRD returns INFO nonzero on example i, LMAX(1)=i
55*> If DHSEQR returns INFO nonzero on example i, LMAX(2)=i
56*> If DTRSEN 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 DGEHRD returned INFO nonzero
63*> NINFO(2) = No. of times DHSEQR returned INFO nonzero
64*> NINFO(3) = No. of times DTRSEN 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 double_eig
88*
89* =====================================================================
90 SUBROUTINE dget38( 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 DOUBLE PRECISION ZERO, ONE, TWO
108 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
109 DOUBLE PRECISION EPSIN
110 parameter( epsin = 5.9605d-8 )
111 INTEGER LDT, LWORK
112 parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
113 INTEGER LIWORK
114 parameter( liwork = ldt*ldt )
115* ..
116* .. Local Scalars ..
117 INTEGER I, INFO, ISCL, ITMP, J, KMIN, M, N, NDIM
118 DOUBLE PRECISION BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN,
119 $ SMLNUM, STMP, TNRM, TOL, TOLIN, V, VIMIN, VMAX,
120 $ VMUL, VRMIN
121* ..
122* .. Local Arrays ..
123 LOGICAL SELECT( LDT )
124 INTEGER IPNT( LDT ), ISELEC( LDT ), IWORK( LIWORK )
125 DOUBLE PRECISION Q( LDT, LDT ), QSAV( LDT, LDT ),
126 $ QTMP( LDT, LDT ), RESULT( 2 ), T( LDT, LDT ),
127 $ TMP( LDT, LDT ), TSAV( LDT, LDT ),
128 $ TSAV1( LDT, LDT ), TTMP( LDT, LDT ), VAL( 3 ),
129 $ WI( LDT ), WITMP( LDT ), WORK( LWORK ),
130 $ WR( LDT ), WRTMP( LDT )
131* ..
132* .. External Functions ..
133 DOUBLE PRECISION DLAMCH, DLANGE
134 EXTERNAL dlamch, dlange
135* ..
136* .. External Subroutines ..
137 EXTERNAL dcopy, dgehrd, dhseqr, dhst01, dlabad, dlacpy,
138 $ dorghr, dscal, dtrsen
139* ..
140* .. Intrinsic Functions ..
141 INTRINSIC dble, max, sqrt
142* ..
143* .. Executable Statements ..
144*
145 eps = dlamch( 'P' )
146 smlnum = dlamch( 'S' ) / eps
147 bignum = one / smlnum
148 CALL dlabad( smlnum, bignum )
149*
150* EPSIN = 2**(-24) = precision to which input data computed
151*
152 eps = max( eps, epsin )
153 rmax( 1 ) = zero
154 rmax( 2 ) = zero
155 rmax( 3 ) = zero
156 lmax( 1 ) = 0
157 lmax( 2 ) = 0
158 lmax( 3 ) = 0
159 knt = 0
160 ninfo( 1 ) = 0
161 ninfo( 2 ) = 0
162 ninfo( 3 ) = 0
163*
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
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 = dlange( 'M', n, n, tmp, ldt, work )
183 DO 160 iscl = 1, 3
184*
185* Scale input matrix
186*
187 knt = knt + 1
188 CALL dlacpy( 'F', n, n, tmp, ldt, t, ldt )
189 vmul = val( iscl )
190 DO 30 i = 1, n
191 CALL dscal( n, vmul, t( 1, i ), 1 )
192 30 CONTINUE
193 IF( tnrm.EQ.zero )
194 $ vmul = one
195 CALL dlacpy( 'F', n, n, t, ldt, tsav, ldt )
196*
197* Compute Schur form
198*
199 CALL dgehrd( 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 160
205 END IF
206*
207* Generate orthogonal matrix
208*
209 CALL dlacpy( 'l', N, N, T, LDT, Q, LDT )
210 CALL DORGHR( N, 1, N, Q, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
211 $ INFO )
212*
213* Compute Schur form
214*
215 CALL DHSEQR( 's', 'v', N, 1, N, T, LDT, WR, WI, Q, LDT, WORK,
216 $ LWORK, INFO )
217.NE. IF( INFO0 ) THEN
218 LMAX( 2 ) = KNT
219 NINFO( 2 ) = NINFO( 2 ) + 1
220 GO TO 160
221 END IF
222*
223* Sort, select eigenvalues
224*
225 DO 40 I = 1, N
226 IPNT( I ) = I
227 SELECT( I ) = .FALSE.
228 40 CONTINUE
229 CALL DCOPY( N, WR, 1, WRTMP, 1 )
230 CALL DCOPY( N, WI, 1, WITMP, 1 )
231 DO 60 I = 1, N - 1
232 KMIN = I
233 VRMIN = WRTMP( I )
234 VIMIN = WITMP( I )
235 DO 50 J = I + 1, N
236.LT. IF( WRTMP( J )VRMIN ) THEN
237 KMIN = J
238 VRMIN = WRTMP( J )
239 VIMIN = WITMP( J )
240 END IF
241 50 CONTINUE
242 WRTMP( KMIN ) = WRTMP( I )
243 WITMP( KMIN ) = WITMP( I )
244 WRTMP( I ) = VRMIN
245 WITMP( I ) = VIMIN
246 ITMP = IPNT( I )
247 IPNT( I ) = IPNT( KMIN )
248 IPNT( KMIN ) = ITMP
249 60 CONTINUE
250 DO 70 I = 1, NDIM
251 SELECT( IPNT( ISELEC( I ) ) ) = .TRUE.
252 70 CONTINUE
253*
254* Compute condition numbers
255*
256 CALL DLACPY( 'f', N, N, Q, LDT, QSAV, LDT )
257 CALL DLACPY( 'f', N, N, T, LDT, TSAV1, LDT )
258 CALL DTRSEN( 'b', 'v', SELECT, N, T, LDT, Q, LDT, WRTMP, WITMP,
259 $ M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
260.NE. IF( INFO0 ) THEN
261 LMAX( 3 ) = KNT
262 NINFO( 3 ) = NINFO( 3 ) + 1
263 GO TO 160
264 END IF
265 SEPTMP = SEP / VMUL
266 STMP = S
267*
268* Compute residuals
269*
270 CALL DHST01( N, 1, N, TSAV, LDT, T, LDT, Q, LDT, WORK, LWORK,
271 $ RESULT )
272 VMAX = MAX( RESULT( 1 ), RESULT( 2 ) )
273.GT. IF( VMAXRMAX( 1 ) ) THEN
274 RMAX( 1 ) = VMAX
275.EQ. IF( NINFO( 1 )0 )
276 $ LMAX( 1 ) = KNT
277 END IF
278*
279* Compare condition number for eigenvalue cluster
280* taking its condition number into account
281*
282 V = MAX( TWO*DBLE( N )*EPS*TNRM, SMLNUM )
283.EQ. IF( TNRMZERO )
284 $ V = ONE
285.GT. IF( VSEPTMP ) THEN
286 TOL = ONE
287 ELSE
288 TOL = V / SEPTMP
289 END IF
290.GT. IF( VSEPIN ) THEN
291 TOLIN = ONE
292 ELSE
293 TOLIN = V / SEPIN
294 END IF
295 TOL = MAX( TOL, SMLNUM / EPS )
296 TOLIN = MAX( TOLIN, SMLNUM / EPS )
297.GT. IF( EPS*( SIN-TOLIN )STMP+TOL ) THEN
298 VMAX = ONE / EPS
299.GT. ELSE IF( SIN-TOLINSTMP+TOL ) THEN
300 VMAX = ( SIN-TOLIN ) / ( STMP+TOL )
301.LT. ELSE IF( SIN+TOLINEPS*( STMP-TOL ) ) THEN
302 VMAX = ONE / EPS
303.LT. ELSE IF( SIN+TOLINSTMP-TOL ) THEN
304 VMAX = ( STMP-TOL ) / ( SIN+TOLIN )
305 ELSE
306 VMAX = ONE
307 END IF
308.GT. IF( VMAXRMAX( 2 ) ) THEN
309 RMAX( 2 ) = VMAX
310.EQ. IF( NINFO( 2 )0 )
311 $ LMAX( 2 ) = KNT
312 END IF
313*
314* Compare condition numbers for invariant subspace
315* taking its condition number into account
316*
317.GT. IF( VSEPTMP*STMP ) THEN
318 TOL = SEPTMP
319 ELSE
320 TOL = V / STMP
321 END IF
322.GT. IF( VSEPIN*SIN ) THEN
323 TOLIN = SEPIN
324 ELSE
325 TOLIN = V / SIN
326 END IF
327 TOL = MAX( TOL, SMLNUM / EPS )
328 TOLIN = MAX( TOLIN, SMLNUM / EPS )
329.GT. IF( EPS*( SEPIN-TOLIN )SEPTMP+TOL ) THEN
330 VMAX = ONE / EPS
331.GT. ELSE IF( SEPIN-TOLINSEPTMP+TOL ) THEN
332 VMAX = ( SEPIN-TOLIN ) / ( SEPTMP+TOL )
333.LT. ELSE IF( SEPIN+TOLINEPS*( SEPTMP-TOL ) ) THEN
334 VMAX = ONE / EPS
335.LT. ELSE IF( SEPIN+TOLINSEPTMP-TOL ) THEN
336 VMAX = ( SEPTMP-TOL ) / ( SEPIN+TOLIN )
337 ELSE
338 VMAX = ONE
339 END IF
340.GT. IF( VMAXRMAX( 2 ) ) THEN
341 RMAX( 2 ) = VMAX
342.EQ. IF( NINFO( 2 )0 )
343 $ LMAX( 2 ) = KNT
344 END IF
345*
346* Compare condition number for eigenvalue cluster
347* without taking its condition number into account
348*
349.LE..AND..LE. IF( SINDBLE( 2*N )*EPS STMPDBLE( 2*N )*EPS ) THEN
350 VMAX = ONE
351.GT. ELSE IF( EPS*SINSTMP ) THEN
352 VMAX = ONE / EPS
353.GT. ELSE IF( SINSTMP ) THEN
354 VMAX = SIN / STMP
355.LT. ELSE IF( SINEPS*STMP ) THEN
356 VMAX = ONE / EPS
357.LT. ELSE IF( SINSTMP ) THEN
358 VMAX = STMP / SIN
359 ELSE
360 VMAX = ONE
361 END IF
362.GT. IF( VMAXRMAX( 3 ) ) THEN
363 RMAX( 3 ) = VMAX
364.EQ. IF( NINFO( 3 )0 )
365 $ LMAX( 3 ) = KNT
366 END IF
367*
368* Compare condition numbers for invariant subspace
369* without taking its condition number into account
370*
371.LE..AND..LE. IF( SEPINV SEPTMPV ) THEN
372 VMAX = ONE
373.GT. ELSE IF( EPS*SEPINSEPTMP ) THEN
374 VMAX = ONE / EPS
375.GT. ELSE IF( SEPINSEPTMP ) THEN
376 VMAX = SEPIN / SEPTMP
377.LT. ELSE IF( SEPINEPS*SEPTMP ) THEN
378 VMAX = ONE / EPS
379.LT. ELSE IF( SEPINSEPTMP ) THEN
380 VMAX = SEPTMP / SEPIN
381 ELSE
382 VMAX = ONE
383 END IF
384.GT. IF( VMAXRMAX( 3 ) ) THEN
385 RMAX( 3 ) = VMAX
386.EQ. IF( NINFO( 3 )0 )
387 $ LMAX( 3 ) = KNT
388 END IF
389*
390* Compute eigenvalue condition number only and compare
391* Update Q
392*
393 VMAX = ZERO
394 CALL DLACPY( 'f', N, N, TSAV1, LDT, TTMP, LDT )
395 CALL DLACPY( 'f', N, N, QSAV, LDT, QTMP, LDT )
396 SEPTMP = -ONE
397 STMP = -ONE
398 CALL DTRSEN( 'e', 'v', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
399 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
400 $ LIWORK, INFO )
401.NE. IF( INFO0 ) THEN
402 LMAX( 3 ) = KNT
403 NINFO( 3 ) = NINFO( 3 ) + 1
404 GO TO 160
405 END IF
406.NE. IF( SSTMP )
407 $ VMAX = ONE / EPS
408.NE. IF( -ONESEPTMP )
409 $ VMAX = ONE / EPS
410 DO 90 I = 1, N
411 DO 80 J = 1, N
412.NE. IF( TTMP( I, J )T( I, J ) )
413 $ VMAX = ONE / EPS
414.NE. IF( QTMP( I, J )Q( I, J ) )
415 $ VMAX = ONE / EPS
416 80 CONTINUE
417 90 CONTINUE
418*
419* Compute invariant subspace condition number only and compare
420* Update Q
421*
422 CALL DLACPY( 'f', N, N, TSAV1, LDT, TTMP, LDT )
423 CALL DLACPY( 'f', N, N, QSAV, LDT, QTMP, LDT )
424 SEPTMP = -ONE
425 STMP = -ONE
426 CALL DTRSEN( 'v', 'v', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
427 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
428 $ LIWORK, INFO )
429.NE. IF( INFO0 ) THEN
430 LMAX( 3 ) = KNT
431 NINFO( 3 ) = NINFO( 3 ) + 1
432 GO TO 160
433 END IF
434.NE. IF( -ONESTMP )
435 $ VMAX = ONE / EPS
436.NE. IF( SEPSEPTMP )
437 $ VMAX = ONE / EPS
438 DO 110 I = 1, N
439 DO 100 J = 1, N
440.NE. IF( TTMP( I, J )T( I, J ) )
441 $ VMAX = ONE / EPS
442.NE. IF( QTMP( I, J )Q( I, J ) )
443 $ VMAX = ONE / EPS
444 100 CONTINUE
445 110 CONTINUE
446*
447* Compute eigenvalue condition number only and compare
448* Do not update Q
449*
450 CALL DLACPY( 'f', N, N, TSAV1, LDT, TTMP, LDT )
451 CALL DLACPY( 'f', N, N, QSAV, LDT, QTMP, LDT )
452 SEPTMP = -ONE
453 STMP = -ONE
454 CALL DTRSEN( 'e', 'n', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
455 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
456 $ LIWORK, INFO )
457.NE. IF( INFO0 ) THEN
458 LMAX( 3 ) = KNT
459 NINFO( 3 ) = NINFO( 3 ) + 1
460 GO TO 160
461 END IF
462.NE. IF( SSTMP )
463 $ VMAX = ONE / EPS
464.NE. IF( -ONESEPTMP )
465 $ VMAX = ONE / EPS
466 DO 130 I = 1, N
467 DO 120 J = 1, N
468.NE. IF( TTMP( I, J )T( I, J ) )
469 $ VMAX = ONE / EPS
470.NE. IF( QTMP( I, J )QSAV( I, J ) )
471 $ VMAX = ONE / EPS
472 120 CONTINUE
473 130 CONTINUE
474*
475* Compute invariant subspace condition number only and compare
476* Do not update Q
477*
478 CALL DLACPY( 'f', N, N, TSAV1, LDT, TTMP, LDT )
479 CALL DLACPY( 'f', N, N, QSAV, LDT, QTMP, LDT )
480 SEPTMP = -ONE
481 STMP = -ONE
482 CALL DTRSEN( 'v', 'n', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
483 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
484 $ LIWORK, INFO )
485.NE. IF( INFO0 ) THEN
486 LMAX( 3 ) = KNT
487 NINFO( 3 ) = NINFO( 3 ) + 1
488 GO TO 160
489 END IF
490.NE. IF( -ONESTMP )
491 $ VMAX = ONE / EPS
492.NE. IF( SEPSEPTMP )
493 $ VMAX = ONE / EPS
494 DO 150 I = 1, N
495 DO 140 J = 1, N
496.NE. IF( TTMP( I, J )T( I, J ) )
497 $ VMAX = ONE / EPS
498.NE. IF( QTMP( I, J )QSAV( I, J ) )
499 $ VMAX = ONE / EPS
500 140 CONTINUE
501 150 CONTINUE
502.GT. IF( VMAXRMAX( 1 ) ) THEN
503 RMAX( 1 ) = VMAX
504.EQ. IF( NINFO( 1 )0 )
505 $ LMAX( 1 ) = KNT
506 END IF
507 160 CONTINUE
508 GO TO 10
509*
510* End of DGET38
511*
512 END
dlabad
subroutine dlabad(small, large)
DLABAD
Definition dlabad.f:74
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dgehrd
subroutine dgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
DGEHRD
Definition dgehrd.f:167
dhseqr
subroutine dhseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
DHSEQR
Definition dhseqr.f:316
dtrsen
subroutine dtrsen(job, compq, select, n, t, ldt, q, ldq, wr, wi, m, s, sep, work, lwork, iwork, liwork, info)
DTRSEN
Definition dtrsen.f:313
dorghr
subroutine dorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
DORGHR
Definition dorghr.f:126
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dhst01
subroutine dhst01(n, ilo, ihi, a, lda, h, ldh, q, ldq, work, lwork, result)
DHST01
Definition dhst01.f:134
dget38
subroutine dget38(rmax, lmax, ninfo, knt, nin)
DGET38
Definition dget38.f:91
max
#define max(a, b)
Definition macros.h:21