Subversion Repositories openrisc

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