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.```
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 *
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
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
Definition: zgehrd.f:169
subroutine zhst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RWORK, RESULT)
ZHST01
Definition: zhst01.f:142
subroutine ztrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, SEP, WORK, LWORK, INFO)
ZTRSEN
Definition: ztrsen.f:266
subroutine zunghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGHR
Definition: zunghr.f:128