From ribbens@ribbens.cs.vt.edu Thu Jun 20 09:55:21 1991
Return-Path: <ribbens@ribbens.cs.vt.edu>
Received: from ribbens.cs.vt.edu by CS.UTK.EDU with SMTP (5.61++/2.5.1s-UTK)
	id AA25151; Thu, 20 Jun 91 09:54:44 -0400
Subject: new sched.shar
To: dongarra@cs.utk.edu
Date: Thu, 20 Jun 91 9:55:50 EDT
X-Mailer: ELM [version 2.2 PL16]
From: ribbens@ribbens.cs.vt.edu (Cal Ribbens)
Message-Id: <677426153.AA10726@ribbens.cs.vt.edu>
Content-Type: text
Content-Length: 138299
Status: R

#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	README
#	cputm.c
#	d_and_c.f
#	data
#	example.f
#	ftsubs.f
#	ftsubs.graph.f
#	maindp.f
#	make1
#	makefile
#	newevdp0.f
#	oldtest
#	putq.c
#	second.f
#	stateig.f
#	statses.f
#	stseswait.f
#	stuffspawn.f
#	testrun
# This archive created: Thu Jun 20 09:39:32 1991
export PATH; PATH=/bin:$PATH
if test -f 'README'
then
	echo shar: will not over-write existing file "'README'"
else
cat << \SHAR_EOF > 'README'
To build the library needed to use schedule type: 

"make sched" or "make graph"

     make sched will produce the schedule library in file sched.a .
           This is the standard libarary for schedule. All examples
           begining with 'x' use this library.

     make graph will produce the trace version of schedule in file graph.a .
           This version produces output which can be viewed using the
           Schedule Trace Facility. All examples begining with 'g' use
           this library.

To run the examples type:

"make xtest; xtest" or "make gtest; gtest"

     This will make and run a test program which solves a triangular
     system of equations using the spawning capabilities of schedule.

"make xdandc; xdandc" or "make gdandc; gdandc"

     This will make and run a test program which implements a recursive
     divide and conquer technique using the spawning capabilities of schedule.

"make xeig; xeig" or "make geig; geig"

     This will make and run a test program which computes the eigenvalues
     of a symmetric tridiagonal eigenvalue problem using a divide and 
     conquer technique.

"make xwait; xwait" or "make gwait; gwait"

     This will make and run a test program which solves the same
     problem as above except that it uses the schedule wait command.

"make xexample; xexample" or "make gexample; gexample"

     This will make and run a test program which computes the
     inner product of two vectors.

A complete make and test can be carried out by doing:

testrun >& testout
diff oldtest testout (The only differences should be the times.)
SHAR_EOF
fi # end of overwriting check
if test -f 'cputm.c'
then
	echo shar: will not over-write existing file "'cputm.c'"
else
cat << \SHAR_EOF > 'cputm.c'
#include <sys/time.h>
#include <sys/resource.h>
#include <stdio.h>

long cputm_()
{
  long seconds, microseconds, milliseconds;
  struct rusage buffer;

  getrusage(RUSAGE_SELF, &buffer);
  seconds = buffer.ru_utime.tv_sec;
  microseconds = buffer.ru_utime.tv_usec;
  milliseconds = 1000*seconds + microseconds/1000;
  return(milliseconds);
}
SHAR_EOF
fi # end of overwriting check
if test -f 'd_and_c.f'
then
	echo shar: will not over-write existing file "'d_and_c.f'"
else
cat << \SHAR_EOF > 'd_and_c.f'
      integer a(256),klevl(8),myid(256)
      external top
      write(6,*) ' input nprocs nlevls '
      read (5,*) nprocs,nlevls
      do 50 j = 1,8
         klevl(j) = j
   50 continue
      call sched(nprocs,top,a,nlevls,klevl,myid)
      do 100 j = 1,2**nlevls-1
         write(6,*) a(j)
  100 continue
      stop
      end
      subroutine top(a,nlevls,klevl,myid)
      integer a(*),klevl(*),myid(*)
      external split
c      write(6,*) ' from top ' , a(1)
      jobtag = 1
      icango =0
      nchks = 0
      myid(1) = 1
      a(1) = 1
c
      call dep(jobtag,icango,nchks,mychkn)
      call putq(jobtag,split,myid,a,nlevls,klevl)
c
      return
      end
      subroutine split(myid,a,nlevls,klevl)
      integer a(*),klevl(*),myid(*),rnode
      external clone
c      write(6,*) ' from split ',a(1)
c
      if (klevl(1) .ge. nlevls) return
c
      lnode = 2*a(1)
      rnode = lnode + 1
      indx = lnode - a(1) + 1
      a(indx) = lnode
      a(indx+1) = rnode
      mytag = myid(a(1))
c
      call nxtag(mytag,jobtag)
      myid(lnode) = jobtag
      call spawn(mytag,jobtag,clone,myid,a(indx),nlevls,klevl(2))
c
      call nxtag(mytag,jobtag)
      myid(rnode) = jobtag
      call spawn(mytag,jobtag,clone,myid,a(indx+1),nlevls,klevl(2))
c
      return
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'data'
then
	echo shar: will not over-write existing file "'data'"
else
cat << \SHAR_EOF > 'data'
6
4
SHAR_EOF
fi # end of overwriting check
if test -f 'example.f'
then
	echo shar: will not over-write existing file "'example.f'"
else
cat << \SHAR_EOF > 'example.f'
      program main
      integer n, k
c
      external parprd
c
      real a(1000), b(1000), temp(50), sigma
      read (5,*) nprocs
      n = 100
      k = 10
c
      do 100 j = 1, n
        a(j) = j
        b(j) = 1
100   continue
c
      call sched(nprocs, parprd, n, k, a, b, temp, sigma)
c
      write (6,*) ' sigma = ', sigma
      stop
      end
c
c
      subroutine parprd(n, k, a, b, temp, sigma)
c
      integer n, k
      real a(*), b(*), temp(*), sigma
c
      integer m1, m2, index, j, jobtag, icango, ncheks, mychkn(2)
c
      external inprod, addup
      save m1, m2
c
      m1 = n/k
      index = 1
      do 200 j = 1, k - 1
        jobtag = j
        icango = 0
        ncheks = 1
        mychkn(1) = k + 1
        call dep(jobtag, icango, ncheks, mychkn)
        call putq(jobtag, inprod, m1, a(index), b(index), temp(j))
        index = index + m1
200   continue
c
      m2 = n - index + 1
      jobtag = k
      icango = 0
      ncheks = 1
      mychkn(1) = k + 1
      call dep(jobtag, icango, ncheks, mychkn)
      call putq(jobtag, inprod, m2, a(index), b(index), temp(k))
      index = index + m1
c
      jobtag = k + 1
      icango = k
      ncheks = 0
      call dep(jobtag, icango, ncheks, mychkn)
      call putq(jobtag, addup, k, sigma, temp)
c
      return
      end
c
c
      subroutine inprod(m, a, b, sigma)
      integer m
      real a(*), b(*), sigma
      sigma = 0.0
      do 100 j = 1, m
        sigma = sigma + a(j)*b(j)
100   continue
      return
      end
c
c
      subroutine addup(k, sigma, temp)
      integer k
      real sigma, temp(*)
      sigma = 0.0
      do 100 j = 1, k
        sigma = sigma + temp(j)
100   continue
      return
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'ftsubs.f'
then
	echo shar: will not over-write existing file "'ftsubs.f'"
else
cat << \SHAR_EOF > 'ftsubs.f'
      subroutine chekin(jobtag)
      integer jobtag
c***********************************************************************
c
c     this subroutine reports unit of computation labeled by 
c     jobtag has completed to all dependent nodes.  these nodes are 
c     recorded in parmq(i,jobtag)  where 6 .le. i .le. nchks+5
c     checkin consists of decrementing the value in each of these
c     locations by 1.  each of these is done in a critical section
c     protected by qlock(ichek)
c
c     if the value in parmq(2,ichek) is 0 where ichek is a process
c     dependent upon this one then ichek is placed on the readyq
c     by entering the critical section protected by rtlock.  the
c     pointer rtail to the tail of the readyq is  incremented
c     unless task done is to be recorded.  task done is placed on
c     the ready q and the pointer rtail left in place if nchks .eq. 0
c     is found. 
c
c     see the common block description in libopn for more detail.
c     
c***********************************************************************
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
c
c
c     common block description:
c
c     a complete common block description is given in the routine libopn
c
c****************************************************************************
c
c     check to see if this process has completed.  if not do not checkin
c
c     first ask if any kids spawned by jobtag
c
      if (parmq(4,jobtag) .ne. 0 .or. parmq(5,jobtag) .ne. 0 ) then
c
c          either kids have been spawned or ientry has been referenced
c          indicating reentry is required
c          find out how many are waiting to complete
c
         if (parmq(4,jobtag) .ne. 0) then
            call s_lock(qlock(jobtag))
               parmq(2,jobtag) = parmq(2,jobtag) - parmq(4,jobtag)
            call s_unlock(qlock(jobtag))
         endif
c
c          reset number of kids
c
         parmq(4,jobtag) = 0
c
c          update the number of times this procedure has been
c          entered
c
         parmq(1,jobtag) = parmq(1,jobtag) + 1
c
c          return without checkin if all the kids have not 
c          checked in to jobtag yet or if a reentry is required.  
c          process jobtag is not done in either case.
c
         if (parmq(2,jobtag) .ne. 0) return
c
c        if ientry has been called but jobtag is not waiting
c        on any kids then jobtag is placed back on the readyq
c
         if (parmq(5,jobtag) .ne. 0) then    
            iwrkr = mod(jobtag,mxces) + 1
            call s_lock(trlock(iwrkr))
               readyq(rtail(iwrkr),iwrkr) = jobtag
               rtail(iwrkr) = rtail(iwrkr) + 1
            call s_unlock(trlock(iwrkr))
            return
         endif
      endif
c
c     the process has completed so chekin proceeds
c
      nchks = parmq(3,jobtag)
c
c     if this is the final process (indicated by nchks .eq. 0) then
c     record task done.  do not advance the tail so task done sequence
c     is set.  all subsequent gtprb queries will leave rhead .eq. rtail
c     with readyq(rhead) .eq. done.
c
      if (nchks .eq. 0) then 
         do 20 iwrkr = 1,mxces
            call s_lock(trlock(iwrkr))
               readyq(rtail(iwrkr),iwrkr) = done
            call s_unlock(trlock(iwrkr))
   20    continue
         return
      endif
      do 50 j = 6,nchks+5
         mychek = parmq(j,jobtag)   
c                  
c     get unique access to the checkin node mychek
c     and checkin by decrementing the appropriate counter
c
         mchkgo = 1
         call s_lock(qlock(mychek))
            parmq(2,mychek) = parmq(2,mychek) - 1
            mchkgo = parmq(2,mychek)
         call s_unlock(qlock(mychek))
c
c     place mychek on readyq if parmq(2,mychek) is 0
c
         if (mchkgo .eq. 0 ) then
            iwrkr = mod(mychek,mxces) + 1
            call s_lock(trlock(iwrkr))
               readyq(rtail(iwrkr),iwrkr) = mychek
               rtail(iwrkr) = rtail(iwrkr) + 1
            call s_unlock(trlock(iwrkr))
         endif
   50 continue
      return
c
c     last card of chekin
c
      end
      subroutine dep(jobtag,icango,nchks,mychkn)
      integer jobtag,icango,nchks,mychkn(*)
c*************************************************************************
c
c  warning - this routine may only be used by driver in a static definition
c            of the data dependencies in the task.
c
c
c      usage
c      subroutine xxx(<parms>)
c      external yyy
c       .
c       .
c       .
c           call dep(jobtag,icango,nchks,mychkn)
c           call putq(jobtag,yyy,<parms2>)
c       .
c       .
c       .
c
c     this subroutine puts data dependencies for problem on the queue.  
c     no synchronization is necessary because each index of a column of 
c     parmq is associated with a jobtag specified by the user and 
c     associated with a unique unit of computation.  the arguments of 
c     dep specify a the data dependencies associated with the unit of  
c     computation labeled by jobtag and are placed in a column of parmq 
c     to the menue specified below.
c
c     jobtag is an integer specifying a unique schedulable process
c
c
c     icango is a positive integer specifying how many processes must check in
c            to this process before it can be placed on the readyq.
c
c     nchks is the number of processes that depend upon the completion of 
c           this process.
c
c     mychkn is an integer array specifying the jobtags of the processes
c            which depend upon completion of this process.
c
c*************************************************************************
c
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
c
c
c     common block description:
c
c     for a complete common block description see the subroutine libopn
c
c
c     place process jobtag on the problem queue
c     no synchronization required to update qtail since
c     only one program work executes this code.
c
      if ((jobtag .le. 0 .or. jobtag .gt. mxprcs) .or. 
     &     icango .lt. 0 .or. nchks .lt. 0) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      incorrect specification of dependencies '
         write(6,*) '      all DEP parameters must be non-negative'
         write(6,*) ' input was '
         write(6,*) '      jobtag ',jobtag ,'.... must be postitive '
         write(6,*) '              but less than ',mxprcs
         write(6,*) '      icango ',icango
         write(6,*) '      nchks  ',nchks
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         stop
      endif
c
      qtail = qtail + 1
      next = jobtag
      parmq(1,next) = 1
      parmq(2,next) = icango
      parmq(3,next) = nchks
      parmq(4,next) = 0
c
c     check to see that exactly one node has ncheks set to 0
c
       if (nchks .eq. 0 .and. done .eq. 0) then
           done = -2
       else
          if (nchks .eq. 0) done = 0
       endif
c
c     specify identifiers of processes which depend on this one
c     if there are too many abort
c
      if (nchks .gt. nslots - 5) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      attempt to place too many dependencies '
         write(6,*) '      on chekin list during call to dep '
         write(6,*) '      with jobtag ',jobtag
         write(6,*) ' '
         write(6,*) '      user tried to place ',nchks ,' dependencies '
         write(6,*) '      the maximum number is ',nslots - 5
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         stop
      endif
c
      do 50 j = 1,nchks
       parmq(j+5,next) = mychkn(j)
        if (mychkn(j) .le. 0) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      incorrect specification of dependencies '
         write(6,*) '      all mychkn entries must be positive'
         write(6,*) ' input was '
         write(6,*) '      mychkn(',j,') = ',mychkn(j)
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         stop
        endif
   50 continue
c
      return
c
c     last card of dep
c
      end
      integer function gtprb(id,jobtag)
c**************************************************************************
c
c     this function gets unique access to the head of the readyq
c     pointed to by   id    and then claims the pointer to the next 
c     schedulable process if there is one and returns with a nonzero 
c     value when there is a process to schedule. if there are no entries
c     in the readyq indexed by   id   then the remaning ready ques are 
c     polled in a round robin manner until schedulable process is found
c     or task done is recorded. if task done has been recorded the value 
c     zero is returned in gtprb.  if a nonzero value is returned in gtprb, 
c     the integer jobtag will contain the identifier of the unit of 
c     computation that is to be executed.
c     
c     input parameter
c
c        id  an integer specifying which readyq to access first
c            for work to do.
c
c     output parameters
c
c       jobtag an integer containing the next process to be executed
c
c       gtprb  the return value of this integer function is:
c       
c               0 if task done has been posted
c
c               nonzero if a schedulable process has been claimed.	
c
c
c***************************************************************************
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
c
c     common block description:
c
c     for a complete common block description see the routine libopn
c
c
      nspins = 0
      iwrkr = id
   10 continue
      mhead = -1
c
c     gain access to head of readyq.  if task done has not been recorded
c     then increment the head of the readyq. otherwise the head pointer
c     is left fixed so the next active process will receive task done.
c
      call s_lock(hrlock(iwrkr))
         if (rhead(iwrkr) .lt. rtail(iwrkr)) then
            mhead = rhead(iwrkr)
            rhead(iwrkr) = rhead(iwrkr) + 1
         endif
      call s_unlock(hrlock(iwrkr))
c
      if (mhead .gt. 0) then
c
c        there was a work unit on the readyq
c
         jobtag = readyq(mhead,iwrkr)
c
         if (jobtag .ne. done) then
c
c           record the subroutine call identifier in gtprb and return
c           the process identifier in jobtag.
c
            gtprb = parmq(1,jobtag)
            if (gtprb .gt. 1 .and. parmq(5,jobtag) .eq. 0) gtprb = -1
c
         else
c
c           task done has been indicated. request a return from subroutine work
c           by returning the value 0 in gtprb.
c
            gtprb = 0
c
         endif
      else 
c
         jobtag = readyq(rhead(iwrkr),iwrkr)
         if (jobtag .eq. done) then
c
c           task done has been posted
c
            gtprb = 0
c
         else
c
c           there was not any work on the readyq
c
            iwrkr = mod((iwrkr+1),mxces)  
            if (iwrkr .eq. 0) iwrkr = mxces
            nspins = nspins + 1
            if (mod(nspins,mxces) .eq. 0) call nops
            go to 10
c
         endif
      endif
      return
c
c     last card of gtprb
c
      end
      subroutine libopn(nproc)
      integer nproc
c************************************************************************
c
c     this routine sets locks and initializes variables
c     and then spawns nproc generic work routines.
c
c     nproc is a positive integer.  care should be taken to 
c           match nproc to the number of physical processors 
c           available.
c
c************************************************************************
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
c
      integer ispace(mxces), itag(mxces)
c
c     common block description:
c
c     common/qdata/
c     
c     parmq is a two dimensional integer array.  each column of 
c           this array represents a schedulable process.  a process is
c           identified by its jobtag which corresponds to a unique 
c           column of parmq.  a column of parmq has the following 
c           entries
c
c                   parmq(1,jobtag) = nentries
c                                     a nonzero integer. if  process jobtag
c                                     is on the readyq then this integer
c                                     is equal to the one plus number of times 
c                                     process jobtag has been entered.
c                                     thus when work executes this process
c                                     the integer is equal to the number
c                                     of times the process has been entered.
c      
c                   parmq(2,jobtag) = icango 
c                                     an integer specifying the number
c                                     of processes that must check in
c                                     before this process may scheduled
c                                     (ie. be placed on the ready queue)
c
c                   parmq(3,jobtag) = nchks
c                                     an integer specifying the number
c                                     of processes that this process 
c                                     must checkin to.  identifiers of
c                                     these processes are recorded below.
c                                     if nchks .eq. 0 then completion of 
c                                     this process signifies completion of
c                                     task.
c
c                   parmq(4,jobtag) = the number of kids spawned by this
c                                     process.  if this value is zero
c                                     then this process has not spawned
c                                     any subprocesses.
c
c                   parmq(5,jobtag) = entry_flag
c                                     has the value 1 if ientry was called
c                                     has the value 0  otherwise      
c
c                   parmq(6:6+nchks,jobtag) is reserved for identifiers of the
c                                                  nchks
c                                     processes that must wait for completion
c                                     of this process before they can execute.
c
c
c     phead pointer to head of parmq
c
c     intspn pointer to first spawned process.  all jobtags
c            with values greater than or equal to intspn will 
c            be spawned processes.
c
c     readyq a two dimensional integer array that holds the jobtags of those
c            processes that are ready to execute.  the k-th column of
c            this array serves as a readyq for the k-th work routine.
c            on executing gtprb, the k-th work  routine will look for work
c            in the k-th readyq first and then the others (round robin).
c            if readyq(j,id) .eq. done has been set then a return from 
c            subroutine work(*,*,id)  is indicated.
c
c     rhead  an integer array.  the i-th entry of rhead is a pointer to the 
c            head of the i-th column of readyq
c
c     rtail  an integer array.  the i-th entry of rtail is a pointer to the 
c            tail of the i-th column of readyq
c
c
c     common/sync/
c
c     qlock is an integer array of locks.  there is one lock for each 
c           column of parmq.  the purpose of this lock is to ensure
c           unique access to a column of parmq during the checkin operation.
c
c     hrlock is an integer lock.  the purpose of this lock is to ensure 
c            unique access to the pointer rhead to the head of the readyq.
c
c     trlock is an integer lock.  the purpose of this lock is to ensure 
c            unique access to the pointer rtail to the tail of the readyq.
c
c     done   is a unique non positive integer set in libopn to indicate
c            task done.
c
      if (nproc .gt. mxces-1) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      user asking for more physical processors'
         write(6,*) '      than are currently available on this system '
         write(6,*) '      the maximum avilable is ', mxces
         write(6,*) ' '
         stop
      endif
c
      ierr = m_set_procs(nproc)
c
c     if all goes well in the calls to dep (one and only one process
c     has ncheks equal to zero), done will be set to zero before the
c     processes begin execution.
c
      done = -1
c
c     set qlocks off
c     initialize readyq(*) = -1  to set done sequence 
c     initialize reentry indicator in parmq(5,*)
c
      do 50 j = 1,mxces
         call s_init_lock(hrlock(j))
         call s_init_lock(trlock(j))
         rhead(j)  = 1
         rtail(j)  = 1
         do 20 i = 1,iprcs
            readyq(i,j) = -1
   20    continue
   50 continue
c
      do 100 j = 1,mxprcs
         call s_init_lock(qlock(j))
         parmq(5,j) = 0
  100 continue
c
c     set readyq locks off
c
      call s_init_lock(qtlock)
c
c     initialize queue pointers
c
      phead = 1
      intspn = 1
      qtail = 2
      next = 1 
c
c     set lock on pointer to head of readyq so 
c     no process may start until all process and data dependencies
c     have been specified by the user supplied routine driver.
c
      do 150 j = 1,mxces
         call s_lock(hrlock(j))
  150 continue
c
c     now spawn virtual processors.  these generic work routines will
c     assume the identity of any schedulable process specified by driver.
c
C$DOACROSS share(ispace)
      do 200 j = 1,nproc
         call work(j,ispace(j))
  200 continue
      return
c
c     last card of libopn
c
      end
      subroutine nxtag(mypar,jobtag)
      integer mypar,jobtag
c*************************************************************************
c
c
c     this subroutine puts data dependencies for problem on the queue.  
c     no synchronization
c     is necessary because each index of a column of parmq is associated
c     with a jobtag specified by the user and associated with a unique
c     schedulable process.  the arguments of putq specify a process and are
c     placed in a column of jobq according to the menue specified in the
c     common block description given below.
c
c     jobtag is an integer specifying a unique schedulable process
c
c
c     icango is a positive integer specifying how many processes must check in
c            to this process before it can be placed on the readyq.
c
c     nchks is the number of processes that depend upon the completion of 
c           this process.
c
c     mchkin is an integer array specifying the jobtags of the processes
c            which depend upon completion of this process.
c
c*************************************************************************
c
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
c
c
c     common block description:
c
c     for a complete common block description see the subroutine libopn
c
c
c
c     place this process on the next slot in the problem queue
c
      call s_lock(qtlock)
         next = qtail
         qtail = qtail + 1
      call s_unlock(qtlock)
c
      if ( next .gt. mxprcs) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      user attempt to create too many processes'
         write(6,*) '      through dynamic spawning '
         write(6,*) '      the maximum allowed is ',mxprcs
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         stop
      endif
c
      jobtag = next
      parmq(1,next) = 1
      parmq(2,next) = 0
      parmq(3,next) = 1     
      parmq(6,next) = mypar
c
c     update the icango counter of the parent process
c     by adding 2 to parmq(2,mypar)... prevents race condition.
c     add 1 to the number of kids spawned by parent mypar
c
         call s_lock(qlock(mypar))
            parmq(2,mypar) = parmq(2,mypar) + 2
            parmq(4,mypar) = parmq(4,mypar) + 1
         call s_unlock(qlock(mypar))
c
c     set number of kids spawned by next to zero
c
         parmq(4,next) = 0
c
      return
c
c     last card of nxtag
c
      end
      subroutine start2
c
c     this routine allows parallel processing to start after user supplied
c     driver has completed by unlocking the head of the readyq
c
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
      logical nostrt
c
c
c     for common block description see subroutine libopn.
c
      if (done .ne. 0) then
         write(6,*) '*************SCHED USER ERROR********************'
         if (done .eq. -1 ) then
            write(6,*) '      no process has set nchks  equal to 0 '
         else
            write(6,*) '      more than one process has set nchks to 0 '
         endif
         write(6,*) '      SCHEDULE will not be able to terminate job'
         write(6,*) '      correctly '
         write(6,*) ' '
         write(6,*) '      check subroutine passed to initial call to'
         write(6,*) '      to see that at exactly one call to DEP  has '
         write(6,*) '      set nchks = 0 '
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         stop
      endif
c
      nostrt = .true.
      do 100 iwrkr = 1,mxces
         if (rhead(iwrkr) .ne. rtail(iwrkr)) nostrt = .false.
  100 continue
      if (nostrt) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      no process had an intitial icango of 0 '
         write(6,*) '      SCHEDULE could not begin '
         write(6,*) ' '
         write(6,*) '      check subroutine passed to initial call to'
         write(6,*) '      to see that at least one call to DEP  has '
         write(6,*) '      set icango = 0 '
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         stop
      endif
c
      intspn = qtail
      do 200 iwrkr = 1,mxces
         call s_unlock(hrlock(iwrkr))
  200 continue
c
      return
c
c     last card of start2
c
      end
      subroutine place(jobtag)
      integer jobtag
c*************************************************************************
c
c
c      this subroutine places a problem on the readyq 
c
c     jobtag is an integer specifying a unique schedulable process
c
c
c     icango is a positive integer specifying how many processes must check in
c            to this process before it can be placed on the readyq.
c
c     nchks is the number of processes that depend upon the completion of 
c           this process.
c
c     mchkin is an integer array specifying the jobtags of the processes
c            which depend upon completion of this process.
c
c*************************************************************************
c
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
c
c
c     common block description:
c
c     for a complete common block description see the subroutine libopn
c
c
c
c     place this process on readyq if icango is 0
c     when icango .eq. 0 this process does not depend on any
c     others.
c
      icango = parmq(2,jobtag)
      if (icango .eq. 0 ) then
         iwrkr = mod(jobtag,mxces) + 1
         call s_lock(trlock(iwrkr))
            readyq(rtail(iwrkr),iwrkr) = jobtag
            rtail(iwrkr) = rtail(iwrkr) + 1
         call s_unlock(trlock(iwrkr))
      endif
c
c     last card of place
c
      return
      end
      integer function ientry(jobtag,nentrs)
c
      integer jobtag
c*****************************************************************************
c
c     this routine will allow process jobtag to continue after
c     spawned processes have all checked in.  it should only be called if 
c     processes have been spawned by jobtag through the use of 
c     the subroutine spawn.  
c
c          go to (1000,2000,...,N000), ientry(jobtag,N)
c     1000 continue
c            .	
c            .	
c            .	
c          do 10 j = 1,nproc
c                 .
c                 . (set parameters to define spawned process)
c                 .
c             call nxtag(myid,jobtag)
c             call spawn(myid,jobtag,subname,<parms>)
c      10  continue
c          return
c     2000 continue
c            .
c            .
c            .
c          return
c     N000 continue
c           <statements>
c          return
c          end
c
c          this subroutine returns the number of times process jobtag
c          has been entered.  if that number is equal to the total
c          number nentrs of expected reentries then parmq(5,jobtag)
c          is set to 0 indicating no more reentries required.
c
c*****************************************************************************
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
c
c     report the entry point where process jobtag should resume
c     computation
c
      if (nentrs .lt. 2) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      user call to IENTRY  with number of    '
         write(6,*) '      labels in nentrs set less than  2 '
         write(6,*) '      from process ',jobtag
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         stop
      endif
c
      ientry = parmq(1,jobtag) 
      if (ientry .lt. nentrs) then
         parmq(5,jobtag) = nentrs
      else
         parmq(5,jobtag) = 0
      endif
c
      return
c
c     last card of ientry
c
      end
      logical function wait(jobtag,ienter)
c
      integer jobtag,ienter
c*****************************************************************************
c
c     this routine will allow process jobtag to continue after
c     spawned processes have all checked in.  it should only be called if 
c     processes have been spawned by jobtag through the use of 
c     the subroutine spawn.  this routine must be used in conjunction with
c     subroutine prtspn.  the required syntax is 
c
c          go to (1000,...,L000,...,N000), ientry(mytag,N)
c     1000 continue
c            .	
c            .	
c            .	
c          do 100 j = 1,nproc
c                 .
c                 . (set parameters to define spawned process)
c                 .
c             call nxtag(mytag,jobtag)
c             call spawn(mytag,jobtag,subname,<parms>)
c      100 continue
c          label = L
c          if (wait(jobtag,label)) return
c     L000 continue
c            .
c            .
c            .
c
c     if this subroutine returns a value of .true. then the calling process
c     jobtag should issue a return.  if a value of .false. is returned then
c     the calling process jobtag should resume execution at the 
c     statement immediately following the reference to wait (ie. at L000 in
c     the example above.  a return value .true. indicates that some spawned
c     process has not yet completed and checked in.  a return value .false.
c     indicates all spawned processes have checked in.
c
c*****************************************************************************
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
c
c        
c     check the icango counter to see if all spawned processes (kids) 
c     have checked in.
c
      icango = 1
      call s_lock(qlock(jobtag))
         icango = parmq(2,jobtag) - parmq(4,jobtag)
      call s_unlock(qlock(jobtag))
c
      if (icango .eq. 0) then
c
c        all kids are done ... dont wait (ie return false)
c
         wait = .false.
c
c        record re_entry label where computation is to 
c        resume after wait is complete 
c
         parmq(1,jobtag) = ienter
c
         if (ienter .gt. parmq(5,jobtag)) then
            write(6,*) '*************SCHED USER ERROR*****************'
            write(6,*) '      executing SCHEDULE function WAIT '
            write(6,*) '      return label larger than the maximum '
            write(6,*) '      specified by user in call to ientry  '
            write(6,*) '      from process jobtag ',jobtag
            write(6,*) ' '
            write(6,*) '      the maximum reentry number is '
            write(6,*) '      ', parmq(5,jobtag)
            write(6,*) ' '
            write(6,*) ' EXECUTION TERMINATED BY SCHED '
            stop
         endif
c
c        set last re_entry indication (parmq(5,jobtag) = 0)
c        if this reentry point corresponds to last one
c        (recorded in parmq(5,jobtag) during call to ientry)
c
         if (ienter .eq. parmq(5,jobtag)) parmq(5,jobtag) = 0
c
      else
c
c        kids are not done 
c
         wait = .true.
c
c        a checkin will be made so set the number of 
c        entries to return label ienter - 1 to get
c        correct entry point after checkin
c
         parmq(1,jobtag) = ienter - 1
c
      endif
c
      return
c
c     last card of wait
c
      end
c
      subroutine lckasn(ilock)
      integer*1 ilock
      call s_init_lock(ilock)
      return
      end
c
      subroutine lockon(ilock)
      integer*1 ilock
      call s_lock(ilock)
      return
      end
c
      subroutine lockoff(ilock)
      integer*1 ilock
      call s_unlock(ilock)
      return
      end
c
      subroutine nops
      j = 1
      return
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'ftsubs.graph.f'
then
	echo shar: will not over-write existing file "'ftsubs.graph.f'"
else
cat << \SHAR_EOF > 'ftsubs.graph.f'
      subroutine chekin(jobtag)
      integer jobtag
c***********************************************************************
c
c     this subroutine reports unit of computation labeled by 
c     jobtag has completed to all dependent nodes.  these nodes are 
c     recorded in parmq(i,jobtag)  where 6 .le. i .le. nchks+5
c     checkin consists of decrementing the value in each of these
c     locations by 1.  each of these is done in a critical section
c     protected by qlock(ichek)
c
c     if the value in parmq(2,ichek) is 0 where ichek is a process
c     dependent upon this one then ichek is placed on the readyq
c     by entering the critical section protected by rtlock.  the
c     pointer rtail to the tail of the readyq is  incremented
c     unless task done is to be recorded.  task done is placed on
c     the ready q and the pointer rtail left in place if nchks .eq. 0
c     is found. 
c
c     see the common block description in libopn for more detail.
c     
c***********************************************************************
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      parameter (nbuffr = 500)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
      integer endgrf
      real igraph
      common /gphout/ endgrf,igraph(nslots,nbuffr)
c
c
c     common block description:
c
c     a complete common block description is given in the routine libopn
c
c****************************************************************************
c
c     check to see if this process has completed.  if not do not checkin
c
c     first ask if any kids spawned by jobtag
c
      if (parmq(4,jobtag) .ne. 0 .or. parmq(5,jobtag) .ne. 0 ) then
c
c          either kids have been spawned or ientry has been referenced
c          indicating reentry is required
c          find out how many are waiting to complete
c
         if (parmq(4,jobtag) .ne. 0) then
            call s_lock(qlock(jobtag))
               parmq(2,jobtag) = parmq(2,jobtag) - parmq(4,jobtag)
            call s_unlock(qlock(jobtag))
         endif
c
c          reset number of kids
c
         parmq(4,jobtag) = 0
c
c          update the number of times this procedure has been
c          entered
c
         parmq(1,jobtag) = parmq(1,jobtag) + 1
c
c          return without checkin if all the kids have not 
c          checked in to jobtag yet or if a reentry is required.  
c          process jobtag is not done in either case.
c
         if (parmq(2,jobtag) .ne. 0) return
c
c        if ientry has been called but jobtag is not waiting
c        on any kids then jobtag is placed back on the readyq
c
         if (parmq(5,jobtag) .ne. 0) then    
            iwrkr = mod(jobtag,mxces) + 1
            call s_lock(trlock(iwrkr))
               readyq(rtail(iwrkr),iwrkr) = jobtag
               rtail(iwrkr) = rtail(iwrkr) + 1
            call s_unlock(trlock(iwrkr))
            return
         endif
      endif
c
c     the process has completed so chekin proceeds
c
      call s_lock(qlock(mxprcs))
            if (endgrf .gt. nbuffr) call dump(endgrf,igraph)
            insrt = endgrf
            endgrf = endgrf + 1
      call s_unlock(qlock(mxprcs))
            if (jobtag .ge. intspn) then
               igraph(1,insrt) = 5
               igraph(2,insrt) = parmq(6,jobtag)
               igraph(3,insrt) = jobtag
               igraph(4,insrt) = second(foo)
            else
               igraph(1,insrt) = 2
               igraph(2,insrt) = jobtag
               igraph(3,insrt) = second(foo)
            endif
c
      nchks = parmq(3,jobtag)
c
c     if this is the final process (indicated by nchks .eq. 0) then
c     record task done.  do not advance the tail so task done sequence
c     is set.  all subsequent gtprb queries will leave rhead .eq. rtail
c     with readyq(rhead) .eq. done.
c
      if (nchks .eq. 0) then 
         do 20 iwrkr = 1,mxces
            call s_lock(trlock(iwrkr))
               readyq(rtail(iwrkr),iwrkr) = done
            call s_unlock(trlock(iwrkr))
   20    continue
         return
      endif
      do 50 j = 6,nchks+5
         mychek = parmq(j,jobtag)   
c                  
c     get unique access to the checkin node mychek
c     and checkin by decrementing the appropriate counter
c
         mchkgo = 1
         call s_lock(qlock(mychek))
            parmq(2,mychek) = parmq(2,mychek) - 1
            mchkgo = parmq(2,mychek)
         call s_unlock(qlock(mychek))
c
c     place mychek on readyq if parmq(2,mychek) is 0
c
         if (mchkgo .eq. 0 ) then
            iwrkr = mod(mychek,mxces) + 1
            call s_lock(trlock(iwrkr))
               readyq(rtail(iwrkr),iwrkr) = mychek
               rtail(iwrkr) = rtail(iwrkr) + 1
            call s_unlock(trlock(iwrkr))
         endif
   50 continue
      return
c
c     last card of chekin
c
      end
      subroutine dep(jobtag,icango,nchks,mychkn)
      integer jobtag,icango,nchks,mychkn(*)
c*************************************************************************
c
c  warning - this routine may only be used by driver in a static definition
c            of the data dependencies in the task.
c
c
c      usage
c      subroutine xxx(<parms>)
c      external yyy
c       .
c       .
c       .
c           call dep(jobtag,icango,nchks,mychkn)
c           call putq(jobtag,yyy,<parms2>)
c       .
c       .
c       .
c
c     this subroutine puts data dependencies for problem on the queue.  
c     no synchronization is necessary because each index of a column of 
c     parmq is associated with a jobtag specified by the user and 
c     associated with a unique unit of computation.  the arguments of 
c     dep specify a the data dependencies associated with the unit of  
c     computation labeled by jobtag and are placed in a column of parmq 
c     to the menue specified below.
c
c     jobtag is an integer specifying a unique schedulable process
c
c
c     icango is a positive integer specifying how many processes must check in
c            to this process before it can be placed on the readyq.
c
c     nchks is the number of processes that depend upon the completion of 
c           this process.
c
c     mychkn is an integer array specifying the jobtags of the processes
c            which depend upon completion of this process.
c
c*************************************************************************
c
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      parameter (nbuffr = 500)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
      integer endgrf
      real igraph
      common /gphout/ endgrf,igraph(nslots,nbuffr)
c
c
c     common block description:
c
c     for a complete common block description see the subroutine libopn
c
c
c     place process jobtag on the problem queue
c     no synchronization required to update qtail since
c     only one program work executes this code.
c
      if ((jobtag .le. 0 .or. jobtag .gt. mxprcs) .or. 
     &     icango .lt. 0 .or. nchks .lt. 0) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      incorrect specification of dependencies '
         write(6,*) '      all DEP parameters must be non-negative'
         write(6,*) ' input was '
         write(6,*) '      jobtag ',jobtag ,'.... must be postitive '
         write(6,*) '              but less than ',mxprcs
         write(6,*) '      icango ',icango
         write(6,*) '      nchks  ',nchks
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         call dump(endgrf,igraph)
         stop
      endif
c
      qtail = qtail + 1
      next = jobtag
      parmq(1,next) = 1
      parmq(2,next) = icango
      parmq(3,next) = nchks
      parmq(4,next) = 0
c
c     check to see that exactly one node has ncheks set to 0
c
       if (nchks .eq. 0 .and. done .eq. 0) then
           done = -2
       else
          if (nchks .eq. 0) done = 0
       endif
c
c     specify identifiers of processes which depend on this one
c     if there are too many abort
c
      if (nchks .gt. nslots - 5) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      attempt to place too many dependencies '
         write(6,*) '      on chekin list during call to dep '
         write(6,*) '      with jobtag ',jobtag
         write(6,*) ' '
         write(6,*) '      user tried to place ',nchks ,' dependencies '
         write(6,*) '      the maximum number is ',nslots - 5
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         call dump(endgrf,igraph)
         stop
      endif
c
      do 50 j = 1,nchks
       parmq(j+5,next) = mychkn(j)
        if (mychkn(j) .le. 0) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      incorrect specification of dependencies '
         write(6,*) '      all mychkn entries must be positive'
         write(6,*) ' input was '
         write(6,*) '      mychkn(',j,') = ',mychkn(j)
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         call dump(endgrf,igraph)
         stop
        endif
   50 continue
c
      call s_lock(qlock(mxprcs))
            if (endgrf .gt. nbuffr) call dump(endgrf,igraph)
            insrt = endgrf
            endgrf = endgrf + 1
      call s_unlock(qlock(mxprcs))
      igraph(1,insrt) = 0
      igraph(2,insrt) = jobtag
      igraph(3,insrt) = icango
      igraph(4,insrt) = nchks
      do 9001 jnsrt = 6,nchks + 5
         igraph(jnsrt-1,insrt) = parmq(jnsrt,next)
 9001 continue
c
      return
c
c     last card of dep
c
      end
      integer function gtprb(id,jobtag)
c**************************************************************************
c
c     this function gets unique access to the head of the readyq
c     pointed to by   id    and then claims the pointer to the next 
c     schedulable process if there is one and returns with a nonzero 
c     value when there is a process to schedule. if there are no entries
c     in the readyq indexed by   id   then the remaning ready ques are 
c     polled in a round robin manner until schedulable process is found
c     or task done is recorded. if task done has been recorded the value 
c     zero is returned in gtprb.  if a nonzero value is returned in gtprb, 
c     the integer jobtag will contain the identifier of the unit of 
c     computation that is to be executed.
c     
c     input parameter
c
c        id  an integer specifying which readyq to access first
c            for work to do.
c
c     output parameters
c
c       jobtag an integer containing the next process to be executed
c
c       gtprb  the return value of this integer function is:
c       
c               0 if task done has been posted
c
c               nonzero if a schedulable process has been claimed.	
c
c
c***************************************************************************
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      parameter (nbuffr = 500)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
      integer endgrf
      real igraph
      common /gphout/ endgrf,igraph(nslots,nbuffr)
c
c     common block description:
c
c     for a complete common block description see the routine libopn
c
c
      nspins = 0
      iwrkr = id
   10 continue
      mhead = -1
c
c     gain access to head of readyq.  if task done has not been recorded
c     then increment the head of the readyq. otherwise the head pointer
c     is left fixed so the next active process will receive task done.
c
      call s_lock(hrlock(iwrkr))
         if (rhead(iwrkr) .lt. rtail(iwrkr)) then
            mhead = rhead(iwrkr)
            rhead(iwrkr) = rhead(iwrkr) + 1
         endif
      call s_unlock(hrlock(iwrkr))
c
      if (mhead .gt. 0) then
c
c        there was a work unit on the readyq
c
         jobtag = readyq(mhead,iwrkr)
         call s_lock(qlock(mxprcs))
            if (endgrf .gt. nbuffr) call dump(endgrf,igraph)
            insrt = endgrf
            endgrf = endgrf + 1
         call s_unlock(qlock(mxprcs))
         if (jobtag .ge. intspn) then
            igraph(1,insrt) = 4
            igraph(2,insrt) = parmq(6,jobtag)
            igraph(3,insrt) = jobtag
            igraph(4,insrt) = second(foo)
         else
            igraph(1,insrt) = 1
            igraph(2,insrt) = jobtag
            igraph(3,insrt) = second(foo)
         endif
c
         if (jobtag .ne. done) then
c
c           record the subroutine call identifier in gtprb and return
c           the process identifier in jobtag.
c
            gtprb = parmq(1,jobtag)
            if (gtprb .gt. 1 .and. parmq(5,jobtag) .eq. 0) gtprb = -1
c
         else
c
c           task done has been indicated. request a return from subroutine work
c           by returning the value 0 in gtprb.
c
            gtprb = 0
c
         endif
      else 
c
         jobtag = readyq(rhead(iwrkr),iwrkr)
         if (jobtag .eq. done) then
c
c           task done has been posted
c
            gtprb = 0
c
         else
c
c           there was not any work on the readyq
c
            iwrkr = mod((iwrkr+1),mxces)  
            if (iwrkr .eq. 0) iwrkr = mxces
            nspins = nspins + 1
            if (mod(nspins,mxces) .eq. 0) call nops
            go to 10
c
         endif
      endif
      return
c
c     last card of gtprb
c
      end
      subroutine libopn(nproc)
      integer nproc
c************************************************************************
c
c     this routine sets locks and initializes variables
c     and then spawns nproc generic work routines.
c
c     nproc is a positive integer.  care should be taken to 
c           match nproc to the number of physical processors 
c           available.
c
c************************************************************************
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      parameter (nbuffr = 500)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
      integer endgrf
      real igraph
      common /gphout/ endgrf,igraph(nslots,nbuffr)
c
      integer ispace(mxces), itag(mxces)
c
c     common block description:
c
c     common/qdata/
c     
c     parmq is a two dimensional integer array.  each column of 
c           this array represents a schedulable process.  a process is
c           identified by its jobtag which corresponds to a unique 
c           column of parmq.  a column of parmq has the following 
c           entries
c
c                   parmq(1,jobtag) = nentries
c                                     a nonzero integer. if  process jobtag
c                                     is on the readyq then this integer
c                                     is equal to the one plus number of times 
c                                     process jobtag has been entered.
c                                     thus when work executes this process
c                                     the integer is equal to the number
c                                     of times the process has been entered.
c      
c                   parmq(2,jobtag) = icango 
c                                     an integer specifying the number
c                                     of processes that must check in
c                                     before this process may scheduled
c                                     (ie. be placed on the ready queue)
c
c                   parmq(3,jobtag) = nchks
c                                     an integer specifying the number
c                                     of processes that this process 
c                                     must checkin to.  identifiers of
c                                     these processes are recorded below.
c                                     if nchks .eq. 0 then completion of 
c                                     this process signifies completion of
c                                     task.
c
c                   parmq(4,jobtag) = the number of kids spawned by this
c                                     process.  if this value is zero
c                                     then this process has not spawned
c                                     any subprocesses.
c
c                   parmq(5,jobtag) = entry_flag
c                                     has the value 1 if ientry was called
c                                     has the value 0  otherwise      
c
c                   parmq(6:6+nchks,jobtag) is reserved for identifiers of the
c                                                  nchks
c                                     processes that must wait for completion
c                                     of this process before they can execute.
c
c
c     phead pointer to head of parmq
c
c     intspn pointer to first spawned process.  all jobtags
c            with values greater than or equal to intspn will 
c            be spawned processes.
c
c     readyq a two dimensional integer array that holds the jobtags of those
c            processes that are ready to execute.  the k-th column of
c            this array serves as a readyq for the k-th work routine.
c            on executing gtprb, the k-th work  routine will look for work
c            in the k-th readyq first and then the others (round robin).
c            if readyq(j,id) .eq. done has been set then a return from 
c            subroutine work(*,*,id)  is indicated.
c
c     rhead  an integer array.  the i-th entry of rhead is a pointer to the 
c            head of the i-th column of readyq
c
c     rtail  an integer array.  the i-th entry of rtail is a pointer to the 
c            tail of the i-th column of readyq
c
c
c     common/sync/
c
c     qlock is an integer array of locks.  there is one lock for each 
c           column of parmq.  the purpose of this lock is to ensure
c           unique access to a column of parmq during the checkin operation.
c
c     hrlock is an integer lock.  the purpose of this lock is to ensure 
c            unique access to the pointer rhead to the head of the readyq.
c
c     trlock is an integer lock.  the purpose of this lock is to ensure 
c            unique access to the pointer rtail to the tail of the readyq.
c
c     done   is a unique non positive integer set in libopn to indicate
c            task done.
c
c    common /gphout/
c
c    endgrf  is an integer pointing to the next available
c            slot in igraph
c
c    igraph  is a two dimensional integer array
c            used as a buffer for graphics output
c            each column of igraph records an event.
c
      if (nproc .gt. mxces-1) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      user asking for more physical processors'
         write(6,*) '      than are currently available on this system '
         write(6,*) '      the maximum avilable is ', mxces
         write(6,*) ' '
         call dump(endgrf,igraph)
         stop
      endif
c
      ierr = m_set_procs(nproc)
c
c     if all goes well in the calls to dep (one and only one process
c     has ncheks equal to zero), done will be set to zero before the
c     processes begin execution.
c
      done = -1
c
c     set qlocks off
c     initialize readyq(*) = -1  to set done sequence 
c     initialize reentry indicator in parmq(5,*)
c
      do 50 j = 1,mxces
         call s_init_lock(hrlock(j))
         call s_init_lock(trlock(j))
         rhead(j)  = 1
         rtail(j)  = 1
         do 20 i = 1,iprcs
            readyq(i,j) = -1
   20    continue
   50 continue
c
      do 100 j = 1,mxprcs
         call s_init_lock(qlock(j))
         parmq(5,j) = 0
  100 continue
c
c     set readyq locks off
c
      call s_init_lock(qtlock)
c
c     initialize queue pointers
c
      phead = 1
      intspn = 1
      qtail = 2
      next = 1 
      endgrf = 1
      open( file='trace.graph',unit=3)
c
c     set lock on pointer to head of readyq so 
c     no process may start until all process and data dependencies
c     have been specified by the user supplied routine driver.
c
      do 150 j = 1,mxces
         call s_lock(hrlock(j))
  150 continue
c
c     now spawn virtual processors.  these generic work routines will
c     assume the identity of any schedulable process specified by driver.
c
C$DOACROSS share(ispace)
      do 200 j = 1,nproc
         call work(j,ispace(j))
  200 continue
      call dump(endgrf,igraph)
      return
c
c     last card of libopn
c
      end
      subroutine nxtag(mypar,jobtag)
      integer mypar,jobtag
c*************************************************************************
c
c
c     this subroutine puts data dependencies for problem on the queue.  
c     no synchronization
c     is necessary because each index of a column of parmq is associated
c     with a jobtag specified by the user and associated with a unique
c     schedulable process.  the arguments of putq specify a process and are
c     placed in a column of jobq according to the menue specified in the
c     common block description given below.
c
c     jobtag is an integer specifying a unique schedulable process
c
c
c     icango is a positive integer specifying how many processes must check in
c            to this process before it can be placed on the readyq.
c
c     nchks is the number of processes that depend upon the completion of 
c           this process.
c
c     mchkin is an integer array specifying the jobtags of the processes
c            which depend upon completion of this process.
c
c*************************************************************************
c
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      parameter (nbuffr = 500)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
      integer endgrf
      real igraph
      common /gphout/ endgrf,igraph(nslots,nbuffr)
c
c
c     common block description:
c
c     for a complete common block description see the subroutine libopn
c
c
c
c     place this process on the next slot in the problem queue
c
      call s_lock(qtlock)
         next = qtail
         qtail = qtail + 1
      call s_unlock(qtlock)
c
      call s_lock(qlock(mxprcs))
         if (endgrf .gt. nbuffr) call dump(endgrf,igraph)
         insrt = endgrf
         endgrf = endgrf + 1
      call s_unlock(qlock(mxprcs))
      igraph(1,insrt) = 3
      igraph(2,insrt) = mypar
      igraph(3,insrt) = next
c
      if ( next .gt. mxprcs) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      user attempt to create too many processes'
         write(6,*) '      through dynamic spawning '
         write(6,*) '      the maximum allowed is ',mxprcs
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         call dump(endgrf,igraph)
         stop
      endif
c
      jobtag = next
      parmq(1,next) = 1
      parmq(2,next) = 0
      parmq(3,next) = 1     
      parmq(6,next) = mypar
c
c     update the icango counter of the parent process
c     by adding 2 to parmq(2,mypar)... prevents race condition.
c     add 1 to the number of kids spawned by parent mypar
c
         call s_lock(qlock(mypar))
            parmq(2,mypar) = parmq(2,mypar) + 2
            parmq(4,mypar) = parmq(4,mypar) + 1
         call s_unlock(qlock(mypar))
c
c     set number of kids spawned by next to zero
c
         parmq(4,next) = 0
c
      return
c
c     last card of nxtag
c
      end
      subroutine start2
c
c     this routine allows parallel processing to start after user supplied
c     driver has completed by unlocking the head of the readyq
c
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      parameter (nbuffr = 500)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
      logical nostrt
      integer endgrf
      real igraph
      common /gphout/ endgrf,igraph(nslots,nbuffr)
c
c
c     for common block description see subroutine libopn.
c
      if (done .ne. 0) then
         write(6,*) '*************SCHED USER ERROR********************'
         if (done .eq. -1 ) then
            write(6,*) '      no process has set nchks  equal to 0 '
         else
            write(6,*) '      more than one process has set nchks to 0 '
         endif
         write(6,*) '      SCHEDULE will not be able to terminate job'
         write(6,*) '      correctly '
         write(6,*) ' '
         write(6,*) '      check subroutine passed to initial call to'
         write(6,*) '      to see that at exactly one call to DEP  has '
         write(6,*) '      set nchks = 0 '
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         call dump(endgrf,igraph)
         stop
      endif
c
      nostrt = .true.
      do 100 iwrkr = 1,mxces
         if (rhead(iwrkr) .ne. rtail(iwrkr)) nostrt = .false.
  100 continue
      if (nostrt) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      no process had an intitial icango of 0 '
         write(6,*) '      SCHEDULE could not begin '
         write(6,*) ' '
         write(6,*) '      check subroutine passed to initial call to'
         write(6,*) '      to see that at least one call to DEP  has '
         write(6,*) '      set icango = 0 '
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         call dump(endgrf,igraph)
         stop
      endif
c
      intspn = qtail
      do 200 iwrkr = 1,mxces
         call s_unlock(hrlock(iwrkr))
  200 continue
c
      return
c
c     last card of start2
c
      end
      subroutine place(jobtag)
      integer jobtag
c*************************************************************************
c
c
c      this subroutine places a problem on the readyq 
c
c     jobtag is an integer specifying a unique schedulable process
c
c
c     icango is a positive integer specifying how many processes must check in
c            to this process before it can be placed on the readyq.
c
c     nchks is the number of processes that depend upon the completion of 
c           this process.
c
c     mchkin is an integer array specifying the jobtags of the processes
c            which depend upon completion of this process.
c
c*************************************************************************
c
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      parameter (nbuffr = 500)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
      integer endgrf
      real igraph
      common /gphout/ endgrf,igraph(nslots,nbuffr)
c
c
c     common block description:
c
c     for a complete common block description see the subroutine libopn
c
c
c
c     place this process on readyq if icango is 0
c     when icango .eq. 0 this process does not depend on any
c     others.
c
      icango = parmq(2,jobtag)
      if (icango .eq. 0 ) then
         iwrkr = mod(jobtag,mxces) + 1
         call s_lock(trlock(iwrkr))
            readyq(rtail(iwrkr),iwrkr) = jobtag
            rtail(iwrkr) = rtail(iwrkr) + 1
         call s_unlock(trlock(iwrkr))
      endif
c
c     last card of place
c
      return
      end
      integer function ientry(jobtag,nentrs)
c
      integer jobtag
c*****************************************************************************
c
c     this routine will allow process jobtag to continue after
c     spawned processes have all checked in.  it should only be called if 
c     processes have been spawned by jobtag through the use of 
c     the subroutine spawn.  
c
c          go to (1000,2000,...,N000), ientry(jobtag,N)
c     1000 continue
c            .	
c            .	
c            .	
c          do 10 j = 1,nproc
c                 .
c                 . (set parameters to define spawned process)
c                 .
c             call nxtag(myid,jobtag)
c             call spawn(myid,jobtag,subname,<parms>)
c      10  continue
c          return
c     2000 continue
c            .
c            .
c            .
c          return
c     N000 continue
c           <statements>
c          return
c          end
c
c          this subroutine returns the number of times process jobtag
c          has been entered.  if that number is equal to the total
c          number nentrs of expected reentries then parmq(5,jobtag)
c          is set to 0 indicating no more reentries required.
c
c*****************************************************************************
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      parameter (nbuffr = 500)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
      integer endgrf
      real igraph
      common /gphout/ endgrf,igraph(nslots,nbuffr)
c
c     report the entry point where process jobtag should resume
c     computation
c
      if (nentrs .lt. 2) then
         write(6,*) '*************SCHED USER ERROR********************'
         write(6,*) '      user call to IENTRY  with number of    '
         write(6,*) '      labels in nentrs set less than  2 '
         write(6,*) '      from process ',jobtag
         write(6,*) ' '
         write(6,*) ' EXECUTION TERMINATED BY SCHED '
         call dump(endgrf,igraph)
         stop
      endif
c
      ientry = parmq(1,jobtag) 
      if (ientry .lt. nentrs) then
         parmq(5,jobtag) = nentrs
      else
         parmq(5,jobtag) = 0
      endif
c
      return
c
c     last card of ientry
c
      end
      subroutine dump(endgrf,igraph)
      parameter (nslots = 30,nbuffr = 500)
      integer endgrf
      real igraph(nslots,nbuffr)
      integer ievent(nslots)
c
c     this routine writes graphics output to a file
c     and resets endgrf to 1
c
c      do 100 j = 1,endgrf
c         write(6,15)  (igraph(i,j),i = 1,nslots)
c 100  continue
c  15  format('dump: igraph',30f5.1)
      do 300 j = 1,endgrf-1
         do 302 i = 1,nslots
            ievent(i) = igraph(i,j)
  302    continue
         if( ievent(1) .eq. 0 )
     $            write(3,301) (ievent(i),i=1,ievent(4)+4)
         if( ievent(1) .eq. 1 )
     $               write(3,303) (ievent(i),i=1,2),igraph(3,j)
         if( ievent(1) .eq. 2 )
     $               write(3,303) (ievent(i),i=1,2),igraph(3,j)
         if( ievent(1) .eq. 3 )
     $               write(3,301) (ievent(i),i=1,3)
         if( ievent(1) .eq. 4 )
     $               write(3,304) (ievent(i),i=1,3),igraph(4,j)
         if( ievent(1) .eq. 5 )
     $               write(3,304) (ievent(i),i=1,3),igraph(4,j)
  301    format(14i8)
  303    format(2i8,1pe16.8)
  304    format(3i8,1pe16.8)
  300 continue
c
      endgrf = 1
c
      return
c
c     last card of dump
c
      end
      logical function wait(jobtag,ienter)
c
      integer jobtag,ienter
c*****************************************************************************
c
c     this routine will allow process jobtag to continue after
c     spawned processes have all checked in.  it should only be called if 
c     processes have been spawned by jobtag through the use of 
c     the subroutine spawn.  this routine must be used in conjunction with
c     subroutine prtspn.  the required syntax is 
c
c          go to (1000,...,L000,...,N000), ientry(mytag,N)
c     1000 continue
c            .	
c            .	
c            .	
c          do 100 j = 1,nproc
c                 .
c                 . (set parameters to define spawned process)
c                 .
c             call nxtag(mytag,jobtag)
c             call spawn(mytag,jobtag,subname,<parms>)
c      100 continue
c          label = L
c          if (wait(jobtag,label)) return
c     L000 continue
c            .
c            .
c            .
c
c     if this subroutine returns a value of .true. then the calling process
c     jobtag should issue a return.  if a value of .false. is returned then
c     the calling process jobtag should resume execution at the 
c     statement immediately following the reference to wait (ie. at L000 in
c     the example above.  a return value .true. indicates that some spawned
c     process has not yet completed and checked in.  a return value .false.
c     indicates all spawned processes have checked in.
c
c*****************************************************************************
      parameter (mxprcs = 1000,iprcs = 120,mxces = 24,nslots = 30)
      parameter (nbuffr = 500)
      integer parmq,readyq,phead,intspn,rhead,rtail,done,qtail
      integer*1 qlock,hrlock,trlock,qtlock
      common /qdata/ parmq(nslots,mxprcs),phead,intspn,
     &               readyq(iprcs,mxces),rhead(mxces),rtail(mxces),
     &               qtail 
      common /qsync/ qlock(mxprcs),hrlock(mxces),trlock(mxces),
     &               done,qtlock
      integer endgrf
      real igraph
      common /gphout/ endgrf,igraph(nslots,nbuffr)
c
c        
c     check the icango counter to see if all spawned processes (kids) 
c     have checked in.
c
      icango = 1
      call s_lock(qlock(jobtag))
         icango = parmq(2,jobtag) - parmq(4,jobtag)
      call s_unlock(qlock(jobtag))
c
      if (icango .eq. 0) then
c
c        all kids are done ... dont wait (ie return false)
c
         wait = .false.
c
c        record re_entry label where computation is to 
c        resume after wait is complete 
c
         parmq(1,jobtag) = ienter
c
         if (ienter .gt. parmq(5,jobtag)) then
            write(6,*) '*************SCHED USER ERROR*****************'
            write(6,*) '      executing SCHEDULE function WAIT '
            write(6,*) '      return label larger than the maximum '
            write(6,*) '      specified by user in call to ientry  '
            write(6,*) '      from process jobtag ',jobtag
            write(6,*) ' '
            write(6,*) '      the maximum reentry number is '
            write(6,*) '      ', parmq(5,jobtag)
            write(6,*) ' '
            write(6,*) ' EXECUTION TERMINATED BY SCHED '
         call dump(endgrf,igraph)
            stop
         endif
c
c        set last re_entry indication (parmq(5,jobtag) = 0)
c        if this reentry point corresponds to last one
c        (recorded in parmq(5,jobtag) during call to ientry)
c
         if (ienter .eq. parmq(5,jobtag)) parmq(5,jobtag) = 0
c
      else
c
c        kids are not done 
c
         wait = .true.
c
c        a checkin will be made so set the number of 
c        entries to return label ienter - 1 to get
c        correct entry point after checkin
c
         parmq(1,jobtag) = ienter - 1
c
      endif
c
      return
c
c     last card of wait
c
      end
c
      subroutine lckasn(ilock)
      integer*1 ilock
      call s_init_lock(ilock)
      return
      end
c
      subroutine lockon(ilock)
      integer*1 ilock
      call s_lock(ilock)
      return
      end
c
      subroutine lockoff(ilock)
      integer*1 ilock
      call s_unlock(ilock)
      return
      end
c
      subroutine nops
      j = 1
      return
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'maindp.f'
then
	echo shar: will not over-write existing file "'maindp.f'"
else
cat << \SHAR_EOF > 'maindp.f'
      double precision esave(500),dsave(500)
      double precision s,t,t2,tnorm,res
      integer icase
c
      double precision d(500),e(500),q(500,500)
      integer n,ldq
      common /prbdef/n,ldq,q,d,e
      common/prfpms/ksect,kgran
      double precision enorm
      external treeql
c
      ldq = 500
c
      write(6,*) ' input nproc ... the number of processors '
      read(5,*) nproc
      do 9999 n = 100,100
      ksect = n/10
      kgran = n/20
      write(6,*) 'nprocs = ',nproc
      write(6,*) 'ksect = ',ksect, ' kgran = ',kgran
      write(6,*)'=============================='
      write(6,*)' n = ',n
      do 9998 icase = 1,1    
      write(6,*)'++++++++++++++++++++++++++++++'
      if( icase .eq. 1 ) write(6,*)' twos on diagonal'
      if( icase .eq. 2 ) write(6,*)' random numbers on diagonal'
      if( icase .eq. 3 ) write(6,*)' glued wilks eps = 1.d-8   '
      if( icase .eq. 4 ) write(6,*)' glued wilks eps = 1.d-14  '
      nsect = 2**ksect
      nd2 = n/2
      go to (112,113,114,115), icase
  112 do 13 i = 1,n
         d(i) = 2.0
         e(i) = 1
   13 continue
      go to 445
  113 do 14 i = 1,n
         d(i) = rand(foo)
         e(i) = rand(foo)
   14 continue
      go to 445
  114 continue
      eps = 1.e-8
      do 151 inum = 1,num
      do 15 i = 1,21
         d((inum-1)*21+i) = iabs(11-i)
         e((inum-1)*21+i) = 1
   15 continue
  151 continue
      go to 444
  115 continue
      eps = 1.e-14
      do 161 inum = 1,num
      do 16 i = 1,n
         d((inum-1)*21+i) = iabs(11-i)
         e((inum-1)*21+i) = 1
   16 continue
  161 continue
  444 continue
      do 443 inum = 1,num
         e(inum*21+1) = eps
  443 continue
  445 continue
      do 10 j = 1,n
         do 5 i = 1,n
              q(i,j) = 0.0
    5    continue
         q(j,j) = 1.0
         dsave(j) = d(j)
         esave(j) = e(j)
   10 continue
c
c     the test problem has been defined now comes the numerical
c     refinements that will avoid cancellation
c
      t1 = second(gtime)
c
         call sched(nproc,treeql,n,ldq,q,d,e,ifail)         
c
      t2 = second(gtime) - t1
      write(6,*) ' time for sesupd ',t2
      if( ifail .gt. 1 ) write(6,*)' deflate from sesupd',ifail
      tnorm = 0.0d0
      res = 0.0d0
      do 530 j = 1,n
         e(1) =   dsave(1)*q(1,j) +
     $          esave(2)*q(2,j) - d(j)*q(1,j)
      do 400 i = 2,n-1
         e(i) = esave(i)*q(i-1,j) + dsave(i)*q(i,j) +
     $          esave(i+1)*q(i+1,j) - d(j)*q(i,j)
  400 continue
         e(n) = esave(n)*q(n-1,j) + dsave(n)*q(n,j)
     $          - d(j)*q(n,j)
      t = enorm(n,e)
      res = max(res,t)
         do 520 i = 1,n
            s = 0.0
            do 510 k = 1,n
               s = s + q(k,i)*q(k,j)
  510       continue
            if (i .eq. j) s = s - 1.0
            t =  abs(s)
         tnorm = max(tnorm,t)
  520    continue
  530 continue
      write(6,*)' the residual for the updated values and vectors',res
      write(6,*)' the upd norm of q*q sup t  is',tnorm
      write (6,*) ' spectrum [ ',d(1) ,',',d(n),']'
 9998 continue
 9999 continue
      stop
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'make1'
then
	echo shar: will not over-write existing file "'make1'"
else
cat << \SHAR_EOF > 'make1'
FILES  = ftsubs.o putq.o second.o cputm.o
FILES2 = ftsubs.graph.o putq.o second.o cputm.o

sched :	$(FILES)
	rm -f sched.a; ar q sched.a $(FILES); ranlib sched.a

graph :	$(FILES2)
	rm -f graph.a; ar q graph.a $(FILES2); ranlib graph.a

.f.o : ;	fortran -g -c -mp -auto $*.f
.c.o : ;	cc -g -c $*.c


SHAR_EOF
fi # end of overwriting check
if test -f 'makefile'
then
	echo shar: will not over-write existing file "'makefile'"
else
cat << \SHAR_EOF > 'makefile'
FILES = stuffspawn.o
xtest :	$(FILES)
	fortran -g -mp $(FILES) sched.a -o xtest

gtest :	$(FILES)
	fortran -g -mp $(FILES) graph.a -o gtest


FILES2 = d_and_c.o
xdandc : $(FILES2)
	fortran -g -mp $(FILES2) sched.a -o xdandc

gdandc : $(FILES2)
	fortran -g -mp $(FILES2) graph.a -o gdandc


FILES4 = maindp.o newevdp0.o stateig.o statses.o 
xeig :	$(FILES4)
	fortran -g -mp $(FILES4) sched.a -o xeig

geig :	$(FILES4)
	fortran -g -mp $(FILES4) graph.a -o geig


FILES5 = maindp.o newevdp0.o stateig.o stseswait.o 
xwait :	$(FILES5)
	fortran -g -mp $(FILES5) sched.a -o xwait

gwait :	$(FILES5)
	fortran -g -mp $(FILES5) graph.a -o gwait


FILES6 = example.o
xexample : $(FILES6)
	fortran -g -mp $(FILES6) sched.a -o xexample

gexample : $(FILES6)
	fortran -g -mp $(FILES6) graph.a -o gexample

sched: 
	make -f make1 sched
graph: 
	make -f make1 graph

clean:
	/bin/rm -f *.o *.fpp *.dbg *~

.f.o : ;	fortran -g -c -mp -auto $*.f
SHAR_EOF
fi # end of overwriting check
if test -f 'newevdp0.f'
then
	echo shar: will not over-write existing file "'newevdp0.f'"
else
cat << \SHAR_EOF > 'newevdp0.f'
      subroutine evupd(n,i,d,z,delta,rho,dlam,ifail)
      implicit real*8 (a-h,o-z)
CVD$G NOCONCUR
c***********************************************************************
c     this subroutine computes the updated eigenvalues of a
c     rank one modification to a symmetric matrix.  it is assumed
c     that the eigenvalues are in the array d, and that
c
c        d(i) .ne. d(j)  for  i .ne. j
c
c     it is also assumed that the eigenvalues are in increasing
c     order and that the value of rho is positive.  this is
c     arranged by the calling subroutine sesupd, and is no loss
c     in generality.
c     it is also assumed that the values in the array z are
c     the squares of the components of the updatingvector.
c
c
c     the method consists of approximating the rational functions
c
c        rho = sum(z(j)/((d(j)-d(i))/rho - lamda): j = i+1,n)
c
c        phi =sum(z(j)/(d(j)-d(i))/rho - lamda): j = 1,i)
c
c     by simple interpolating rational functions.  this avoids
c     the need for safeguarding by bisection since
c     the convergence is monotone, and quadratic from any starting
c     point that is greater than zero but less than the solution
c
c
c     input variables...
c        n      the length of all arrays
c
c        i      the i - th  eigenvalue is computed
c
c        d      the original eigenvalues.  it is assumed that they are
c               in order, d(i) .lt. d(j)  for i .lt. j.
c
c        z      this array of length n contains the squares of
c               of the components of the updating vector
c
c        delta  this array of length n contains (d(j) - lamda(i))/rho
c               in its j - th component.  these values will be used
c               to update the eigenvectors in sesupd.
c
c        rho    this is the scalar in the the symmetric updating
c               formula.
c
c     dlam      this scalar will contain the value of the i - th
c               updated eigenvalue on return
c
c        ifail  this integer variable indicates failure of the
c               updating process with value 1,  and success
c               with value 0.
c
c
c
c**************************************************************************
      integer i,n,im1,ip1,ip2,niter
      dimension d(n),z(n), delta(n)
      double precision  d,z,rho,delta,zero,one,two,tsave,
     1                 del,phi,dphi,psi,dpsi,lambda,oldlam,
     2                 t,temp,a,b,d1,w,eps,eps1,eta,dlam,dmax
         double precision epslon,enorm
         zero = 0.0d0
         one  = 1.0d0
         two  = 2.0d0
c
c     eps is machine precision
c     eta is the relative accuracy requirement on the roots.
c
c         these values should be adjusted for the particular machine.
c
         eps = epslon(one)
         eta = eps*8.
c
         im1 = i - 1
         ip1 = i + 1
         ip2 = i + 2
         niter = 1
         lambda = zero
         oldlam = zero
         del = d(i)
         do 100 j = 1,n
            delta(j) = (d(j) - del)/rho
  100    continue
         dmax = max(abs(d(1)), abs(d(n)))
         dmax = abs(rho) + dmax
c
c     calculate initial guess
c
         if (i .lt. n)
     *   then
            del = d(ip1)
         else
            del = d(n) + rho
         endif
         a = zero
         do 200 j = 1,im1
            a = a + rho*z(j)/(d(j) - del)
  200    continue
         b = zero
         do 220 j = ip2,n
            b = b + rho*z(j)/(d(j) - del)
  220    continue
         a = a + (b + one)
         if (ip1 .gt. n)
     *   then
            t = z(i)/abs(a)
         else
            t = a*delta(ip1)
            b = t + z(i) + z(ip1)
            if (b .ge. zero)
     *      then
               t = two*z(i)*delta(ip1)/
     *             (b + sqrt(abs(b*b - 4.0d0*t*z(i))))
            else
               t = (b - sqrt(b*b - 4.0d0*t*z(i)))/(two*a)
            endif
            t = t/two
         endif
c
c     test to see that the initial guess is not too close to endpoint
c
      if (ip1 .le. n .and. t .ge. .9*delta(ip1))
     *                              t = .9*delta(ip1)
c
c     update the values of the array delta
c
  250 continue
         tsave = abs(t)
         do 300 j = 1,n
             delta(j) = delta(j) - t
  300    continue
         lambda = lambda + t
         dlam = d(i) + rho*lambda
c
c     evaluate psi and the derivative dpsi
c
         dpsi = zero
         psi  = zero
         do 400 j = 1,i
            t = z(j)/delta(j)
            psi = psi + t
            dpsi = dpsi + t/delta(j)
  400    continue
c
c     evaluate phi and the derivative dphi
c
         dphi = zero
         phi  = zero
         if (i .eq. n) go to 600
         do 500 j = ip1,n
            t = z(j)/delta(j)
            phi = phi + t
            dphi = dphi + t/delta(j)
  500    continue
c
c     test for convergence
c        return if the test is satisfied
c
  600    continue
         w = one + phi + psi
         eps1 = eps*dmax*sqrt(dphi + dpsi)
         if ((abs(w) .le. eps1) .and. (abs(lambda - oldlam) .le.
     1                eta*oldlam) ) then
            return
         endif
c
c     return with ifail = -1 if convergence has not ocurred
c     within 45 iterations.
c
         if (niter .lt. 45) go to 650
            ifail = -1
            return
  650    continue
         niter = niter + 1
         oldlam = lambda
 
c
c     calculate the new step
c
         if (i .ne. n) go to 700
c
c           otherwise
c           calculate the step for the special case  i = n
c
            t = (w*psi)/dpsi
            go to 250
  700    continue
         del = delta(ip1)
         temp = psi/dpsi
         d1 = one + phi - del*dphi
         a = (del*(one + phi) + psi*temp)/d1 + temp
         b = (two*temp*del*w)/d1
         t = sqrt(abs(a*a - two*b))
         t = b/(a + t)
         if (t .lt. -eps) t = -tsave/two
      go to 250
c
c     last card of evupd
c
      end
      subroutine deflat(j,jlam,jdel,indx,n,ldq,q,d,z,tol1)
      implicit real*8 (a-h,o-z)
CVD$G NOCONCUR
c
c
c     rotate z components to 0
c     assumes norm(z) .eq. 1
c
      integer indx(*),jdel,j,jlam,n,ldq
      double precision q(ldq,*),d(*),z(*),gamma,sigma,tau,t1,t2,
     *                 tol,tol1
c
           tol = tol1
           gamma = z(jlam)     
           sigma = z(j)        
           tau = pythag(gamma,sigma)
           jdel = -1
           if (tau .gt. 0.0d0) 
     *     then
              delta = d(j) - d(jlam)
              gamma = gamma/tau
              sigma = sigma/tau
              if (abs(delta*gamma*sigma) .le. tol) 
     *        then
                 if (abs(gamma) .le. abs(sigma))
     *           then
                   z(j) = tau
                    z(jlam) = 0.0d0
                    jjlam = indx(jlam)
                    jj = indx(j)
                    do 100 i = 1,n
                       t1 = q(i,jjlam)
                       t2 = q(i,jj)
                       q(i,jjlam) = -t1*sigma + t2*gamma
                       q(i,jj) = t1*gamma + t2*sigma
  100               continue
                    temp = d(jlam)*sigma**2 + d(j)*gamma**2
                    d(j) = d(jlam)*gamma**2 + d(j)*sigma**2
                    d(jlam) = temp
                    jdel = jlam
                    jlam = j
                 else
                    z(jlam) = tau    
                    z(j) = 0.0d0
                    jjlam = indx(jlam)
                    jj = indx(j)
                    do 200 i = 1,n
                       t1 = q(i,jjlam)
                       t2 = q(i,jj)
                       q(i,jjlam) =  t1*gamma + t2*sigma
                       q(i,jj) = -t1*sigma + t2*gamma
  200               continue
                    temp = d(jlam)*gamma**2 + d(j)*sigma**2
                    d(j) = d(jlam)*sigma**2 + d(j)*gamma**2
                    d(jlam) = temp
                    jdel = j
                 endif
              endif
            endif
c
      return
c
c     last card of deflat
      end
      function rand(foo)
CVD$G NOCONCUR
      data init/1365/
      init = mod(3125*init,65536)
      rand = (init - 32768.)/(2.**15)
      return
      end
      double precision function enorm(n,x)
CVD$G NOCONCUR
      integer n
      double precision x(n),temp,tmax
      tmax = 0.0d0
      do 100 j = 1,n
         temp = max(tmax,abs(x(j)))
         if (temp .gt. tmax) tmax = temp
  100 continue
      if (tmax .eq. 0.0d0) then       
         enorm = tmax
      else
         temp = 0.0d0
         do 200 j = 1,n
            temp = temp + (x(j)/tmax)**2
  200    continue
         enorm = tmax*sqrt(temp)
      endif
      return
 
      end
      subroutine tql2(nm,n,d,e,z,ierr)
      implicit real*8 (a-h,o-z)
CVD$G NOCONCUR
c
      integer i,j,k,l,m,n,ii,l1,l2,nm,mml,ierr
      double precision d(n),e(n),z(nm,n)
      double precision c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2,pythag
c
c     this subroutine is a translation of the algol procedure tql2,
c     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
c     wilkinson.
c     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
c
c     this subroutine finds the eigenvalues and eigenvectors
c     of a symmetric tridiagonal matrix by the ql method.
c     the eigenvectors of a full symmetric matrix can also
c     be found if  tred2  has been used to reduce this
c     full matrix to tridiagonal form.
c
c     on input
c
c        nm must be set to the row dimension of two-dimensional
c          array parameters as declared in the calling program
c          dimension statement.
c
c        n is the order of the matrix.
c
c        d contains the diagonal elements of the input matrix.
c
c        e contains the subdiagonal elements of the input matrix
c          in its last n-1 positions.  e(1) is arbitrary.
c
c        z contains the transformation matrix produced in the
c          reduction by  tred2, if performed.  if the eigenvectors
c          of the tridiagonal matrix are desired, z must contain
c          the identity matrix.
c
c      on output
c
c        d contains the eigenvalues in ascending order.  if an
c          error exit is made, the eigenvalues are correct but
c          unordered for indices 1,2,...,ierr-1.
c
c        e has been destroyed.
c
c        z contains orthonormal eigenvectors of the symmetric
c          tridiagonal (or full) matrix.  if an error exit is made,
c          z contains the eigenvectors associated with the stored
c          eigenvalues.
c
c        ierr is set to
c          zero       for normal return,
c          j          if the j-th eigenvalue has not been
c                     determined after 30 iterations.
c
c     calls pythag for  sqrt(a*a + b*b) .
c
c     questions and comments should be directed to burton s. garbow,
c     mathematics and computer science div, argonne national laboratory
c
c     this version dated august 1983.
c
c     ------------------------------------------------------------------
c
      ierr = 0
      if (n .eq. 1) go to 1001
c
      do 100 i = 2, n
  100 e(i-1) = e(i)
c
      f = 0.0d0
      tst1 = 0.0d0
      e(n) = 0.0d0
c
      do 240 l = 1, n
         j = 0
         h = abs(d(l)) + abs(e(l))
         if (tst1 .lt. h) tst1 = h
c     .......... look for small sub-diagonal element ..........
         do 110 m = l, n
            tst2 = tst1 + abs(e(m))
            if (tst2 .eq. tst1) go to 120
c     .......... e(n) is always zero, so there is no exit
c                through the bottom of the loop ..........
  110    continue
c
  120    if (m .eq. l) go to 220
  130    if (j .eq. 30) go to 1000
         j = j + 1
c     .......... form shift ..........
         l1 = l + 1
         l2 = l1 + 1
         g = d(l)
         p = (d(l1) - g) / (2.0d0 * e(l))
         r = pythag(p,1.0d0)
         d(l) = e(l) / (p + sign(r,p))
         d(l1) = e(l) * (p + sign(r,p))
         dl1 = d(l1)
         h = g - d(l)
         if (l2 .gt. n) go to 145
c
         do 140 i = l2, n
  140    d(i) = d(i) - h
c
  145    f = f + h
c     .......... ql transformation ..........
         p = d(m)
         c = 1.0d0
         c2 = c
         el1 = e(l1)
         s = 0.0d0
         mml = m - l
c     .......... for i=m-1 step -1 until l do -- ..........
         do 200 ii = 1, mml
            c3 = c2
            c2 = c
            s2 = s
            i = m - ii
            g = c * e(i)
            h = c * p
            r = pythag(p,e(i))
            e(i+1) = s * r
            s = e(i) / r
            c = p / r
            p = c * d(i) - s * g
            d(i+1) = h + s * (c * g + s * d(i))
c     .......... form vector ..........
            do 180 k = 1, n
               h = z(k,i+1)
               z(k,i+1) = s * z(k,i) + c * h
               z(k,i) = c * z(k,i) - s * h
  180       continue
c
  200    continue
c
         p = -s * s2 * c3 * el1 * e(l) / dl1
         e(l) = s * p
         d(l) = c * p
         tst2 = tst1 + abs(e(l))
         if (tst2 .gt. tst1) go to 130
  220    d(l) = d(l) + f
  240 continue
c     .......... order eigenvalues and eigenvectors ..........
      do 300 ii = 2, n
         i = ii - 1
         k = i
         p = d(i)
c
         do 260 j = ii, n
            if (d(j) .ge. p) go to 260
            k = j
            p = d(j)
  260    continue
c
         if (k .eq. i) go to 300
         d(k) = d(i)
         d(i) = p
c
         do 280 j = 1, n
            p = z(j,i)
            z(j,i) = z(j,k)
            z(j,k) = p
  280    continue
c
  300 continue
c
      go to 1001
c     .......... set error -- no convergence to an
c                eigenvalue after 30 iterations ..........
 1000 ierr = l
 1001 return
      end
      double precision function pythag(a,b)
      double precision a,b
c
c     finds sqrt(a**2+b**2) without overflow or destructive underflow
c
      double precision p,r,s,t,u
      p = max(abs(a),abs(b))
      if (p .eq. 0.0d0) go to 20
      r = (min(abs(a),abs(b))/p)**2
   10 continue
         t = 4.0d0 + r
         if (t .eq. 4.0d0) go to 20
         s = r/t
         u = 1.0d0 + 2.0d0*s
         p = u*p
         r = (s/u)**2 * r
      go to 10
   20 pythag = p
      return
      end
      double precision function epslon (x)
      double precision x
c
c     estimate unit roundoff in quantities of size x.
c
      double precision a,b,c,eps
c
c     this program should function properly on all systems
c     satisfying the following two assumptions,
c        1.  the base used in representing floating point
c            numbers is not a power of three.
c        2.  the quantity  a  in statement 10 is represented to
c            the accuracy used in floating point variables
c            that are stored in memory.
c     the statement number 10 and the go to 10 are intended to
c     force optimizing compilers to generate code satisfying
c     assumption 2.
c     under these assumptions, it should be true that,
c            a  is not exactly equal to four-thirds,
c            b  has a zero for its last bit or digit,
c            c  is not exactly equal to one,
c            eps  measures the separation of 1.0 from
c                 the next larger floating point number.
c     the developers of eispack would appreciate being informed
c     about any systems where these assumptions do not hold.
c
c     this version dated 4/6/83.
c
      a = 4.0d0/3.0d0
   10 b = a - 1.0d0
      c = b + b + b
      eps = abs(c-1.0d0)
      if (eps .eq. 0.0d0) go to 10
      epslon = eps*abs(x)
      epslon = 1.0d-16
      return
      end
      subroutine evdrv(k,kstart,kstop,ldq,n,q,d,rho,z,dlamda,q2,
     *                     w,indxp,indx,ksnrho,ifail)
      implicit real*8 (a-h,o-z)
CVD$G NOCONCUR
c
c
      integer k,kstart,kstop,ldq,n,ifail
      integer indx(n),indxp(n),ksnrho
      logical rhoge0
      double precision 
     *     q(ldq,n),d(n),z(n),q2(ldq,n),dlamda(n),
     *     w(n),rho
c
c     local variables
c
      double precision x(300),delta(300)
      double precision t,s,dlam 
      double precision enorm
      rhoge0 = .true.
      if (ksnrho .lt. 0) rhoge0 = .false.
c
         do 2700 j = kstart,kstop
            i = j
            if (.not. rhoge0) i = k - i + 1
            call evupd(k,i,dlamda,w,delta,rho,dlam,ifail)
c
c        if the zero finder failed the computation is terminated
c
            if (ifail .eq. 1) return
c
            if (rhoge0)
     *      then
              jp = indxp(j)
            else
              dlam = -dlam
              jp = indxp(k-j+1)
            endif
c
            d(jp) = dlam
c
c        compute the updated eigenvectors
c
               do 2200 i1 = 1,n
                  x(i1) = 0.0d0 
 2200          continue
            do 2500 jj = 1,k
               j3 = jj
               if (.not. rhoge0) j3 = k - jj + 1
                  jp = indxp(j3)
                  jjpp = indx(jp)
                  s = z(jp)/delta(j3)
c                 qq(jj,j) = s 
                  do 2300 i1 = 1,n
                     x(i1) = x(i1) + q(i1,jjpp)*s
 2300             continue
 2500       continue
            t =  enorm(n,x)
            j3 = j
            if (.not. rhoge0) j3 = k - j + 1
            jp = indxp(j3)
            do 2600 i1 = 1,n
               q2(i1,jp) = x(i1)/t
 2600       continue
 2700    continue
c
      return
c
c     last card of evdrv
c
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'oldtest'
then
	echo shar: will not over-write existing file "'oldtest'"
else
cat << \SHAR_EOF > 'oldtest'
	make -f make1 sched
	fortran -g -c -mp -auto ftsubs.f
	cc -g -c putq.c
	fortran -g -c -mp -auto second.f
	cc -g -c cputm.c
	rm -f sched.a; ar q sched.a ftsubs.o putq.o second.o cputm.o; ranlib sched.a
ar: creating sched.a
	make -f make1 graph
	fortran -g -c -mp -auto ftsubs.graph.f
	rm -f graph.a; ar q graph.a ftsubs.graph.o putq.o second.o cputm.o; ranlib graph.a
ar: creating graph.a
	fortran -g -c -mp -auto stuffspawn.f
	fortran -g -mp stuffspawn.o sched.a -o xtest
  input nprocs 
  time =    4.170000     nprocs =            6
   1.000000000000000    
  -2.000000000000000    
  -6.000000000000000    
  -8.000000000000000    
  -10.00000000000000    
  -12.00000000000000    
  -14.00000000000000    
  -16.00000000000000    
   2.000000000000000    
  -3.000000000000000    
  -8.000000000000000    
  -10.00000000000000    
  -12.00000000000000    
  -14.00000000000000    
  -16.00000000000000    
   3.000000000000000    
  -4.000000000000000    
  -10.00000000000000    
  -12.00000000000000    
  -14.00000000000000    
  -16.00000000000000    
   4.000000000000000    
  -5.000000000000000    
  -12.00000000000000    
  -14.00000000000000    
  -16.00000000000000    
   5.000000000000000    
  -6.000000000000000    
  -14.00000000000000    
  -16.00000000000000    
   6.000000000000000    
  -7.000000000000000    
  -16.00000000000000    
   7.000000000000000    
  -8.000000000000000    
   8.000000000000000    
FORTRAN STOP
	fortran -g -mp stuffspawn.o graph.a -o gtest
  input nprocs 
  time =    4.950000     nprocs =            6
   1.000000000000000    
  -2.000000000000000    
  -6.000000000000000    
  -8.000000000000000    
  -10.00000000000000    
  -12.00000000000000    
  -14.00000000000000    
  -16.00000000000000    
   2.000000000000000    
  -3.000000000000000    
  -8.000000000000000    
  -10.00000000000000    
  -12.00000000000000    
  -14.00000000000000    
  -16.00000000000000    
   3.000000000000000    
  -4.000000000000000    
  -10.00000000000000    
  -12.00000000000000    
  -14.00000000000000    
  -16.00000000000000    
   4.000000000000000    
  -5.000000000000000    
  -12.00000000000000    
  -14.00000000000000    
  -16.00000000000000    
   5.000000000000000    
  -6.000000000000000    
  -14.00000000000000    
  -16.00000000000000    
   6.000000000000000    
  -7.000000000000000    
  -16.00000000000000    
   7.000000000000000    
  -8.000000000000000    
   8.000000000000000    
FORTRAN STOP
	fortran -g -c -mp -auto d_and_c.f
	fortran -g -mp d_and_c.o sched.a -o xdandc
  input nprocs nlevls 
           1
           2
           3
           4
           5
           6
           7
           8
           9
          10
          11
          12
          13
          14
          15
FORTRAN STOP
	fortran -g -mp d_and_c.o graph.a -o gdandc
  input nprocs nlevls 
           1
           2
           3
           4
           5
           6
           7
           8
           9
          10
          11
          12
          13
          14
          15
FORTRAN STOP
	fortran -g -c -mp -auto maindp.f
	fortran -g -c -mp -auto newevdp0.f
	fortran -g -c -mp -auto stateig.f
	fortran -g -c -mp -auto statses.f
	fortran -g -mp maindp.o newevdp0.o stateig.o statses.o  sched.a -o xeig
  input nproc ... the number of processors 
 nprocs =            6
 ksect =           10 kgran =            5
 ==============================
  n =          100
 ++++++++++++++++++++++++++++++
  twos on diagonal
  time for sesupd    3.430000264197588    
  the residual for the updated values and vectors  3.6584310490749194E-15
  the upd norm of q*q sup t  is  1.5525097483634620E-15
  spectrum [   3.8688057328107606E-03,   3.999032564583977    ]
FORTRAN STOP
	fortran -g -mp maindp.o newevdp0.o stateig.o statses.o  graph.a -o geig
  input nproc ... the number of processors 
 nprocs =            6
 ksect =           10 kgran =            5
 ==============================
  n =          100
 ++++++++++++++++++++++++++++++
  twos on diagonal
  time for sesupd    3.830000072717667    
  the residual for the updated values and vectors  3.6584310490749194E-15
  the upd norm of q*q sup t  is  1.5525097483634620E-15
  spectrum [   3.8688057328107606E-03,   3.999032564583977    ]
FORTRAN STOP
	fortran -g -c -mp -auto stseswait.f
	fortran -g -mp maindp.o newevdp0.o stateig.o stseswait.o  sched.a -o xwait
  input nproc ... the number of processors 
 nprocs =            6
 ksect =           10 kgran =            5
 ==============================
  n =          100
 ++++++++++++++++++++++++++++++
  twos on diagonal
  time for sesupd    3.470000218600035    
  the residual for the updated values and vectors  3.6584310490749194E-15
  the upd norm of q*q sup t  is  1.5525097483634620E-15
  spectrum [   3.8688057328107606E-03,   3.999032564583977    ]
FORTRAN STOP
	fortran -g -mp maindp.o newevdp0.o stateig.o stseswait.o  graph.a -o gwait
  input nproc ... the number of processors 
 nprocs =            6
 ksect =           10 kgran =            5
 ==============================
  n =          100
 ++++++++++++++++++++++++++++++
  twos on diagonal
  time for sesupd    3.890000294893980    
  the residual for the updated values and vectors  3.6584310490749194E-15
  the upd norm of q*q sup t  is  1.5525097483634620E-15
  spectrum [   3.8688057328107606E-03,   3.999032564583977    ]
FORTRAN STOP
	fortran -g -c -mp -auto example.f
	fortran -g -mp example.o sched.a -o xexample
  sigma =    5050.000    
FORTRAN STOP
	fortran -g -mp example.o graph.a -o gexample
  sigma =    5050.000    
FORTRAN STOP
SHAR_EOF
fi # end of overwriting check
if test -f 'putq.c'
then
	echo shar: will not over-write existing file "'putq.c'"
else
cat << \SHAR_EOF > 'putq.c'
#include <stdio.h>

#define MAXPARMS	20

struct parms
{
/*  int   static_link; */
  int   (*subname)();
  long  *parms[MAXPARMS];
};
shared struct parms indx[1000]; 
sched_(nprocs,parms)
int *nprocs;
struct parms	parms;
/*  
    this procedure obtains nprocs physical processors devoted
    to the the execution of the parallel program indicated through parms
    which is a structure whose first entry is a subroutine name and whose
    remaining entries are parameters appearing in the calling sequence
    of that subroutine.
*/
{
        int libopn_();
	bcopy(&parms, &indx[0], sizeof(struct parms));
/*      
           the subroutine name and prameter list have been copied and 
           placed in a special slot on the parmq      
          
           then libopn is invoked to initialize pointers, grab physical
           processors and begin the computation
*/
        libopn_(nprocs);
	return(0);
}
putq_(jobtag,parms)
int *jobtag;
struct parms	parms;
/*  
    this procedure puts the descriptor of a schedulable process <jobtag>
    onto the problem queue.  this process will be scheduled for execution
    when its data dependencies have been satisfied (indicated by icango==0).
    the argument parms is a structure whose first entry is a subroutine name 
    and whose remaining entries are parameters appearing in the calling sequence
    of that subroutine.
*/
{
        register int j;
        int place_();
        j = *jobtag;
	bcopy(&parms, &indx[j], sizeof(struct parms));
/*
        first the parms block is copied into the slot pointed to by 
        by jobtag and then this descriptor is placed on the problem 
        queue
*/ 
        place_(jobtag);
	return(0);
}
spawn_(parent,jobtag,parms)
int *parent,*jobtag;
struct parms	parms;
/*  
    this procedure puts the descriptor of a schedulable process <jobtag>
    onto the problem queue.  this process will be scheduled for execution
    when its data dependencies have been satisfied (indicated by icango==0).
    the argument parms is a structure whose first entry is a subroutine name 
    and whose remaining entries are parameters appearing in the calling sequence
    of that subroutine.
    
    the action of this procedure differs from putq in that the user does not
    assign jobtags or data dependencies.  a parent may spawn any number of 
    children but these child processes only report to the parent.
*/
{
        register int j,i;
        int place_(),clone_();
        j = *jobtag;
        i = *parent;
	bcopy(&parms, &indx[j], sizeof(struct parms));
/*
        first the parms block is copied into the slot pointed to by 
        by jobtag and then this descriptor is placed on the problem 
        queue
*/ 
        if (indx[j].subname == clone_) indx[j].subname = indx[i].subname;
/*
        here we ask if this is a recursive spawning.  if so the name
        clone has been called instead of subname so we replace the name
        clone by subname.
*/
        place_(jobtag);
	return(0);
}
clone_()
{
/*
        this is a dummy routine to satisfy unresolved external
*/
        return(0);
}
work_(id,jobtag)
int *id,*jobtag;
{
        int start2_(),gtprb_();
        register int j,myjob;

        j = *id;
        if (j == 1) 
/*
        the worker whose id is 1 will execute the subroutine passed to 
        sched.  this subroutine executes the static data dependency graph.
        this graph must have at least one node.
*/
        {
	   indx[0].subname(indx[0].parms[0], indx[0].parms[1],
			       indx[0].parms[2], indx[0].parms[3],
			       indx[0].parms[4], indx[0].parms[5],
			       indx[0].parms[6], indx[0].parms[7],
			       indx[0].parms[8], indx[0].parms[9],
			       indx[0].parms[10], indx[0].parms[11],
			       indx[0].parms[12], indx[0].parms[13],
			       indx[0].parms[14], indx[0].parms[15],
			       indx[0].parms[16], indx[0].parms[17],
			       indx[0].parms[18], indx[0].parms[19]);
              start2_();
        }
         myjob = gtprb_(id,jobtag);
         while (myjob != 0) 
         {
           j = *jobtag;
           if (myjob <= -1 )
           {
/*
              reenter... simple spawning was done
              all kids completed and no reentry
              is required.  this indicates
              jobtag is all done and checkin can proceed.
*/
              chekin_(jobtag);
              myjob = gtprb_(id,jobtag);
           }
           else
           {
/*
               call subname(<parms>)..........
*/
	       indx[j].subname(indx[j].parms[0], indx[j].parms[1],
			       indx[j].parms[2], indx[j].parms[3],
			       indx[j].parms[4], indx[j].parms[5],
			       indx[j].parms[6], indx[j].parms[7],
			       indx[j].parms[8], indx[j].parms[9],
			       indx[j].parms[10], indx[j].parms[11],
			       indx[j].parms[12], indx[j].parms[13],
			       indx[j].parms[14], indx[j].parms[15],
			       indx[j].parms[16], indx[j].parms[17],
			       indx[j].parms[18], indx[j].parms[19]);
             chekin_(jobtag); 
             myjob = gtprb_(id,jobtag);
           }
          }
          return(0);
}
SHAR_EOF
fi # end of overwriting check
if test -f 'second.f'
then
	echo shar: will not over-write existing file "'second.f'"
else
cat << \SHAR_EOF > 'second.f'
$SYSTEM
      real function second(t)
c
c     this routine will gather the user time for a process.
c     it has resolution of a millisecond
c     and uses the unix system call getrusage.
c     see the unix manual for details
c     returns time in seconds.
c
      integer cputm
c
      itime = cputm()
      second = float(itime)/1000.0
c
c     this statement is here to bump the time by a bit
c     incase no the interval was too small.
c
c      second = second + second*1.0e-6
c
      return
      end

SHAR_EOF
fi # end of overwriting check
if test -f 'stateig.f'
then
	echo shar: will not over-write existing file "'stateig.f'"
else
cat << \SHAR_EOF > 'stateig.f'
      subroutine treeql(n,ldq,q,d,e,ifail)         
      implicit real*8 (a-h,o-z)
CVD$G NOCONCUR
c                                                                 *
      integer n,ldq
      double precision d(*),e(*),q(ldq,*)
c
      double precision z(500),dlamda(500),x(500),q2(500,500),w(500),
     *                 delta(500),rho(300)
      integer indx(500),indxp(500),nn(3,300),kwork(50,300),ifail
      common /wspace/z,dlamda,x,q2,w,delta,rho,indx,indxp,nn,kwork
      common/prfpms/nsmall,kgran
      external sesupd,tql2
c
c     this subroutine splits a problem in two parts
c
      nn(1,1) = n
      nn(2,1) = 1 
      nn(3,1) = 1
      ifail = 0
      nsize = n
      nlevl = -1
 1111 continue
      nlevl = nlevl + 1
      if (nsize .gt. 2*nsmall) then
c
c        the problem is large enough to split one more level
c        define the splits here in the 100 loop and then
c        place a call to sesupd on the queue to glue the results
c        together
c
         nsize = nsize/2
         do 100 id = 2**nlevl,2**(nlevl+1) - 1
c
            idl = 2*id
            idr = idl + 1
            n1 = nn(1,id)/2
c
            nn(1,idl) = n1 
            nn(2,idl) = idl
            nn(3,idl) = nn(3,id)
c
            nn(1,idr) =  nn(1,id) - n1   
            nn(2,idr) = idr      
            nn(3,idr) = nn(3,id) + n1
c
            isplt1 = nn(3,idr) 
            isplt  = isplt1 - 1
            n1p1 = n1 + 1
            rho(id) = e(isplt1)
            alpha = d(isplt)
            alphap = d(isplt1)
            sigma = 1.0d0
c
            if (sign(1.0d0,alpha)*sign(1.0d0,alphap) .ge. 0.0d0)
     *      then
               if (alpha .lt. 0.0d0 .or. alphap .lt. 0.0d0)
     *         then
                  if(rho(id) .ne. 0.0d0)
     *            then
                    rho(id) = - rho(id)                      
                    sigma = -sigma
                  endif
               endif
            endif
            kwork(50,id) = sigma
            d(isplt) = d(isplt) - rho(id)         
            d(isplt1) = d(isplt1) - rho(id)         
c
c     the problem is large enough to split into two parts
c
            isplt = nn(3,id) 
c
            jobtag = id
            icango = 2
            nchks = 1
            if (jobtag .eq. 1) nchks = 0
            list = id/2
            call dep(jobtag,icango,nchks,list)    
            call putq(jobtag,sesupd,nn(2,id),ldq,nn(1,id),
     *               q(isplt,isplt),d(isplt),rho(id),z(isplt),
     *               x(isplt),dlamda(isplt),q2(isplt,isplt),
     *               delta(isplt),w(isplt),indxp(isplt),
     *               indx(isplt),kwork(1,id),ifail)
c
  100    continue
         go to 1111
      else
c
c       we have reached the lowest level
c       granularity now warrants call to tql2
c
         do 200 id = 2**nlevl,2**(nlevl+1) - 1
c
            isplt = nn(3,id) 
c
            jobtag = id
            icango = 0
            ncheks = 1
            list = id/2
c
            call dep(jobtag,icango,ncheks,list)
            call putq(jobtag,tql2,ldq,nn(1,id),d(isplt),e(isplt),
     *                  q(isplt,isplt),ifail) 

c

  200    continue
      endif
c
      return
c
c     last card of split
c
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'statses.f'
then
	echo shar: will not over-write existing file "'statses.f'"
else
cat << \SHAR_EOF > 'statses.f'
      subroutine sesupd(myid,ldq,n,q,d,rho,z,x,dlamda,q2,delta,w,indxp,
     *                                          indx,kwork,ifail)
      implicit real*8 (a-h,o-z)
CVD$G NOCONCUR
c***********************************************************************
c
c
c     this subroutine will compute the updated eigensystem of a
c     of a symmetric matrix after modification by a rank one
c     symmetric matrix.
c
c     a = qdq' + rho*z*z'
c
c     it is assumed that the eigenvectors of the original matrix
c     are stored in q, and the eigenvalues are in d.
c     the algorithm consists of three stages...
c
c
c        the first stage constists of deflating the size of
c        the problem when there multiple eigenvalues or if there
c        zero of the vector q'z.  for each such ocurrence the dimension
c        is reduced by one.
c
c        the second stage consists of calculating the updated
c        eigenvalues of the reduced problem.  this requires a call
c        to the zero finding routine evupd.
c
c        the final stage consists of computing the updated eigenvectors
c        directly using the updated eigenvalue.
c
c
c     the algorithm requires o(n**2) operations to update the
c     eigenvectors, but n**3 + o(n**2) to update the eigenvectors.
c
c
c     input variables...
c
c        n      the dimension of the problem.  q is n x n
c
c        q      an  n x n  matrix that contains the eigenvectors of
c               the original matrix on input and the updated
c               eigenvectors on output.
c
c        d      a vector of length n. the original eigenvalues are
c               contained in d on input.  the updated eigenvectors
c               are contained on output.
c
c        rho    a scalar
c
c        z      a vector of length n.  on input this vector
c               containes the updating vector.  the contents of z
c               are destroyed during the updating process.
c
c        x      a working array of length n.
c
c        dlamda a working array of length n
c
c        q2     a working array of dimension  n x n
c
c        delta  a working array of length n
c
c        w      a working array of length n.
c
c        indx   an integer array of length n.
c
c        ifail  this integer variable indicates failure of the
c               updating process with value 1,  and success
c               with value 0.
c
c     called subroutines...
c
c        evupd  a subroutine for calculating the updated eigenvalues.
c
c
c***********************************************************************
      integer ldq,n,ifail
      integer indx(n),indxp(n),kwork(*)
      integer ksize(20),kstart(20),kstop(20),kbins,msize,krem,prcsys
      integer kgran,ksect
      common/prfpms/ksect,kgran
      logical rhoge0
      double precision 
     *       q(ldq,n),d(n),z(n),x(n),dlamda(n),q2(ldq,n),delta(n)
     *       ,w(n),rho
      double precision eps,zero,one,two,s,t,evsprd,dmax,dlam
      double precision epslon,enorm
      external evdrv
c
c
c     eps is machine precision
c
      ifail = 0
      zero = 0.0d0
      one = 1.0d0
      two = 2.0d0
      eps = epslon(one)
      rhoge0 = (rho .ge. zero)
      go to (1111,2222),ientry(myid,2)
 1111 continue
c
c     compute q(transpose)*z
c
         n1 = n/2
         n1p1 = n1 + 1
         do 100 j = 1,n1
            z(j) = q(n1,j)
  100    continue
         sigma = kwork(50)
         if (sigma .gt. 0.0) then
            do 200 j = n1p1,n
               z(j) = q(n1p1,j)
  200       continue
         else
            do 201 j = n1p1,n
               z(j) = -q(n1p1,j)
  201       continue
         endif
c
c     normalize z so that norm(z) = 1
c
      t = enorm(n,z)
      do 300 j = 1,n
         z(j) = z(j)/t
         indx(j) = j
  300 continue
      rho = rho*t*t
c
c     calculate the allowable deflation tolerence  
c
        tol = (1.0d2)*eps*max(abs(d(1)),abs(d(n)))
         if (abs(rho) .le. tol) return
c
c
c     order the eigenvalues
c
      nm1 = n - 1
      do 600 j = 1,nm1
         t = d(1)
         s = z(1)
         inx = indx(1)
         k = n - j + 1
         do 500 i = 2,k
            im1 = i - 1
            if (d(i) .lt. t) go to 400
               t = d(i)
               s = z(i)
               inx = indx(i)
               go to 500
  400       continue
            d(im1) = d(i)
            d(i) = t
            z(im1) = z(i)
            z(i) = s
            indx(im1) = indx(i)
            indx(i) = inx
  500    continue
  600 continue
c
c     if there multiple eigenvalues then the problem deflates.
c     here the number of equal eigenvalues are found.
c     then an elementary reflector is computed to rotate the
c     corresponding eigensubspace so that certain components of
c     z are zero in this new basis.
c
      if (rhoge0) 
     *then
         k = 0
         k2 = n + 1
         do 670 j = 1,n
            if(rho*abs(z(j)) .le. tol) 
     *      then
c
c              deflate due to small z component
c
               k2 = k2 - 1
               indxp(k2) = j
               if (j .eq. n) go to 1111
            else
               jlam = j
               go to 700
            endif
  670    continue
  700    continue
            j = j + 1
            if (j .gt. n) go to 800
            if(rho*abs(z(j)) .le. tol) 
c           if(rho*abs(z(j)) .le. 1.0d-10) 
     *      then
c
c              deflate due to small z component
c
               k2 = k2 - 1
               indxp(k2) = j
            else
               call deflat(j,jlam,jdef,indx,n,ldq,q,d,z,tol)
               if (jdef .le. 0) 
     *         then
                  k = k + 1
                  w(k) = z(jlam)**2    
                  dlamda(k) = d(jlam)
                  indxp(k) = jlam
                  jlam = j
               else
                  k2 = k2 - 1
                  indxp(k2) = jdef
               endif
            endif 
            go to 700
  800    continue
c
c        record the last eigenvalue
c
            k = k + 1
            w(k) = z(jlam)**2    
            dlamda(k) = d(jlam)
            indxp(k) = jlam
      else
         jlam = n
         k = n + 1
         k2 = 0     
         do 671 j = n,1,-1
            if(abs(rho*z(j)) .le. tol) 
c           if(abs(rho*z(j)) .le. 1.0d-10) 
     *      then
c
c              deflate due to small z component
c
               k2 = k2 + 1
               indxp(n-k2+1) = j
               if (j .eq. 1) go to 1111
            else
               jlam = j
               go to 701
            endif
  671    continue
  701    continue
            j = j - 1
            if (j .lt. 1) go to 801
            if(abs(rho*z(j)) .le. tol) 
     *      then
c
c              deflate due to small z component
c
               k2 = k2 + 1
               indxp(n-k2+1) = j
            else
               call deflat(j,jlam,jdef,indx,n,ldq,q,d,z,tol)
               if (jdef .le. 0) 
     *         then
                  k = k - 1
                  w(n-k+1) = z(jlam)**2    
                  dlamda(n-k+1) = -d(jlam)
                  indxp(n-k+1) = jlam
                  jlam = j
               else
                  k2 = k2 + 1
                  indxp(n-k2+1) = jdef
               endif
            endif 
            go to 701
  801    continue
c
c        record the last eigenvalue
c
         k = k - 1
         w(n-k+1) = z(jlam)**2    
         dlamda(n-k+1) = -d(jlam)
         indxp(n-k+1) = jlam
         k = n - k + 1
         rho = -rho
      endif
c
c     compute the updated eigenvalues of the deflated problem
c
c
      kbins = k/kgran + 1
      kbins = min(kbins,20)
      krem = mod(k,kbins)
      msize = (k - krem)/kbins
      do 900 j = 1,kbins
         ksize(j) = msize
  900 continue
      do 910 j = 1,krem
         ksize(j) = ksize(j) + 1
  910 continue
      kstart(1) = 1
      kstop(1) = ksize(1)
      kwork(49) = k
      kwork(48) = -1
      if (rhoge0) kwork(48) = 1
      if ( kbins .gt. 1 ) then
c
         do 920 j = 2,kbins    
               kstart(j) = kstop(j-1) + 1
               kstop(j) = kstop(j-1) + ksize(j)
               kwork(j) = kstart(j)
               kwork(kbins+j) = kstop(j)
               call nxtag(myid,jtemp)
               call spawn(myid,jtemp,evdrv,kwork(49),
     *                    kwork(j),kwork(kbins+j),
     *                    ldq,n,q,d,rho,z,dlamda,q2,
     *                    w,indxp,indx,kwork(48),ifail)
  920    continue
         call evdrv(k,kstart(1),kstop(1),ldq,n,q,d,rho,z,
     *              dlamda,q2,w,indxp,indx,kwork(48),ifail)
c
c        kbins-1  processes were spawned 
c        place a join point here and resume when all spawned 
c        processes have checked in
c
c
      else
c
c        there was not sufficient granularity to warrant spawning
c        of additional processes to compute roots
c
            kstop(1) = k
         call evdrv(k,kstart(1),kstop(1),ldq,n,q,d,rho,z,
     *              dlamda,q2,w,indxp,indx,kwork(48),ifail)
      endif
      return
c
 2222 continue
c
c     store the updated eigenvectors back into q
c
      k = kwork(49) 
      do 1100 j = k+1,n
         jp = indxp(j)
         jjpp = indx(jp)
         do 1000 i = 1,n
            q2(i,jp) = q(i,jjpp)
 1000    continue
 1100 continue
      do 1300 j = 1,n
         do 1200 i = 1,n
            q(i,j) = q2(i,j)
 1200    continue
 1300 continue
c
      return
c
c     last card of sesupd
c
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'stseswait.f'
then
	echo shar: will not over-write existing file "'stseswait.f'"
else
cat << \SHAR_EOF > 'stseswait.f'
      subroutine sesupd(myid,ldq,n,q,d,rho,z,x,dlamda,q2,delta,w,indxp,
     *                                          indx,kwork,ifail)
      implicit real*8 (a-h,o-z)
CVD$G NOCONCUR
c***********************************************************************
c
c
c     this subroutine will compute the updated eigensystem of a
c     of a symmetric matrix after modification by a rank one
c     symmetric matrix.
c
c     a = qdq' + rho*z*z'
c
c     it is assumed that the eigenvectors of the original matrix
c     are stored in q, and the eigenvalues are in d.
c     the algorithm consists of three stages...
c
c
c        the first stage constists of deflating the size of
c        the problem when there multiple eigenvalues or if there
c        zero of the vector q'z.  for each such ocurrence the dimension
c        is reduced by one.
c
c        the second stage consists of calculating the updated
c        eigenvalues of the reduced problem.  this requires a call
c        to the zero finding routine evupd.
c
c        the final stage consists of computing the updated eigenvectors
c        directly using the updated eigenvalue.
c
c
c     the algorithm requires o(n**2) operations to update the
c     eigenvectors, but n**3 + o(n**2) to update the eigenvectors.
c
c
c     input variables...
c
c        n      the dimension of the problem.  q is n x n
c
c        q      an  n x n  matrix that contains the eigenvectors of
c               the original matrix on input and the updated
c               eigenvectors on output.
c
c        d      a vector of length n. the original eigenvalues are
c               contained in d on input.  the updated eigenvectors
c               are contained on output.
c
c        rho    a scalar
c
c        z      a vector of length n.  on input this vector
c               containes the updating vector.  the contents of z
c               are destroyed during the updating process.
c
c        x      a working array of length n.
c
c        dlamda a working array of length n
c
c        q2     a working array of dimension  n x n
c
c        delta  a working array of length n
c
c        w      a working array of length n.
c
c        indx   an integer array of length n.
c
c        ifail  this integer variable indicates failure of the
c               updating process with value 1,  and success
c               with value 0.
c
c     called subroutines...
c
c        evupd  a subroutine for calculating the updated eigenvalues.
c
c
c***********************************************************************
      integer ldq,n,ifail
      integer indx(n),indxp(n),kwork(*)
      integer ksize(20),kstart(20),kstop(20),kbins,msize,krem,prcsys
      integer kgran,ksect
      common/prfpms/ksect,kgran
      common/wsync/iwrite
      logical rhoge0,wait
      double precision eps,zero,one,two,s,t,evsprd,dmax,dlam
      double precision epslon,enorm
      double precision 
     *       q(ldq,n),d(n),z(n),x(n),dlamda(n),q2(ldq,n),delta(n)
     *       ,w(n),rho
      external evdrv
c
c
c     eps is machine precision
c
c     write(6,*) ' in ses ldq n rho ifail  ',ldq,n,rho,ifail
c     write(6,*) 'fses myid rho q d z x ',myid,rho,q(1,1),d(1),z(1),x(1)
      ifail = 0
      zero = 0.0d0
      one = 1.0d0
      two = 2.0d0
      eps = epslon(one)
      rhoge0 = (rho .ge. zero)
      go to (1111,2222),ientry(myid,2)
 1111 continue
c
c     compute q(transpose)*z
c
         n1 = n/2
         n1p1 = n1 + 1
         do 100 j = 1,n1
            z(j) = q(n1,j)
  100    continue
         sigma = kwork(50)
c        write(6,*) ' sigma at 2 ',sigma
         if (sigma .gt. 0.0) then
            do 200 j = n1p1,n
               z(j) = q(n1p1,j)
  200       continue
         else
            do 201 j = n1p1,n
               z(j) = -q(n1p1,j)
  201       continue
         endif
c
c     normalize z so that norm(z) = 1
c
      t = enorm(n,z)
c     call lockon(iwrite)
c           write(6,*) ' n t z1 zn ',n,t,z(1),z(n)
c     call lockoff(iwrite)
      do 300 j = 1,n
         z(j) = z(j)/t
         indx(j) = j
  300 continue
      rho = rho*t*t
c
c     calculate the allowable deflation tolerence  
c
        tol = (1.0d2)*eps*max(abs(d(1)),abs(d(n)))
         if (abs(rho) .le. tol) return
c
c
c     order the eigenvalues
c
c      write(6,*) ' before order evs ',d(1),d(n)
      nm1 = n - 1
      do 600 j = 1,nm1
         t = d(1)
         s = z(1)
         inx = indx(1)
         k = n - j + 1
         do 500 i = 2,k
            im1 = i - 1
            if (d(i) .lt. t) go to 400
               t = d(i)
               s = z(i)
               inx = indx(i)
               go to 500
  400       continue
            d(im1) = d(i)
            d(i) = t
            z(im1) = z(i)
            z(i) = s
            indx(im1) = indx(i)
            indx(i) = inx
  500    continue
  600 continue
c      write(6,*) ' after order evs ',d(1),d(n)
c
c     if there multiple eigenvalues then the problem deflates.
c     here the number of equal eigenvalues are found.
c     then an elementary reflector is computed to rotate the
c     corresponding eigensubspace so that certain components of
c     z are zero in this new basis.
c
      if (rhoge0) 
     *then
         k = 0
         k2 = n + 1
         do 670 j = 1,n
            if(rho*abs(z(j)) .le. tol) 
     *      then
c
c              deflate due to small z component
c
               k2 = k2 - 1
               indxp(k2) = j
               if (j .eq. n) go to 1111
            else
               jlam = j
               go to 700
            endif
  670    continue
  700    continue
            j = j + 1
            if (j .gt. n) go to 800
            if(rho*abs(z(j)) .le. tol) 
c           if(rho*abs(z(j)) .le. 1.0d-10) 
     *      then
c
c              deflate due to small z component
c
               k2 = k2 - 1
               indxp(k2) = j
            else
               call deflat(j,jlam,jdef,indx,n,ldq,q,d,z,tol)
               if (jdef .le. 0) 
     *         then
                  k = k + 1
                  w(k) = z(jlam)**2    
                  dlamda(k) = d(jlam)
                  indxp(k) = jlam
                  jlam = j
               else
                  k2 = k2 - 1
                  indxp(k2) = jdef
               endif
            endif 
            go to 700
  800    continue
c
c        record the last eigenvalue
c
            k = k + 1
            w(k) = z(jlam)**2    
            dlamda(k) = d(jlam)
            indxp(k) = jlam
      else
         jlam = n
         k = n + 1
         k2 = 0     
         do 671 j = n,1,-1
            if(abs(rho*z(j)) .le. tol) 
c           if(abs(rho*z(j)) .le. 1.0d-10) 
     *      then
c
c              deflate due to small z component
c
               k2 = k2 + 1
               indxp(n-k2+1) = j
               if (j .eq. 1) go to 1111
            else
               jlam = j
               go to 701
            endif
  671    continue
  701    continue
            j = j - 1
            if (j .lt. 1) go to 801
            if(abs(rho*z(j)) .le. tol) 
     *      then
c
c              deflate due to small z component
c
               k2 = k2 + 1
               indxp(n-k2+1) = j
            else
               call deflat(j,jlam,jdef,indx,n,ldq,q,d,z,tol)
               if (jdef .le. 0) 
     *         then
                  k = k - 1
                  w(n-k+1) = z(jlam)**2    
                  dlamda(n-k+1) = -d(jlam)
                  indxp(n-k+1) = jlam
                  jlam = j
               else
                  k2 = k2 + 1
                  indxp(n-k2+1) = jdef
               endif
            endif 
            go to 701
  801    continue
c
c        record the last eigenvalue
c
         k = k - 1
         w(n-k+1) = z(jlam)**2    
         dlamda(n-k+1) = -d(jlam)
         indxp(n-k+1) = jlam
         k = n - k + 1
         rho = -rho
      endif
c     write(6,*) ' begin a problem ------------------------------- '
c     write(6,*) ' '
c     write(6,*) k,rho
c     do 6666 jj = 1,k
c        write(6,*) z(indxp(jj)),w(jj),dlamda(jj)
c6666 continue
c     write(6,*) ' begin a problem ------------------------------- '
c
c
c     compute the updated eigenvalues of the deflated problem
c
c
      kbins = k/kgran + 1
      kbins = min(kbins,20)
      krem = mod(k,kbins)
      msize = (k - krem)/kbins
      do 900 j = 1,kbins
         ksize(j) = msize
  900 continue
      do 910 j = 1,krem
         ksize(j) = ksize(j) + 1
  910 continue
c        call lockon(iwrite)
c        write(6,*) ' from ses w ',w(1),w(2),w(k)
c        call lockoff(iwrite)
      kstart(1) = 1
      kstop(1) = ksize(1)
      kwork(49) = k
      kwork(48) = -1
      if (rhoge0) kwork(48) = 1
c     write(6,*) ' rhoge0 kwork  ',rhoge0,kwork(48)
c
      if ( kbins .gt. 1 ) then
c     write(6,*) ' about to spawn',kbins,' evdrvs from ses ',myid
c
         do 920 j = 2,kbins    
               kstart(j) = kstop(j-1) + 1
               kstop(j) = kstop(j-1) + ksize(j)
               kwork(j) = kstart(j)
               kwork(kbins+j) = kstop(j)
c              write(6,*) ' bef spawn evdrv ',kwork(49),
c    *                  kwork(j),kwork(kbins+j),kwork(48)
c        call lockon(iwrite)
c        write(6,*) ' about to spawn evdrv ',k,kstart(j),kstop(j)
c        call lockoff(iwrite)
               call nxtag(myid,jtemp)
               call spawn(myid,jtemp,evdrv,kwork(49),
     *                    kwork(j),kwork(kbins+j),
     *                    ldq,n,q,d,rho,z,dlamda,q2,
     *                    w,indxp,indx,kwork(48),ifail)
c               write(6,*) ' tag assigned ',jtemp
  920    continue
c        call lockon(iwrite)
c        write(6,*) ' about to call evdrv ',k,kstart(1),kstop(1)
c        call lockoff(iwrite)
         call evdrv(k,kstart(1),kstop(1),ldq,n,q,d,rho,z,
     *              dlamda,q2,w,indxp,indx,kwork(48),ifail)
c
c        kbins-1  processes were spawned 
c        place a join point here and resume when all spawned 
c        processes have checked in
c
c
      else
c
c        there was not sufficient granularity to warrant spawning
c        of additional processes to compute roots
c
            kstop(1) = k
c     write(6,*) 'bef rho q d z x ',rho,q(1,1),d(1),z(1),x(1)
c        call lockon(iwrite)
c        write(6,*) ' about to call evdrv ',k,kstart(1),kstop(1)
c        call lockoff(iwrite)
         call evdrv(k,kstart(1),kstop(1),ldq,n,q,d,rho,z,
     *              dlamda,q2,w,indxp,indx,kwork(48),ifail)
      endif
      if( wait(myid,2) ) return
c
 2222 continue
c
c     store the updated eigenvectors back into q
c
      k = kwork(49) 
      do 1100 j = k+1,n
         jp = indxp(j)
         jjpp = indx(jp)
         do 1000 i = 1,n
            q2(i,jp) = q(i,jjpp)
 1000    continue
 1100 continue
      do 1300 j = 1,n
         do 1200 i = 1,n
            q(i,j) = q2(i,j)
 1200    continue
 1300 continue
c     call lockon(iwrite)
c        write(6,*) ' after ret in ses ',k,n
c           write(6,*) ' k n z1 zn ',k,n,q(1,1),q(n,1)
c     call lockoff(iwrite)
c     ifail = n-k
c
      return
c
c     last card of sesupd
c
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'stuffspawn.f'
then
	echo shar: will not over-write existing file "'stuffspawn.f'"
else
cat << \SHAR_EOF > 'stuffspawn.f'
      double precision a,b
      common /prbdef/ a(1000),b(100),itmp(10),jtmp(10),myid(10)
      EXTERNAL PARALG
      nblks = 8
      do 10 j = 1,10
         itmp(j) = j
         jtmp(j) = j
   10 continue
      write(6,*) ' input nprocs '
      read (5,*) nprocs
      t1 = second(foo)
c
      CALL SCHED(nprocs,paralg,nblks,a,b,itmp,jtmp,myid)
c
      t2 = second(foo)
      write(6,*) ' time = ',t2-t1,' nprocs = ',nprocs
      nn = 36
      do 100 j = 1,nn
      write(6,*) a(j)
  100 continue 
      stop
      end
c
      subroutine paralg(n,a,b,itmp,jtmp,myid)
      integer itmp(*),jtmp(*),myid(*)
      double precision a(*),b(*)
      integer mychkn(1)
      EXTERNAL STUFF1
c
c     this is the driver for filling a packed triangular matrix with 
c     i on the i-th diagonal and -(i+j) in the (i,j) off diagonal position
c
      icount = 1
      do 200 j = 1,n-1
c
c
c        the j-th diagonal waits for the diagonal above to complete
c
c        the j-th diagonal completion will allow 
c        the (j+1)-st diagonal to start
c
c
            jobtag = j
            icango =  1
            if (jobtag .eq. 1) icango = 0
            nchks = 1
            mychkn(1) = j+1
c
c        we just set up data dependencies and are ready to put
c        this process on the queue
c
            CALL DEP(jobtag,icango,nchks,mychkn)
            myid(jobtag) = jobtag          
            CALL PUTQ(jobtag,stuff1,myid(jobtag),n,
     &                        a(icount),jtmp(j),itmp(1))
c
c        when the data dependencies for process jobtag are satisfied
c        the following call will be made
c
c          call  stuff1(myid....,itmp(1))
c
         icount = icount + (n-j+1)
  200 continue
c
         icango = 1
         nchks = 0
         mychkn(1) = n+1
c
         CALL DEP(n,icango,nchks,mychkn)
         myid(n) = n
         CALL PUTQ(n,stuff1,myid(n),n,a(icount),jtmp(n),itmp(1))
c
      return
      end
c
      subroutine stuff1(myid,n,a,j,itmp)
      double precision a(*)
      integer myid,n,j,itmp(*)
      logical wait
      EXTERNAL STUFF2,STUFF3
c
c         write(6,*) ' enter stuff1 ',myid,j
c
c         write(6,*) ' enter stuff1 ',ientry(myid),myid
         go to (1111,2222,3333),ientry(myid,3)
 1111    continue
         ii = 2
         do 100 i = 2,n  
c
            CALL NXTAG(myid,jdummy)
c               write(6,*) ' about to spawn id jd ',myid,jdummy
            CALL SPAWN(myid,jdummy,stuff2,a(ii + 1),itmp(i),itmp(j))
c
c           this spawns a process that will execute a call to stuff2
c           and report completion  to process MYID
c
            ii = ii + 1
  100    continue
         call stuff2(a(2),itmp(1),itmp(j))
         if (wait(myid,2)) return
c
c        return to help out and then return here (at label 2222) 
c        on the next reentry
c
 2222    continue
         ii = 2
         do 200 i = 2,n  
c
            CALL NXTAG(myid,jdummy)
c               write(6,*) ' about to spawn 3 id jd ',myid,jdummy
            CALL SPAWN(myid,jdummy,stuff3,a(ii + 1),itmp(i),itmp(j))
c
c           this spawns a process that will execute a call to stuff2
c           and report completion  to process MYID
c
            ii = ii + 1
  200    continue
         call stuff2(a(2),itmp(1),itmp(j))
         if (wait(myid,3)) return
 3333    continue
         a(1) = j
c
      return
      end
      subroutine stuff2(a,i,j)
      double precision a(*),one
      one = 1.0d0
      a(1) = 0.0d0
c         write(6,*) ' enter stuff2 ',i,j
         do 100 kk = 1,10000
            a(1) = a(1) + one
  100    continue
         a(1) = -(i+j)
      return
      end
      subroutine stuff3(a,i,j)
      double precision a(*),one,save
      one = 1.0d0
      save = a(1) 
c         write(6,*) ' enter stuff3 ',i,j
         do 100 kk = 1,10000
            a(1) = a(1) + one
  100    continue
         a(1) = -(i+j) + save
      return
      end
c
SHAR_EOF
fi # end of overwriting check
if test -f 'testrun'
then
	echo shar: will not over-write existing file "'testrun'"
else
cat << \SHAR_EOF > 'testrun'
make sched; make graph
make xtest; xtest < data
make gtest; gtest < data
make xdandc; xdandc < data
make gdandc; gdandc < data
make xeig; xeig < data
make geig; geig < data
make xwait; xwait < data
make gwait; gwait < data
make xexample; xexample < data
make gexample; gexample < data
SHAR_EOF
chmod +x 'testrun'
fi # end of overwriting check
#	End of shell archive
exit 0