LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zchkgk.f
Go to the documentation of this file.
1*> \brief \b ZCHKGK
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 ZCHKGK( NIN, NOUT )
12*
13* .. Scalar Arguments ..
14* INTEGER NIN, NOUT
15* ..
16*
17*
18*> \par Purpose:
19* =============
20*>
21*> \verbatim
22*>
23*> ZCHKGK tests ZGGBAK, 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 complex16_eig
51*
52* =====================================================================
53 SUBROUTINE zchkgk( 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, LRWORK
69 parameter( lde = 50, ldf = 50, ldwork = 50,
70 $ lrwork = 6*50 )
71 DOUBLE PRECISION ZERO
72 parameter( zero = 0.0d+0 )
73 COMPLEX*16 CZERO, CONE
74 parameter( czero = ( 0.0d+0, 0.0d+0 ),
75 $ cone = ( 1.0d+0, 0.0d+0 ) )
76* ..
77* .. Local Scalars ..
78 INTEGER I, IHI, ILO, INFO, J, KNT, M, N, NINFO
79 DOUBLE PRECISION ANORM, BNORM, EPS, RMAX, VMAX
80 COMPLEX*16 CDUM
81* ..
82* .. Local Arrays ..
83 INTEGER LMAX( 4 )
84 DOUBLE PRECISION LSCALE( LDA ), RSCALE( LDA ), RWORK( LRWORK )
85 COMPLEX*16 A( LDA, LDA ), AF( LDA, LDA ), B( LDB, LDB ),
86 $ BF( LDB, LDB ), E( LDE, LDE ), F( LDF, LDF ),
87 $ VL( LDVL, LDVL ), VLF( LDVL, LDVL ),
88 $ VR( LDVR, LDVR ), VRF( LDVR, LDVR ),
89 $ WORK( LDWORK, LDWORK )
90* ..
91* .. External Functions ..
92 DOUBLE PRECISION DLAMCH, ZLANGE
93 EXTERNAL dlamch, zlange
94* ..
95* .. External Subroutines ..
96 EXTERNAL zgemm, zggbak, zggbal, zlacpy
97* ..
98* .. Intrinsic Functions ..
99 INTRINSIC abs, dble, dimag, max
100* ..
101* .. Statement Functions ..
102 DOUBLE PRECISION CABS1
103* ..
104* .. Statement Function definitions ..
105 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
106* ..
107* .. Executable Statements ..
108*
109 lmax( 1 ) = 0
110 lmax( 2 ) = 0
111 lmax( 3 ) = 0
112 lmax( 4 ) = 0
113 ninfo = 0
114 knt = 0
115 rmax = zero
116*
117 eps = dlamch( 'Precision' )
118*
119 10 CONTINUE
120 READ( nin, fmt = * )n, m
121 IF( n.EQ.0 )
122 $ GO TO 100
123*
124 DO 20 i = 1, n
125 READ( nin, fmt = * )( a( i, j ), j = 1, n )
126 20 CONTINUE
127*
128 DO 30 i = 1, n
129 READ( nin, fmt = * )( b( i, j ), j = 1, n )
130 30 CONTINUE
131*
132 DO 40 i = 1, n
133 READ( nin, fmt = * )( vl( i, j ), j = 1, m )
134 40 CONTINUE
135*
136 DO 50 i = 1, n
137 READ( nin, fmt = * )( vr( i, j ), j = 1, m )
138 50 CONTINUE
139*
140 knt = knt + 1
141*
142 anorm = zlange( 'M', n, n, a, lda, rwork )
143 bnorm = zlange( 'M', n, n, b, ldb, rwork )
144*
145 CALL zlacpy( 'FULL', n, n, a, lda, af, lda )
146 CALL zlacpy( 'FULL', n, n, b, ldb, bf, ldb )
147*
148 CALL zggbal( 'B', n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
149 $ rwork, info )
150 IF( info.NE.0 ) THEN
151 ninfo = ninfo + 1
152 lmax( 1 ) = knt
153 END IF
154*
155 CALL zlacpy( 'FULL', n, m, vl, ldvl, vlf, ldvl )
156 CALL zlacpy( 'FULL', n, m, vr, ldvr, vrf, ldvr )
157*
158 CALL zggbak( 'B', 'L', n, ilo, ihi, lscale, rscale, m, vl, ldvl,
159 $ info )
160 IF( info.NE.0 ) THEN
161 ninfo = ninfo + 1
162 lmax( 2 ) = knt
163 END IF
164*
165 CALL zggbak( 'B', 'R', n, ilo, ihi, lscale, rscale, m, vr, ldvr,
166 $ info )
167 IF( info.NE.0 ) THEN
168 ninfo = ninfo + 1
169 lmax( 3 ) = knt
170 END IF
171*
172* Test of ZGGBAK
173*
174* Check tilde(VL)'*A*tilde(VR) - VL'*tilde(A)*VR
175* where tilde(A) denotes the transformed matrix.
176*
177 CALL zgemm( 'N', 'N', n, m, n, cone, af, lda, vr, ldvr, czero,
178 $ work, ldwork )
179 CALL zgemm( 'C', 'N', m, m, n, cone, vl, ldvl, work, ldwork,
180 $ czero, e, lde )
181*
182 CALL zgemm( 'N', 'N', n, m, n, cone, a, lda, vrf, ldvr, czero,
183 $ work, ldwork )
184 CALL zgemm( 'C', 'N', m, m, n, cone, vlf, ldvl, work, ldwork,
185 $ czero, f, ldf )
186*
187 vmax = zero
188 DO 70 j = 1, m
189 DO 60 i = 1, m
190 vmax = max( vmax, cabs1( e( i, j )-f( i, j ) ) )
191 60 CONTINUE
192 70 CONTINUE
193 vmax = vmax / ( eps*max( anorm, bnorm ) )
194 IF( vmax.GT.rmax ) THEN
195 lmax( 4 ) = knt
196 rmax = vmax
197 END IF
198*
199* Check tilde(VL)'*B*tilde(VR) - VL'*tilde(B)*VR
200*
201 CALL zgemm( 'N', 'N', n, m, n, cone, bf, ldb, vr, ldvr, czero,
202 $ work, ldwork )
203 CALL zgemm( 'C', 'N', m, m, n, cone, vl, ldvl, work, ldwork,
204 $ czero, e, lde )
205*
206 CALL zgemm( 'n', 'n', n, m, n, cone, b, ldb, vrf, ldvr, czero,
207 $ work, ldwork )
208 CALL zgemm( 'C', 'N', m, m, n, cone, vlf, ldvl, work, ldwork,
209 $ czero, f, ldf )
210*
211 vmax = zero
212 DO 90 j = 1, m
213 DO 80 i = 1, m
214 vmax = max( vmax, cabs1( e( i, j )-f( i, j ) ) )
215 80 CONTINUE
216 90 CONTINUE
217 vmax = vmax / ( eps*max( anorm, bnorm ) )
218 IF( vmax.GT.rmax ) THEN
219 lmax( 4 ) = knt
220 rmax = vmax
221 END IF
222*
223 GO TO 10
224*
225 100 CONTINUE
226*
227 WRITE( nout, fmt = 9999 )
228 9999 FORMAT( 1x, '.. test output of ZGGBAK .. ' )
229*
230 WRITE( nout, fmt = 9998 )rmax
231 9998 FORMAT( ' value of largest test error =', d12.3 )
232 WRITE( nout, fmt = 9997 )lmax( 1 )
233 9997 FORMAT( ' example number where ZGGBAL info is not 0 =', i4 )
234 WRITE( nout, fmt = 9996 )lmax( 2 )
235 9996 FORMAT( ' example number where ZGGBAK(L) info is not 0 =', i4 )
236 WRITE( nout, fmt = 9995 )lmax( 3 )
237 9995 FORMAT( ' example number where ZGGBAK(R) info is not 0 =', i4 )
238 WRITE( nout, fmt = 9994 )lmax( 4 )
239 9994 FORMAT( ' example number having largest error =', i4 )
240 WRITE( nout, fmt = 9992 )ninfo
241 9992 FORMAT( ' number of examples where info is not 0 =', i4 )
242 WRITE( nout, fmt = 9991 )knt
243 9991 FORMAT( ' total number of examples tested =', i4 )
244*
245 RETURN
246*
247* End of ZCHKGK
248*
249 END
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
subroutine zggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
ZGGBAK
Definition zggbak.f:148
subroutine zggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
ZGGBAL
Definition zggbal.f:177
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
subroutine zchkgk(nin, nout)
ZCHKGK
Definition zchkgk.f:54