84
85
86
87
88
89
90 INTEGER KNT, LMAX, NIN, NINFO
91 DOUBLE PRECISION RMAX
92
93
94
95
96
97 INTEGER LDT
98 parameter( ldt = 10 )
99 DOUBLE PRECISION ZERO, ONE, TWO
100 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
101 DOUBLE PRECISION LARGE
102 parameter( large = 1.0d6 )
103 COMPLEX*16 CONE
104 parameter( cone = 1.0d0 )
105
106
107 CHARACTER TRANA, TRANB
108 INTEGER I, IMLA, IMLAD, IMLB, IMLC, INFO, ISGN, ITRANA,
109 $ ITRANB, J, M, N
110 DOUBLE PRECISION BIGNUM, EPS, RES, RES1, SCALE, SMLNUM, TNRM,
111 $ XNRM
112 COMPLEX*16 RMUL
113
114
115 DOUBLE PRECISION DUM( 1 ), VM1( 3 ), VM2( 3 )
116 COMPLEX*16 A( LDT, LDT ), ATMP( LDT, LDT ), B( LDT, LDT ),
117 $ BTMP( LDT, LDT ), C( LDT, LDT ),
118 $ CSAV( LDT, LDT ), CTMP( LDT, LDT )
119
120
121 DOUBLE PRECISION DLAMCH, ZLANGE
123
124
126
127
128 INTRINSIC abs, dble, max, sqrt
129
130
131
132
133
135 smlnum =
dlamch(
'S' ) / eps
136 bignum = one / smlnum
137 CALL dlabad( smlnum, bignum )
138
139
140
141 vm1( 1 ) = sqrt( smlnum )
142 vm1( 2 ) = one
143 vm1( 3 ) = large
144 vm2( 1 ) = one
145 vm2( 2 ) = one + two*eps
146 vm2( 3 ) = two
147
148 knt = 0
149 ninfo = 0
150 lmax = 0
151 rmax = zero
152
153
154
155 10 CONTINUE
156 READ( nin, fmt = * )m, n
157 IF( n.EQ.0 )
158 $ RETURN
159 DO 20 i = 1, m
160 READ( nin, fmt = * )( atmp( i, j ), j = 1, m )
161 20 CONTINUE
162 DO 30 i = 1, n
163 READ( nin, fmt = * )( btmp( i, j ), j = 1, n )
164 30 CONTINUE
165 DO 40 i = 1, m
166 READ( nin, fmt = * )( ctmp( i, j ), j = 1, n )
167 40 CONTINUE
168 DO 170 imla = 1, 3
169 DO 160 imlad = 1, 3
170 DO 150 imlb = 1, 3
171 DO 140 imlc = 1, 3
172 DO 130 itrana = 1, 2
173 DO 120 itranb = 1, 2
174 DO 110 isgn = -1, 1, 2
175 IF( itrana.EQ.1 )
176 $ trana = 'N'
177 IF( itrana.EQ.2 )
178 $ trana = 'C'
179 IF( itranb.EQ.1 )
180 $ tranb = 'N'
181 IF( itranb.EQ.2 )
182 $ tranb = 'C'
183 tnrm = zero
184 DO 60 i = 1, m
185 DO 50 j = 1, m
186 a( i, j ) = atmp( i, j )*vm1( imla )
187 tnrm = max( tnrm, abs( a( i, j ) ) )
188 50 CONTINUE
189 a( i, i ) = a( i, i )*vm2( imlad )
190 tnrm = max( tnrm, abs( a( i, i ) ) )
191 60 CONTINUE
192 DO 80 i = 1, n
193 DO 70 j = 1, n
194 b( i, j ) = btmp( i, j )*vm1( imlb )
195 tnrm = max( tnrm, abs( b( i, j ) ) )
196 70 CONTINUE
197 80 CONTINUE
198 IF( tnrm.EQ.zero )
199 $ tnrm = one
200 DO 100 i = 1, m
201 DO 90 j = 1, n
202 c( i, j ) = ctmp( i, j )*vm1( imlc )
203 csav( i, j ) = c( i, j )
204 90 CONTINUE
205 100 CONTINUE
206 knt = knt + 1
207 CALL ztrsyl( trana, tranb, isgn, m, n, a,
208 $ ldt, b, ldt, c, ldt, scale,
209 $ info )
210 IF( info.NE.0 )
211 $ ninfo = ninfo + 1
212 xnrm =
zlange(
'M', m, n, c, ldt, dum )
213 rmul = cone
214 IF( xnrm.GT.one .AND. tnrm.GT.one ) THEN
215 IF( xnrm.GT.bignum / tnrm ) THEN
216 rmul = max( xnrm, tnrm )
217 rmul = cone / rmul
218 END IF
219 END IF
220 CALL zgemm( trana,
'N', m, n, m, rmul, a,
221 $ ldt, c, ldt, -scale*rmul, csav,
222 $ ldt )
223 CALL zgemm(
'N', tranb, m, n, n,
224 $ dble( isgn )*rmul, c, ldt, b,
225 $ ldt, cone, csav, ldt )
226 res1 =
zlange(
'M', m, n, csav, ldt, dum )
227 res = res1 / max( smlnum, smlnum*xnrm,
228 $ ( ( abs( rmul )*tnrm )*eps )*xnrm )
229 IF( res.GT.rmax ) THEN
230 lmax = knt
231 rmax = res
232 END IF
233 110 CONTINUE
234 120 CONTINUE
235 130 CONTINUE
236 140 CONTINUE
237 150 CONTINUE
238 160 CONTINUE
239 170 CONTINUE
240 GO TO 10
241
242
243
double precision function dlamch(CMACH)
DLAMCH
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
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 ...
subroutine ztrsyl(TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO)
ZTRSYL