1*> \brief \b DGET38
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 DGET38( 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*> DGET38 tests DTRSEN, 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 DHST01 or comparing
41*>                    different calls to DTRSEN
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 DGEHRD returns INFO nonzero on example i, LMAX(1)=i
55*>          If DHSEQR returns INFO nonzero on example i, LMAX(2)=i
56*>          If DTRSEN 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 DGEHRD returned INFO nonzero
63*>          NINFO(2) = No. of times DHSEQR returned INFO nonzero
64*>          NINFO(3) = No. of times DTRSEN 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 double_eig
88*
89*  =====================================================================
90      SUBROUTINE DGET38( 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      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*
512      END
513