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 REAL ZERO
72 parameter( zero = 0.0e+0 )
73 COMPLEX CZERO, CONE
74 parameter( czero = ( 0.0e+0, 0.0e+0 ),
75 $ cone = ( 1.0e+0, 0.0e+0 ) )
76
77
78 INTEGER I, IHI, ILO, INFO, J, KNT, M, N, NINFO
79 REAL ANORM, BNORM, EPS, RMAX, VMAX
80 COMPLEX CDUM
81
82
83 INTEGER LMAX( 4 )
84 REAL LSCALE( LDA ), RSCALE( LDA ), RWORK( LRWORK )
85 COMPLEX 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 REAL CLANGE, SLAMCH
94
95
97
98
99 INTRINSIC abs, aimag, max, real
100
101
102 REAL CABS1
103
104
105 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( 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 =
slamch(
'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 =
clange(
'M', n, n, a, lda, rwork )
143 bnorm =
clange(
'M', n, n, b, ldb, rwork )
144
145 CALL clacpy(
'FULL', n, n, a, lda, af, lda )
146 CALL clacpy(
'FULL', n, n, b, ldb, bf, ldb )
147
148 CALL cggbal(
'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 clacpy(
'FULL', n, m, vl, ldvl, vlf, ldvl )
156 CALL clacpy(
'FULL', n, m, vr, ldvr, vrf, ldvr )
157
158 CALL cggbak(
'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 cggbak(
'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 cgemm(
'N',
'N', n, m, n, cone, af, lda, vr, ldvr, czero,
178 $ work, ldwork )
179 CALL cgemm(
'C',
'N', m, m, n, cone, vl, ldvl, work, ldwork,
181
182 CALL cgemm(
'N',
'N', n, m, n, cone, a, lda, vrf, ldvr, czero,
183 $ work, ldwork )
184 CALL cgemm(
'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 cgemm(
'N',
'N', n, m, n, cone, bf, ldb, vr, ldvr, czero,
202 $ work, ldwork )
203 CALL cgemm(
'C',
'N', m, m, n, cone, vl, ldvl, work, ldwork,
205
206 CALL cgemm(
'n',
'n', n, m, n, cone, b, ldb, vrf, ldvr, czero,
207 $ work, ldwork )
208 CALL cgemm(
'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 CGGBAK .. ' )
229
230 WRITE( nout, fmt = 9998 )rmax
231 9998 FORMAT( ' value of largest test error =', e12.3 )
232 WRITE( nout, fmt = 9997 )lmax( 1 )
233 9997 FORMAT( ' example number where CGGBAL info is not 0 =', i4 )
234 WRITE( nout, fmt = 9996 )lmax( 2 )
235 9996 FORMAT( ' example number where CGGBAK(L) info is not 0 =', i4 )
236 WRITE( nout, fmt = 9995 )lmax( 3 )
237 9995 FORMAT( ' example number where CGGBAK(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 cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
subroutine cggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
CGGBAK
subroutine cggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
CGGBAL
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...