1*> \brief \b DCHKGK
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 DCHKGK( NIN, NOUT )
12*
13*       .. Scalar Arguments ..
14*       INTEGER            NIN, NOUT
15*       ..
16*
17*
18*> \par Purpose:
19*  =============
20*>
21*> \verbatim
22*>
23*> DCHKGK tests DGGBAK, a routine for backward balancing  of
24*> a matrix pair (A, B).
25*> \endverbatim
26*
27*  Arguments:
28*  ==========
29*
30*> \param[in] NIN
31*> \verbatim
32*>          NIN is INTEGER
33*>          The logical unit number for input.  NIN > 0.
34*> \endverbatim
35*>
36*> \param[in] NOUT
37*> \verbatim
38*>          NOUT is INTEGER
39*>          The logical unit number for output.  NOUT > 0.
40*> \endverbatim
41*
42*  Authors:
43*  ========
44*
45*> \author Univ. of Tennessee
46*> \author Univ. of California Berkeley
47*> \author Univ. of Colorado Denver
48*> \author NAG Ltd.
49*
50*> \date December 2016
51*
52*> \ingroup double_eig
53*
54*  =====================================================================
55      SUBROUTINE DCHKGK( NIN, NOUT )
56*
57*  -- LAPACK test routine (version 3.7.0) --
58*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
59*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
60*     December 2016
61*
62*     .. Scalar Arguments ..
63      INTEGER            NIN, NOUT
64*     ..
65*
66*  =====================================================================
67*
68*     .. Parameters ..
69      INTEGER            LDA, LDB, LDVL, LDVR
70      PARAMETER          ( LDA = 50, LDB = 50, LDVL = 50, LDVR = 50 )
71      INTEGER            LDE, LDF, LDWORK
72      PARAMETER          ( LDE = 50, LDF = 50, LDWORK = 50 )
73      DOUBLE PRECISION   ZERO, ONE
74      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
75*     ..
76*     .. Local Scalars ..
77      INTEGER            I, IHI, ILO, INFO, J, KNT, M, N, NINFO
78      DOUBLE PRECISION   ANORM, BNORM, EPS, RMAX, VMAX
79*     ..
80*     .. Local Arrays ..
81      INTEGER            LMAX( 4 )
82      DOUBLE PRECISION   A( LDA, LDA ), AF( LDA, LDA ), B( LDB, LDB ),
83     $                   BF( LDB, LDB ), E( LDE, LDE ), F( LDF, LDF ),
84     $                   LSCALE( LDA ), RSCALE( LDA ), VL( LDVL, LDVL ),
85     $                   VLF( LDVL, LDVL ), VR( LDVR, LDVR ),
86     $                   VRF( LDVR, LDVR ), WORK( LDWORK, LDWORK )
87*     ..
88*     .. External Functions ..
89      DOUBLE PRECISION   DLAMCH, DLANGE
90      EXTERNAL           DLAMCH, DLANGE
91*     ..
92*     .. External Subroutines ..
93      EXTERNAL           DGEMM, DGGBAK, DGGBAL, DLACPY
94*     ..
95*     .. Intrinsic Functions ..
96      INTRINSIC          ABS, MAX
97*     ..
98*     .. Executable Statements ..
99*
100*     Initialization
101*
102      LMAX( 1 ) = 0
103      LMAX( 2 ) = 0
104      LMAX( 3 ) = 0
105      LMAX( 4 ) = 0
106      NINFO = 0
107      KNT = 0
108      RMAX = ZERO
109*
110      EPS = DLAMCH( 'Precision' )
111*
112   10 CONTINUE
113      READ( NIN, FMT = * )N, M
114      IF( N.EQ.0 )
115     $   GO TO 100
116*
117      DO 20 I = 1, N
118         READ( NIN, FMT = * )( A( I, J ), J = 1, N )
119   20 CONTINUE
120*
121      DO 30 I = 1, N
122         READ( NIN, FMT = * )( B( I, J ), J = 1, N )
123   30 CONTINUE
124*
125      DO 40 I = 1, N
126         READ( NIN, FMT = * )( VL( I, J ), J = 1, M )
127   40 CONTINUE
128*
129      DO 50 I = 1, N
130         READ( NIN, FMT = * )( VR( I, J ), J = 1, M )
131   50 CONTINUE
132*
133      KNT = KNT + 1
134*
135      ANORM = DLANGE( 'M', N, N, A, LDA, WORK )
136      BNORM = DLANGE( 'M', N, N, B, LDB, WORK )
137*
138      CALL DLACPY( 'FULL', N, N, A, LDA, AF, LDA )
139      CALL DLACPY( 'FULL', N, N, B, LDB, BF, LDB )
140*
141      CALL DGGBAL( 'B', N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
142     $             WORK, INFO )
143      IF( INFO.NE.0 ) THEN
144         NINFO = NINFO + 1
145         LMAX( 1 ) = KNT
146      END IF
147*
148      CALL DLACPY( 'FULL', N, M, VL, LDVL, VLF, LDVL )
149      CALL DLACPY( 'FULL', N, M, VR, LDVR, VRF, LDVR )
150*
151      CALL DGGBAK( 'B', 'L', N, ILO, IHI, LSCALE, RSCALE, M, VL, LDVL,
152     $             INFO )
153      IF( INFO.NE.0 ) THEN
154         NINFO = NINFO + 1
155         LMAX( 2 ) = KNT
156      END IF
157*
158      CALL DGGBAK( 'B', 'R', N, ILO, IHI, LSCALE, RSCALE, M, VR, LDVR,
159     $             INFO )
160      IF( INFO.NE.0 ) THEN
161         NINFO = NINFO + 1
162         LMAX( 3 ) = KNT
163      END IF
164*
165*     Test of DGGBAK
166*
167*     Check tilde(VL)'*A*tilde(VR) - VL'*tilde(A)*VR
168*     where tilde(A) denotes the transformed matrix.
169*
170      CALL DGEMM( 'N', 'N', N, M, N, ONE, AF, LDA, VR, LDVR, ZERO, WORK,
171     $            LDWORK )
172      CALL DGEMM( 'T', 'N', M, M, N, ONE, VL, LDVL, WORK, LDWORK, ZERO,
173     $            E, LDE )
174*
175      CALL DGEMM( 'N', 'N', N, M, N, ONE, A, LDA, VRF, LDVR, ZERO, WORK,
176     $            LDWORK )
177      CALL DGEMM( 'T', 'N', M, M, N, ONE, VLF, LDVL, WORK, LDWORK, ZERO,
178     $            F, LDF )
179*
180      VMAX = ZERO
181      DO 70 J = 1, M
182         DO 60 I = 1, M
183            VMAX = MAX( VMAX, ABS( E( I, J )-F( I, J ) ) )
184   60    CONTINUE
185   70 CONTINUE
186      VMAX = VMAX / ( EPS*MAX( ANORM, BNORM ) )
187      IF( VMAX.GT.RMAX ) THEN
188         LMAX( 4 ) = KNT
189         RMAX = VMAX
190      END IF
191*
192*     Check tilde(VL)'*B*tilde(VR) - VL'*tilde(B)*VR
193*
194      CALL DGEMM( 'N', 'N', N, M, N, ONE, BF, LDB, VR, LDVR, ZERO, WORK,
195     $            LDWORK )
196      CALL DGEMM( 'T', 'N', M, M, N, ONE, VL, LDVL, WORK, LDWORK, ZERO,
197     $            E, LDE )
198*
199      CALL DGEMM( 'N', 'N', N, M, N, ONE, B, LDB, VRF, LDVR, ZERO, WORK,
200     $            LDWORK )
201      CALL DGEMM( 'T', 'N', M, M, N, ONE, VLF, LDVL, WORK, LDWORK, ZERO,
202     $            F, LDF )
203*
204      VMAX = ZERO
205      DO 90 J = 1, M
206         DO 80 I = 1, M
207            VMAX = MAX( VMAX, ABS( E( I, J )-F( I, J ) ) )
208   80    CONTINUE
209   90 CONTINUE
210      VMAX = VMAX / ( EPS*MAX( ANORM, BNORM ) )
211      IF( VMAX.GT.RMAX ) THEN
212         LMAX( 4 ) = KNT
213         RMAX = VMAX
214      END IF
215*
216      GO TO 10
217*
218  100 CONTINUE
219*
220      WRITE( NOUT, FMT = 9999 )
221 9999 FORMAT( 1X, '.. test output of DGGBAK .. ' )
222*
223      WRITE( NOUT, FMT = 9998 )RMAX
224 9998 FORMAT( ' value of largest test error                  =', D12.3 )
225      WRITE( NOUT, FMT = 9997 )LMAX( 1 )
226 9997 FORMAT( ' example number where DGGBAL info is not 0    =', I4 )
227      WRITE( NOUT, FMT = 9996 )LMAX( 2 )
228 9996 FORMAT( ' example number where DGGBAK(L) info is not 0 =', I4 )
229      WRITE( NOUT, FMT = 9995 )LMAX( 3 )
230 9995 FORMAT( ' example number where DGGBAK(R) info is not 0 =', I4 )
231      WRITE( NOUT, FMT = 9994 )LMAX( 4 )
232 9994 FORMAT( ' example number having largest error          =', I4 )
233      WRITE( NOUT, FMT = 9993 )NINFO
234 9993 FORMAT( ' number of examples where info is not 0       =', I4 )
235      WRITE( NOUT, FMT = 9992 )KNT
236 9992 FORMAT( ' total number of examples tested              =', I4 )
237*
238      RETURN
239*
240*     End of DCHKGK
241*
242      END
243