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

Definition at line 93 of file dget38.f.

93 *
94 * -- LAPACK test routine (version 3.4.0) --
95 * -- LAPACK is a software package provided by Univ. of Tennessee, --
96 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
97 * November 2011
98 *
99 * .. Scalar Arguments ..
100  INTEGER knt, nin
101 * ..
102 * .. Array Arguments ..
103  INTEGER lmax( 3 ), ninfo( 3 )
104  DOUBLE PRECISION rmax( 3 )
105 * ..
106 *
107 * =====================================================================
108 *
109 * .. Parameters ..
110  DOUBLE PRECISION zero, one, two
111  parameter ( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
112  DOUBLE PRECISION epsin
113  parameter ( epsin = 5.9605d-8 )
114  INTEGER ldt, lwork
115  parameter ( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
116  INTEGER liwork
117  parameter ( liwork = ldt*ldt )
118 * ..
119 * .. Local Scalars ..
120  INTEGER i, info, iscl, itmp, j, kmin, m, n, ndim
121  DOUBLE PRECISION bignum, eps, s, sep, sepin, septmp, sin,
122  $ smlnum, stmp, tnrm, tol, tolin, v, vimin, vmax,
123  $ vmul, vrmin
124 * ..
125 * .. Local Arrays ..
126  LOGICAL select( ldt )
127  INTEGER ipnt( ldt ), iselec( ldt ), iwork( liwork )
128  DOUBLE PRECISION q( ldt, ldt ), qsav( ldt, ldt ),
129  $ qtmp( ldt, ldt ), result( 2 ), t( ldt, ldt ),
130  $ tmp( ldt, ldt ), tsav( ldt, ldt ),
131  $ tsav1( ldt, ldt ), ttmp( ldt, ldt ), val( 3 ),
132  $ wi( ldt ), witmp( ldt ), work( lwork ),
133  $ wr( ldt ), wrtmp( ldt )
134 * ..
135 * .. External Functions ..
136  DOUBLE PRECISION dlamch, dlange
137  EXTERNAL dlamch, dlange
138 * ..
139 * .. External Subroutines ..
140  EXTERNAL dcopy, dgehrd, dhseqr, dhst01, dlabad, dlacpy,
141  $ dorghr, dscal, dtrsen
142 * ..
143 * .. Intrinsic Functions ..
144  INTRINSIC dble, max, sqrt
145 * ..
146 * .. Executable Statements ..
147 *
148  eps = dlamch( 'P' )
149  smlnum = dlamch( 'S' ) / eps
150  bignum = one / smlnum
151  CALL dlabad( smlnum, bignum )
152 *
153 * EPSIN = 2**(-24) = precision to which input data computed
154 *
155  eps = max( eps, epsin )
156  rmax( 1 ) = zero
157  rmax( 2 ) = zero
158  rmax( 3 ) = zero
159  lmax( 1 ) = 0
160  lmax( 2 ) = 0
161  lmax( 3 ) = 0
162  knt = 0
163  ninfo( 1 ) = 0
164  ninfo( 2 ) = 0
165  ninfo( 3 ) = 0
166 *
167  val( 1 ) = sqrt( smlnum )
168  val( 2 ) = one
169  val( 3 ) = sqrt( sqrt( bignum ) )
170 *
171 * Read input data until N=0. Assume input eigenvalues are sorted
172 * lexicographically (increasing by real part, then decreasing by
173 * imaginary part)
174 *
175  10 CONTINUE
176  READ( nin, fmt = * )n, ndim
177  IF( n.EQ.0 )
178  $ RETURN
179  READ( nin, fmt = * )( iselec( i ), i = 1, ndim )
180  DO 20 i = 1, n
181  READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
182  20 CONTINUE
183  READ( nin, fmt = * )sin, sepin
184 *
185  tnrm = dlange( 'M', n, n, tmp, ldt, work )
186  DO 160 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 30 i = 1, n
194  CALL dscal( n, vmul, t( 1, i ), 1 )
195  30 CONTINUE
196  IF( tnrm.EQ.zero )
197  $ vmul = one
198  CALL dlacpy( 'F', n, n, t, ldt, tsav, ldt )
199 *
200 * Compute Schur form
201 *
202  CALL dgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
203  $ info )
204  IF( info.NE.0 ) THEN
205  lmax( 1 ) = knt
206  ninfo( 1 ) = ninfo( 1 ) + 1
207  GO TO 160
208  END IF
209 *
210 * Generate orthogonal matrix
211 *
212  CALL dlacpy( 'L', n, n, t, ldt, q, ldt )
213  CALL dorghr( n, 1, n, q, ldt, work( 1 ), work( n+1 ), lwork-n,
214  $ info )
215 *
216 * Compute Schur form
217 *
218  CALL dhseqr( 'S', 'V', n, 1, n, t, ldt, wr, wi, q, ldt, work,
219  $ lwork, info )
220  IF( info.NE.0 ) THEN
221  lmax( 2 ) = knt
222  ninfo( 2 ) = ninfo( 2 ) + 1
223  GO TO 160
224  END IF
225 *
226 * Sort, select eigenvalues
227 *
228  DO 40 i = 1, n
229  ipnt( i ) = i
230  SELECT( i ) = .false.
231  40 CONTINUE
232  CALL dcopy( n, wr, 1, wrtmp, 1 )
233  CALL dcopy( n, wi, 1, witmp, 1 )
234  DO 60 i = 1, n - 1
235  kmin = i
236  vrmin = wrtmp( i )
237  vimin = witmp( i )
238  DO 50 j = i + 1, n
239  IF( wrtmp( j ).LT.vrmin ) THEN
240  kmin = j
241  vrmin = wrtmp( j )
242  vimin = witmp( j )
243  END IF
244  50 CONTINUE
245  wrtmp( kmin ) = wrtmp( i )
246  witmp( kmin ) = witmp( i )
247  wrtmp( i ) = vrmin
248  witmp( i ) = vimin
249  itmp = ipnt( i )
250  ipnt( i ) = ipnt( kmin )
251  ipnt( kmin ) = itmp
252  60 CONTINUE
253  DO 70 i = 1, ndim
254  SELECT( ipnt( iselec( i ) ) ) = .true.
255  70 CONTINUE
256 *
257 * Compute condition numbers
258 *
259  CALL dlacpy( 'F', n, n, q, ldt, qsav, ldt )
260  CALL dlacpy( 'F', n, n, t, ldt, tsav1, ldt )
261  CALL dtrsen( 'B', 'V', SELECT, n, t, ldt, q, ldt, wrtmp, witmp,
262  $ m, s, sep, work, lwork, iwork, liwork, info )
263  IF( info.NE.0 ) THEN
264  lmax( 3 ) = knt
265  ninfo( 3 ) = ninfo( 3 ) + 1
266  GO TO 160
267  END IF
268  septmp = sep / vmul
269  stmp = s
270 *
271 * Compute residuals
272 *
273  CALL dhst01( n, 1, n, tsav, ldt, t, ldt, q, ldt, work, lwork,
274  $ result )
275  vmax = max( result( 1 ), result( 2 ) )
276  IF( vmax.GT.rmax( 1 ) ) THEN
277  rmax( 1 ) = vmax
278  IF( ninfo( 1 ).EQ.0 )
279  $ lmax( 1 ) = knt
280  END IF
281 *
282 * Compare condition number for eigenvalue cluster
283 * taking its condition number into account
284 *
285  v = max( two*dble( n )*eps*tnrm, smlnum )
286  IF( tnrm.EQ.zero )
287  $ v = one
288  IF( v.GT.septmp ) THEN
289  tol = one
290  ELSE
291  tol = v / septmp
292  END IF
293  IF( v.GT.sepin ) THEN
294  tolin = one
295  ELSE
296  tolin = v / sepin
297  END IF
298  tol = max( tol, smlnum / eps )
299  tolin = max( tolin, smlnum / eps )
300  IF( eps*( sin-tolin ).GT.stmp+tol ) THEN
301  vmax = one / eps
302  ELSE IF( sin-tolin.GT.stmp+tol ) THEN
303  vmax = ( sin-tolin ) / ( stmp+tol )
304  ELSE IF( sin+tolin.LT.eps*( stmp-tol ) ) THEN
305  vmax = one / eps
306  ELSE IF( sin+tolin.LT.stmp-tol ) THEN
307  vmax = ( stmp-tol ) / ( sin+tolin )
308  ELSE
309  vmax = one
310  END IF
311  IF( vmax.GT.rmax( 2 ) ) THEN
312  rmax( 2 ) = vmax
313  IF( ninfo( 2 ).EQ.0 )
314  $ lmax( 2 ) = knt
315  END IF
316 *
317 * Compare condition numbers for invariant subspace
318 * taking its condition number into account
319 *
320  IF( v.GT.septmp*stmp ) THEN
321  tol = septmp
322  ELSE
323  tol = v / stmp
324  END IF
325  IF( v.GT.sepin*sin ) THEN
326  tolin = sepin
327  ELSE
328  tolin = v / sin
329  END IF
330  tol = max( tol, smlnum / eps )
331  tolin = max( tolin, smlnum / eps )
332  IF( eps*( sepin-tolin ).GT.septmp+tol ) THEN
333  vmax = one / eps
334  ELSE IF( sepin-tolin.GT.septmp+tol ) THEN
335  vmax = ( sepin-tolin ) / ( septmp+tol )
336  ELSE IF( sepin+tolin.LT.eps*( septmp-tol ) ) THEN
337  vmax = one / eps
338  ELSE IF( sepin+tolin.LT.septmp-tol ) THEN
339  vmax = ( septmp-tol ) / ( sepin+tolin )
340  ELSE
341  vmax = one
342  END IF
343  IF( vmax.GT.rmax( 2 ) ) THEN
344  rmax( 2 ) = vmax
345  IF( ninfo( 2 ).EQ.0 )
346  $ lmax( 2 ) = knt
347  END IF
348 *
349 * Compare condition number for eigenvalue cluster
350 * without taking its condition number into account
351 *
352  IF( sin.LE.dble( 2*n )*eps .AND. stmp.LE.dble( 2*n )*eps ) THEN
353  vmax = one
354  ELSE IF( eps*sin.GT.stmp ) THEN
355  vmax = one / eps
356  ELSE IF( sin.GT.stmp ) THEN
357  vmax = sin / stmp
358  ELSE IF( sin.LT.eps*stmp ) THEN
359  vmax = one / eps
360  ELSE IF( sin.LT.stmp ) THEN
361  vmax = stmp / sin
362  ELSE
363  vmax = one
364  END IF
365  IF( vmax.GT.rmax( 3 ) ) THEN
366  rmax( 3 ) = vmax
367  IF( ninfo( 3 ).EQ.0 )
368  $ lmax( 3 ) = knt
369  END IF
370 *
371 * Compare condition numbers for invariant subspace
372 * without taking its condition number into account
373 *
374  IF( sepin.LE.v .AND. septmp.LE.v ) THEN
375  vmax = one
376  ELSE IF( eps*sepin.GT.septmp ) THEN
377  vmax = one / eps
378  ELSE IF( sepin.GT.septmp ) THEN
379  vmax = sepin / septmp
380  ELSE IF( sepin.LT.eps*septmp ) THEN
381  vmax = one / eps
382  ELSE IF( sepin.LT.septmp ) THEN
383  vmax = septmp / sepin
384  ELSE
385  vmax = one
386  END IF
387  IF( vmax.GT.rmax( 3 ) ) THEN
388  rmax( 3 ) = vmax
389  IF( ninfo( 3 ).EQ.0 )
390  $ lmax( 3 ) = knt
391  END IF
392 *
393 * Compute eigenvalue condition number only and compare
394 * Update Q
395 *
396  vmax = zero
397  CALL dlacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
398  CALL dlacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
399  septmp = -one
400  stmp = -one
401  CALL dtrsen( 'E', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
402  $ witmp, m, stmp, septmp, work, lwork, iwork,
403  $ liwork, info )
404  IF( info.NE.0 ) THEN
405  lmax( 3 ) = knt
406  ninfo( 3 ) = ninfo( 3 ) + 1
407  GO TO 160
408  END IF
409  IF( s.NE.stmp )
410  $ vmax = one / eps
411  IF( -one.NE.septmp )
412  $ vmax = one / eps
413  DO 90 i = 1, n
414  DO 80 j = 1, n
415  IF( ttmp( i, j ).NE.t( i, j ) )
416  $ vmax = one / eps
417  IF( qtmp( i, j ).NE.q( i, j ) )
418  $ vmax = one / eps
419  80 CONTINUE
420  90 CONTINUE
421 *
422 * Compute invariant subspace condition number only and compare
423 * Update Q
424 *
425  CALL dlacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
426  CALL dlacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
427  septmp = -one
428  stmp = -one
429  CALL dtrsen( 'V', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
430  $ witmp, m, stmp, septmp, work, lwork, iwork,
431  $ liwork, info )
432  IF( info.NE.0 ) THEN
433  lmax( 3 ) = knt
434  ninfo( 3 ) = ninfo( 3 ) + 1
435  GO TO 160
436  END IF
437  IF( -one.NE.stmp )
438  $ vmax = one / eps
439  IF( sep.NE.septmp )
440  $ vmax = one / eps
441  DO 110 i = 1, n
442  DO 100 j = 1, n
443  IF( ttmp( i, j ).NE.t( i, j ) )
444  $ vmax = one / eps
445  IF( qtmp( i, j ).NE.q( i, j ) )
446  $ vmax = one / eps
447  100 CONTINUE
448  110 CONTINUE
449 *
450 * Compute eigenvalue condition number only and compare
451 * Do not update Q
452 *
453  CALL dlacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
454  CALL dlacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
455  septmp = -one
456  stmp = -one
457  CALL dtrsen( 'E', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
458  $ witmp, m, stmp, septmp, work, lwork, iwork,
459  $ liwork, info )
460  IF( info.NE.0 ) THEN
461  lmax( 3 ) = knt
462  ninfo( 3 ) = ninfo( 3 ) + 1
463  GO TO 160
464  END IF
465  IF( s.NE.stmp )
466  $ vmax = one / eps
467  IF( -one.NE.septmp )
468  $ vmax = one / eps
469  DO 130 i = 1, n
470  DO 120 j = 1, n
471  IF( ttmp( i, j ).NE.t( i, j ) )
472  $ vmax = one / eps
473  IF( qtmp( i, j ).NE.qsav( i, j ) )
474  $ vmax = one / eps
475  120 CONTINUE
476  130 CONTINUE
477 *
478 * Compute invariant subspace condition number only and compare
479 * Do not update Q
480 *
481  CALL dlacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
482  CALL dlacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
483  septmp = -one
484  stmp = -one
485  CALL dtrsen( 'V', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wrtmp,
486  $ witmp, m, stmp, septmp, work, lwork, iwork,
487  $ liwork, info )
488  IF( info.NE.0 ) THEN
489  lmax( 3 ) = knt
490  ninfo( 3 ) = ninfo( 3 ) + 1
491  GO TO 160
492  END IF
493  IF( -one.NE.stmp )
494  $ vmax = one / eps
495  IF( sep.NE.septmp )
496  $ vmax = one / eps
497  DO 150 i = 1, n
498  DO 140 j = 1, n
499  IF( ttmp( i, j ).NE.t( i, j ) )
500  $ vmax = one / eps
501  IF( qtmp( i, j ).NE.qsav( i, j ) )
502  $ vmax = one / eps
503  140 CONTINUE
504  150 CONTINUE
505  IF( vmax.GT.rmax( 1 ) ) THEN
506  rmax( 1 ) = vmax
507  IF( ninfo( 1 ).EQ.0 )
508  $ lmax( 1 ) = knt
509  END IF
510  160 CONTINUE
511  GO TO 10
512 *
513 * End of DGET38
514 *
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
dhst01
subroutine dhst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
DHST01
Definition: dhst01.f:136
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
dorghr
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
Definition: dorghr.f:128
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
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:315
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: