integer ind(2001),ind2(2001)
      integer nn(10),pivrow
      real x(2001),y(2001)
      real a(2001),b(2001)
      real m(2001,5)
      real t1
      real time(10,20)
      logical lbig
c
c     numbn is the number of different n's to be timed
c
      data numbn /7/
c
c     nrout is the number of different routines to be timed
c
      data nrout /11/
c
c     ntimes is the number times to average over
c
      data ntimes /500/
c
c size of ldm must be at least twice nn(numbn)
      ldm = 2001
c      call link(" o=indird,unit6=(o,create,text)//")
      write(6,*)' Times will be averaged over ',ntimes,' calls.'
c
c     compute the overhead of the clock
c
      time0 = second(foo)
      do 1 j = 1, ntimes
         call dummy(j,x,y,a,b,m,ldm,ind,ind2,t1)
    1 continue
      time0 = second(foo) - time0
c
c main loop .. initialized every pass to avoid overflow.
c
      do 190 in = 1,numbn
         call init(nn,x,y,a,b,m,ldm,ind,ind2,t1)
         n = nn(in)
c
c     sparse saxpy
c
c         write(6,*)'n',n
c         write(6,*)' s-saxpy'
         time1 = second(foo)
         do 11 j = 1, ntimes
            k = 1
            do 10 i = 1, n
               x(ind(i)) = x(ind(i)) + a(k)*t1
               k=k+1
   10       continue
            call dummy(j,x,y,a,b,m,ldm,ind,ind2,t1)
   11    continue
         time1 = second(foo) - time1 - time0
         time(in,1) = time1/ntimes
         time2 = second(foo)
         do 21 j = 1, ntimes
            k = 1
            do 20 i = 1, n
               x(k) = x(k) + a(ind(i))*t1
               k=k+1
   20       continue
            call dummy(j,x,y,a,b,m,ldm,ind,ind2,t1)
   21    continue
         time2 = second(foo) - time2 - time0
         time(in,2) = time2/ntimes
c
c     sparse sdot
c
c         write(6,*)' s-sdot'
         time3 = second(foo)
         do 31 j = 1, ntimes
            t1 = 0.0
            k = 1
            do 30 i = 1, n
               t1 = t1 + a(k)*x(ind(i))
               k = k + 1
   30       continue
            call dummy(j,x,y,a,b,m,ind,ind2,t1)
   31    continue
         time3 = second(foo) - time3 - time0
         time(in,3) = time3/ntimes
c
c     normal saxpy like loop
c
c         write(6,*)' saxpy'
         time4 = second(foo)
         do 41 j = 1, ntimes
            do 40 i = 1,n
               a(i) = a(i) + t1*y(i)
   40       continue
            call dummy(j,x,y,a,b,m,ind,ind2,t1)
   41    continue
         time4 = second(foo) - time4 - time0
         time(in,4) = time4/ntimes
c
c     normal sdot like loop
c
c         write(6,*)' sdot'
         time5 = second(foo)
         do 51 j = 1, ntimes
            t1 = 0.0
            do 50 i = 1,n
               t1 = t1 + a(i)*b(i)
   50       continue
            call dummy(j,x,y,a,b,m,ind,ind2,t1)
   51    continue
         time5 = second(foo) - time5 - time0
         time(in,5) = time5/ntimes
c
c
c     gather loop
c
c         write(6,*)' gather'
         time6 = second(foo)
         do61 j = 1, ntimes
            do 60 i = 1,n
               a(i) = x(ind(i))
   60       continue
            call dummy(j,x,y,a,b,m,ind,ind2,t1)
   61    continue
         time6 = second(foo) - time6 - time0
         time(in,6) = time6/ntimes
c
c     scatter loop
c
c         write(6,*)' scatter'
         time7 = second(foo)
         do 71 j = 1, ntimes
            do 70 i = 1,n
               x(ind(i)) = a(i)
   70       continue
            call dummy(j,x,y,a,b,m,ind,ind2,t1)
   71    continue
         time7 = second(foo) - time7 - time0
         time(in,7) = time7/ntimes
c
c     more complicated sparse saxpy
c
c         write(6,*)' c-s-saxpy'
         time8 = second(foo)
         do 81 j = 1, ntimes
            ibase = 0
            jj = 3
            do 80 l = 1, n
               ioffs = ibase + l + 3
               m(ind(l),jj) = m(ind(l),jj) + a(ioffs)*t1
   80       continue
            call dummy(j,x,y,a,b,m,ind,ind2,t1)
   81    continue
         time8 = second(foo) - time8 - time0
         time(in,8) = time8/ntimes
c
c more complicated sparse sdot
c
c         write(6,*)' c-s-sdot'
         time9 = second(foo)
         do 101 j = 1, ntimes
            t1 = 0.0
            jj = 3
            inci = 0
            do 100 i = 1, n
               ioff = inci + i
               t1 = t1 + a(ioff)*m(ind(i),jj)
  100       continue
            call dummy(j,x,y,a,b,m,ind,ind2,t1)
  101    continue
         time9 = second(foo)- time9 - time0
         time(in,9) = time9/ntimes
c
c     double indirect address
c
c         write(6,*)' 2-ind-add'
         time10 = second(foo)
         do 121 j = 1, ntimes
            k = 1
            do 120 ii = 1, n
               idx = ii + ind2(ind(ii))
               a(idx) = a(idx) + a(k)
               k = k + 1
  120       continue
            call dummy(j,x,y,a,b,m,ind,ind2,t1)
  121    continue
         time10 = second(foo) - time10 - time0
         time(in,10) = time10/ntimes
c
c     more complicated double indirect address.   ... ma30 inner-loop
c
c         write(6,*)' c-2-ind-add'
         time12 = second(foo)
         do 132 j = 1, ntimes
            call set(n,j,lbig,ijpos,idrop,iop,big,tol,ind,ind2)
            call reset(n,ind)
            call dummy(j,x,y,a,b,m,ind,ind2,t1)
  132    continue
         time12 = second(foo) - time12
         time11 = second(foo)
         do 131 j = 1, ntimes
            call set(n,j,lbig,ijpos,idrop,iop,big,tol,ind,ind2)
            do 130 ii = 1, n
               jj = ind(ii)
               if( ind2(jj) .gt. 0 ) go to 130
                  iop= iop + 1
                  pivrow = ijpos - ind2(jj)
                  x(ii) = x(ii) + t1*x(pivrow)
                  if( lbig ) big = max1(abs(x(ii)),big)
                  if( abs(x(ii)) .lt. tol ) idrop = idrop + 1
                  ind(pivrow) = -ind(pivrow)
  130       continue
            call reset(n,ind)
            call dummy(j,x,y,a,b,m,ind,ind2,t1)
  131    continue
c           write(6,*) ' n,iop',n,iop
         time11 = second(foo) - time11 - time12
         time(in,11) = time11/ntimes
  190 continue
c
c     compute r(infinity)
c
      do 200 i = 1,nrout
         if(  time(numbn,i) .ne. time(numbn-1,i) ) then
            time(numbn+1,i) = ( nn(numbn) - nn(numbn-1) ) /
     $                        ( time(numbn,i) - time(numbn-1,i) )
         else
            time(numbn+1,i) = 0.0
         endif
  200 continue
c
c     compute n sub 1/2 and scale r(infinity)
c
      do 300 i = 1,nrout
         time(numbn+2,i) = time(2,i)*time(numbn+1,i)-nn(2)
         time(numbn+1,i) = 2.e-6*time(numbn+1,i)
  300 continue
         time(numbn+1,6) = time(numbn+1,6)/2.0
         time(numbn+1,7) = time(numbn+1,7)/2.0
         time(numbn+1,11) = time(numbn+1,11)/2.0
      write(6,*)'length s-saxpy    s-saxpy-2    s-sdot       saxpy'
      call outp(1,4,nn,time,numbn)
      write(6,*) ' '
      write(6,*) 'length   sdot      gather     scatter     c-s-saxpy'
      call outp(5,8,nn,time,numbn)
      write(6,*) ' '
      write(6,*) 'length c-s-sdot  2-ind-add  c-2-ind-add'
      call outp(9,11,nn,time,numbn)
      stop
      end
      subroutine outp(ist,iend,nn,time,numbn)
      real time(10,20)
      integer nn(10)
      do 401 i = 1,numbn
         write(6,400)nn(i),(time(i,j),j=ist,iend)
  400    format(i5,(8e12.4))
  401 continue
      write(6,410)(time(numbn+1,j),j=ist,iend)
  410 format(/'r inf',(5f12.4))
      write(6,420)(time(numbn+2,j),j=ist,iend)
  420 format('nhalf',(5f12.4))
      return
      end
      subroutine init(nn,x,y,a,b,m,ldm,ind,ind2,t1)
      real x(2001),y(2001),a(2001),b(2001),m(ldm,5),t1
      integer nn(10),ind(2001),ind2(2001)
c
      nn(1) = 5
      nn(2) = 10
      nn(3) = 20
      nn(4) = 30
      nn(5) = 100
      nn(6) = 500
      nn(7) = 1000
c
c     initialize the world
c
      do 10 i = 1,nn(7)
         n = nn(7)
         x(i) = 1.0/float(i)
         y(i) = 1.0/float(i)
         a(i) = 1.0/float(i)
         b(i) = 1.0/float(i)
         ind(i) = nn(7) - i + 1
         ind2(i) = i
         do 20 j = 1,5
            m(i,j) = float(i) + float(n-j+1)/float(n)
   20    continue
   10 continue
      t1 = .5
      return
      end
      subroutine dummy(j,x,y,a,b,m,ldm,ind,ind2,t1)
c
      real x(2001),y(2001),a(2001),b(2001),m(ldm,5),t1
      integer ind(2001),ind2(2001)
c
      return
      end
      subroutine set(n,j,lbig,ijpos,idrop,iop,big,tol,ind,ind2)
      integer ind(*),ind2(*)
      logical lbig
            lbig = .true.
            ijpos = n
            idrop = 0
            iop = 0
            big = 0.0
            tol = 1.0e-8
            if (j.gt.1) go to 500
            do 100 k=1,n,2
              ind2(k) = - iabs(ind2(k))
 100        continue
            do 200 k=1,n
              ind(k) = n-k+1
 200        continue
            do 300 k=n+1,2*n
              ind(k) = 2*n-k+1
 300        continue
 500  return
      end
      subroutine reset(n,ind)
      integer ind(*)
        do 100 k=n+1,2*n
          ind(k) = iabs(ind(k))
 100    continue
      return
      end

.