LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dhst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
DHST01
Definition: dhst01.f:134
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
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
Definition: dgehrd.f:167
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
Definition: dhseqr.f:316
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
Definition: dorghr.f:126
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: