summaryrefslogtreecommitdiff
path: root/gcc/testsuite/gfortran.dg/g77/980310-3.f
blob: 565602378593e13b1e43068b8675576ff30e913c (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
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