Subversion Repositories scarts

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