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
138
139
140 vm1( 1 ) = sqrt( smlnum )
141 vm1( 2 ) = one
142 vm1( 3 ) = large
143 vm2( 1 ) = one
144 vm2( 2 ) = one + two*eps
145 vm2( 3 ) = two
146
147 knt = 0
148 ninfo = 0
149 lmax = 0
150 rmax = zero
151
152
153
154 10 CONTINUE
155 READ( nin, fmt = * )m, n
156 IF( n.EQ.0 )
157 $ RETURN
158 DO 20 i = 1, m
159 READ( nin, fmt = * )( atmp( i, j ), j = 1, m )
160 20 CONTINUE
161 DO 30 i = 1, n
162 READ( nin, fmt = * )( btmp( i, j ), j = 1, n )
163 30 CONTINUE
164 DO 40 i = 1, m
165 READ( nin, fmt = * )( ctmp( i, j ), j = 1, n )
166 40 CONTINUE
167 DO 170 imla = 1, 3
168 DO 160 imlad = 1, 3
169 DO 150 imlb = 1, 3
170 DO 140 imlc = 1, 3
171 DO 130 itrana = 1, 2
172 DO 120 itranb = 1, 2
173 DO 110 isgn = -1, 1, 2
174 IF( itrana.EQ.1 )
175 $ trana = 'N'
176 IF( itrana.EQ.2 )
177 $ trana = 'C'
178 IF( itranb.EQ.1 )
179 $ tranb = 'N'
180 IF( itranb.EQ.2 )
181 $ tranb = 'C'
182 tnrm = zero
183 DO 60 i = 1, m
184 DO 50 j = 1, m
185 a( i, j ) = atmp( i, j )*vm1( imla )
186 tnrm = max( tnrm, abs( a( i, j ) ) )
187 50 CONTINUE
188 a( i, i ) = a( i, i )*vm2( imlad )
189 tnrm = max( tnrm, abs( a( i, i ) ) )
190 60 CONTINUE
191 DO 80 i = 1, n
192 DO 70 j = 1, n
193 b( i, j ) = btmp( i, j )*vm1( imlb )
194 tnrm = max( tnrm, abs( b( i, j ) ) )
195 70 CONTINUE
196 80 CONTINUE
197 IF( tnrm.EQ.zero )
198 $ tnrm = one
199 DO 100 i = 1, m
200 DO 90 j = 1, n
201 c( i, j ) = ctmp( i, j )*vm1( imlc )
202 csav( i, j ) = c( i, j )
203 90 CONTINUE
204 100 CONTINUE
205 knt = knt + 1
206 CALL ztrsyl( trana, tranb, isgn, m, n, a,
207 $ ldt, b, ldt, c, ldt, scale,
208 $ info )
209 IF( info.NE.0 )
210 $ ninfo = ninfo + 1
211 xnrm =
zlange(
'M', m, n, c, ldt, dum )
212 rmul = cone
213 IF( xnrm.GT.one .AND. tnrm.GT.one ) THEN
214 IF( xnrm.GT.bignum / tnrm ) THEN
215 rmul = max( xnrm, tnrm )
216 rmul = cone / rmul
217 END IF
218 END IF
219 CALL zgemm( trana,
'N', m, n, m, rmul, a,
220 $ ldt, c, ldt, -scale*rmul, csav,
221 $ ldt )
222 CALL zgemm(
'N', tranb, m, n, n,
223 $ dble( isgn )*rmul, c, ldt, b,
224 $ ldt, cone, csav, ldt )
225 res1 =
zlange(
'M', m, n, csav, ldt, dum )
226 res = res1 / max( smlnum, smlnum*xnrm,
227 $ ( ( abs( rmul )*tnrm )*eps )*xnrm )
228 IF( res.GT.rmax ) THEN
229 lmax = knt
230 rmax = res
231 END IF
232 110 CONTINUE
233 120 CONTINUE
234 130 CONTINUE
235 140 CONTINUE
236 150 CONTINUE
237 160 CONTINUE
238 170 CONTINUE
239 GO TO 10
240
241
242
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
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 ...
subroutine ztrsyl(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, info)
ZTRSYL