LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011

Definition at line 92 of file dget37.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: