LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dget38()

subroutine dget38 ( double precision, dimension( 3 )  RMAX,
integer, dimension( 3 )  LMAX,
integer, dimension( 3 )  NINFO,
integer  KNT,
integer  NIN 
)

DGET38

Purpose:
 DGET38 tests DTRSEN, a routine for estimating condition numbers of a
 cluster of eigenvalues and/or its associated right invariant subspace

 The test matrices are read from a file with logical unit number NIN.
Parameters
[out]RMAX
          RMAX is DOUBLE PRECISION array, dimension (3)
          Values of the largest test ratios.
          RMAX(1) = largest residuals from DHST01 or comparing
                    different calls to DTRSEN
          RMAX(2) = largest error in reciprocal condition
                    numbers taking their conditioning into account
          RMAX(3) = largest error in reciprocal condition
                    numbers not taking their conditioning into
                    account (may be larger than RMAX(2))
[out]LMAX
          LMAX is INTEGER array, dimension (3)
          LMAX(i) is example number where largest test ratio
          RMAX(i) is achieved. Also:
          If DGEHRD returns INFO nonzero on example i, LMAX(1)=i
          If DHSEQR returns INFO nonzero on example i, LMAX(2)=i
          If DTRSEN returns INFO nonzero on example i, LMAX(3)=i
[out]NINFO
          NINFO is INTEGER array, dimension (3)
          NINFO(1) = No. of times DGEHRD returned INFO nonzero
          NINFO(2) = No. of times DHSEQR returned INFO nonzero
          NINFO(3) = No. of times DTRSEN returned INFO nonzero
[out]KNT
          KNT is INTEGER
          Total number of examples tested.
[in]NIN
          NIN is INTEGER
          Input logical unit number.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 90 of file dget38.f.

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 IF( info.NE.0 ) 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 IF( wrtmp( j ).LT.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 IF( info.NE.0 ) 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 IF( vmax.GT.rmax( 1 ) ) THEN
274 rmax( 1 ) = vmax
275 IF( ninfo( 1 ).EQ.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 IF( tnrm.EQ.zero )
284 $ v = one
285 IF( v.GT.septmp ) THEN
286 tol = one
287 ELSE
288 tol = v / septmp
289 END IF
290 IF( v.GT.sepin ) 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 IF( eps*( sin-tolin ).GT.stmp+tol ) THEN
298 vmax = one / eps
299 ELSE IF( sin-tolin.GT.stmp+tol ) THEN
300 vmax = ( sin-tolin ) / ( stmp+tol )
301 ELSE IF( sin+tolin.LT.eps*( stmp-tol ) ) THEN
302 vmax = one / eps
303 ELSE IF( sin+tolin.LT.stmp-tol ) THEN
304 vmax = ( stmp-tol ) / ( sin+tolin )
305 ELSE
306 vmax = one
307 END IF
308 IF( vmax.GT.rmax( 2 ) ) THEN
309 rmax( 2 ) = vmax
310 IF( ninfo( 2 ).EQ.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 IF( v.GT.septmp*stmp ) THEN
318 tol = septmp
319 ELSE
320 tol = v / stmp
321 END IF
322 IF( v.GT.sepin*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 IF( eps*( sepin-tolin ).GT.septmp+tol ) THEN
330 vmax = one / eps
331 ELSE IF( sepin-tolin.GT.septmp+tol ) THEN
332 vmax = ( sepin-tolin ) / ( septmp+tol )
333 ELSE IF( sepin+tolin.LT.eps*( septmp-tol ) ) THEN
334 vmax = one / eps
335 ELSE IF( sepin+tolin.LT.septmp-tol ) THEN
336 vmax = ( septmp-tol ) / ( sepin+tolin )
337 ELSE
338 vmax = one
339 END IF
340 IF( vmax.GT.rmax( 2 ) ) THEN
341 rmax( 2 ) = vmax
342 IF( ninfo( 2 ).EQ.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 IF( sin.LE.dble( 2*n )*eps .AND. stmp.LE.dble( 2*n )*eps ) THEN
350 vmax = one
351 ELSE IF( eps*sin.GT.stmp ) THEN
352 vmax = one / eps
353 ELSE IF( sin.GT.stmp ) THEN
354 vmax = sin / stmp
355 ELSE IF( sin.LT.eps*stmp ) THEN
356 vmax = one / eps
357 ELSE IF( sin.LT.stmp ) THEN
358 vmax = stmp / sin
359 ELSE
360 vmax = one
361 END IF
362 IF( vmax.GT.rmax( 3 ) ) THEN
363 rmax( 3 ) = vmax
364 IF( ninfo( 3 ).EQ.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 IF( sepin.LE.v .AND. septmp.LE.v ) THEN
372 vmax = one
373 ELSE IF( eps*sepin.GT.septmp ) THEN
374 vmax = one / eps
375 ELSE IF( sepin.GT.septmp ) THEN
376 vmax = sepin / septmp
377 ELSE IF( sepin.LT.eps*septmp ) THEN
378 vmax = one / eps
379 ELSE IF( sepin.LT.septmp ) THEN
380 vmax = septmp / sepin
381 ELSE
382 vmax = one
383 END IF
384 IF( vmax.GT.rmax( 3 ) ) THEN
385 rmax( 3 ) = vmax
386 IF( ninfo( 3 ).EQ.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 IF( info.NE.0 ) THEN
402 lmax( 3 ) = knt
403 ninfo( 3 ) = ninfo( 3 ) + 1
404 GO TO 160
405 END IF
406 IF( s.NE.stmp )
407 $ vmax = one / eps
408 IF( -one.NE.septmp )
409 $ vmax = one / eps
410 DO 90 i = 1, n
411 DO 80 j = 1, n
412 IF( ttmp( i, j ).NE.t( i, j ) )
413 $ vmax = one / eps
414 IF( qtmp( i, j ).NE.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 IF( info.NE.0 ) THEN
430 lmax( 3 ) = knt
431 ninfo( 3 ) = ninfo( 3 ) + 1
432 GO TO 160
433 END IF
434 IF( -one.NE.stmp )
435 $ vmax = one / eps
436 IF( sep.NE.septmp )
437 $ vmax = one / eps
438 DO 110 i = 1, n
439 DO 100 j = 1, n
440 IF( ttmp( i, j ).NE.t( i, j ) )
441 $ vmax = one / eps
442 IF( qtmp( i, j ).NE.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 IF( info.NE.0 ) THEN
458 lmax( 3 ) = knt
459 ninfo( 3 ) = ninfo( 3 ) + 1
460 GO TO 160
461 END IF
462 IF( s.NE.stmp )
463 $ vmax = one / eps
464 IF( -one.NE.septmp )
465 $ vmax = one / eps
466 DO 130 i = 1, n
467 DO 120 j = 1, n
468 IF( ttmp( i, j ).NE.t( i, j ) )
469 $ vmax = one / eps
470 IF( qtmp( i, j ).NE.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 IF( info.NE.0 ) THEN
486 lmax( 3 ) = knt
487 ninfo( 3 ) = ninfo( 3 ) + 1
488 GO TO 160
489 END IF
490 IF( -one.NE.stmp )
491 $ vmax = one / eps
492 IF( sep.NE.septmp )
493 $ vmax = one / eps
494 DO 150 i = 1, n
495 DO 140 j = 1, n
496 IF( ttmp( i, j ).NE.t( i, j ) )
497 $ vmax = one / eps
498 IF( qtmp( i, j ).NE.qsav( i, j ) )
499 $ vmax = one / eps
500 140 CONTINUE
501 150 CONTINUE
502 IF( vmax.GT.rmax( 1 ) ) THEN
503 rmax( 1 ) = vmax
504 IF( ninfo( 1 ).EQ.0 )
505 $ lmax( 1 ) = knt
506 END IF
507 160 CONTINUE
508 GO TO 10
509*
510* End of DGET38
511*
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
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
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
dhst01
subroutine dhst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
DHST01
Definition: dhst01.f:134
dlange
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:114
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
dorghr
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
Definition: dorghr.f:126
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
Here is the call graph for this function:
Here is the caller graph for this function: