1c { dg-do compile }
2c
3c This demonstrates a problem with g77 and pic on x86 where
4c egcs 1.0.1 and earlier will generate bogus assembler output.
5c unfortunately, gas accepts the bogus acssembler output and
6c generates code that almost works.
7c
8
9
10C Date: Wed, 17 Dec 1997 23:20:29 +0000
11C From: Joao Cardoso <jcardoso@inescn.pt>
12C To: egcs-bugs@cygnus.com
13C Subject: egcs-1.0 f77 bug on OSR5
14C When trying to compile the Fortran file that I enclose bellow,
15C I got an assembler error:
16C
17C ./g77 -B./ -fpic -O -c scaleg.f
18C /usr/tmp/cca002D8.s:123:syntax error at (
19C
20C ./g77 -B./ -fpic -O0 -c scaleg.f
21C /usr/tmp/cca002EW.s:246:invalid operand combination: leal
22C
23C Compiling without the -fpic flag runs OK.
24
25      subroutine scaleg (n,ma,a,mb,b,low,igh,cscale,cperm,wk)
26c
27c     *****parameters:
28      integer igh,low,ma,mb,n
29      double precision a(ma,n),b(mb,n),cperm(n),cscale(n),wk(n,6)
30c
31c     *****local variables:
32      integer i,ir,it,j,jc,kount,nr,nrp2
33      double precision alpha,basl,beta,cmax,coef,coef2,coef5,cor,
34     *                 ew,ewc,fi,fj,gamma,pgamma,sum,t,ta,tb,tc
35c
36c     *****fortran functions:
37      double precision dabs, dlog10, dsign
38c     float
39c
40c     *****subroutines called:
41c     none
42c
43c     ---------------------------------------------------------------
44c
45c     *****purpose:
46c     scales the matrices a and b in the generalized eigenvalue
47c     problem a*x = (lambda)*b*x such that the magnitudes of the
48c     elements of the submatrices of a and b (as specified by low
49c     and igh) are close to unity in the least squares sense.
50c     ref.:  ward, r. c., balancing the generalized eigenvalue
51c     problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981,
52c     141-152.
53c
54c     *****parameter description:
55c
56c     on input:
57c
58c       ma,mb   integer
59c               row dimensions of the arrays containing matrices
60c               a and b respectively, as declared in the main calling
61c               program dimension statement;
62c
63c       n       integer
64c               order of the matrices a and b;
65c
66c       a       real(ma,n)
67c               contains the a matrix of the generalized eigenproblem
68c               defined above;
69c
70c       b       real(mb,n)
71c               contains the b matrix of the generalized eigenproblem
72c               defined above;
73c
74c       low     integer
75c               specifies the beginning -1 for the rows and
76c               columns of a and b to be scaled;
77c
78c       igh     integer
79c               specifies the ending -1 for the rows and columns
80c               of a and b to be scaled;
81c
82c       cperm   real(n)
83c               work array.  only locations low through igh are
84c               referenced and altered by this subroutine;
85c
86c       wk      real(n,6)
87c               work array that must contain at least 6*n locations.
88c               only locations low through igh, n+low through n+igh,
89c               ..., 5*n+low through 5*n+igh are referenced and
90c               altered by this subroutine.
91c
92c     on output:
93c
94c       a,b     contain the scaled a and b matrices;
95c
96c       cscale  real(n)
97c               contains in its low through igh locations the integer
98c               exponents of 2 used for the column scaling factors.
99c               the other locations are not referenced;
100c
101c       wk      contains in its low through igh locations the integer
102c               exponents of 2 used for the row scaling factors.
103c
104c     *****algorithm notes:
105c     none.
106c
107c     *****history:
108c     written by r. c. ward.......
109c     modified 8/86 by bobby bodenheimer so that if
110c       sum = 0 (corresponding to the case where the matrix
111c       doesn't need to be scaled) the routine returns.
112c
113c     ---------------------------------------------------------------
114c
115      if (low .eq. igh) go to 410
116      do 210 i = low,igh
117         wk(i,1) = 0.0d0
118         wk(i,2) = 0.0d0
119         wk(i,3) = 0.0d0
120         wk(i,4) = 0.0d0
121         wk(i,5) = 0.0d0
122         wk(i,6) = 0.0d0
123         cscale(i) = 0.0d0
124         cperm(i) = 0.0d0
125  210 continue
126c
127c     compute right side vector in resulting linear equations
128c
129      basl = dlog10(2.0d0)
130      do 240 i = low,igh
131         do 240 j = low,igh ! { dg-warning "Obsolescent feature: Shared DO termination" }
132            tb = b(i,j)
133            ta = a(i,j)
134            if (ta .eq. 0.0d0) go to 220
135            ta = dlog10(dabs(ta)) / basl
136  220       continue
137            if (tb .eq. 0.0d0) go to 230
138            tb = dlog10(dabs(tb)) / basl
139  230       continue
140            wk(i,5) = wk(i,5) - ta - tb
141            wk(j,6) = wk(j,6) - ta - tb
142  240 continue
143      nr = igh-low+1
144      coef = 1.0d0/float(2*nr)
145      coef2 = coef*coef
146      coef5 = 0.5d0*coef2
147      nrp2 = nr+2
148      beta = 0.0d0
149      it = 1
150c
151c     start generalized conjugate gradient iteration
152c
153  250 continue
154      ew = 0.0d0
155      ewc = 0.0d0
156      gamma = 0.0d0
157      do 260 i = low,igh
158         gamma = gamma + wk(i,5)*wk(i,5) + wk(i,6)*wk(i,6)
159         ew = ew + wk(i,5)
160         ewc = ewc + wk(i,6)
161  260 continue
162      gamma = coef*gamma - coef2*(ew**2 + ewc**2)
163     +        - coef5*(ew - ewc)**2
164      if (it .ne. 1) beta = gamma / pgamma
165      t = coef5*(ewc - 3.0d0*ew)
166      tc = coef5*(ew - 3.0d0*ewc)
167      do 270 i = low,igh
168         wk(i,2) = beta*wk(i,2) + coef*wk(i,5) + t
169         cperm(i) = beta*cperm(i) + coef*wk(i,6) + tc
170  270 continue
171c
172c     apply matrix to vector
173c
174      do 300 i = low,igh
175         kount = 0
176         sum = 0.0d0
177         do 290 j = low,igh
178            if (a(i,j) .eq. 0.0d0) go to 280
179            kount = kount+1
180            sum = sum + cperm(j)
181  280       continue
182            if (b(i,j) .eq. 0.0d0) go to 290
183            kount = kount+1
184            sum = sum + cperm(j)
185  290    continue
186         wk(i,3) = float(kount)*wk(i,2) + sum
187  300 continue
188      do 330 j = low,igh
189         kount = 0
190         sum = 0.0d0
191         do 320 i = low,igh
192            if (a(i,j) .eq. 0.0d0) go to 310
193            kount = kount+1
194            sum = sum + wk(i,2)
195  310       continue
196            if (b(i,j) .eq. 0.0d0) go to 320
197            kount = kount+1
198            sum = sum + wk(i,2)
199  320    continue
200         wk(j,4) = float(kount)*cperm(j) + sum
201  330 continue
202      sum = 0.0d0
203      do 340 i = low,igh
204         sum = sum + wk(i,2)*wk(i,3) + cperm(i)*wk(i,4)
205  340 continue
206      if(sum.eq.0.0d0) return
207      alpha = gamma / sum
208c
209c     determine correction to current iterate
210c
211      cmax = 0.0d0
212      do 350 i = low,igh
213         cor = alpha * wk(i,2)
214         if (dabs(cor) .gt. cmax) cmax = dabs(cor)
215         wk(i,1) = wk(i,1) + cor
216         cor = alpha * cperm(i)
217         if (dabs(cor) .gt. cmax) cmax = dabs(cor)
218         cscale(i) = cscale(i) + cor
219  350 continue
220      if (cmax .lt. 0.5d0) go to 370
221      do 360 i = low,igh
222         wk(i,5) = wk(i,5) - alpha*wk(i,3)
223         wk(i,6) = wk(i,6) - alpha*wk(i,4)
224  360 continue
225      pgamma = gamma
226      it = it+1
227      if (it .le. nrp2) go to 250
228c
229c     end generalized conjugate gradient iteration
230c
231  370 continue
232      do 380 i = low,igh
233         ir = wk(i,1) + dsign(0.5d0,wk(i,1))
234         wk(i,1) = ir
235         jc = cscale(i) + dsign(0.5d0,cscale(i))
236         cscale(i) = jc
237  380 continue
238c
239c     scale a and b
240c
241      do 400 i = 1,igh
242         ir = wk(i,1)
243         fi = 2.0d0**ir
244         if (i .lt. low) fi = 1.0d0
245         do 400 j =low,n ! { dg-warning "Obsolescent feature: Shared DO termination" }
246            jc = cscale(j)
247            fj = 2.0d0**jc
248            if (j .le. igh) go to 390
249            if (i .lt. low) go to 400
250            fj = 1.0d0
251  390       continue
252            a(i,j) = a(i,j)*fi*fj
253            b(i,j) = b(i,j)*fi*fj
254  400 continue
255  410 continue
256      return
257c
258c     last line of scaleg
259c
260      end
261