93 INTEGER knt, lmax, nin, ninfo
101 parameter ( ldt = 10 )
103 parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
105 parameter ( large = 1.0e6 )
107 parameter ( cone = 1.0e0 )
110 CHARACTER trana, tranb
111 INTEGER i, imla, imlad, imlb, imlc, info, isgn, itrana,
113 REAL bignum, eps, res, res1, scale, smlnum, tnrm,
118 REAL dum( 1 ), vm1( 3 ), vm2( 3 )
119 COMPLEX a( ldt, ldt ), atmp( ldt, ldt ), b( ldt, ldt ),
120 $ btmp( ldt, ldt ), c( ldt, ldt ),
121 $ csav( ldt, ldt ), ctmp( ldt, ldt )
131 INTRINSIC abs, max,
REAL, sqrt
138 smlnum =
slamch(
'S' ) / eps
139 bignum = one / smlnum
140 CALL slabad( smlnum, bignum )
144 vm1( 1 ) = sqrt( smlnum )
148 vm2( 2 ) = one + two*eps
159 READ( nin, fmt = * )m, n
163 READ( nin, fmt = * )( atmp( i, j ), j = 1, m )
166 READ( nin, fmt = * )( btmp( i, j ), j = 1, n )
169 READ( nin, fmt = * )( ctmp( i, j ), j = 1, n )
177 DO 110 isgn = -1, 1, 2
189 a( i, j ) = atmp( i, j )*vm1( imla )
190 tnrm = max( tnrm, abs( a( i, j ) ) )
192 a( i, i ) = a( i, i )*vm2( imlad )
193 tnrm = max( tnrm, abs( a( i, i ) ) )
197 b( i, j ) = btmp( i, j )*vm1( imlb )
198 tnrm = max( tnrm, abs( b( i, j ) ) )
205 c( i, j ) = ctmp( i, j )*vm1( imlc )
206 csav( i, j ) = c( i, j )
210 CALL ctrsyl( trana, tranb, isgn, m, n, a,
211 $ ldt, b, ldt, c, ldt, scale,
215 xnrm =
clange(
'M', m, n, c, ldt, dum )
217 IF( xnrm.GT.one .AND. tnrm.GT.one )
THEN
218 IF( xnrm.GT.bignum / tnrm )
THEN
219 rmul = max( xnrm, tnrm )
223 CALL cgemm( trana,
'N', m, n, m, rmul, a,
224 $ ldt, c, ldt, -scale*rmul, csav,
226 CALL cgemm(
'N', tranb, m, n, n,
227 $
REAL( isgn )*rmul, c, ldt, b,
228 $ ldt, cone, csav, ldt )
229 res1 =
clange(
'M', m, n, csav, ldt, dum )
230 res = res1 / max( smlnum, smlnum*xnrm,
231 $ ( ( abs( rmul )*tnrm )*eps )*xnrm )
232 IF( res.GT.rmax )
THEN
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine ctrsyl(TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO)
CTRSYL
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 ...
real function slamch(CMACH)
SLAMCH
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM