gcc/testsuite/g77.f-torture/compile/980310-3.f

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