261 lines
7.5 KiB
Fortran
261 lines
7.5 KiB
Fortran
c { dg-do compile }
|
|
c
|
|
c This demonstrates a problem with g77 and pic on x86 where
|
|
c egcs 1.0.1 and earlier will generate bogus assembler output.
|
|
c unfortunately, gas accepts the bogus acssembler output and
|
|
c generates code that almost works.
|
|
c
|
|
|
|
|
|
C Date: Wed, 17 Dec 1997 23:20:29 +0000
|
|
C From: Joao Cardoso <jcardoso@inescn.pt>
|
|
C To: egcs-bugs@cygnus.com
|
|
C Subject: egcs-1.0 f77 bug on OSR5
|
|
C When trying to compile the Fortran file that I enclose bellow,
|
|
C I got an assembler error:
|
|
C
|
|
C ./g77 -B./ -fpic -O -c scaleg.f
|
|
C /usr/tmp/cca002D8.s:123:syntax error at (
|
|
C
|
|
C ./g77 -B./ -fpic -O0 -c scaleg.f
|
|
C /usr/tmp/cca002EW.s:246:invalid operand combination: leal
|
|
C
|
|
C Compiling without the -fpic flag runs OK.
|
|
|
|
subroutine scaleg (n,ma,a,mb,b,low,igh,cscale,cperm,wk)
|
|
c
|
|
c *****parameters:
|
|
integer igh,low,ma,mb,n
|
|
double precision a(ma,n),b(mb,n),cperm(n),cscale(n),wk(n,6)
|
|
c
|
|
c *****local variables:
|
|
integer i,ir,it,j,jc,kount,nr,nrp2
|
|
double precision alpha,basl,beta,cmax,coef,coef2,coef5,cor,
|
|
* ew,ewc,fi,fj,gamma,pgamma,sum,t,ta,tb,tc
|
|
c
|
|
c *****fortran functions:
|
|
double precision dabs, dlog10, dsign
|
|
c float
|
|
c
|
|
c *****subroutines called:
|
|
c none
|
|
c
|
|
c ---------------------------------------------------------------
|
|
c
|
|
c *****purpose:
|
|
c scales the matrices a and b in the generalized eigenvalue
|
|
c problem a*x = (lambda)*b*x such that the magnitudes of the
|
|
c elements of the submatrices of a and b (as specified by low
|
|
c and igh) are close to unity in the least squares sense.
|
|
c ref.: ward, r. c., balancing the generalized eigenvalue
|
|
c problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981,
|
|
c 141-152.
|
|
c
|
|
c *****parameter description:
|
|
c
|
|
c on input:
|
|
c
|
|
c ma,mb integer
|
|
c row dimensions of the arrays containing matrices
|
|
c a and b respectively, as declared in the main calling
|
|
c program dimension statement;
|
|
c
|
|
c n integer
|
|
c order of the matrices a and b;
|
|
c
|
|
c a real(ma,n)
|
|
c contains the a matrix of the generalized eigenproblem
|
|
c defined above;
|
|
c
|
|
c b real(mb,n)
|
|
c contains the b matrix of the generalized eigenproblem
|
|
c defined above;
|
|
c
|
|
c low integer
|
|
c specifies the beginning -1 for the rows and
|
|
c columns of a and b to be scaled;
|
|
c
|
|
c igh integer
|
|
c specifies the ending -1 for the rows and columns
|
|
c of a and b to be scaled;
|
|
c
|
|
c cperm real(n)
|
|
c work array. only locations low through igh are
|
|
c referenced and altered by this subroutine;
|
|
c
|
|
c wk real(n,6)
|
|
c work array that must contain at least 6*n locations.
|
|
c only locations low through igh, n+low through n+igh,
|
|
c ..., 5*n+low through 5*n+igh are referenced and
|
|
c altered by this subroutine.
|
|
c
|
|
c on output:
|
|
c
|
|
c a,b contain the scaled a and b matrices;
|
|
c
|
|
c cscale real(n)
|
|
c contains in its low through igh locations the integer
|
|
c exponents of 2 used for the column scaling factors.
|
|
c the other locations are not referenced;
|
|
c
|
|
c wk contains in its low through igh locations the integer
|
|
c exponents of 2 used for the row scaling factors.
|
|
c
|
|
c *****algorithm notes:
|
|
c none.
|
|
c
|
|
c *****history:
|
|
c written by r. c. ward.......
|
|
c modified 8/86 by bobby bodenheimer so that if
|
|
c sum = 0 (corresponding to the case where the matrix
|
|
c doesn't need to be scaled) the routine returns.
|
|
c
|
|
c ---------------------------------------------------------------
|
|
c
|
|
if (low .eq. igh) go to 410
|
|
do 210 i = low,igh
|
|
wk(i,1) = 0.0d0
|
|
wk(i,2) = 0.0d0
|
|
wk(i,3) = 0.0d0
|
|
wk(i,4) = 0.0d0
|
|
wk(i,5) = 0.0d0
|
|
wk(i,6) = 0.0d0
|
|
cscale(i) = 0.0d0
|
|
cperm(i) = 0.0d0
|
|
210 continue
|
|
c
|
|
c compute right side vector in resulting linear equations
|
|
c
|
|
basl = dlog10(2.0d0)
|
|
do 240 i = low,igh
|
|
do 240 j = low,igh
|
|
tb = b(i,j)
|
|
ta = a(i,j)
|
|
if (ta .eq. 0.0d0) go to 220
|
|
ta = dlog10(dabs(ta)) / basl
|
|
220 continue
|
|
if (tb .eq. 0.0d0) go to 230
|
|
tb = dlog10(dabs(tb)) / basl
|
|
230 continue
|
|
wk(i,5) = wk(i,5) - ta - tb
|
|
wk(j,6) = wk(j,6) - ta - tb
|
|
240 continue
|
|
nr = igh-low+1
|
|
coef = 1.0d0/float(2*nr)
|
|
coef2 = coef*coef
|
|
coef5 = 0.5d0*coef2
|
|
nrp2 = nr+2
|
|
beta = 0.0d0
|
|
it = 1
|
|
c
|
|
c start generalized conjugate gradient iteration
|
|
c
|
|
250 continue
|
|
ew = 0.0d0
|
|
ewc = 0.0d0
|
|
gamma = 0.0d0
|
|
do 260 i = low,igh
|
|
gamma = gamma + wk(i,5)*wk(i,5) + wk(i,6)*wk(i,6)
|
|
ew = ew + wk(i,5)
|
|
ewc = ewc + wk(i,6)
|
|
260 continue
|
|
gamma = coef*gamma - coef2*(ew**2 + ewc**2)
|
|
+ - coef5*(ew - ewc)**2
|
|
if (it .ne. 1) beta = gamma / pgamma
|
|
t = coef5*(ewc - 3.0d0*ew)
|
|
tc = coef5*(ew - 3.0d0*ewc)
|
|
do 270 i = low,igh
|
|
wk(i,2) = beta*wk(i,2) + coef*wk(i,5) + t
|
|
cperm(i) = beta*cperm(i) + coef*wk(i,6) + tc
|
|
270 continue
|
|
c
|
|
c apply matrix to vector
|
|
c
|
|
do 300 i = low,igh
|
|
kount = 0
|
|
sum = 0.0d0
|
|
do 290 j = low,igh
|
|
if (a(i,j) .eq. 0.0d0) go to 280
|
|
kount = kount+1
|
|
sum = sum + cperm(j)
|
|
280 continue
|
|
if (b(i,j) .eq. 0.0d0) go to 290
|
|
kount = kount+1
|
|
sum = sum + cperm(j)
|
|
290 continue
|
|
wk(i,3) = float(kount)*wk(i,2) + sum
|
|
300 continue
|
|
do 330 j = low,igh
|
|
kount = 0
|
|
sum = 0.0d0
|
|
do 320 i = low,igh
|
|
if (a(i,j) .eq. 0.0d0) go to 310
|
|
kount = kount+1
|
|
sum = sum + wk(i,2)
|
|
310 continue
|
|
if (b(i,j) .eq. 0.0d0) go to 320
|
|
kount = kount+1
|
|
sum = sum + wk(i,2)
|
|
320 continue
|
|
wk(j,4) = float(kount)*cperm(j) + sum
|
|
330 continue
|
|
sum = 0.0d0
|
|
do 340 i = low,igh
|
|
sum = sum + wk(i,2)*wk(i,3) + cperm(i)*wk(i,4)
|
|
340 continue
|
|
if(sum.eq.0.0d0) return
|
|
alpha = gamma / sum
|
|
c
|
|
c determine correction to current iterate
|
|
c
|
|
cmax = 0.0d0
|
|
do 350 i = low,igh
|
|
cor = alpha * wk(i,2)
|
|
if (dabs(cor) .gt. cmax) cmax = dabs(cor)
|
|
wk(i,1) = wk(i,1) + cor
|
|
cor = alpha * cperm(i)
|
|
if (dabs(cor) .gt. cmax) cmax = dabs(cor)
|
|
cscale(i) = cscale(i) + cor
|
|
350 continue
|
|
if (cmax .lt. 0.5d0) go to 370
|
|
do 360 i = low,igh
|
|
wk(i,5) = wk(i,5) - alpha*wk(i,3)
|
|
wk(i,6) = wk(i,6) - alpha*wk(i,4)
|
|
360 continue
|
|
pgamma = gamma
|
|
it = it+1
|
|
if (it .le. nrp2) go to 250
|
|
c
|
|
c end generalized conjugate gradient iteration
|
|
c
|
|
370 continue
|
|
do 380 i = low,igh
|
|
ir = wk(i,1) + dsign(0.5d0,wk(i,1))
|
|
wk(i,1) = ir
|
|
jc = cscale(i) + dsign(0.5d0,cscale(i))
|
|
cscale(i) = jc
|
|
380 continue
|
|
c
|
|
c scale a and b
|
|
c
|
|
do 400 i = 1,igh
|
|
ir = wk(i,1)
|
|
fi = 2.0d0**ir
|
|
if (i .lt. low) fi = 1.0d0
|
|
do 400 j =low,n
|
|
jc = cscale(j)
|
|
fj = 2.0d0**jc
|
|
if (j .le. igh) go to 390
|
|
if (i .lt. low) go to 400
|
|
fj = 1.0d0
|
|
390 continue
|
|
a(i,j) = a(i,j)*fi*fj
|
|
b(i,j) = b(i,j)*fi*fj
|
|
400 continue
|
|
410 continue
|
|
return
|
|
c
|
|
c last line of scaleg
|
|
c
|
|
end
|