LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dchkgk.f
Go to the documentation of this file.
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*> \ingroup double_eig
51*
52* =====================================================================
53 SUBROUTINE dchkgk( NIN, NOUT )
54*
55* -- LAPACK test routine --
56* -- LAPACK is a software package provided by Univ. of Tennessee, --
57* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
58*
59* .. Scalar Arguments ..
60 INTEGER NIN, NOUT
61* ..
62*
63* =====================================================================
64*
65* .. Parameters ..
66 INTEGER LDA, LDB, LDVL, LDVR
67 parameter( lda = 50, ldb = 50, ldvl = 50, ldvr = 50 )
68 INTEGER LDE, LDF, LDWORK
69 parameter( lde = 50, ldf = 50, ldwork = 50 )
70 DOUBLE PRECISION ZERO, ONE
71 parameter( zero = 0.0d+0, one = 1.0d+0 )
72* ..
73* .. Local Scalars ..
74 INTEGER I, IHI, ILO, INFO, J, KNT, M, N, NINFO
75 DOUBLE PRECISION ANORM, BNORM, EPS, RMAX, VMAX
76* ..
77* .. Local Arrays ..
78 INTEGER LMAX( 4 )
79 DOUBLE PRECISION A( LDA, LDA ), AF( LDA, LDA ), B( LDB, LDB ),
80 $ BF( LDB, LDB ), E( LDE, LDE ), F( LDF, LDF ),
81 $ LSCALE( LDA ), RSCALE( LDA ), VL( LDVL, LDVL ),
82 $ VLF( LDVL, LDVL ), VR( LDVR, LDVR ),
83 $ VRF( LDVR, LDVR ), WORK( LDWORK, LDWORK )
84* ..
85* .. External Functions ..
86 DOUBLE PRECISION DLAMCH, DLANGE
87 EXTERNAL dlamch, dlange
88* ..
89* .. External Subroutines ..
90 EXTERNAL dgemm, dggbak, dggbal, dlacpy
91* ..
92* .. Intrinsic Functions ..
93 INTRINSIC abs, max
94* ..
95* .. Executable Statements ..
96*
97* Initialization
98*
99 lmax( 1 ) = 0
100 lmax( 2 ) = 0
101 lmax( 3 ) = 0
102 lmax( 4 ) = 0
103 ninfo = 0
104 knt = 0
105 rmax = zero
106*
107 eps = dlamch( 'Precision' )
108*
109 10 CONTINUE
110 READ( nin, fmt = * )n, m
111 IF( n.EQ.0 )
112 $ GO TO 100
113*
114 DO 20 i = 1, n
115 READ( nin, fmt = * )( a( i, j ), j = 1, n )
116 20 CONTINUE
117*
118 DO 30 i = 1, n
119 READ( nin, fmt = * )( b( i, j ), j = 1, n )
120 30 CONTINUE
121*
122 DO 40 i = 1, n
123 READ( nin, fmt = * )( vl( i, j ), j = 1, m )
124 40 CONTINUE
125*
126 DO 50 i = 1, n
127 READ( nin, fmt = * )( vr( i, j ), j = 1, m )
128 50 CONTINUE
129*
130 knt = knt + 1
131*
132 anorm = dlange( 'M', n, n, a, lda, work )
133 bnorm = dlange( 'M', n, n, b, ldb, work )
134*
135 CALL dlacpy( 'FULL', n, n, a, lda, af, lda )
136 CALL dlacpy( 'FULL', n, n, b, ldb, bf, ldb )
137*
138 CALL dggbal( 'B', n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
139 $ work, info )
140 IF( info.NE.0 ) THEN
141 ninfo = ninfo + 1
142 lmax( 1 ) = knt
143 END IF
144*
145 CALL dlacpy( 'FULL', n, m, vl, ldvl, vlf, ldvl )
146 CALL dlacpy( 'FULL', n, m, vr, ldvr, vrf, ldvr )
147*
148 CALL dggbak( 'B', 'L', n, ilo, ihi, lscale, rscale, m, vl, ldvl,
149 $ info )
150 IF( info.NE.0 ) THEN
151 ninfo = ninfo + 1
152 lmax( 2 ) = knt
153 END IF
154*
155 CALL dggbak( 'B', 'R', n, ilo, ihi, lscale, rscale, m, vr, ldvr,
156 $ info )
157 IF( info.NE.0 ) THEN
158 ninfo = ninfo + 1
159 lmax( 3 ) = knt
160 END IF
161*
162* Test of DGGBAK
163*
164* Check tilde(VL)'*A*tilde(VR) - VL'*tilde(A)*VR
165* where tilde(A) denotes the transformed matrix.
166*
167 CALL dgemm( 'N', 'N', n, m, n, one, af, lda, vr, ldvr, zero, work,
168 $ ldwork )
169 CALL dgemm( 'T', 'N', m, m, n, one, vl, ldvl, work, ldwork, zero,
170 $ e, lde )
171*
172 CALL dgemm( 'N', 'N', n, m, n, one, a, lda, vrf, ldvr, zero, work,
173 $ ldwork )
174 CALL dgemm( 'T', 'N', m, m, n, one, vlf, ldvl, work, ldwork, zero,
175 $ f, ldf )
176*
177 vmax = zero
178 DO 70 j = 1, m
179 DO 60 i = 1, m
180 vmax = max( vmax, abs( e( i, j )-f( i, j ) ) )
181 60 CONTINUE
182 70 CONTINUE
183 vmax = vmax / ( eps*max( anorm, bnorm ) )
184 IF( vmax.GT.rmax ) THEN
185 lmax( 4 ) = knt
186 rmax = vmax
187 END IF
188*
189* Check tilde(VL)'*B*tilde(VR) - VL'*tilde(B)*VR
190*
191 CALL dgemm( 'N', 'N', n, m, n, one, bf, ldb, vr, ldvr, zero, work,
192 $ ldwork )
193 CALL dgemm( 'T', 'N', m, m, n, one, vl, ldvl, work, ldwork, zero,
194 $ e, lde )
195*
196 CALL dgemm( 'N', 'N', n, m, n, one, b, ldb, vrf, ldvr, zero, work,
197 $ ldwork )
198 CALL dgemm( 'T', 'N', m, m, n, one, vlf, ldvl, work, ldwork, zero,
199 $ f, ldf )
200*
201 vmax = zero
202 DO 90 j = 1, m
203 DO 80 i = 1, m
204 vmax = max( vmax, abs( e( i, j )-f( i, j ) ) )
205 80 CONTINUE
206 90 CONTINUE
207 vmax = vmax / ( eps*max( anorm, bnorm ) )
208 IF( vmax.GT.rmax ) THEN
209 lmax( 4 ) = knt
210 rmax = vmax
211 END IF
212*
213 GO TO 10
214*
215 100 CONTINUE
216*
217 WRITE( nout, fmt = 9999 )
218 9999 FORMAT( 1x, '.. test output of DGGBAK .. ' )
219*
220 WRITE( nout, fmt = 9998 )rmax
221 9998 FORMAT( ' value of largest test error =', d12.3 )
222 WRITE( nout, fmt = 9997 )lmax( 1 )
223 9997 FORMAT( ' example number where DGGBAL info is not 0 =', i4 )
224 WRITE( nout, fmt = 9996 )lmax( 2 )
225 9996 FORMAT( ' example number where DGGBAK(L) info is not 0 =', i4 )
226 WRITE( nout, fmt = 9995 )lmax( 3 )
227 9995 FORMAT( ' example number where DGGBAK(R) info is not 0 =', i4 )
228 WRITE( nout, fmt = 9994 )lmax( 4 )
229 9994 FORMAT( ' example number having largest error =', i4 )
230 WRITE( nout, fmt = 9993 )ninfo
231 9993 FORMAT( ' number of examples where info is not 0 =', i4 )
232 WRITE( nout, fmt = 9992 )knt
233 9992 FORMAT( ' total number of examples tested =', i4 )
234*
235 RETURN
236*
237* End of DCHKGK
238*
239 END
subroutine dchkgk(nin, nout)
DCHKGK
Definition dchkgk.f:54
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
DGGBAK
Definition dggbak.f:147
subroutine dggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
DGGBAL
Definition dggbal.f:177
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103