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

◆ dget37()

subroutine dget37 ( double precision, dimension( 3 )  rmax,
integer, dimension( 3 )  lmax,
integer, dimension( 3 )  ninfo,
integer  knt,
integer  nin 
)

DGET37

Purpose:
 DGET37 tests DTRSNA, a routine for estimating condition numbers of
 eigenvalues and/or right eigenvectors of a matrix.

 The test matrices are read from a file with logical unit number NIN.
Parameters
[out]RMAX
          RMAX is DOUBLE PRECISION array, dimension (3)
          Value of the largest test ratio.
          RMAX(1) = largest ratio comparing different calls to DTRSNA
          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 DTRSNA 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 DTRSNA 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 89 of file dget37.f.

90*
91* -- LAPACK test routine --
92* -- LAPACK is a software package provided by Univ. of Tennessee, --
93* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
94*
95* .. Scalar Arguments ..
96 INTEGER KNT, NIN
97* ..
98* .. Array Arguments ..
99 INTEGER LMAX( 3 ), NINFO( 3 )
100 DOUBLE PRECISION RMAX( 3 )
101* ..
102*
103* =====================================================================
104*
105* .. Parameters ..
106 DOUBLE PRECISION ZERO, ONE, TWO
107 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
108 DOUBLE PRECISION EPSIN
109 parameter( epsin = 5.9605d-8 )
110 INTEGER LDT, LWORK
111 parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
112* ..
113* .. Local Scalars ..
114 INTEGER I, ICMP, IFND, INFO, ISCL, J, KMIN, M, N
115 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
116 $ VIMIN, VMAX, VMUL, VRMIN
117* ..
118* .. Local Arrays ..
119 LOGICAL SELECT( LDT )
120 INTEGER IWORK( 2*LDT ), LCMP( 3 )
121 DOUBLE PRECISION DUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
122 $ S( LDT ), SEP( LDT ), SEPIN( LDT ),
123 $ SEPTMP( LDT ), SIN( LDT ), STMP( LDT ),
124 $ T( LDT, LDT ), TMP( LDT, LDT ), VAL( 3 ),
125 $ WI( LDT ), WIIN( LDT ), WITMP( LDT ),
126 $ WORK( LWORK ), WR( LDT ), WRIN( LDT ),
127 $ WRTMP( LDT )
128* ..
129* .. External Functions ..
130 DOUBLE PRECISION DLAMCH, DLANGE
131 EXTERNAL dlamch, dlange
132* ..
133* .. External Subroutines ..
134 EXTERNAL dcopy, dgehrd, dhseqr, dlacpy, dscal, dtrevc,
135 $ dtrsna
136* ..
137* .. Intrinsic Functions ..
138 INTRINSIC dble, max, sqrt
139* ..
140* .. Executable Statements ..
141*
142 eps = dlamch( 'P' )
143 smlnum = dlamch( 'S' ) / eps
144 bignum = one / smlnum
145*
146* EPSIN = 2**(-24) = precision to which input data computed
147*
148 eps = max( eps, epsin )
149 rmax( 1 ) = zero
150 rmax( 2 ) = zero
151 rmax( 3 ) = zero
152 lmax( 1 ) = 0
153 lmax( 2 ) = 0
154 lmax( 3 ) = 0
155 knt = 0
156 ninfo( 1 ) = 0
157 ninfo( 2 ) = 0
158 ninfo( 3 ) = 0
159*
160 val( 1 ) = sqrt( smlnum )
161 val( 2 ) = one
162 val( 3 ) = sqrt( bignum )
163*
164* Read input data until N=0. Assume input eigenvalues are sorted
165* lexicographically (increasing by real part, then decreasing by
166* imaginary part)
167*
168 10 CONTINUE
169 READ( nin, fmt = * )n
170 IF( n.EQ.0 )
171 $ RETURN
172 DO 20 i = 1, n
173 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
174 20 CONTINUE
175 DO 30 i = 1, n
176 READ( nin, fmt = * )wrin( i ), wiin( i ), sin( i ), sepin( i )
177 30 CONTINUE
178 tnrm = dlange( 'M', n, n, tmp, ldt, work )
179*
180* Begin test
181*
182 DO 240 iscl = 1, 3
183*
184* Scale input matrix
185*
186 knt = knt + 1
187 CALL dlacpy( 'F', n, n, tmp, ldt, t, ldt )
188 vmul = val( iscl )
189 DO 40 i = 1, n
190 CALL dscal( n, vmul, t( 1, i ), 1 )
191 40 CONTINUE
192 IF( tnrm.EQ.zero )
193 $ vmul = one
194*
195* Compute eigenvalues and eigenvectors
196*
197 CALL dgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
198 $ info )
199 IF( info.NE.0 ) THEN
200 lmax( 1 ) = knt
201 ninfo( 1 ) = ninfo( 1 ) + 1
202 GO TO 240
203 END IF
204 DO 60 j = 1, n - 2
205 DO 50 i = j + 2, n
206 t( i, j ) = zero
207 50 CONTINUE
208 60 CONTINUE
209*
210* Compute Schur form
211*
212 CALL dhseqr( 'S', 'N', n, 1, n, t, ldt, wr, wi, dum, 1, work,
213 $ lwork, info )
214 IF( info.NE.0 ) THEN
215 lmax( 2 ) = knt
216 ninfo( 2 ) = ninfo( 2 ) + 1
217 GO TO 240
218 END IF
219*
220* Compute eigenvectors
221*
222 CALL dtrevc( 'Both', 'All', SELECT, n, t, ldt, le, ldt, re,
223 $ ldt, n, m, work, info )
224*
225* Compute condition numbers
226*
227 CALL dtrsna( 'Both', 'All', SELECT, n, t, ldt, le, ldt, re,
228 $ ldt, s, sep, n, m, work, n, iwork, info )
229 IF( info.NE.0 ) THEN
230 lmax( 3 ) = knt
231 ninfo( 3 ) = ninfo( 3 ) + 1
232 GO TO 240
233 END IF
234*
235* Sort eigenvalues and condition numbers lexicographically
236* to compare with inputs
237*
238 CALL dcopy( n, wr, 1, wrtmp, 1 )
239 CALL dcopy( n, wi, 1, witmp, 1 )
240 CALL dcopy( n, s, 1, stmp, 1 )
241 CALL dcopy( n, sep, 1, septmp, 1 )
242 CALL dscal( n, one / vmul, septmp, 1 )
243 DO 80 i = 1, n - 1
244 kmin = i
245 vrmin = wrtmp( i )
246 vimin = witmp( i )
247 DO 70 j = i + 1, n
248 IF( wrtmp( j ).LT.vrmin ) THEN
249 kmin = j
250 vrmin = wrtmp( j )
251 vimin = witmp( j )
252 END IF
253 70 CONTINUE
254 wrtmp( kmin ) = wrtmp( i )
255 witmp( kmin ) = witmp( i )
256 wrtmp( i ) = vrmin
257 witmp( i ) = vimin
258 vrmin = stmp( kmin )
259 stmp( kmin ) = stmp( i )
260 stmp( i ) = vrmin
261 vrmin = septmp( kmin )
262 septmp( kmin ) = septmp( i )
263 septmp( i ) = vrmin
264 80 CONTINUE
265*
266* Compare condition numbers for eigenvalues
267* taking their condition numbers into account
268*
269 v = max( two*dble( n )*eps*tnrm, smlnum )
270 IF( tnrm.EQ.zero )
271 $ v = one
272 DO 90 i = 1, n
273 IF( v.GT.septmp( i ) ) THEN
274 tol = one
275 ELSE
276 tol = v / septmp( i )
277 END IF
278 IF( v.GT.sepin( i ) ) THEN
279 tolin = one
280 ELSE
281 tolin = v / sepin( i )
282 END IF
283 tol = max( tol, smlnum / eps )
284 tolin = max( tolin, smlnum / eps )
285 IF( eps*( sin( i )-tolin ).GT.stmp( i )+tol ) THEN
286 vmax = one / eps
287 ELSE IF( sin( i )-tolin.GT.stmp( i )+tol ) THEN
288 vmax = ( sin( i )-tolin ) / ( stmp( i )+tol )
289 ELSE IF( sin( i )+tolin.LT.eps*( stmp( i )-tol ) ) THEN
290 vmax = one / eps
291 ELSE IF( sin( i )+tolin.LT.stmp( i )-tol ) THEN
292 vmax = ( stmp( i )-tol ) / ( sin( i )+tolin )
293 ELSE
294 vmax = one
295 END IF
296 IF( vmax.GT.rmax( 2 ) ) THEN
297 rmax( 2 ) = vmax
298 IF( ninfo( 2 ).EQ.0 )
299 $ lmax( 2 ) = knt
300 END IF
301 90 CONTINUE
302*
303* Compare condition numbers for eigenvectors
304* taking their condition numbers into account
305*
306 DO 100 i = 1, n
307 IF( v.GT.septmp( i )*stmp( i ) ) THEN
308 tol = septmp( i )
309 ELSE
310 tol = v / stmp( i )
311 END IF
312 IF( v.GT.sepin( i )*sin( i ) ) THEN
313 tolin = sepin( i )
314 ELSE
315 tolin = v / sin( i )
316 END IF
317 tol = max( tol, smlnum / eps )
318 tolin = max( tolin, smlnum / eps )
319 IF( eps*( sepin( i )-tolin ).GT.septmp( i )+tol ) THEN
320 vmax = one / eps
321 ELSE IF( sepin( i )-tolin.GT.septmp( i )+tol ) THEN
322 vmax = ( sepin( i )-tolin ) / ( septmp( i )+tol )
323 ELSE IF( sepin( i )+tolin.LT.eps*( septmp( i )-tol ) ) THEN
324 vmax = one / eps
325 ELSE IF( sepin( i )+tolin.LT.septmp( i )-tol ) THEN
326 vmax = ( septmp( i )-tol ) / ( sepin( i )+tolin )
327 ELSE
328 vmax = one
329 END IF
330 IF( vmax.GT.rmax( 2 ) ) THEN
331 rmax( 2 ) = vmax
332 IF( ninfo( 2 ).EQ.0 )
333 $ lmax( 2 ) = knt
334 END IF
335 100 CONTINUE
336*
337* Compare condition numbers for eigenvalues
338* without taking their condition numbers into account
339*
340 DO 110 i = 1, n
341 IF( sin( i ).LE.dble( 2*n )*eps .AND. stmp( i ).LE.
342 $ dble( 2*n )*eps ) THEN
343 vmax = one
344 ELSE IF( eps*sin( i ).GT.stmp( i ) ) THEN
345 vmax = one / eps
346 ELSE IF( sin( i ).GT.stmp( i ) ) THEN
347 vmax = sin( i ) / stmp( i )
348 ELSE IF( sin( i ).LT.eps*stmp( i ) ) THEN
349 vmax = one / eps
350 ELSE IF( sin( i ).LT.stmp( i ) ) THEN
351 vmax = stmp( i ) / sin( i )
352 ELSE
353 vmax = one
354 END IF
355 IF( vmax.GT.rmax( 3 ) ) THEN
356 rmax( 3 ) = vmax
357 IF( ninfo( 3 ).EQ.0 )
358 $ lmax( 3 ) = knt
359 END IF
360 110 CONTINUE
361*
362* Compare condition numbers for eigenvectors
363* without taking their condition numbers into account
364*
365 DO 120 i = 1, n
366 IF( sepin( i ).LE.v .AND. septmp( i ).LE.v ) THEN
367 vmax = one
368 ELSE IF( eps*sepin( i ).GT.septmp( i ) ) THEN
369 vmax = one / eps
370 ELSE IF( sepin( i ).GT.septmp( i ) ) THEN
371 vmax = sepin( i ) / septmp( i )
372 ELSE IF( sepin( i ).LT.eps*septmp( i ) ) THEN
373 vmax = one / eps
374 ELSE IF( sepin( i ).LT.septmp( i ) ) THEN
375 vmax = septmp( i ) / sepin( i )
376 ELSE
377 vmax = one
378 END IF
379 IF( vmax.GT.rmax( 3 ) ) THEN
380 rmax( 3 ) = vmax
381 IF( ninfo( 3 ).EQ.0 )
382 $ lmax( 3 ) = knt
383 END IF
384 120 CONTINUE
385*
386* Compute eigenvalue condition numbers only and compare
387*
388 vmax = zero
389 dum( 1 ) = -one
390 CALL dcopy( n, dum, 0, stmp, 1 )
391 CALL dcopy( n, dum, 0, septmp, 1 )
392 CALL dtrsna( 'Eigcond', 'All', SELECT, n, t, ldt, le, ldt, re,
393 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
394 IF( info.NE.0 ) THEN
395 lmax( 3 ) = knt
396 ninfo( 3 ) = ninfo( 3 ) + 1
397 GO TO 240
398 END IF
399 DO 130 i = 1, n
400 IF( stmp( i ).NE.s( i ) )
401 $ vmax = one / eps
402 IF( septmp( i ).NE.dum( 1 ) )
403 $ vmax = one / eps
404 130 CONTINUE
405*
406* Compute eigenvector condition numbers only and compare
407*
408 CALL dcopy( n, dum, 0, stmp, 1 )
409 CALL dcopy( n, dum, 0, septmp, 1 )
410 CALL dtrsna( 'Veccond', 'All', SELECT, n, t, ldt, le, ldt, re,
411 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
412 IF( info.NE.0 ) THEN
413 lmax( 3 ) = knt
414 ninfo( 3 ) = ninfo( 3 ) + 1
415 GO TO 240
416 END IF
417 DO 140 i = 1, n
418 IF( stmp( i ).NE.dum( 1 ) )
419 $ vmax = one / eps
420 IF( septmp( i ).NE.sep( i ) )
421 $ vmax = one / eps
422 140 CONTINUE
423*
424* Compute all condition numbers using SELECT and compare
425*
426 DO 150 i = 1, n
427 SELECT( i ) = .true.
428 150 CONTINUE
429 CALL dcopy( n, dum, 0, stmp, 1 )
430 CALL dcopy( n, dum, 0, septmp, 1 )
431 CALL dtrsna( 'Bothcond', 'Some', SELECT, n, t, ldt, le, ldt,
432 $ re, ldt, stmp, septmp, n, m, work, n, iwork,
433 $ info )
434 IF( info.NE.0 ) THEN
435 lmax( 3 ) = knt
436 ninfo( 3 ) = ninfo( 3 ) + 1
437 GO TO 240
438 END IF
439 DO 160 i = 1, n
440 IF( septmp( i ).NE.sep( i ) )
441 $ vmax = one / eps
442 IF( stmp( i ).NE.s( i ) )
443 $ vmax = one / eps
444 160 CONTINUE
445*
446* Compute eigenvalue condition numbers using SELECT and compare
447*
448 CALL dcopy( n, dum, 0, stmp, 1 )
449 CALL dcopy( n, dum, 0, septmp, 1 )
450 CALL dtrsna( 'Eigcond', 'Some', SELECT, n, t, ldt, le, ldt, re,
451 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
452 IF( info.NE.0 ) THEN
453 lmax( 3 ) = knt
454 ninfo( 3 ) = ninfo( 3 ) + 1
455 GO TO 240
456 END IF
457 DO 170 i = 1, n
458 IF( stmp( i ).NE.s( i ) )
459 $ vmax = one / eps
460 IF( septmp( i ).NE.dum( 1 ) )
461 $ vmax = one / eps
462 170 CONTINUE
463*
464* Compute eigenvector condition numbers using SELECT and compare
465*
466 CALL dcopy( n, dum, 0, stmp, 1 )
467 CALL dcopy( n, dum, 0, septmp, 1 )
468 CALL dtrsna( 'Veccond', 'Some', SELECT, n, t, ldt, le, ldt, re,
469 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
470 IF( info.NE.0 ) THEN
471 lmax( 3 ) = knt
472 ninfo( 3 ) = ninfo( 3 ) + 1
473 GO TO 240
474 END IF
475 DO 180 i = 1, n
476 IF( stmp( i ).NE.dum( 1 ) )
477 $ vmax = one / eps
478 IF( septmp( i ).NE.sep( i ) )
479 $ vmax = one / eps
480 180 CONTINUE
481 IF( vmax.GT.rmax( 1 ) ) THEN
482 rmax( 1 ) = vmax
483 IF( ninfo( 1 ).EQ.0 )
484 $ lmax( 1 ) = knt
485 END IF
486*
487* Select first real and first complex eigenvalue
488*
489 IF( wi( 1 ).EQ.zero ) THEN
490 lcmp( 1 ) = 1
491 ifnd = 0
492 DO 190 i = 2, n
493 IF( ifnd.EQ.1 .OR. wi( i ).EQ.zero ) THEN
494 SELECT( i ) = .false.
495 ELSE
496 ifnd = 1
497 lcmp( 2 ) = i
498 lcmp( 3 ) = i + 1
499 CALL dcopy( n, re( 1, i ), 1, re( 1, 2 ), 1 )
500 CALL dcopy( n, re( 1, i+1 ), 1, re( 1, 3 ), 1 )
501 CALL dcopy( n, le( 1, i ), 1, le( 1, 2 ), 1 )
502 CALL dcopy( n, le( 1, i+1 ), 1, le( 1, 3 ), 1 )
503 END IF
504 190 CONTINUE
505 IF( ifnd.EQ.0 ) THEN
506 icmp = 1
507 ELSE
508 icmp = 3
509 END IF
510 ELSE
511 lcmp( 1 ) = 1
512 lcmp( 2 ) = 2
513 ifnd = 0
514 DO 200 i = 3, n
515 IF( ifnd.EQ.1 .OR. wi( i ).NE.zero ) THEN
516 SELECT( i ) = .false.
517 ELSE
518 lcmp( 3 ) = i
519 ifnd = 1
520 CALL dcopy( n, re( 1, i ), 1, re( 1, 3 ), 1 )
521 CALL dcopy( n, le( 1, i ), 1, le( 1, 3 ), 1 )
522 END IF
523 200 CONTINUE
524 IF( ifnd.EQ.0 ) THEN
525 icmp = 2
526 ELSE
527 icmp = 3
528 END IF
529 END IF
530*
531* Compute all selected condition numbers
532*
533 CALL dcopy( icmp, dum, 0, stmp, 1 )
534 CALL dcopy( icmp, dum, 0, septmp, 1 )
535 CALL dtrsna( 'Bothcond', 'Some', SELECT, n, t, ldt, le, ldt,
536 $ re, ldt, stmp, septmp, n, m, work, n, iwork,
537 $ info )
538 IF( info.NE.0 ) THEN
539 lmax( 3 ) = knt
540 ninfo( 3 ) = ninfo( 3 ) + 1
541 GO TO 240
542 END IF
543 DO 210 i = 1, icmp
544 j = lcmp( i )
545 IF( septmp( i ).NE.sep( j ) )
546 $ vmax = one / eps
547 IF( stmp( i ).NE.s( j ) )
548 $ vmax = one / eps
549 210 CONTINUE
550*
551* Compute selected eigenvalue condition numbers
552*
553 CALL dcopy( icmp, dum, 0, stmp, 1 )
554 CALL dcopy( icmp, dum, 0, septmp, 1 )
555 CALL dtrsna( 'Eigcond', 'Some', SELECT, n, t, ldt, le, ldt, re,
556 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
557 IF( info.NE.0 ) THEN
558 lmax( 3 ) = knt
559 ninfo( 3 ) = ninfo( 3 ) + 1
560 GO TO 240
561 END IF
562 DO 220 i = 1, icmp
563 j = lcmp( i )
564 IF( stmp( i ).NE.s( j ) )
565 $ vmax = one / eps
566 IF( septmp( i ).NE.dum( 1 ) )
567 $ vmax = one / eps
568 220 CONTINUE
569*
570* Compute selected eigenvector condition numbers
571*
572 CALL dcopy( icmp, dum, 0, stmp, 1 )
573 CALL dcopy( icmp, dum, 0, septmp, 1 )
574 CALL dtrsna( 'Veccond', 'Some', SELECT, n, t, ldt, le, ldt, re,
575 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
576 IF( info.NE.0 ) THEN
577 lmax( 3 ) = knt
578 ninfo( 3 ) = ninfo( 3 ) + 1
579 GO TO 240
580 END IF
581 DO 230 i = 1, icmp
582 j = lcmp( i )
583 IF( stmp( i ).NE.dum( 1 ) )
584 $ vmax = one / eps
585 IF( septmp( i ).NE.sep( j ) )
586 $ vmax = one / eps
587 230 CONTINUE
588 IF( vmax.GT.rmax( 1 ) ) THEN
589 rmax( 1 ) = vmax
590 IF( ninfo( 1 ).EQ.0 )
591 $ lmax( 1 ) = knt
592 END IF
593 240 CONTINUE
594 GO TO 10
595*
596* End of DGET37
597*
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
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
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
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dtrevc
subroutine dtrevc(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, info)
DTREVC
Definition dtrevc.f:222
dtrsna
subroutine dtrsna(job, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, s, sep, mm, m, work, ldwork, iwork, info)
DTRSNA
Definition dtrsna.f:265
Here is the call graph for this function:
Here is the caller graph for this function: