1*> \brief \b ZGET38
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*  Definition:
9*  ===========
10*
11*       SUBROUTINE ZGET38( RMAX, LMAX, NINFO, KNT, NIN )
12*
13*       .. Scalar Arguments ..
14*       INTEGER            KNT, NIN
15*       ..
16*       .. Array Arguments ..
17*       INTEGER            LMAX( 3 ), NINFO( 3 )
18*       DOUBLE PRECISION   RMAX( 3 )
19*       ..
20*
21*
22*> \par Purpose:
23*  =============
24*>
25*> \verbatim
26*>
27*> ZGET38 tests ZTRSEN, a routine for estimating condition numbers of a
28*> cluster of eigenvalues and/or its associated right invariant subspace
29*>
30*> The test matrices are read from a file with logical unit number NIN.
31*> \endverbatim
32*
33*  Arguments:
34*  ==========
35*
36*> \param[out] RMAX
37*> \verbatim
38*>          RMAX is DOUBLE PRECISION array, dimension (3)
39*>          Values of the largest test ratios.
40*>          RMAX(1) = largest residuals from ZHST01 or comparing
41*>                    different calls to ZTRSEN
42*>          RMAX(2) = largest error in reciprocal condition
43*>                    numbers taking their conditioning into account
44*>          RMAX(3) = largest error in reciprocal condition
45*>                    numbers not taking their conditioning into
46*>                    account (may be larger than RMAX(2))
47*> \endverbatim
48*>
49*> \param[out] LMAX
50*> \verbatim
51*>          LMAX is INTEGER array, dimension (3)
52*>          LMAX(i) is example number where largest test ratio
53*>          RMAX(i) is achieved. Also:
54*>          If ZGEHRD returns INFO nonzero on example i, LMAX(1)=i
55*>          If ZHSEQR returns INFO nonzero on example i, LMAX(2)=i
56*>          If ZTRSEN returns INFO nonzero on example i, LMAX(3)=i
57*> \endverbatim
58*>
59*> \param[out] NINFO
60*> \verbatim
61*>          NINFO is INTEGER array, dimension (3)
62*>          NINFO(1) = No. of times ZGEHRD returned INFO nonzero
63*>          NINFO(2) = No. of times ZHSEQR returned INFO nonzero
64*>          NINFO(3) = No. of times ZTRSEN returned INFO nonzero
65*> \endverbatim
66*>
67*> \param[out] KNT
68*> \verbatim
69*>          KNT is INTEGER
70*>          Total number of examples tested.
71*> \endverbatim
72*>
73*> \param[in] NIN
74*> \verbatim
75*>          NIN is INTEGER
76*>          Input logical unit number.
77*> \endverbatim
78*
79*  Authors:
80*  ========
81*
82*> \author Univ. of Tennessee
83*> \author Univ. of California Berkeley
84*> \author Univ. of Colorado Denver
85*> \author NAG Ltd.
86*
87*> \ingroup complex16_eig
88*
89*  =====================================================================
90      SUBROUTINE ZGET38( RMAX, LMAX, NINFO, KNT, NIN )
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      INTEGER            LDT, LWORK
108      PARAMETER          ( LDT = 20, LWORK = 2*LDT*( 10+LDT ) )
109      DOUBLE PRECISION   ZERO, ONE, TWO
110      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
111      DOUBLE PRECISION   EPSIN
112      PARAMETER          ( EPSIN = 5.9605D-8 )
113      COMPLEX*16         CZERO
114      PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
115*     ..
116*     .. Local Scalars ..
117      INTEGER            I, INFO, ISCL, ISRT, ITMP, J, KMIN, M, N, NDIM
118      DOUBLE PRECISION   BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN,
119     $                   SMLNUM, STMP, TNRM, TOL, TOLIN, V, VMAX, VMIN,
120     $                   VMUL
121*     ..
122*     .. Local Arrays ..
123      LOGICAL            SELECT( LDT )
124      INTEGER            IPNT( LDT ), ISELEC( LDT )
125      DOUBLE PRECISION   RESULT( 2 ), RWORK( LDT ), VAL( 3 ),
126     $                   WSRT( LDT )
127      COMPLEX*16         Q( LDT, LDT ), QSAV( LDT, LDT ),
128     $                   QTMP( LDT, LDT ), T( LDT, LDT ),
129     $                   TMP( LDT, LDT ), TSAV( LDT, LDT ),
130     $                   TSAV1( LDT, LDT ), TTMP( LDT, LDT ), W( LDT ),
131     $                   WORK( LWORK ), WTMP( LDT )
132*     ..
133*     .. External Functions ..
134      DOUBLE PRECISION   DLAMCH, ZLANGE
135      EXTERNAL           DLAMCH, ZLANGE
136*     ..
137*     .. External Subroutines ..
138      EXTERNAL           DLABAD, ZDSCAL, ZGEHRD, ZHSEQR, ZHST01, ZLACPY,
139     $                   ZTRSEN, ZUNGHR
140*     ..
141*     .. Intrinsic Functions ..
142      INTRINSIC          DBLE, DIMAG, MAX, SQRT
143*     ..
144*     .. Executable Statements ..
145*
146      EPS = DLAMCH( 'P' )
147      SMLNUM = DLAMCH( 'S' ) / EPS
148      BIGNUM = ONE / SMLNUM
149      CALL DLABAD( SMLNUM, BIGNUM )
150*
151*     EPSIN = 2**(-24) = precision to which input data computed
152*
153      EPS = MAX( EPS, EPSIN )
154      RMAX( 1 ) = ZERO
155      RMAX( 2 ) = ZERO
156      RMAX( 3 ) = ZERO
157      LMAX( 1 ) = 0
158      LMAX( 2 ) = 0
159      LMAX( 3 ) = 0
160      KNT = 0
161      NINFO( 1 ) = 0
162      NINFO( 2 ) = 0
163      NINFO( 3 ) = 0
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, ISRT
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 = ZLANGE( 'M', N, N, TMP, LDT, RWORK )
183      DO 200 ISCL = 1, 3
184*
185*        Scale input matrix
186*
187         KNT = KNT + 1
188         CALL ZLACPY( 'F', N, N, TMP, LDT, T, LDT )
189         VMUL = VAL( ISCL )
190         DO 30 I = 1, N
191            CALL ZDSCAL( N, VMUL, T( 1, I ), 1 )
192   30    CONTINUE
193         IF( TNRM.EQ.ZERO )
194     $      VMUL = ONE
195         CALL ZLACPY( 'F', N, N, T, LDT, TSAV, LDT )
196*
197*        Compute Schur form
198*
199         CALL ZGEHRD( 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 200
205         END IF
206*
207*        Generate unitary matrix
208*
209         CALL ZLACPY( 'L', N, N, T, LDT, Q, LDT )
210         CALL ZUNGHR( N, 1, N, Q, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
211     $                INFO )
212*
213*        Compute Schur form
214*
215         DO 50 J = 1, N - 2
216            DO 40 I = J + 2, N
217               T( I, J ) = CZERO
218   40       CONTINUE
219   50    CONTINUE
220         CALL ZHSEQR( 'S', 'V', N, 1, N, T, LDT, W, Q, LDT, WORK, LWORK,
221     $                INFO )
222         IF( INFO.NE.0 ) THEN
223            LMAX( 2 ) = KNT
224            NINFO( 2 ) = NINFO( 2 ) + 1
225            GO TO 200
226         END IF
227*
228*        Sort, select eigenvalues
229*
230         DO 60 I = 1, N
231            IPNT( I ) = I
232            SELECT( I ) = .FALSE.
233   60    CONTINUE
234         IF( ISRT.EQ.0 ) THEN
235            DO 70 I = 1, N
236               WSRT( I ) = DBLE( W( I ) )
237   70       CONTINUE
238         ELSE
239            DO 80 I = 1, N
240               WSRT( I ) = DIMAG( W( I ) )
241   80       CONTINUE
242         END IF
243         DO 100 I = 1, N - 1
244            KMIN = I
245            VMIN = WSRT( I )
246            DO 90 J = I + 1, N
247               IF( WSRT( J ).LT.VMIN ) THEN
248                  KMIN = J
249                  VMIN = WSRT( J )
250               END IF
251   90       CONTINUE
252            WSRT( KMIN ) = WSRT( I )
253            WSRT( I ) = VMIN
254            ITMP = IPNT( I )
255            IPNT( I ) = IPNT( KMIN )
256            IPNT( KMIN ) = ITMP
257  100    CONTINUE
258         DO 110 I = 1, NDIM
259            SELECT( IPNT( ISELEC( I ) ) ) = .TRUE.
260  110    CONTINUE
261*
262*        Compute condition numbers
263*
264         CALL ZLACPY( 'F', N, N, Q, LDT, QSAV, LDT )
265         CALL ZLACPY( 'F', N, N, T, LDT, TSAV1, LDT )
266         CALL ZTRSEN( 'B', 'V', SELECT, N, T, LDT, Q, LDT, WTMP, M, S,
267     $                SEP, WORK, LWORK, INFO )
268         IF( INFO.NE.0 ) THEN
269            LMAX( 3 ) = KNT
270            NINFO( 3 ) = NINFO( 3 ) + 1
271            GO TO 200
272         END IF
273         SEPTMP = SEP / VMUL
274         STMP = S
275*
276*        Compute residuals
277*
278         CALL ZHST01( N, 1, N, TSAV, LDT, T, LDT, Q, LDT, WORK, LWORK,
279     $                RWORK, RESULT )
280         VMAX = MAX( RESULT( 1 ), RESULT( 2 ) )
281         IF( VMAX.GT.RMAX( 1 ) ) THEN
282            RMAX( 1 ) = VMAX
283            IF( NINFO( 1 ).EQ.0 )
284     $         LMAX( 1 ) = KNT
285         END IF
286*
287*        Compare condition number for eigenvalue cluster
288*        taking its condition number into account
289*
290         V = MAX( TWO*DBLE( N )*EPS*TNRM, SMLNUM )
291         IF( TNRM.EQ.ZERO )
292     $      V = ONE
293         IF( V.GT.SEPTMP ) THEN
294            TOL = ONE
295         ELSE
296            TOL = V / SEPTMP
297         END IF
298         IF( V.GT.SEPIN ) THEN
299            TOLIN = ONE
300         ELSE
301            TOLIN = V / SEPIN
302         END IF
303         TOL = MAX( TOL, SMLNUM / EPS )
304         TOLIN = MAX( TOLIN, SMLNUM / EPS )
305         IF( EPS*( SIN-TOLIN ).GT.STMP+TOL ) THEN
306            VMAX = ONE / EPS
307         ELSE IF( SIN-TOLIN.GT.STMP+TOL ) THEN
308            VMAX = ( SIN-TOLIN ) / ( STMP+TOL )
309         ELSE IF( SIN+TOLIN.LT.EPS*( STMP-TOL ) ) THEN
310            VMAX = ONE / EPS
311         ELSE IF( SIN+TOLIN.LT.STMP-TOL ) THEN
312            VMAX = ( STMP-TOL ) / ( SIN+TOLIN )
313         ELSE
314            VMAX = ONE
315         END IF
316         IF( VMAX.GT.RMAX( 2 ) ) THEN
317            RMAX( 2 ) = VMAX
318            IF( NINFO( 2 ).EQ.0 )
319     $         LMAX( 2 ) = KNT
320         END IF
321*
322*        Compare condition numbers for invariant subspace
323*        taking its condition number into account
324*
325         IF( V.GT.SEPTMP*STMP ) THEN
326            TOL = SEPTMP
327         ELSE
328            TOL = V / STMP
329         END IF
330         IF( V.GT.SEPIN*SIN ) THEN
331            TOLIN = SEPIN
332         ELSE
333            TOLIN = V / SIN
334         END IF
335         TOL = MAX( TOL, SMLNUM / EPS )
336         TOLIN = MAX( TOLIN, SMLNUM / EPS )
337         IF( EPS*( SEPIN-TOLIN ).GT.SEPTMP+TOL ) THEN
338            VMAX = ONE / EPS
339         ELSE IF( SEPIN-TOLIN.GT.SEPTMP+TOL ) THEN
340            VMAX = ( SEPIN-TOLIN ) / ( SEPTMP+TOL )
341         ELSE IF( SEPIN+TOLIN.LT.EPS*( SEPTMP-TOL ) ) THEN
342            VMAX = ONE / EPS
343         ELSE IF( SEPIN+TOLIN.LT.SEPTMP-TOL ) THEN
344            VMAX = ( SEPTMP-TOL ) / ( SEPIN+TOLIN )
345         ELSE
346            VMAX = ONE
347         END IF
348         IF( VMAX.GT.RMAX( 2 ) ) THEN
349            RMAX( 2 ) = VMAX
350            IF( NINFO( 2 ).EQ.0 )
351     $         LMAX( 2 ) = KNT
352         END IF
353*
354*        Compare condition number for eigenvalue cluster
355*        without taking its condition number into account
356*
357         IF( SIN.LE.DBLE( 2*N )*EPS .AND. STMP.LE.DBLE( 2*N )*EPS ) THEN
358            VMAX = ONE
359         ELSE IF( EPS*SIN.GT.STMP ) THEN
360            VMAX = ONE / EPS
361         ELSE IF( SIN.GT.STMP ) THEN
362            VMAX = SIN / STMP
363         ELSE IF( SIN.LT.EPS*STMP ) THEN
364            VMAX = ONE / EPS
365         ELSE IF( SIN.LT.STMP ) THEN
366            VMAX = STMP / SIN
367         ELSE
368            VMAX = ONE
369         END IF
370         IF( VMAX.GT.RMAX( 3 ) ) THEN
371            RMAX( 3 ) = VMAX
372            IF( NINFO( 3 ).EQ.0 )
373     $         LMAX( 3 ) = KNT
374         END IF
375*
376*        Compare condition numbers for invariant subspace
377*        without taking its condition number into account
378*
379         IF( SEPIN.LE.V .AND. SEPTMP.LE.V ) THEN
380            VMAX = ONE
381         ELSE IF( EPS*SEPIN.GT.SEPTMP ) THEN
382            VMAX = ONE / EPS
383         ELSE IF( SEPIN.GT.SEPTMP ) THEN
384            VMAX = SEPIN / SEPTMP
385         ELSE IF( SEPIN.LT.EPS*SEPTMP ) THEN
386            VMAX = ONE / EPS
387         ELSE IF( SEPIN.LT.SEPTMP ) THEN
388            VMAX = SEPTMP / SEPIN
389         ELSE
390            VMAX = ONE
391         END IF
392         IF( VMAX.GT.RMAX( 3 ) ) THEN
393            RMAX( 3 ) = VMAX
394            IF( NINFO( 3 ).EQ.0 )
395     $         LMAX( 3 ) = KNT
396         END IF
397*
398*        Compute eigenvalue condition number only and compare
399*        Update Q
400*
401         VMAX = ZERO
402         CALL ZLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
403         CALL ZLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
404         SEPTMP = -ONE
405         STMP = -ONE
406         CALL ZTRSEN( 'E', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP,
407     $                M, STMP, SEPTMP, WORK, LWORK, INFO )
408         IF( INFO.NE.0 ) THEN
409            LMAX( 3 ) = KNT
410            NINFO( 3 ) = NINFO( 3 ) + 1
411            GO TO 200
412         END IF
413         IF( S.NE.STMP )
414     $      VMAX = ONE / EPS
415         IF( -ONE.NE.SEPTMP )
416     $      VMAX = ONE / EPS
417         DO 130 I = 1, N
418            DO 120 J = 1, N
419               IF( TTMP( I, J ).NE.T( I, J ) )
420     $            VMAX = ONE / EPS
421               IF( QTMP( I, J ).NE.Q( I, J ) )
422     $            VMAX = ONE / EPS
423  120       CONTINUE
424  130    CONTINUE
425*
426*        Compute invariant subspace condition number only and compare
427*        Update Q
428*
429         CALL ZLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
430         CALL ZLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
431         SEPTMP = -ONE
432         STMP = -ONE
433         CALL ZTRSEN( 'V', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP,
434     $                M, STMP, SEPTMP, WORK, LWORK, INFO )
435         IF( INFO.NE.0 ) THEN
436            LMAX( 3 ) = KNT
437            NINFO( 3 ) = NINFO( 3 ) + 1
438            GO TO 200
439         END IF
440         IF( -ONE.NE.STMP )
441     $      VMAX = ONE / EPS
442         IF( SEP.NE.SEPTMP )
443     $      VMAX = ONE / EPS
444         DO 150 I = 1, N
445            DO 140 J = 1, N
446               IF( TTMP( I, J ).NE.T( I, J ) )
447     $            VMAX = ONE / EPS
448               IF( QTMP( I, J ).NE.Q( I, J ) )
449     $            VMAX = ONE / EPS
450  140       CONTINUE
451  150    CONTINUE
452*
453*        Compute eigenvalue condition number only and compare
454*        Do not update Q
455*
456         CALL ZLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
457         CALL ZLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
458         SEPTMP = -ONE
459         STMP = -ONE
460         CALL ZTRSEN( 'E', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP,
461     $                M, STMP, SEPTMP, WORK, LWORK, INFO )
462         IF( INFO.NE.0 ) THEN
463            LMAX( 3 ) = KNT
464            NINFO( 3 ) = NINFO( 3 ) + 1
465            GO TO 200
466         END IF
467         IF( S.NE.STMP )
468     $      VMAX = ONE / EPS
469         IF( -ONE.NE.SEPTMP )
470     $      VMAX = ONE / EPS
471         DO 170 I = 1, N
472            DO 160 J = 1, N
473               IF( TTMP( I, J ).NE.T( I, J ) )
474     $            VMAX = ONE / EPS
475               IF( QTMP( I, J ).NE.QSAV( I, J ) )
476     $            VMAX = ONE / EPS
477  160       CONTINUE
478  170    CONTINUE
479*
480*        Compute invariant subspace condition number only and compare
481*        Do not update Q
482*
483         CALL ZLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
484         CALL ZLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
485         SEPTMP = -ONE
486         STMP = -ONE
487         CALL ZTRSEN( 'V', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP,
488     $                M, STMP, SEPTMP, WORK, LWORK, INFO )
489         IF( INFO.NE.0 ) THEN
490            LMAX( 3 ) = KNT
491            NINFO( 3 ) = NINFO( 3 ) + 1
492            GO TO 200
493         END IF
494         IF( -ONE.NE.STMP )
495     $      VMAX = ONE / EPS
496         IF( SEP.NE.SEPTMP )
497     $      VMAX = ONE / EPS
498         DO 190 I = 1, N
499            DO 180 J = 1, N
500               IF( TTMP( I, J ).NE.T( I, J ) )
501     $            VMAX = ONE / EPS
502               IF( QTMP( I, J ).NE.QSAV( I, J ) )
503     $            VMAX = ONE / EPS
504  180       CONTINUE
505  190    CONTINUE
506         IF( VMAX.GT.RMAX( 1 ) ) THEN
507            RMAX( 1 ) = VMAX
508            IF( NINFO( 1 ).EQ.0 )
509     $         LMAX( 1 ) = KNT
510         END IF
511  200 CONTINUE
512      GO TO 10
513*
514*     End of ZGET38
515*
516      END
517