54
55
56
57
58
59
60 INTEGER NIN, NOUT
61
62
63
64
65
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
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
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
92 DOUBLE PRECISION DLAMCH, ZLANGE
94
95
97
98
99 INTRINSIC abs, dble, dimag, max
100
101
102 DOUBLE PRECISION CABS1
103
104
105 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
106
107
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
173
174
175
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,
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
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,
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
248
logical function lde(ri, rj, lr)
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
subroutine zggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
ZGGBAK
subroutine zggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
ZGGBAL
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...