LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zget38 ( double precision, dimension( 3 )  RMAX,
integer, dimension( 3 )  LMAX,
integer, dimension( 3 )  NINFO,
integer  KNT,
integer  NIN 
)

ZGET38

Purpose: 
 ZGET38 tests ZTRSEN, 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 ZHST01 or comparing
                    different calls to ZTRSEN
          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 ZGEHRD returns INFO nonzero on example i, LMAX(1)=i
          If ZHSEQR returns INFO nonzero on example i, LMAX(2)=i
          If ZTRSEN returns INFO nonzero on example i, LMAX(3)=i
[out]NINFO
          NINFO is INTEGER array, dimension (3)
          NINFO(1) = No. of times ZGEHRD returned INFO nonzero
          NINFO(2) = No. of times ZHSEQR returned INFO nonzero
          NINFO(3) = No. of times ZTRSEN 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 zget38.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  INTEGER ldt, lwork
111  parameter ( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
112  DOUBLE PRECISION zero, one, two
113  parameter ( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
114  DOUBLE PRECISION epsin
115  parameter ( epsin = 5.9605d-8 )
116  COMPLEX*16 czero
117  parameter ( czero = ( 0.0d+0, 0.0d+0 ) )
118 * ..
119 * .. Local Scalars ..
120  INTEGER i, info, iscl, isrt, itmp, j, kmin, m, n, ndim
121  DOUBLE PRECISION bignum, eps, s, sep, sepin, septmp, sin,
122  $ smlnum, stmp, tnrm, tol, tolin, v, vmax, vmin,
123  $ vmul
124 * ..
125 * .. Local Arrays ..
126  LOGICAL select( ldt )
127  INTEGER ipnt( ldt ), iselec( ldt )
128  DOUBLE PRECISION result( 2 ), rwork( ldt ), val( 3 ),
129  $ wsrt( ldt )
130  COMPLEX*16 q( ldt, ldt ), qsav( ldt, ldt ),
131  $ qtmp( ldt, ldt ), t( ldt, ldt ),
132  $ tmp( ldt, ldt ), tsav( ldt, ldt ),
133  $ tsav1( ldt, ldt ), ttmp( ldt, ldt ), w( ldt ),
134  $ work( lwork ), wtmp( ldt )
135 * ..
136 * .. External Functions ..
137  DOUBLE PRECISION dlamch, zlange
138  EXTERNAL dlamch, zlange
139 * ..
140 * .. External Subroutines ..
141  EXTERNAL dlabad, zdscal, zgehrd, zhseqr, zhst01, zlacpy,
142  $ ztrsen, zunghr
143 * ..
144 * .. Intrinsic Functions ..
145  INTRINSIC dble, dimag, max, sqrt
146 * ..
147 * .. Executable Statements ..
148 *
149  eps = dlamch( 'P' )
150  smlnum = dlamch( 'S' ) / eps
151  bignum = one / smlnum
152  CALL dlabad( smlnum, bignum )
153 *
154 * EPSIN = 2**(-24) = precision to which input data computed
155 *
156  eps = max( eps, epsin )
157  rmax( 1 ) = zero
158  rmax( 2 ) = zero
159  rmax( 3 ) = zero
160  lmax( 1 ) = 0
161  lmax( 2 ) = 0
162  lmax( 3 ) = 0
163  knt = 0
164  ninfo( 1 ) = 0
165  ninfo( 2 ) = 0
166  ninfo( 3 ) = 0
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, isrt
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 = zlange( 'M', n, n, tmp, ldt, rwork )
186  DO 200 iscl = 1, 3
187 *
188 * Scale input matrix
189 *
190  knt = knt + 1
191  CALL zlacpy( 'F', n, n, tmp, ldt, t, ldt )
192  vmul = val( iscl )
193  DO 30 i = 1, n
194  CALL zdscal( n, vmul, t( 1, i ), 1 )
195  30 CONTINUE
196  IF( tnrm.EQ.zero )
197  $ vmul = one
198  CALL zlacpy( 'F', n, n, t, ldt, tsav, ldt )
199 *
200 * Compute Schur form
201 *
202  CALL zgehrd( 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 200
208  END IF
209 *
210 * Generate unitary matrix
211 *
212  CALL zlacpy( 'L', n, n, t, ldt, q, ldt )
213  CALL zunghr( n, 1, n, q, ldt, work( 1 ), work( n+1 ), lwork-n,
214  $ info )
215 *
216 * Compute Schur form
217 *
218  DO 50 j = 1, n - 2
219  DO 40 i = j + 2, n
220  t( i, j ) = czero
221  40 CONTINUE
222  50 CONTINUE
223  CALL zhseqr( 'S', 'V', n, 1, n, t, ldt, w, q, ldt, work, lwork,
224  $ info )
225  IF( info.NE.0 ) THEN
226  lmax( 2 ) = knt
227  ninfo( 2 ) = ninfo( 2 ) + 1
228  GO TO 200
229  END IF
230 *
231 * Sort, select eigenvalues
232 *
233  DO 60 i = 1, n
234  ipnt( i ) = i
235  SELECT( i ) = .false.
236  60 CONTINUE
237  IF( isrt.EQ.0 ) THEN
238  DO 70 i = 1, n
239  wsrt( i ) = dble( w( i ) )
240  70 CONTINUE
241  ELSE
242  DO 80 i = 1, n
243  wsrt( i ) = dimag( w( i ) )
244  80 CONTINUE
245  END IF
246  DO 100 i = 1, n - 1
247  kmin = i
248  vmin = wsrt( i )
249  DO 90 j = i + 1, n
250  IF( wsrt( j ).LT.vmin ) THEN
251  kmin = j
252  vmin = wsrt( j )
253  END IF
254  90 CONTINUE
255  wsrt( kmin ) = wsrt( i )
256  wsrt( i ) = vmin
257  itmp = ipnt( i )
258  ipnt( i ) = ipnt( kmin )
259  ipnt( kmin ) = itmp
260  100 CONTINUE
261  DO 110 i = 1, ndim
262  SELECT( ipnt( iselec( i ) ) ) = .true.
263  110 CONTINUE
264 *
265 * Compute condition numbers
266 *
267  CALL zlacpy( 'F', n, n, q, ldt, qsav, ldt )
268  CALL zlacpy( 'F', n, n, t, ldt, tsav1, ldt )
269  CALL ztrsen( 'B', 'V', SELECT, n, t, ldt, q, ldt, wtmp, m, s,
270  $ sep, work, lwork, info )
271  IF( info.NE.0 ) THEN
272  lmax( 3 ) = knt
273  ninfo( 3 ) = ninfo( 3 ) + 1
274  GO TO 200
275  END IF
276  septmp = sep / vmul
277  stmp = s
278 *
279 * Compute residuals
280 *
281  CALL zhst01( n, 1, n, tsav, ldt, t, ldt, q, ldt, work, lwork,
282  $ rwork, result )
283  vmax = max( result( 1 ), result( 2 ) )
284  IF( vmax.GT.rmax( 1 ) ) THEN
285  rmax( 1 ) = vmax
286  IF( ninfo( 1 ).EQ.0 )
287  $ lmax( 1 ) = knt
288  END IF
289 *
290 * Compare condition number for eigenvalue cluster
291 * taking its condition number into account
292 *
293  v = max( two*dble( n )*eps*tnrm, smlnum )
294  IF( tnrm.EQ.zero )
295  $ v = one
296  IF( v.GT.septmp ) THEN
297  tol = one
298  ELSE
299  tol = v / septmp
300  END IF
301  IF( v.GT.sepin ) THEN
302  tolin = one
303  ELSE
304  tolin = v / sepin
305  END IF
306  tol = max( tol, smlnum / eps )
307  tolin = max( tolin, smlnum / eps )
308  IF( eps*( sin-tolin ).GT.stmp+tol ) THEN
309  vmax = one / eps
310  ELSE IF( sin-tolin.GT.stmp+tol ) THEN
311  vmax = ( sin-tolin ) / ( stmp+tol )
312  ELSE IF( sin+tolin.LT.eps*( stmp-tol ) ) THEN
313  vmax = one / eps
314  ELSE IF( sin+tolin.LT.stmp-tol ) THEN
315  vmax = ( stmp-tol ) / ( sin+tolin )
316  ELSE
317  vmax = one
318  END IF
319  IF( vmax.GT.rmax( 2 ) ) THEN
320  rmax( 2 ) = vmax
321  IF( ninfo( 2 ).EQ.0 )
322  $ lmax( 2 ) = knt
323  END IF
324 *
325 * Compare condition numbers for invariant subspace
326 * taking its condition number into account
327 *
328  IF( v.GT.septmp*stmp ) THEN
329  tol = septmp
330  ELSE
331  tol = v / stmp
332  END IF
333  IF( v.GT.sepin*sin ) THEN
334  tolin = sepin
335  ELSE
336  tolin = v / sin
337  END IF
338  tol = max( tol, smlnum / eps )
339  tolin = max( tolin, smlnum / eps )
340  IF( eps*( sepin-tolin ).GT.septmp+tol ) THEN
341  vmax = one / eps
342  ELSE IF( sepin-tolin.GT.septmp+tol ) THEN
343  vmax = ( sepin-tolin ) / ( septmp+tol )
344  ELSE IF( sepin+tolin.LT.eps*( septmp-tol ) ) THEN
345  vmax = one / eps
346  ELSE IF( sepin+tolin.LT.septmp-tol ) THEN
347  vmax = ( septmp-tol ) / ( sepin+tolin )
348  ELSE
349  vmax = one
350  END IF
351  IF( vmax.GT.rmax( 2 ) ) THEN
352  rmax( 2 ) = vmax
353  IF( ninfo( 2 ).EQ.0 )
354  $ lmax( 2 ) = knt
355  END IF
356 *
357 * Compare condition number for eigenvalue cluster
358 * without taking its condition number into account
359 *
360  IF( sin.LE.dble( 2*n )*eps .AND. stmp.LE.dble( 2*n )*eps ) THEN
361  vmax = one
362  ELSE IF( eps*sin.GT.stmp ) THEN
363  vmax = one / eps
364  ELSE IF( sin.GT.stmp ) THEN
365  vmax = sin / stmp
366  ELSE IF( sin.LT.eps*stmp ) THEN
367  vmax = one / eps
368  ELSE IF( sin.LT.stmp ) THEN
369  vmax = stmp / sin
370  ELSE
371  vmax = one
372  END IF
373  IF( vmax.GT.rmax( 3 ) ) THEN
374  rmax( 3 ) = vmax
375  IF( ninfo( 3 ).EQ.0 )
376  $ lmax( 3 ) = knt
377  END IF
378 *
379 * Compare condition numbers for invariant subspace
380 * without taking its condition number into account
381 *
382  IF( sepin.LE.v .AND. septmp.LE.v ) THEN
383  vmax = one
384  ELSE IF( eps*sepin.GT.septmp ) THEN
385  vmax = one / eps
386  ELSE IF( sepin.GT.septmp ) THEN
387  vmax = sepin / septmp
388  ELSE IF( sepin.LT.eps*septmp ) THEN
389  vmax = one / eps
390  ELSE IF( sepin.LT.septmp ) THEN
391  vmax = septmp / sepin
392  ELSE
393  vmax = one
394  END IF
395  IF( vmax.GT.rmax( 3 ) ) THEN
396  rmax( 3 ) = vmax
397  IF( ninfo( 3 ).EQ.0 )
398  $ lmax( 3 ) = knt
399  END IF
400 *
401 * Compute eigenvalue condition number only and compare
402 * Update Q
403 *
404  vmax = zero
405  CALL zlacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
406  CALL zlacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
407  septmp = -one
408  stmp = -one
409  CALL ztrsen( 'E', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
410  $ m, stmp, septmp, work, lwork, info )
411  IF( info.NE.0 ) THEN
412  lmax( 3 ) = knt
413  ninfo( 3 ) = ninfo( 3 ) + 1
414  GO TO 200
415  END IF
416  IF( s.NE.stmp )
417  $ vmax = one / eps
418  IF( -one.NE.septmp )
419  $ vmax = one / eps
420  DO 130 i = 1, n
421  DO 120 j = 1, n
422  IF( ttmp( i, j ).NE.t( i, j ) )
423  $ vmax = one / eps
424  IF( qtmp( i, j ).NE.q( i, j ) )
425  $ vmax = one / eps
426  120 CONTINUE
427  130 CONTINUE
428 *
429 * Compute invariant subspace condition number only and compare
430 * Update Q
431 *
432  CALL zlacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
433  CALL zlacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
434  septmp = -one
435  stmp = -one
436  CALL ztrsen( 'V', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
437  $ m, stmp, septmp, work, lwork, info )
438  IF( info.NE.0 ) THEN
439  lmax( 3 ) = knt
440  ninfo( 3 ) = ninfo( 3 ) + 1
441  GO TO 200
442  END IF
443  IF( -one.NE.stmp )
444  $ vmax = one / eps
445  IF( sep.NE.septmp )
446  $ vmax = one / eps
447  DO 150 i = 1, n
448  DO 140 j = 1, n
449  IF( ttmp( i, j ).NE.t( i, j ) )
450  $ vmax = one / eps
451  IF( qtmp( i, j ).NE.q( i, j ) )
452  $ vmax = one / eps
453  140 CONTINUE
454  150 CONTINUE
455 *
456 * Compute eigenvalue condition number only and compare
457 * Do not update Q
458 *
459  CALL zlacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
460  CALL zlacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
461  septmp = -one
462  stmp = -one
463  CALL ztrsen( 'E', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
464  $ m, stmp, septmp, work, lwork, info )
465  IF( info.NE.0 ) THEN
466  lmax( 3 ) = knt
467  ninfo( 3 ) = ninfo( 3 ) + 1
468  GO TO 200
469  END IF
470  IF( s.NE.stmp )
471  $ vmax = one / eps
472  IF( -one.NE.septmp )
473  $ vmax = one / eps
474  DO 170 i = 1, n
475  DO 160 j = 1, n
476  IF( ttmp( i, j ).NE.t( i, j ) )
477  $ vmax = one / eps
478  IF( qtmp( i, j ).NE.qsav( i, j ) )
479  $ vmax = one / eps
480  160 CONTINUE
481  170 CONTINUE
482 *
483 * Compute invariant subspace condition number only and compare
484 * Do not update Q
485 *
486  CALL zlacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
487  CALL zlacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
488  septmp = -one
489  stmp = -one
490  CALL ztrsen( 'V', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
491  $ m, stmp, septmp, work, lwork, info )
492  IF( info.NE.0 ) THEN
493  lmax( 3 ) = knt
494  ninfo( 3 ) = ninfo( 3 ) + 1
495  GO TO 200
496  END IF
497  IF( -one.NE.stmp )
498  $ vmax = one / eps
499  IF( sep.NE.septmp )
500  $ vmax = one / eps
501  DO 190 i = 1, n
502  DO 180 j = 1, n
503  IF( ttmp( i, j ).NE.t( i, j ) )
504  $ vmax = one / eps
505  IF( qtmp( i, j ).NE.qsav( i, j ) )
506  $ vmax = one / eps
507  180 CONTINUE
508  190 CONTINUE
509  IF( vmax.GT.rmax( 1 ) ) THEN
510  rmax( 1 ) = vmax
511  IF( ninfo( 1 ).EQ.0 )
512  $ lmax( 1 ) = knt
513  END IF
514  200 CONTINUE
515  GO TO 10
516 *
517 * End of ZGET38
518 *
zlacpy
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
zgehrd
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
Definition: zgehrd.f:169
zhst01
subroutine zhst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RWORK, RESULT)
ZHST01
Definition: zhst01.f:142
ztrsen
subroutine ztrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, SEP, WORK, LWORK, INFO)
ZTRSEN
Definition: ztrsen.f:266
zunghr
subroutine zunghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGHR
Definition: zunghr.f:128
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
zlange
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
zhseqr
subroutine zhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
ZHSEQR
Definition: zhseqr.f:301
zdscal
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: