


		    PNOTICE_gtss.alm                01/13/87  1102.7r w 01/13/87  1102.7        2853


	dec	1			"version 1 structure
	dec	1			"no. of pnotices
	dec	3			"no. of STIs
	dec	100			"lgth of all pnotices + no. of pnotices
          acc       "Copyright (c) 1987 by Massachusetts Institute of
Technology and Honeywell Information Systems, Inc."

	aci	"C1TSLM090000"
	aci	"C2TSLM090000"
	aci	"C3TSLM090000"
	end
   



		    attach.pl1                      10/25/77  1244.0r w 10/21/77  1534.7        3681


attach : proc (istat);

dcl  iox_$attach_name entry (char (*), ptr, char (*), ptr, fixed bin (35));
dcl  fname char (32) varying;
dcl  ptr ptr;
dcl  istat fixed bin (35);
dcl (sysprint, sysin) file;
dcl  null builtin;

	put skip list ("Data filename?");
	get list (fname);
	call iox_$attach_name ("file10", ptr, "vfile_ "||fname, null, istat);
	close file (sysprint), file (sysin);
	return;

     end attach;
   



		    detach.pl1                      10/25/77  1244.0r w 10/21/77  1534.7        2997


detach : proc;

dcl  iox_$find_iocb entry (char(*), ptr, fixed bin(35));
dcl  iox_$close entry (ptr,fixed bin(35));
dcl  iox_$detach_iocb entry (ptr, fixed bin(35));
dcl  ptr ptr;
dcl  code fixed bin(35);

call iox_$find_iocb ("file10", ptr, code);
call iox_$close (ptr,code);
call iox_$detach_iocb (ptr, code);
return;

end detach;
   



		    flat.pl1                        10/25/77  1244.0r w 10/21/77  1534.7        2754


flat:    proc(seed) returns (float bin);
    dcl  random_$uniform   entry (fixed bin (35), float bin (27));
    dcl  s                 fixed bin(35),
         seed              fixed bin (35),
         xret              float bin (27);

    call random_$uniform(seed,xret);
    return(xret);
    end flat;
  



		    smlrpattach.pl1                 10/25/77  1244.0r w 10/21/77  1534.7        3366


smlrpattach : proc (istat, fname);

dcl  iox_$attach_name entry (char(*), ptr, char(*), ptr, fixed bin(35));
dcl  fname char(32) varying;
dcl  ptr ptr;
dcl  istat fixed bin(35);
dcl (sysprint, sysin) file;
dcl  null builtin;

put skip list ("Data filename?");
get list (fname);
call iox_$attach_name ("file10", ptr, "vfile_ "||fname, null, istat);
return;

end smlrpattach;
  



		    tss_fortran.fortran             10/25/77  1244.0r w 10/21/77  1538.0      866133


c  ********normal probability function ******
 function anpf(x1,x2,xmu,sigma)
 1 temp=1.4142136*sigma
 u1=(x1-xmu)/temp
 u2=(x2-xmu)/temp
 if (u1) 300,200,100
 100 if (u2) 200,200,900
 300 if (u2) 900,200,200
 200 anpf=0.5*(errf(abs(u2))*sign(1.0,u2)-errf(abs(u1))*sign(1.0,u1))
 go to 5001
 900 anpf =0.5*(errf(-abs(u1))*sign(1.0,u1)-errf(-abs(u2))*sign(1.0,u2))
 5001 return
 end
c   beta
c *****incomplete beta function approximation*****
cSTART beta
       function beta(ix,bbb,xx,yy)
       isign=ix
       bp=bbb
       x=xx
       y=yy
       if(isign-1)100,200,200
 100   c1=gbeta(x,y,0.,bp)
       c2=gbeta(x,y,0.,1.)
       beta=c1/c2
       return
   200 if(x-1.)232,225,232
   225 if(y-1.)229,228,229
   232 if(y-1.)237,233,237
 237   bet=gbeta(x,y,0.,1.)
       beta=x/(x+y)
       do 250 i=1,5
       c1=gbeta(x,y,0.,beta)
       c2=ftion(x,y,beta)
       beta=beta+(bet*bp-c1)/c2
       if(beta-1.)10,40,40
    10 if(beta)20,20,250
    20 beta =0.
    30 go to 300
    40 beta=1.
       go to 300
   250 continue
       go to 300
   228 beta=bp
       go to 300
   229 beta=1.-(1.-bp)**(1./y)
       go to 300
   233 beta=bp**(1./x)
       go to 300
   300 return
       end
cSTART gbeta
       function gbeta(x,y,a,b)
       dimension u(8),h(8)
       data u/.04750625,.14080178,.22900839,.30893812,
&      .37770220,.43281560,.47228751,.49470047/
       data h/.09472531,.09130171,.08457826,.07479799,
&      .06231449,.04757926,.03112676,.01357623/
       sum=0.0
       do 20 i=1,8
       arg1=((b-a)*u(i)+(b+a)/2.)
       arg2=(a-b)*u(i)+(b+a)/2.
       x1=ftion(x,y,arg1)
       x2=ftion(x,y,arg2)
    20 sum=sum+h(i)*(x1+x2)
       gbeta=(b-a)*sum
       return
       end
cSTART ftion
       function ftion(x,y,z)
       xu=x-1.
       if(z-1.e-37)100,100,20
    20 if(xu*alog(z)-80.)30,200,200
    30 if(xu*alog(z)+80.)300,300,40
    40 ch=xu*alog(z)
       a=z**xu
       go to 50
   100 z=1.e-37
       go to 40
   200 a=1.e37
       go to 50
   300 a=0.
    50 yv=y-1.
       q=1.-z
       if(q-1.e-37)400,400,60
    60 if(yv*alog(q)-80.)70,500,500
    70 if(yv*alog(q)+80.)600,600,80
    80 ch=yv*alog(q)
       b=q**yv
    90 if(a)82,95,82
    82 if(b)85,95,85
    85 if((alog(abs(a))+alog(abs(b))).ge.(-80.))go to 95
       ftion=1.e-36
       go to 96
 95    ftion=a*b
    96 return
   400 q=1.e-37
       go to 80
   500 b=1.e37
       go to 90
   600 b=0.
       go to 90
       end
c   errf   mod dec 6,1972
cSTART errf
 function errf(w)
c ******program to evaluate the error function and its compliment
 double precision a(25),b(30),f
 data a/16443152242714.d-13,-9049760497548.d-13,
& 643570883797.d-13,196418177368.d-13,-1244215694.d-13,
& -9101941905.d-13,-1796219835.d-13,139836786.d-13,
& 164789417.d-13,39009267.d-13,-893145.d-13,-3747896.d-13,
& 1298818.d-13,136773.d-13,77107.d-13,46810.d-13,
& 11844.d-13,-5.d-13,-1384.d-13,-652.d-13,145.d-13,
& 10.d-13,24.d-13,11.d-13,2.d-13/
 m=24
 x=abs(w)
 if(x-.01)1,2,2
 1 xerr=2.0/(3.0*1.77245385)*x*(3.0-x**2)
 go to 6
 2 z=(x-1.0)/(x+1.0)
 do 3 i=1,30

 b(i)=0.0
 3 continue
 do 4 i=1,m
 m1=(m+1)-i
 b(m1)=2.0*z*b(m1+1)-b(m1+2)+a(m1+1)
 4 continue
 f=-b(2)+z*b(1)+.5*a(1)
 xerr=1.0-(1.0/1.77245385)*(exp(-(x**2)))*f
 if(x-.01)6,7,7
 6 cerr=1.0-xerr
 go to 5
 7 cerr=(1.0/1.77245385)*(exp(-(x**2)))*f
 5 if(w)9,8,8
 8 errf=xerr

 go to 13
 9 errf=cerr
 13 return
 end
c linefit
cSTART linefit
 subroutine linefit (in,last,x,y,s,yi)
    dimension x(last),y(last)
       n=last-in+1
       s1=0
       s2=0
       s3=0
       s4=0
       do 10 i=in,last
       s1=s1+x(i)
       s2=s2+y(i)
       s3=s3+x(i)*y(i)
       s4=s4+x(i)*x(i)
 10  continue
       s=(n*s3-s1*s2)/(n*s4-s1*s1)
       yi=(s2*s4-s1*s3)/(n*s4-s1*s1)
       return
    end
c   lineq
c  ***********solution of simultaneous linear equations ***********
c  *************  gaussian elimination  ***************************
cSTART lineq
 subroutine lineq(a,b,naarg,nbarg,idim)
 dimension a(idim,naarg),b(idim,nbarg)
 1 na=naarg
 nb=nbarg
 ns=0
 if(na)10,5001,60
 10 if(nb)20,5001,30
 20 nb=iabs(nb)
 na=iabs(na)
 ns=1
 go to 40
 30 na=iabs(na)
 ns=-1
 40 napl1=na+1
 60 do 291 j1=1,na
 temp=0.
 if(ns)120,122,125
 120 do 121 j2=j1,na
 if(abs(a(j2,j1)).lt.temp) go to 121
 temp=abs(a(j2,j1))
 a(j1,napl1)=j2
 121 continue
 go to 125
 122 do 124 j2=j1,na
 if(abs(a(j2,j1)).lt.temp) go to 124
 temp=abs(a(j2,j1))
 ibig=j2
 124 continue
 go to 126
 125 ibig=a(j1,napl1)
 126 if(ibig.eq.j1) go to 201
 if(ns)135,135,160
c  ********** rearrange rows to place largest absolute **********
c  ********** value in pivot position  **************************
 135 do 141 j2=j1,na
 temp=a(j1,j2)
 a(j1,j2)=a(ibig,j2)
 141 a(ibig,j2)=temp
 160 do 161 j2=1,nb
 temp=b(j1,j2)
 b(j1,j2)=b(ibig,j2)
 161 b(ibig,j2)=temp
 201 if(j1.eq.na) go to 301
 n1=j1+1
c  ******** compute new coefficients in remaining rows  *********
 do 281 j2=n1,na
 if(ns)240,240,250
 240 a(j2,j1) = a(j2,j1)/a(j1,j1)
 do 241 j3=n1,na
 241 a(j2,j3)=a(j2,j3)-a(j2,j1)*a(j1,j3)
 250 do 251 j3=1,nb
 251 b(j2,j3)=b(j2,j3)-a(j2,j1)*b(j1,j3)
 281 continue
 291 continue
c  *****************  obtain solutions  *************************
 301 do 391 j1=1,nb
 b(na,j1)=b(na,j1)/a(na,na)
 if(na.eq.1) go to 391
 n1=na
 321 do 341 j2=n1,na
 341 b(n1-1,j1)=b(n1-1,j1)-b(j2,j1)*a(n1-1,j2)
 b(n1-1,j1)=b(n1-1,j1)/a(n1-1,n1-1)
 n1=n1-1
 if(n1.ne.1) go to 321
 391 continue
 5001 return
 end
c   mtinv
c      *************matrix inversion*************************
cSTART mtinv
 subroutine mtinv(a,nrarg,ncarg,idim,label)
 dimension a(idim,ncarg),label(nrarg)
 1 nr=nrarg
 nc=ncarg
 do 21 j1=1,nr
 21 label(j1)=j1
 do 291 j1=1,nr
c      **************find remaining row containing largest***
c      **************absolute value in pivotal column********
 101 temp=0.0
 do 121 j2=j1,nr
 if(abs(a(j2,j1)).lt.temp) go to 121
 temp=abs(a(j2,j1))
 ibig=j2
 121 continue
 if(ibig.eq.j1)go to 201
c      **************rearrange rows to place largest absolute
c      **************value in pivot position*****************
 do 141 j2=1,nc
 temp=a(j1,j2)
 a(j1,j2)=a(ibig,j2)
 141 a(ibig,j2)=temp
 i=label(j1)
 label(j1)=label(ibig)
 label(ibig)=i
c     ::::compute coefficients in pivotal row::::
 201 temp=a(j1,j1)
 a(j1,j1)=1.0
 do 221 j2=1,nc
 221 a(j1,j2)=a(j1,j2)/temp
c      **************compute coefficients in other rows******
 do 281 j2=1,nr
 if(j2.eq.j1) go to 281
 temp=a(j2,j1)
 a(j2,j1)=0.0
 do 241 j3=1,nc
 241 a(j2,j3)=a(j2,j3)-temp*a(j1,j3)
 281 continue
 291 continue
c      **************interchange columns according to********
c      **************interchanges of rows of original matrix*
 301 n1=nr-1
 do 391 j1=1,n1
 do 321 j2=j1,nr
 if(label(j2).ne.j1) go to 321
 if(j2.eq.j1) go to 391
 go to 341
 321 continue
 341 do 361 j3=1,nr
 temp=a(j3,j1)
 a(j3,j1)=a(j3,j2)
 361 a(j3,j2)=temp
 label(j2)=label(j1)
 391 continue
 5001 return
 end
c   mtmpy -- rev. april 1971
c honeywell time sharing applications
c      **************matrix multiplication*******************
cSTART mtmpy
 subroutine mtmpy(ind,a,b,c,larg,marg,n)
 dimension a(25,999),b(25,999),c(25,999)
 1 l=iabs(larg)
 m=iabs(marg)
 i=ind+1
 go to (101,201,301,401),i
 101 do 121 j1=1,l
 do 121 j2=1,n
 c(j1,j2)=0.0
 do 121 j3=1,m
 if(larg)102,5001,110
 102 if(marg)103,5001,105
c      **************t(matrix) x t(matrix)   t(a(m,l))*t(b(n,
 103 temp=a(j3,j1)*b(j2,j3)
 go to 121
c      **************t(matrix) x matrix      t(a(m,l))*b(m,n)
 105 temp=a(j3,j1)*b(j3,j2)
 go to 121
 110 if(marg)111,5001,115
c      **************matrix x t(matrix)      a(l,m)*t(b(n,m))
 111 temp=a(j1,j3)*b(j2,j3)
 go to 121
c      **************matrix x matrix         a(l,m)*b(m,n)=c(
 115 temp=a(j1,j3)*b(j3,j2)
 121 c(j1,j2)=c(j1,j2)+temp
 go to 5001
c      **************diagonal x diagonal     a(l,1)*b(l,1)=c(
 201 do 221 j1=1,l
 221 c(j1,1)=a(j1,1)*b(j1,1)
 go to 5001
 301 do 321 j1=1,l
 do 321 j2=1,m
 if(marg)310,5001,315
c      **************diagonal x t(matrix)    a(l,1)*t(b(m,l))
 310 temp=a(j1,1)*b(j2,j1)
 go to 321
c      **************diagonal x matrix       a(l,1)*b(l,m)=c(
 315 temp=a(j1,1)*b(j1,j2)
 321 c(j1,j2)=temp
 go to 5001
 401 do 421 j1=1,l
 do 421 j2=1,m
 if(larg)410,5001,415
c      **************t(matrix) x diagonal    t(a(m,l))*b(m,1)
 410 temp=a(j2,j1)*b(j2,1)
 go to 421
c      **************matrix x diagonal       a(l,m)*b(m,1)=c(
 415 temp=a(j1,j2)*b(j2,1)
 421 c(j1,j2)=temp
 5001 return
 end
c   plot  --  july 1972
cSTART plot
      subroutine plot(x,y,ymax,ymin,n,ind,ntot)
      dimension ylabel(5),y(n)
save
      character m72*1(72),iblank*1,markx*1,marki*1,mark*1(9)
      data iblank,markx,marki/" ","x","|"/
      data mark/"*",".","+","(",")","#","$","-",","/
      cwidth = 70.
      if(ind)20,20,1
    1 ncount=0
      icount=0
      scale=(ymax-ymin)/cwidth
      lpo=-ymin/scale+1.5
      if(lpo.eq.1)lpo=2
      if(lpo)2,2,3
    2 lpo=2
    3 do 10 j=1,72
   10 m72(j)=" "
      ylabel(1)=ymin
      ylabel(5)=ymax
      ylabel(3)=(ylabel(1)+ylabel(5))/2.
      ylabel(2)=(ylabel(1)+ylabel(3))/2.
      ylabel(4)=(ylabel(3)+ylabel(5))/2.
      print 115,ylabel
      return
   20 icount=icount+1
      if(x)910,888,910
   888 isign=1
      print 110
      go to 907
  900 print 106
      go to 907
  907 return
  910 do 930 j=1,n
      if(y(j).le.ymin)go to 921
      z=(y(j)-ymin)/scale+1.5
      if(z-1.)921,921,920
  920 if(z-cwidth)928,928,925
  921 l=2
      go to 929
  925 l=72
      go to 929
  928 l=z
      if(l)921,921,929
  929 m72(lpo)=marki
      m72(l)=mark(j)
  930 continue
  931 do 934 j=1,72
      if(m72(j).eq.iblank)go to 934
  933 jj=j
  934 continue
      ncount=ncount+1
      if(ncount-1)900,940,936
  936 if((ncount/10)*10-ncount)938,940,938
  938 print 102,(m72(j),j=2,jj)
      go to 962
  940 m72(lpo)=markx
  950 print 103,x,(m72(j),j=13,jj)
  962 do 970 j=1,72
  970 m72(j)=iblank
      if(icount.eq.ntot.and.isign.ne.1)go to 975
      return
  975 isign=1
      print 110
      print 115,ylabel
  102 format(1x,71a1)
  103 format(1x,1pe11.4,60a1)
  106 format(" input error")
  110 format(2h y,2(16(1h-),1hy),2(17(1h-),1hy))
  115 format(1x,1pe11.4,1pe14.4,1p3e15.4)
      return
      end
c rndnrm---normal random number generator
cSTART rndnrm
 function rndnrm(s)
 2 r1=flat(s);w=flat(s)-.5;x=abs(w)-.2
 4 if(x) 6 ,6,8
 6 if(r1-.9955698)32,32,30
 8 if(x-.17) 9,9,14
 9 if(r1-.59035036) 32,32,10
 10 r1mtl=81.3*x*x*x-.9895+r1
  if(r1mtl)32,32,12
 12 if(r1mtl-.0235)30,30,2
 14 if(x-.24) 16,16,20
 16 if(r1-2.024+8.4*x)18,18,2
 18 if(10.13*x-2.31245+r1)32,32,30
 20 if(r1-.008) 30 ,30,2
 30 z=.250000004-w*w
 rndnrm=w*(1.77079939+.402454407/z)
 t=.119047619*(.5-z)/(z*z)+.523809524
 t=t*exp(-.5*rndnrm*rndnrm)
 if(r1-t)40,40,2
 32 rndnrm=w*(1.77079939+.402454407/(.250000004-w*w))
 40 return
 end
c tdist
cSTART tdist
 function tdist(ix,a,x)
 a1 = a/2.
 a2 = .5
 if(ix .eq. 1) go to 20
 if(x .lt. 1.e-20) go to 10
 y = a/(a+x*x)
 tdist = 1.-beta(0,y,a1,a2)

 5 return
 10 tdist = 0.
 go to 5
 20 y = 1.-x
 fx = beta(1,y,a1,a2)
 tdist = sqrt(a*((1./fx)-1.))
 go to 5
 end
 subroutine rkpb1(exit,temp,x,dx,y,f,n)

c **************runge=kutta integration    ********
c *******calculate derivatines    ***************
c *******store functions and derivatives **********
 dimension temp(999),y(n),f(n)
external exit
c *******initialization ***************************
 101 m=n+1
     temp(2*m+1)=x
     do 121 i=2,m
     ip2m=i+2*m
 121 temp(ip2m)=y(i-1)
     call exit
     temp(3*m+1)=dx
     do 141 i=2,m
     ip3m=i+3*m
 141 temp(ip3m)=f(i-1)
5001 return
     end
 subroutine rkpb2(deriv,temp,x,dx,y,f,n)
 dimension temp(999),y(n),f(n),rkcon(4)
external deriv
 data rkcon/1.,2.,2.,1./
c ************integrate to next point **********
  401 m=n+1
c **************restore previous values
      x=temp(2*m+1)
      do 421 i=2,m
      ip2m=i+2*m
      y(i-1)=temp(ip2m)
      ip3m=ip2m+m
      f(i-1)=temp(ip3m)
 421  temp(i)=0.0
      do 491 j=1,4
c *****************compute ki,s,delts  ***********
      do 441 i=2,m
      ipm =i+m
      temp(ipm)=f(i-1)*dx
 441  temp(i)=temp(i)+temp(ipm)*rkcon(j)
      if (j.eq.4) go to 501
      x=temp(2*m+1)+dx/rkcon(j+1)
      do 461 i=2,m
      ipm=i+m
      ip2m=ipm+m
  461 y(i-1)=temp(ip2m)+temp(ipm)/rkcon(j+1)
      call deriv
  491 continue
  501 x=temp(2*m+1)+dx
      do 521 i=2,m
      ip2m=i+2*m
  521 y(i-1)=temp(ip2m)+temp(i)/6.0
 5001 return
      end
c    gaspiia Modified Feb 1975
c     Start gasp
 subroutine gasp(nset,qset)
 dimension nset(999),qset(999)
 common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
 common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
 character fname*16
double precision stat
equivalence (istat,stat)
external com_err_ (descriptors), attach (descriptors),
& datan (descriptors), montr (descriptors), rmove (descriptors),
& events (descriptors), sumry (descriptors), otput (descriptors),
& error (descriptors), detach
 ncrdr=10
 nprnt=6
 not = 0
 call  attach (istat)
 if (istat) 5656,4545,5656
 5656 call com_err_(stat,"gasp","error in attaching")
 stop
 4545 continue
  1  call datan(nset,qset)
 jevnt = 101
 call montr (nset,qset)
 write (nprnt,403)
 403 format(1h ,3x,24h**intermediate results**//)
 10 call rmove(mfe(1),1,nset,qset)
 tnow = atrib(1)
 jevnt = jtrib(1)
 if(jevnt - 100)13,12,6
 13 i = jevnt
 call events (i,nset,qset)
 if (mstop) 40,8,20
 40 mstop = 0
 if (norpt) 14,22,42
 20 if(tnow-tfin)8,22,22
 22 call sumry(nset,qset)
 call otput (nset,qset)
 42 if(nruns-1)14,9,23
 23 nruns = nruns - 1
 nrun = nrun + 1
 go to 1
 14 call error(93,nset,qset)
 6 call montr(nset,qset)
 go to 10
 12 if(jmnit)14,30,31
 30 jmnit = 1
 go to 10
 31 jmnit = 0
 go to 10
 8 if(jmnit)14,10,32
 32 jtrib(1)=jevnt
 jevnt = 100
 call montr(nset,qset)
 go to 10
   9 end file 10
 call detach
  return
 end
function amax(arg1,arg2)
if (arg1 - arg2) 2,1,1
1 amax = arg1
return
2 amax = arg2
return
end
function amin (arg1,arg2)
if (arg1-arg2) 1,1,2
1 amin = arg1
return
2 amin = arg2
return
end
cSTART colct
 subroutine colct (x,n,nset,qset)
 dimension nset(999),qset(999)
 common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
 common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
external error(descriptors)
 if (n) 2,2,1
 2 call error(90,nset,qset)
 1 if (n- nclct) 3,3,2
 3 suma(n,1) = suma(n,1)+x
 suma(n,2) = suma(n,2)+x*x
 suma(n,3) = suma(n,3)+1.0
 suma(n,4) = amin (suma(n,4),x)
 suma(n,5) = amax (suma(n,5),x)
 return
 end
cSTART datan
 subroutine datan(nset,qset)
 dimension nset(999),qset(999)
 common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
 common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
external error(descriptors), exit, drand(descriptors),set(descriptors),filem(descriptors)
 if (not)23,1,2
 2 nt=nep
 go to (1,5,6,41,42,8,43,299,15,20),nt
 23 call error(95,nset,qset)
 1 not = 1
 nrun = 1
  read(ncrdr,100) name
c	name =	12 alpha chars
 100  format(6a2)
 read (ncrdr,101) nproj,mon,nday,nyr,nruns
c	nproj =	4 chars
c	mon =	2 numeric chars
c	nday =	2 numeric chars
c	nyr =	4 numeric chars
c	nruns =	4 numeric chars
 101 format (v)
 if(nruns) 30,30,5
 30 call exit
 5 read (ncrdr,803) nprms,nhist,nclct,nstat,id,im,noq,mxc,imm
 803 format (v)
 if (nhist) 41,41,6
 6 read (ncrdr,803) (ncels(i),i=1,nhist)
 41 read (ncrdr,803) (krank(i),i=1,noq)
 42 read (ncrdr,803) (inn(i),i=1,noq)
 if (nprms) 23,43,8
 8 do 9 i = 1,nprms
 read (ncrdr,803) (param(i,j),j=1,4)
 9 continue
 43 read (ncrdr, 803) mstop,jclr,norpt,nep,tbeg,tfin,jseed
 if (jseed) 26,26,27
 27 iseed=jseed
 call drand(iseed,rnum)
 tnow = tbeg
 do 142 j=1,noq
 142 qtime(j)=tnow
 26 jmnit = 0
 299 do 300 js = 1,id
 read (ncrdr,803)jq,(jtrib(jk),jk=1,im)
 if(jq) 44,15,320
 44 init=1
 call set(1,nset,qset)
 go to 300
 320 read(ncrdr,803) (atrib(jk),jk=1,imm)
 call filem(jq,nset,qset)
 300 continue
 15 if( jclr )20,20,10
 10 if(nclct)23,110,116
 116 do 18 i = 1,nclct
 do 17 j =1,3
 17 suma(i,j) = 0.
 suma(i,4) = 1.0e20
 18 suma(i,5)= -1.0e20
 110 if (nstat)23,111,117
 117 do 360 i = 1,nstat
 do 370 j = 1,3
 370 ssuma(i,j) = 0.
 ssuma(i,4) = 1.0e20
 360 ssuma(i,5) = -1.0e20
 111 if(nhist)23,20,118
 118 do 380 k = 1,nhist
 do 380 l = 1,mxc
 380 jcels(k,l) = 0
 20 write (nprnt,102) nproj,name,mon,nday,nyr,nrun
 102 format (1h ,1x,22hsimulation project no.,i4,2x,2hby,2x,
& 6a2//,2x,4hdate,i3,1h/,i3,1h/,i5,12x,10hrun number,i5//)
 if(nprms ) 60,60,62
 62 do 64 i=1,nprms
 64 write (nprnt,107) i,(param(i,j),j=1,4)
 107 format(2x,14h parameter no.,i5,4f12.4)
 60 return
 end
subroutine drand(iseed,rnum)
external randu(descriptors)
  call randu (iseed,rnum)
return
end
cSTART error
 subroutine error(j,nset,qset)
 dimension nset(999),qset(999)
 common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
 common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
 write (nprnt,100) j,tnow
 100 format(//2x,16herror exit, type,i3,7h error.//19h file status
& a time,f10.4/)
 write (nprnt,200)
 200 format(18x,4hnset/)
 do 210 i=1,id
 il=(i-1)*mxx+1
 iv=il+mxx-1
 210 write(nprnt,90) i,(nset(ij),ij=il,iv)
 90 format(1x,i5,5x,8i8,/,13x,4i8)
 write (nprnt,202)
 202 format(//18x,4hqset/)
 do 215 i=1,id
 il=(i-1)*imm+1
 iv=il+imm-1
 215 write(nprnt,95) i,(qset(ij),ij=il,iv)
 95 format(1x,i5,4x,5e13.6,/,11x,5e13.6)
 if(nclct) 7,7,8
 8 write (nprnt,98)
 98 format(/11h array suma,/)
 do 110 i=1,nclct
 110 write (nprnt,80) i,(suma(i,k),k=1,5)
 80 format(i10,5f10.4)
 write (nprnt,99)
 7 if(nstat)9,9,10
 10 write (nprnt,97)
 97 format(/12h array ssuma/)
 do 111 i=1,nstat
 111 write (nprnt,80) i,(ssuma(i,k),k=1,5)
 write (nprnt,99)
 9 if(nhist) 11,11,12
 12 write (nprnt,96)
 96 format (/12h array jcels/)
 do 112 i=1,nhist
 ncl=ncels (i)+2
 112 write (nprnt,26) i,(jcels(i,k),k=1,ncl)
 26 format(2x,i3,4x,15i4,/11x,15i4)
 11 nfool = 0
 if (nfool) 3,4,3
 3 return
 4 stop
 99 format(1h )
 end
cSTART filem
 subroutine filem (jq,nset,qset)
 dimension nset(999),qset(999)
 common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
 common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
 external error(descriptors), set(descriptors)
 if (mfa - id ) 2,2,3
 3 write (nprnt,4)
 4 format (//24h overlap set given below/)
 call error (87,nset,qset)
 2 indx = (mfa - 1) * imm
 do 1 i = 1,imm
 indx = indx + 1
 1 qset(indx) = atrib(i)
 indx = (mfa - 1) * mxx
 do 10 i = 1,im
 indx = indx + 1
 10 nset(indx) = jtrib(i)
 call set (jq,nset,qset)
 return
 end
cSTART histo
 subroutine histo (x1,a,w,n)
 common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
 common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
 external exit
 if (n-nhist) 11,11,2
 2 print 250,n
 250 format(19h error in histogram,i4//)
 call exit
 11 if(n)2,2,3
 3 x = x1 - a
 if (x)6,7,7
 6 ic = 1
 go to 8
 7 ic = x/w + 2.
 if (ic - ncels(n) - 1) 8,8,9
 9 ic = ncels(n)+2
 8 jcels(n,ic) = jcels(n,ic) + 1
 return
 end
subroutine montr(nset,qset)
 dimension nset(999),qset(999)
 common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
 common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
 if (jevnt - 101) 9,7,9
 7 print 8
 8 format(44h do you want to see a gasp job storage dump?)
 print 10
 10 format(10h0=no,1=yes)
 read 11,no
 11 format (i1)
 if(no.eq.0)return
 write(nprnt,100) tnow
 100 format(1h1,2x,31h**gasp job storage area dump at,f10.4,
& 2x,12htime units**//)
 write (nprnt,200)
 200 format(5x,4hnset/)
 do 210 i=1,id
 il=(i-1)*mxx+1
 iv=il+mxx-1
 210 write(nprnt,90) i,(nset(ij),ij=il,iv)
 90 format(2x,i5,5x,12i8)
 write (nprnt,202)
 202 format(//5x,4hqset/)
 do 215 i=1,id
 il=(i-1)*imm+1
 iv=il+imm-1
 215 write(nprnt,95) i,(qset(ij),ij=il,iv)
 95 format(2x,i5,4x,5e13.6,/12x,5e13.6)
end file 11
 return
 9 if(mfe(1))3,6,1
 1 if(jmnit-1)5,4,3
 3 write (nprnt,199)
 199 format(///2x,26h error exit,type 99 error.)
endfile 11
 4 indx=mfe(1)
 il=(indx-1)*mxx+1
 iv=il+mxx-1
 write (nprnt,103) tnow,jtrib(1),(nset(i),i=il,iv)
 103 format (/2x,23hcurrent event....time =,f8.2,5x,7hevent =,i7,
&/2x,17hnext event(nset).,(6i8))
 il=(indx-1)*imm+1
 iv=il+imm-1
 write(nprnt,120)(qset(i),i=il,iv)
 120 format(/2x,19hnext event(qset)...,(4e12.4))
 5 return
 6 write (nprnt,104) tnow
 104 format (2x,19h file 1 is empty at,f10.2)
 go to 5
 end
cSTART tmst
 subroutine tmst (x,t,n,nset,qset)
 dimension nset(999),qset(999)
 common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
 common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
 external error (descriptors)
 if (n) 2,2,1
 2 call error(91,nset,qset)
 1 if(n-nstat)3,3,2
 3 tt= t-ssuma(n,1)
 ssuma(n,1) = ssuma(n,1) + tt
 ssuma(n,2) = ssuma(n,2)+x*tt
 ssuma(n,3) = ssuma(n,3)+x*x*tt
 ssuma(n,4) = amin (ssuma(n,4),x)
 ssuma(n,5) = amax (ssuma(n,5),x)
 return
 end
cSTART sumry
 subroutine sumry (nset,qset)
 dimension nset(999),qset(999)
 common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
 common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
 external exit,prntq(descriptors)
 write (nprnt,21)
 21 format (1h ,1x,23h**gasp summary report**/)
 write (nprnt,102) nproj,name,mon,nday,nyr,nrun
 102 format (1x,22hsimulation project no.,i4,2x,2hby,2x,
& 6a2//,1x,4hdate,i3,1h/,i3,1h/,i5,12x,10hrun number,i5//)
 if (nprms) 147,147,146
 146 do 64 i=1,nprms
 64 write (nprnt,107) i,(param(i,j),j=1,4)
 107 format(2x,14h parameter no.,i5,4f12.4)
 147 if(nclct)5,60,66
 5 write (nprnt,199)
 199 format(///2x,26herror exit, type 98 error.)
 call exit
 66 write (nprnt,23)
 23 format (//2x,18h**generated data**/ 2x,4hcode,4x,4hmean,6x,6hs
&.dev.,5x,4hmin.,7x,4hmax.,5x,4hobs./)
 do 2 i=1,nclct
 if(suma(i,3))5,62,61
 62 write (nprnt,63) i
 63 format(2x,i3,10x,18hno values recorded)
 go to 2
 61 xs = suma(i,1)
 xss = suma(i,2)
 xn = suma(i,3)
 avg = xs/xn
 std=(((xn*xss)-(xs*xs))/(xn*(xn-1.0)))**.5
 n = xn
 write (nprnt,24) i,avg,std,suma(i,4),suma(i,5),n
 24 format (2x,i3,4f11.4,i7)
 2 continue
 60 if(nstat)5,67,4
 4 write (nprnt,29)
 29 format (/2x,23h**time generated data**/ 2x,4hcode,4x,4hmean,6x,
&8hstd.dev.,5x,4hmin.,7x,4hmax.,3x,10htotal time/)
 do 6 i = 1,nstat
 if(ssuma(i,1))5,71,72
 71 write (nprnt,63) i
 go to 6
 72 xt = ssuma(i,1)
 xs = ssuma(i,2)
 xss = ssuma(i,3)
 avg = xs/xt
 std = (xss/xt-avg*avg)**.5
 write (nprnt,30) i,avg,std,ssuma(i,4),ssuma(i,5),xt
 30 format (2x,i3,5f11.4)
 6 continue
 67 if(nhist)5,75,9
 9 write (nprnt,25)
 25 format (/2x,37h**generated frequency distributions**/ 2x,4hcod
&e,20x,10hhistograms/)
 do 12 i=1,nhist
 ncl = ncels (i)+2
 12 write (nprnt,26) i,(jcels(i,j),j=1,ncl)
 26 format (2x,i3,5x,11i4/10x,11i4/)
 75 do 15 i = 1,noq
 iimm=i
 15 call prntq(iimm,nset,qset)
 return
 end
cSTART set
 subroutine set (jq,nset,qset)
 dimension nset(999),qset(999)
 common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
 common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
 external error(descriptors)
 if (init-1) 27,28,27
 28 kol = 7777
 kof = 8888
 kle = 9999
 mx = im +1
 mxx=im+2
 maxqs = id * imm
 maxns = id * mxx
 do 2 j = 1,maxqs
 2 qset(j) = 0.0
 do 4 j = 1,maxns
 4 nset(j) = 0
 do 1 i = 1,id
 indx = i * mxx
 nset(indx - 1) = i + 1
 1 nset(indx) = i - 1
 nset(maxns - 1) = kof
 do 3 k = 1,noq

 nq(k)=0
 mlc(k)=0
 mfe(k)=0
 maxnq(k) = 0
 mle(k)=0
 enq(k)=0.0
 vnq(k)=0.0
 3 qtime(k)=tnow

 mfa = 1
 init = 0
 out = 0.0
 return
 27 mfex = mfe(jq)

 knt = 2
 ks = krank(jq)

 ksj = 1
 if (ks - 100) 1020,100,1000
 1000 ksj = 2
 ks = ks - 100
 1020 if (out - 1.0) 8,5,100
 8 indx = mfa * mxx - 1
 nxfa = nset(indx)
 if (inn(jq)-1) 100,7,6
 7 mlex=mle(jq)
 if (mlex) 100,10,11
 10 indx = mfa * mxx
 nset(indx) = kle
 mfe(jq) = mfa
 17 indx = mfa * mxx - 1
 nset(indx) = kol
 mle(jq) = mfa
 14 mfa =nxfa
 indx = nxfa*mxx
 nset(indx) = kle
 xnq = nq(jq)
 enq(jq) = enq(jq)+xnq*(tnow-qtime(jq))
 vnq(jq)=vnq(jq)+xnq*xnq*(tnow-qtime(jq))
 qtime(jq) = tnow
 nq(jq) = nq(jq) + 1
 maxnq(jq)=xmax(maxnq(jq),nq(jq))
 mlc(jq)=mfe(jq)
 return
 11 go to (1100,1120),ksj
 1100 indx1 = (mfa - 1) * imm + ks
 indx2 = (mlex - 1) * imm + ks
 if (qset(indx1) - qset(indx2)) 12,13,13
 1120 indx1 = (mfa - 1) * mxx + ks
 indx2 = (mlex - 1) * mxx + ks
 if (nset(indx1) - nset(indx2)) 12,13,13
 13 indx = mlex * mxx - 1
 msu = nset(indx)
 nset(indx) = mfa
 indx = mfa * mxx
 nset(indx) = mlex
 go to (18,17),knt
 18 indx = mfa * mxx - 1
 nset(indx) = msu
 indx = msu * mxx
 nset(indx) = mfa
 go to 14
 12 knt = 1
 indx = mlex * mxx
 mlex = nset(indx)
 if(mlex-kle) 11,16,11
 16 indx = mfa * mxx
 nset(indx) = kle
 mfe(jq) = mfa
 26 indx = mfa * mxx - 1
 nset(indx) = mfex
 indx = mfex * mxx
 nset(indx) = mfa
 go to 14
 6 if (mfex) 100,10,19
 19 go to (1200,1220),ksj
 1200 indx1 = (mfa - 1) * imm + ks
 indx2 = (mfex - 1) * imm + ks
 if (qset(indx1) - qset(indx2)) 20,21,21
 1220 indx1 = (mfa - 1) * mxx + ks
 indx2 = (mfex - 1) * mxx + ks
 if (nset(indx1) - nset(indx2)) 20,21,21
 20 knt = 1
 mpre = mfex
 indx = mfex * mxx - 1
 mfex = nset(indx)
 if (mfex-kol) 19,24,19
 21 go to (22,16),knt
 22 knt = 2
 24 indx = mfa * mxx
 nset(indx) = mpre
 indx = mpre * mxx - 1
 nset(indx) = mfa
 go to (17,26), knt
 5 out = 0.0
 indx = (mlc(jq) - 1) * imm
 do 32 i=1,imm
 indx = indx + 1
 32 qset(indx) = 0.0
 indx = (mlc(jq) - 1) * mxx
 do 1300 i= 1,im
 indx = indx + 1
 1300 nset(indx) = 0
 indx = mlc(jq) * mxx
 jl = nset(indx - 1)
 jk = nset(indx)
 if (jl- kol) 33,34,33
 33 if (jk- kle) 35,36,35
 35 indx = jk * mxx - 1
 nset(indx) = jl
 indx = jl * mxx
 nset(indx) = jk
 37 indx = mlc(jq) * mxx - 1
 nset(indx) = mfa
 nset(indx+1) = kle
 indx=mfa*mxx
 nset(indx)=mlc(jq)
 mfa = mlc(jq)
 mlc(jq) = mfe(jq)
 xnq = nq(jq)
 enq(jq)=enq(jq)+xnq*(tnow-qtime(jq))
 vnq(jq)=vnq(jq)+xnq*xnq*(tnow-qtime(jq))
 qtime(jq) = tnow
 nq(jq) = nq(jq)-1
 return
 36 indx = jl * mxx
 nset(indx) = kle
 mfe(jq) = jl
 go to 37
 34 if (jk-kle) 38,39,38
 38 indx = jk *mxx - 1
 nset(indx) = kol
 mle(jq) = jk
 go to 37
 39 mfe(jq) = 0
 mle(jq) = 0
 go to 37
 100 call error (88,nset,qset)
 return
 end
cSTART rmove
 subroutine rmove (kcoll,jq,nset,qset)
 dimension nset(999),qset(999),kcoll(4)
 common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
 common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&6(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
 external error(descriptors), set(descriptors)
 kcol = kcoll(1)
 if (kcol) 16,16,2
 16 call error (97,nset,qset)
 2 mlc(jq) = kcol
 indx = (kcol - 1) * imm
 do 3 i = 1,imm
 indx = indx + 1
 3 atrib(i) = qset(indx)
 indx = (kcol - 1) * mxx
 do 10 i =1,im
 indx = indx + 1
 10 jtrib(i) = nset(indx)
 out = 1.
 call set (jq,nset,qset)
 return
 end
cSTART prntq
 subroutine prntq (jq,nset,qset)
 dimension nset(999),qset(999)
 common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
 common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
 write (nprnt,100) jq
 if (tnow - tbeg) 12,12,13
 12 write (nprnt,105)
 105 format(/2x,25h no printout tnow = tbeg //)
 go to 2
 13 xnq=nq(jq)
 x=(enq(jq)+xnq*(tnow-qtime(jq)))/(tnow-tbeg)
 std=((vnq(jq)+xnq*xnq*(tnow-qtime(jq)))/(tnow-tbeg)-x*x)**0.5
 write (nprnt,104) x,std,maxnq(jq)
 write (nprnt,101)

& nsq = 1
 write(nprnt,200)
 200 format(2x,4hnset/)
 230 line = mfe(jq)
 if (line-1) 4,1,1
 4 write (nprnt,102)
 2 return
 1 l1 = line - 1
 go to (202,201),nsq
 202 indx = l1 * mxx
 ib = indx + 1
 ie = indx + mxx
 write (nprnt,106) line, (nset(i),i=ib,ie)
 go to 210
 201 indx = l1 * imm
 ib = indx + 1
 ie = indx + imm
 write (nprnt,103) line, (qset(i),i=ib,ie)
 210 indx = line * mxx - 1
 line = nset(indx)
 if (line-7777) 1,2220,5
 2220 if (nsq-2) 221,2,2
 221 nsq = nsq + 1
 write (nprnt,205)
 205 format (//2x,4hqset/)
 go to 230
 5 write (nprnt,199)
 199 format(///2x,26herror exit, type 94 error.)
 100 format(//2x,24h file printout, file no.,i3)
 101 format (/2x,14h file contents//)
 102 format(/2x,17hthe file is empty//)
 103 format (2x,i5,4x,5e13.6/12x,5e13.6)
 104 format(/2x,27haverage number in file was,f10.4,/2x,9hstd. dev
& 18x,f10.4,/2x,7hmaximum,24x,i4)
 106 format (2x,i5,5x,8i8,/13x,4i8)
 stop
 end
subroutine randu(i,x)
x=flat(i)
return
end
function xmax(iarg1,iarg2)
if (iarg1 - iarg2) 2,2,1
1 xmax = iarg1
return
2 xmax = iarg2
return
end
cSTART events
      subroutine events (ix, nset, qset)
      dimension nset(60), qset(80)
      common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
      common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssum
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
      common xist1,xist2,xisys,tld,tbd,xbuz(2),titem,cbalk,tisys,block
      go to (1,2,2,3),ix
    1 call arrvl (nset, qset)
      return
    2 call endsv (nset, qset)
      return
    3 call endsm (nset, qset)
      return

      end
cSTART arrvl
      subroutine arrvl (nset, qset)
      dimension nset(60), qset(80)
      common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
      common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
      common xist1,xist2,xisys,tld,tbd,xbuz(2),titem,cbalk,tisys,block
      call drand(iseed,rnum)
      atrib(1)= tnow- param(1,1)*alog(rnum)
      jtrib(1) = 1
      call filem (1, nset, qset)
      titem=titem+1.
      if(xist1 - 5.)11,10,10
   10 cbalk = cbalk + 1.
      return
   11 call tmst (xisys, tnow, 1, nset, qset)
      if(xist1)7,8,9
    7 call error (31, nset, qset)
      return
    8 xist1 = xist1 + 1.
      xisys = xisys + 1.
      call tmst (xbuz(1), tnow, 2, nset, qset)
      xbuz(1) = 1.
      call drand(iseed,rnum)
      atrib(1) = tnow - param(2,1)*alog(rnum)
      jtrib(1) = 2
      atrib(3) = tnow
      call filem (1, nset, qset)
      return
    9 xist1 = xist1 + 1.
      xisys = xisys + 1.
      atrib(3)=tnow
      call filem (2, nset, qset)
      return

      end
cSTART endsv
      subroutine endsv (nset, qset)
      dimension nset(60), qset(80)
      common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
      common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
      common xist1,xist2,xisys,tld,tbd,xbuz(2),titem,cbalk,tisys,block
      ii = jtrib(1)
      ii = ii - 1
      go to (1,2),ii
    1 jj = xist2 + 1.
      go to (3,4,5,6),jj
    3 xist1 = xist1 - 1.
      xist2 = xist2 + 1.
      y = tnow + .2
      call tmst (xbuz(2), y, 3, nset, qset)
      xbuz(2) = 1.
      call drand(iseed,rnum)
      atrib(1) = y - param (3,1)*alog(rnum)
      jtrib(1) = 3
      call filem (1, nset, qset)
   10 if(nq(2))7,8,9
    7 call error (41, nset, qset)
      return
    8 call tmst (xbuz(1), tnow, 2, nset, qset)
      xbuz(1) = 0.
      return
    9 call rmove (mfe(2), 2, nset, qset)
      call drand(iseed,rnum)
      atrib(1) = tnow - param (2,1)*alog(rnum)
      jtrib(1) = 2
      call filem (1, nset, qset)
      call tmst (xbuz(1), tnow, 2, nset, qset)
      xbuz(1) = 1.0
      return
    4 xist1 = xist1 - 1.
      xist2 = xist2 + 1.
      atrib(4) = tnow + .1
      call filem (3, nset, qset)
      go to 10
    5 xist1 = xist1 - 1.
      xist2 = xist2 + 1.
      atrib(4) = tnow
      call filem (3, nset, qset)
      go to 10
    6 call tmst (block, tnow, 4, nset, qset)
      block = 100.
      call tmst (xbuz(1), tnow, 2, nset, qset)
      xbuz(1) = 0.
      call filem (4, nset, qset)
      return
    2 call tmst (xisys, tnow, 1, nset, qset)
      xisys = xisys - 1.
      xist2 = xist2 - 1.
      tisys = tnow - atrib(3)
      call colct (tisys, 1, nset, qset)
      call histo(tisys,.5,.5,1)
      tbd = tnow - tld
      tld = tnow
      call colct (tbd, 2, nset, qset)
      call histo(tbd,0.0,.175,2)
      if(nq(3))7,11,12
   11 call tmst (xbuz(2), tnow, 3, nset, qset)
      xbuz(2) = 0.0
      return
   12 call rmove (mfe(3), 3, nset, qset)
      call drand(iseed,rnum)
      atrib(1) = tnow - param (3,1)*alog(rnum)
      jtrib(1) = 3
      call filem (1, nset, qset)
      if(block)13,14,15
   13 call error (51, nset, qset)
   14 return
   15 call tmst (block, tnow, 4, nset, qset)
      block = 0.0
      call rmove (mfe(4), 4, nset, qset)
      go to 5
      end
cSTART otput
      subroutine otput (nset, qset)
      dimension nset(60), qset(80)
      common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
      common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
     common xist1,xist2,xisys,tld,tbd,xbuz(2),titem,cbalk,tisys,block
      write(nprnt,10)  (param(k,1),k=1,3)
   10 format(/1x,28hmean time between arrivals =f4.2 /,1x,33hmean servi
&ce time for station 1 =,f4.2 /,1x,33hmean service time for statio
&n 2 =,f4.2)
      ybalk = cbalk*100./titem
      write (nprnt,20) ybalk,cbalk,titem
   20 format(1x,32hpercent of items subcontracted =,f6.2 /1x,
&31hnumber of items subcontracted =,f5.0 /,1x,13htotal items =,
&f6.0)
      return
      end
cSTART endsm
      subroutine endsm (nset, qset)
      dimension nset(60), qset(80)
      common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
      common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
      common xist1,xist2,xisys,tld,tbd,xbuz(2),titem,cbalk,tisys,block
   12 if(nq(1))7,10,9
    7 call error (61, nset, qset)
    9 call rmove (mfe(1), 1, nset, qset)
      tnow = atrib(1)
      jj = jtrib(1)
      go to (12,13,13),jj
   13 call endsv (nset, qset)
      go to 12
   10 call tmst (xisys  , tnow, 1, nset, qset)
      call tmst (xbuz(1), tnow, 2, nset, qset)
      call tmst (xbuz(2), tnow, 3, nset, qset)
      call tmst (block  , tnow, 4, nset, qset)
      mstop = -1
      norpt = 0
      return
      end
c   ********eigenvalues.eigemvectors of a real symme
c   ********matrix-jacobi-corbato method************
 subroutine eig1(atrix2,atrix1,norow1,eps,amax,ir,norow2,norow3)
 dimension atrix1(norow3,norow1),atrix2(100000),amax(norow1),ir(norow1)
 isubs(i,j) = norows*(i-1)+j-(i*(i-1))/2
 isubd(i,j) = norow2*(j-1)+i
 1000 norows=norow1
 1001 enorow=norows
 1003 if(norow2.ne.1) go to 1007
 1005 mfudge=1
 go to 1010
 1007 mfudge = 0
c   mfudge = 1 if atrix2 is stored in a single dimensioned
c   mfudge = 0 if atrix2 is stored in a double dimensioned
 1010 if(eps)3601,1011,1012
 1011 sqeps = 1.0e-6
 go to 1015
 1012 sqeps = sqrt(eps)
 1015 nocols = norows
 do 1060 i = 1,norows
 do 1050 j=1,nocols
 1050 atrix1(i,j) = 0.0
 amax(i) = 0.
 1060 atrix1(i,i) = 1.0
 nminus = norows - 1
 if(nminus)3601,3601,1080
 1080 do 1090 j=2,norows
 k = j-1
 do 1090 i = 1,k
 if (mfudge) 3601,1082,1083
 1082 loc = isubd(i,j)
 go to 1084
 1083 loc = isubs(i,j)
 1084 if (abs(atrix2(loc))-amax(j)) 1090,1085,1085
c   *****ir-row of max.abs.element for column j ********
c   *****amax-max.abs. element of column j(not including diago
 1085 amax(j) = abs(atrix2(loc))
 ir(j) = i
 1090 continue
 1100 biga=0.
 do 1125 j=2,norows
 if(amax(j)-biga)1125,1120,1120
 1120 biga=amax(j)
 l=ir(j)
 m=j
 1125 continue
c   l = row,m = col of max. element
 if (mfudge) 3601,1126,1127
 1126 lam = isubd(l,m)
 mam = lam-l+m
 lal = isubd(l,l)
 go to 1138
 1127 lam = isubs(l,m)
 lal = lam-m+l
 mam = isubs(m,m)
c   *****convergence test   ******
 1138 if (enorow*abs(atrix2(lam))-sqeps) 3490,1150,1150
 1150 lminus =l-1
 lplust = l+1
 mminus = m-1

 mplust = m+1

c   ****calculations associated with max.element at point (l,m)
c   ****only those elements change that are on row l,m and col
 1200 r = sqrt((atrix2(lal)-atrix2(mam))**2+4.*atrix2(lam)*atrix2(lam))
 t = atrix2(lal)-atrix2(mam)
 1400 cos = sqrt(.5*(1.+abs(t)/r))
 sin = -atrix2(lam)/(r*cos)
 if (t) 1401,1402,1402
 1401 sin = -sin
 1402 sinsin = sin*sin
 coscos = cos*cos
 sincos = sin*cos
 prod = 2.0*atrix2(lam)*sincos
 2000 atrix2(lam) = 0.0
 amax(l) = 0.
 amax(m) = 0.
 2010 if (l .eq. 1) go to 2110
 do 2040 i=1,lminus
 if (mfudge) 3601,2012,2013
 2012 ial = isubd(i,l)
 iam = isubd(i,m)
 go to 2014
 2013 ial = isubs(i,l)
 iam = ial-l+m
 2014 temper = atrix2(ial)
c   *****claculate new elements on col.l and m to row l-1,upda
 atrix2(ial) = atrix2(ial)*cos-atrix2(iam)*sin
 atrix2(iam) = atrix2(iam)*cos+temper*sin
 if (abs(atrix2(ial))-amax(l)) 2030,2015,2015
 2015 amax(l) = abs(atrix2(ial))
 ir(l) = i
 2030 if (abs(atrix2(iam))-amax(m)) 2040,2035,2035
 2035 amax(m) = abs(atrix2(iam))
 ir(m) = i
 2040 continue
c   ****if l .eq. 1 max. element is on row 1 ******
 2110 if (l+1 .eq. m) go to 2210
 do 2140 k=lplust,mminus
 if (mfudge) 3601,2112,2113
 2112 lak = isubd(l,k)
 ltemp = lak -l
 kam = isubd(k,m)
 go to 2114
 2113 lak = isubs(l,k)
 kam = isubs(k,m)
 2114 temper = atrix2(lak)
c   *****calculate new elments across row l and down col m to m-
 atrix2(lak) = atrix2(lak)*cos-atrix2(kam)*sin
 atrix2(kam) = atrix2(kam)*cos+temper*sin
 if (ir(k)-l) 2115,2120,2115
 2115 if (abs(atrix2(lak))-amax(k)) 2130,2117,2117
 2117 amax(k) = abs(atrix2(lak))
 ir(k) = l
 go to 2130
 2120 amax(k)=0.
 km1=k-1
 do 2125 it=1,km1
 if (mfudge) 3601,2121,2122
 2121 itak = ltemp+it
 go to 2123
 2122 itak = isubs(it,k)
 2123 if (abs(atrix2(itak))-amax(k)) 2125,2124,2124
 2124 amax(k) = abs(atrix2(itak))
 ir(k) = it
 2125 continue
 2130 if (abs(atrix2(kam))-amax(m)) 2140,2135,2135
 2135 amax(m) = abs(atrix2(kam))
 ir(m) = k
 2140 continue
c   *****if l+1 =m, no calculations across row l or down column m
 2210 if (m .eq. norows) go to 2310
 do 2240 j=mplust,norows
 if (mfudge) 3601,2212,2213
 2212 laj = isubd(l,j)
 ltemp = laj-l
 maj = ltemp+m
 go to 2214
 2213 laj = isubs(l,j)
 maj= isubs(m,j)
 2214 temper = atrix2(laj)
c   *****calculate new elements across row l and row m
 atrix2(laj) = atrix2(laj)*cos-atrix2(maj)*sin
 atrix2(maj) = atrix2(maj)*cos+temper*sin
 if(ir(j)-l) 2215,2225,2215
 2215 if (ir(j)-m) 2217,2225,2217
 2217 if (abs(atrix2(laj))-amax(j)) 2220,2218,2218
 2218 amax(j) = abs(atrix2(laj))
 ir(j) =l
 2220 if (abs(atrix2(maj))-amax(j)) 2240,2221,2221
 2221 amax(j) = abs(atrix2(maj))
 ir(j) = m
 go to 2240
 2225 amax(j)=0.
 jm1=j-1
 do 2230 it=1,jm1
 if(mfudge) 3601,2226,2227
 2226 itaj=ltemp+it
 go to 2228
 2227 itaj= isubs(it,j)
 2228 if (abs(atrix2(itaj))-amax(j)) 2230,2229,2229
 2229 amax(j) = abs(atrix2(itaj))
 ir(j)=it
 2230 continue
 2240 continue
c   *****if m .eq. norows,no more elements change*****
 2310 temper = atrix2(lal)
c   *****compute new diagonal elements ****
 atrix2(lal)=atrix2(lal)*coscos-prod+atrix2(mam)*sinsin
 atrix2(mam)=atrix2(mam)*coscos+prod+temper*sinsin
 2500 do 2520 i=1,norows
 temper = atrix1(i,l)
c   *****compute vectors, column l and m
 atrix1(i,l) = atrix1(i,l)*cos-atrix1(i,m)*sin
 atrix1(i,m) = atrix1(i,m)*cos+temper*sin
 2520 continue
 go to 1100
 3490 if(mfudge) 3601,3601,3500
 3500 do 3505 i = 2, norows
c   *****store eigenvalues in first n locations of atrix2*****
 iai = isubs(i,i)
 3505 atrix2(i) = atrix2(iai)
 3601 return
 end
c   ampb1
cSTART ampb1
 subroutine ampb1(ind,exit,temp,x,dx,y,f,n,icount,niter,mtst)
c      **************adams-moulton integration***************
 dimension temp(999),y(n),f(n)
external rkpb1(descriptors)
c      **************calculate derivatives*******************
c      **************store functions and derivatives*********
 101 if(ind .ne. 0) go to 111
 icount =1
 go to 302
 111 m=n+1
 do 191 i=1,m
c      **************store previous sets of derivatives******
 do 121 j1=1,icount
 i1=i+(icount+3-j1)*m
 i2=i1+m
 121 temp(i2)=temp(i1)
 191 continue
 201 if(ind.le.0) go to 301
c      **************would like to double dx next time*******
c      **************can not do if icount less than 6********
 if(icount.lt.6) go to 301
 dx=2.0*dx
 do 261 j1=1,3
 do 241 i=1,m
 i1=i+(j1+3)*m
 if(i.eq.1) go to 221
 i2=i+(2*j1+3)*m
 temp(i1)=temp(i2)
 go to 241
 221 temp(i1)=2.0*temp(i1)
 241 continue
 261 continue
 icount=3
c      **************increment count of sets of derivatives**
 301 icount=min0(icount+1,6)
 302 nnn=n
 call rkpb1(exit,temp,x,dx,y,f,nnn)
 ind=-1
 return
 end
cSTART gecos.Wed1752.65

c   ampb2
c      **************integrate to next point*****************
 subroutine ampb2(ind,exit,temp,x,dx,y,f,n,icount,niter,mtst)
 dimension temp(999),y(n),f(n),ak(2)
 external rkpb2(descriptors),exit
 data ak/.92962963,-.070370370/
 m=n+1
 n1=niter+1
c      **************restore previous values*****************
 x=temp(2*m+1)
 do 421 i=2,m
 ip2m=i+2*m
 y(i-1)=temp(ip2m)
 ip3m=ip2m+m
 421 f(i-1)=temp(ip3m)
 if(ind.gt.0) go to 501
 if(icount.ge.4) go to 701
 go to 601
c      **************reduce dx*******************************
 501 dx=dx/2.0

 icount=1
c      **************integration by runge-kutta**************
 601 nnn=n
 call rkpb2(exit,temp,x,dx,y,f,nnn)
 go to 1001
c      **************integration by adams-moulton************
c      **************independent variable********************
 701 x=x+dx
 temp(m+1)=x
c      **************predictor values************************
 801 do 821 i=2,m
 ipm=i+m
 ip2m=ipm+m
 ip3m=ip2m+m
 ip4m=ip3m+m
 ip5m=ip4m+m
 ip6m=ip5m+m
 y(i-1)=temp(ip2m)+dx*(55.0*temp(ip3m)-59.0*temp(ip4m)+
& 37.0*temp(ip5m)-9.0*temp(ip6m))/24.0
 temp(ipm)=y(i-1)
 if(mtst.ne.0.and.icount.gt.4) y(i-1)=y(i-1)+ak(1)/ak(2)*temp
& (i)
 821 continue
c      **************corrector values************************
 901 do 991 j1=1,n1
  call exit
 do 921 i=2,m
 ipm=i+m
 ip2m=ipm+m
 ip3m=ip2m+m
 ip4m=ip3m+m
 ip5m=ip4m+m
 y(i-1)=temp(ip2m)+dx*(9.0*f(i-1)+19.0*temp(ip3m)-5.0*temp(ip
& 4m)+
& temp(ip5m))/24.0
 temp(i)=ak(2)*(y(i-1)-temp(ipm))
 if(mtst.ne.0) y(i-1)=y(i-1)+temp(i)
 921 continue
 991 continue
c      **************restore normal mode*********************
 1001 ind=-1
 5001 return
 end
cSTART gecos.Fri1603.50

c   secant
        subroutine secant(n,ni,cc,fm,y,f,q,z,s,g,x,idim,eval,ierr)
 dimension y(999),f(999),q(999),z(999),s(999),g(idim,999),x(idim,999)
       fm=1.0e+37
       np=n+1
       jpa=0
       nit=ni
    10 y(np)=1.
       f(np)=1.
       z(np)=1.
       bf=abs(y(1))
       do 20 i=2,n
    20 bf=amax1(bf,abs(y(i)))
       do 30 i=1,n
       if(y(i))30,25,30
    25 y(i)=(1.0e-6)*bf
    30 z(i)=y(i)
       l=np
    35 y(l-1)=1.1*y(l-1)
    40 y(l)=z(l)
   45 call eval(f,y)
       af=abs(f(1))
       if(n-1)55,55,47
    47 do 50 i=2,n
    50 af=amax1(af,abs(f(i)))
    55 if(af-fm)60,80,80
    60 fm=af
       if (fm)5000,5000,65
    65 do 70 i=1,np
    70 s(i)=y(i)
    80 if(l)120,120,85
    85 do 90 i=1,np
       x(i,l)=y(i)
    90 g(i,l)=f(i)
       l=l-1
       if(l-1)100,40,35
   100 call mtinv(g,np,np,idim,q)
   105 do 110 i=1,np
       y(i)=0.
       do 110 j=1,np
   110 y(i)=y(i)+x(i,j)*g(j,np)
       if(abs(y(np)-1.)-.01)45,190,190
   120 af=abs(y(1)-z(1))
       z(1)=y(1)
       do 130 i=2,n
       af=amax1(af,abs(y(i)-z(i)))
   130 z(i)=y(i)
       if(af-cc)5000,5000,140
   140 nit=nit-1
       if(nit)192,192,145
   145 do 150 i=1,np
       q(i)=0.
       do 150 j=1,np
   150 q(i)=q(i)+g(i,j)*f(j)
       af=abs(q(1))
       k=1
       do 160 i=2,np
       afi=abs(q(i))
       if(afi-af)160,160,155
   155 af=afi
       k=i
   160 continue
       do 170 i=1,np
   170 g(k,i)=g(k,i)/q(k)
       do 181 i=1,np
       x(i,k)=y(i)
       if(i-k)175,181,175
   175 do 180 j=1,np
       g(i,j)=g(i,j)-q(i)*g(k,j)
   180 continue
   181 continue
       go to 105
   190 ierr=1
       go to 193
   192 ierr=2
   193 if(jpa)195,195,5010
   195 jpa=1
       nit=ni
       do 200 i=1,n
   200 y(i)=s(i)
       go to 10
  5000 ierr=0
  5010 return
       end
 subroutine corrl2
c **************  corrl2  ***************************************
cSTART gecos.Wed1751.30
        dimension x(2000),xr(1000),t(2)
  2     format(///" number of elements in each array")
  3     format(//" input x,y arrays"/)
  4     format(v)
   50   format(v)
        print 2
        read 4, n
        jn=2*n
        print 3
        read 50, (x(i),i=1,jn)
        t(1)=0.;t(2)=0.
        do 20 j=1,2
        do 19 i=1,n
        idum=(j-1)*n+i
        num=0;num1=0;num2=0
        do 10 k=1,n
        idum1=(j-1)*n+k
        if(x(idum)-x(idum1))10,6,5
  6     num1=num1+1
  5     num=num+1
  10    continue
        num=num-num1+1
        if(num1-1)8,16,8
  8     num2=num
        num3=num1-1
        do 15 kk=1,num3
        num2=num2+num+kk
        xnum2=num2
        xnum1=num1
        xnum=xnum2/xnum1
  15    continue
        t(j)=t(j)+(xnum1**3-xnum1)/(12.*xnum1)
        go to 17
  16    xnum=num
    17  xr(i)=xnum
  19     continue
        do 20 i=1,n
        idum=(j-1)*n+i
        x(idum)=xr(i)
   20   continue
        sum=0.
        do 30 i=1,n
        idum2=(i-1)*2
        xr(i)=x(i)-x(n+i)
        sum=sum+(x(i)-x(n+i))**2
  30    continue
        xn=n
        xsqr=(xn**3-xn)/12.-t(1)
        ysqr=(xn**3-xn)/12.-t(2)
        rs=(xsqr+ysqr-sum)/(2.*sqrt(xsqr*ysqr))
        print 25,rs
  25    format(//" The spearman correlation coefficient: rs=",e20.8)
        if(n-10)200,100,100
  100    tt=rs*sqrt((xn-2.)/(1.-rs**2))
        xnn=n-2
        prob=(1.-tdist(0,xnn,tt))/2.
        print 150,prob
  150   format(//" This value of rs is significant at the",
&       e20.8,"  level"///)
  200   stop
        end
 subroutine eigsr
c   Modified Feb 1975
       dimension a(25,25),b(25,25),am(25),ir(25),d(25)
 character iy*4,icy*4
  1 print 2
 2     format(/"          EIGSR"/)
       print 100
 100   format(49h Do you desire user instructions, type yes or no?)
       read 1010,iy
       if("yes" .ne. iy)go to 200
      print 102
  102 format(/" This program finds the Eigenvalues and Eigenvectors of a read")
       print 142
 142   format(47h symmetric matrix by the jacobi-corbato method.//)
       print 150
 150   format(27h The matrix is of the form:/)
     print 103
 103 format(" a11  a12  a1n    where the a(i,j) are real(floating point) and")
     print 106
 106 format(" a21  a22  a2n    n(fixed point) cannot exceed 25:")
 155 format(14h an1  an2  ann//)
       print 155
 160   format(49h Since the matrix is symmetric, only the elements)
 161   format(37h on and above the diagonal are input.//)
       print 160
       print 161
      print 107
 107 format(/" The program types Order ?    the user types the order of the matrix."/)
       print 108
 108 format(" The program types  a(1,1)=   the user types the first row elements")
     print 112
 112 format (" the program types  a(2,2)=   the user types the
& second row starting")
 169   format(48h     with the diagonal element.-etc.-  to a(n,n)/)
       print 169
   print 114
 114 format (" Input is typed in with each value containing a decimal point and
&separated by commas, a carriage return ends the field."/)
   print 115
 115 format(" After a(n,n) is input, the program provides the
& opportunity")
 182   format(46h and instructions for correcting typing errors//)
       print 182
 190   format(15h Now you try it//)
       print 190
 200   print 201
 201   format(7h Order?)
       read 1004, n
  1004 format (v)
       if(n .le. 0)go to 450
       nij=n*(n+1)/2
       nm1=n-1
       do 300 i=1,n
       print 1000,i,i
 300   read 1006, (a(i,j),j=i,n)
 1006  format(v)
 1000  format(3h a(,i2,1h,,i2,3h) =)
     print 118
     print 119
 118 format(" Are any of the above a(i,j) elements typed incorrectly?")
 119 format(" If user wishes to correct an element, type yes,
& otherwise type no")
 305  format(" Any corrections?")
  306  print 305
 307   read 1010, icy
 308   format(27h illegal command, try again)
       if("no" .eq. icy)go to 350
 315   if("yes" .eq. icy)go to 320
 316   print 308
       go to 306
    320 print 120
 120 format (" correct element by typing i subscript(row)")
     print 122
 122 format(" j subscript, value, carriage return")
 325   read 1008, i,j,atemp
 1008 format(v)
 335   if(i-j)340,340,336
    336 print 130
 130 format(" illegal subscript, i greater than j, or i or
& j greater than n. try again")
       go to 325
 340   if(i-n)341,341,336
  341  if(j-n)342,342,336
 342   a(i,j)=atemp
 343   print 305
       read 1010, icy
  1010 format(v)
       if("no"  .eq. icy)go to 350
 345   if("yes" .eq. icy)go to 325
 346   print 308
       go to 343
 350   na=25
       nb=25
       eps=0.
       do 352 i=1,n
       d(i)=a(i,i)
       do 352 j=i,n
 352   a(j,i)=a(i,j)
       call eig1(a,b,n,eps,am,ir,na,nb)
 369   format(//7h-Order=,i3)
 370   format(//2x,16h0     Eigenvalue,i3,16h     Eigenvector/)
 371   format(1pe20.7,1pe19.7)
 372   format(1pe39.7)
       print 369,n
       do 400 i=1,n
       print 370,i
       print 371,a(i,i),b(1,i)
       do 400 j=2,n
  400  print 372,b(j,i)
       do 410 i=1,n
       a(i,i)=d(i)
       do 410 j=i,n
 410   a(i,j)=a(j,i)
       t=0.
       do 430 i=1,n
       do 420 j=1,n
       d(j)=0.
       do 420 k=1,n
 420   d(j)=d(j)+a(j,k)*b(k,i)
       do 430 j=1,n
       if(j-i)422,430,422
 422   ts=0.
       do 425 k=1,n
 425   ts=ts+b(k,j)*d(k)
       t=t+ts*ts
 430   continue
       print 250

 250 format(" The sum of the squares of the off diagonal elements of")
 443 format("    xt a x =",1pe16.7,"   where xt = x-transpose")
       print 443,t
       go to 200
 450   stop
       end
 subroutine gaspsamp
c   gaspsamp --- sample driver routines for gaspiia
cSTART gecos.Wed1738.40
      dimension nset(60), qset(80)
      common id,im,init,jevnt,jmnit,mfa,mstop,mx,mxc,nclct,nhist,
&noq,norpt,not,nprms,nrun,nruns,nstat,out,iseed,tnow,
&tbeg,tfin,mxx,nprnt,ncrdr,nep,vnq(4),imm,maxqs,maxns
      common atrib(10),enq(4),inn(4),jcels(5,22),krank(4),maxnq(4),m
&fe(4),mlc(4),mle(4),ncels(5),nq(4),param(20,4),qtime(4),ssuma
&(10,5),suma(10,5),name(6),nproj,mon,nday,nyr,jclr,jtrib(12)
      common xist1,xist2,xisys,tld,tbd,xbuz(2),titem,cbalk,tisys,block
 external gasp (descriptors)
      xbuz(1)=1.0
      xbuz(2)=1.0
      xist1 = 4.
      xist2 = 1.
      xisys = 5.
      tld = 0.
      titem=0.
      cbalk = 0.
      block = 0.
      call gasp (nset, qset)
 stop
 end
 subroutine kilter
c    Modified Feb 1975
c     kilter - Out of kilter algorithm
 common arcs,i(100),j(100),cost(100),hi(100),lo(100),flow(100)
 dimension pi(100),na(100),nb(100)
 integer a,aok,cok,c,del,e,eps,src,snk
 integer flow,pi,arcs,cost,hi,lo
integer cst,hibnd,bndlo,istat
 character fname*16
 double precision stat
equivalence (istat,stat)
external com_err_ (descriptors), attach (descriptors), detach
 iin=10
  print 999
 999 format(" Out of kilter algo.")
 do 108 jr=1,25
 i (jr)=0
 j (jr)=0
 hi (jr)=0
 lo (jr)=0
 cost(jr)=0
 108 flow(jr)=0
 nodes=1
 call attach (istat)
 if (istat) 3535,4545,3535
 3535 call com_err_(istat,"kilter","error in attach")
stop
 4545 arcs=1
 109 read(iin,1,end=116)iwrd,k,l,cst,hibnd,bndlo
 1 format(v)
 cost(arcs)=cst
 hi(arcs)=hibnd
 lo(arcs)=bndlo
 i(arcs)=k
 j(arcs)=l
 if(l.gt.nodes)nodes=l
 arcs=arcs+1
 go to 109
 116 arcs=arcs-1
 do 5 m=1,arcs
 5 flow(m)=0
 do 11 m=1,nodes
 11 pi(m)=0
 do 161 k=1,100
 161 na(k)=0
 do 10 a=1,arcs
 if(lo(a).gt.hi(a)) go to 39
 10 continue
c set inf to max available integer
 16 inf=999999
 aok=0
c find out of kilter arc
 20 do 21 a=1,arcs
 ia=i(a)
 ja = j(a)
 c=cost(a)+pi( ia )-pi( ja )
 if((flow(a).lt.lo(a)).or.(c.lt.0.and.flow(a).lt.hi(a))) go to 22
 if((flow(a).gt.hi(a)).or.(c.gt.0.and.flow(a).gt.lo(a))) go to 23
 21 continue
c no remaining out of kilter arcs
 go to 38
 22 src=j(a)
 snk=i(a)
 e=+1
 go to 24
 23 src=i(a)
 snk=j(a)
 e=-1
 go to 24
c attempt to bring out of kilter arcs into kilter
 24 if((a.eq.aok).and.(na(src).ne.0)) go to 25
 aok=a
 do 26 n=1,nodes
 na(n)= 0
 26 nb(n)=0
 na(src)=iabs(snk)*e
 nb(src)=iabs(aok)*e
 25 cok=c
 27 lab=0
 do 30 a=1,arcs
 ia=i(a)
 ja = j(a)
 if((na(ia).eq.0.and.na(ja).eq.0).or.(na(ia).ne.0 .and.na(ja).ne.0)
&)goto 30
 c= cost(a)+pi(ia) - pi(ja)
 if(na(ia).eq.0) go to 28
 if(flow(a).ge.hi(a).or.(flow(a).ge.lo(a).and.c.gt.0))go to 30
 na(ja) = i(a)
 nb(ja) = a
 go to 29
 28 if(flow(a).le.lo(a).or.(flow(a).le.hi(a).and.c.lt.0))go to 30
 ia = i(a)
 na(ia) = -j(a)
 nb(ia) = -a
 29 lab = 1
c node labeled, test for breakthru
 if(na(snk).ne.0) go to 33
 30 continue
c no breakthru
 if(lab.ne.0) go to 27
c determine change to pi vector
 del = inf
 do 31 a=1,arcs
 ia = i(a)
 ja =j(a)
 if((na(ia).eq.0.and.na(ja).eq.0).or.(na(ia).ne.0.and.na(ja).ne.0)
&)goto 31
 c=cost(a)+pi(ia)-pi(ja)
 if(na(ja).eq.0.and.flow(a).lt.hi(a)) del= min0(del,c)
 if(na(ja).ne.0.and.flow(a).gt.lo(a)) del= min0(del,-c)
 31 continue
 ccok=cok
 if(del.eq.inf.and.(flow(aok).eq.hi(aok).or.flow(aok).eq.lo(aok)))
&del=abs(ccok)
 if(del.eq.inf) go to 39
c exit, no feasible flow pattern
c change pi vector by computed del
 do 32 n=1,nodes
 if(na(n).eq.0) pi(n)=pi(n)+del
 32 continue
c find another out-of-kilter arc
 go to 20
c breakthru, compute incremental flow
 33 eps=inf
 ni=src
 34 nj=iabs(na(ni))
 a=iabs(nb(ni))
 c=cost(a)-isign(iabs(pi(ni) -pi(nj)),nb(ni))
 if(nb(ni).lt.0) go to 35
 if(c.gt.0.and.flow(a).lt.lo(a)) eps=min0(eps,lo(a)-flow(a))
 if(c.le.0.and.flow(a).lt.hi(a)) eps=min0(eps,hi(a)-flow(a))
 go to 36
 35 if(c.lt.0.and.flow(a).gt.hi(a)) eps=min0(eps,flow(a)-hi(a))
 if(c.ge.0.and.flow(a).gt.lo(a)) eps=min0(eps,flow(a)-lo(a))
 36 ni=nj
 if(ni.ne.src) go to 34
c change flow vector by computed eps
 37 nj=iabs(na(ni))
 a=iabs(nb(ni))
 flow(a)=flow(a)+isign(eps,nb(ni))
 ni=nj
 if(ni.ne.src) go to 37
c find another out of kilter arc
 aok=0
 go to 20
  39 print 995
 995 format(" solution infeasible")
 38 call koutput
 call detach
 stop
end
cSTART koutput
 subroutine koutput
 common arcs,m(100),n(100),cost(100),hi(100),lo(100),flow(100)
 integer arcs,cost,hi,flow
 print 204
 204 format (//"         Network status and minimum cost flow")
 200 print 201
 201 format(1h0,4x,"i",5x,"j",6x,"cost",5x,"u.bnd",5x,"l.bnd",
&6x,"flow")
 do 202 i=1,arcs
 202 print 203,m(i),n(i),cost(i),hi(i),lo(i),flow(i)
 203 format(1h0,i5,i6,4i10)
 return
 end
 subroutine lnprog
c   lnprog -- Modified Feb 1975
 common t,t1,j,a(31,51),ci(31),cj(51),b(31),idc(51),idr(31)
 dimension isen(31),key(4)
 character ll*3,key*4,min*4,isen*1,fname*32
 double precision stat
integer istat
 equivalence (istat,stat)
external com_err_ (descriptors), attach (descriptors), detach
 lpdata=10
 data key,min/" ","<","=",">","min"/
print,"   Multics LP Program"
print, " "
 call attach (istat)
 if (istat) 3535,3434,3535
 3535 call com_err_(istat,"lnprog","error in attach")
 stop
 3434 tol1=0.1e-5
 tol2=0.1e-4
 tol4=10.0e-6
 tol5=1.0e+25
 10 read(lpdata,11,end=999) (idr(i),i=1,10)
 11 format(10a4)
 read(lpdata,950)m,n,ll
 print 12,(idr(i), i=1,10)
 12 format(/10x,10a4/)
 read(lpdata,950)(cj(j),j=1,n)
 print 13,m,n
 13 format(10x,i2," ROWS x ",i2," COLS"/)
 do 1111 i=1,m
 read(lpdata,950)(a(i,j),j=1,n),isen(i),b(i),ci(i)
 1111 continue
 fmin=1.0
 if(ll .eq.min)fmin=-1.0
 40 do 70 j=1,n
 cj(j)=fmin*cj(j)
 70 idc(j)=10*j
 do 100 i=1,m
 ci(i)=fmin*ci(i); ll=isen(i); kr=10*i+1001
 if(ll.ne.key(4))go to 90
 80 kr=kr+2; b(i)=-b(i)
 do 85 j=1,n
 85 a(i,j)=-a(i,j)
 90 if(ll.ne.key(3))go to 100
 95 ci(i)=-tol5; kr=kr+1
 100 idr(i)=kr
 zm=0.0; infs=0
 do 130 i=1,m
 if(zm.le.b(i))go to 130

 zm=b(i)
 130 continue
 if(zm)140,180,180
 140 zm=abs(zm); n=n+1
 do 160 i=1,m
 if(b(i) .lt.0.)go to 150
 a(i,n)=0.0
 go to 160
 150 b(i)=b(i)+zm; a(i,n)=1.0
 160 continue
 m=m+1; idr(m)=10*m+2002; idc(n)=9990
 do 170 j=1,n
 170 a(m,j)=0.0
 a(m,n)=1.0; b(m)=zm; ci(m)=-tol5; cj(n)=0.0
 180 print,"1=print iteration log, 2=suppress"
 read,ksw
 if(ksw-1)190,192,190
 190 ksw=2
 go to 200
 192 print,"   ITR  V IN V OUT      OBJ FUNCT"
 itr=0;k=0;l=0
 200 zm=0.0
 do 210 i=1,m
 210 zm=zm+b(i)*ci(i)
 go to (220,230),ksw
 220 print 900,itr,l/10,k/10,fmin*zm
 itr=itr+1
 230 dj=-tol1; kc=0
 do 250 j=1,n
 if(idc(j) .le.0)go to 250

 t=zj(m)
 if(dj .le. t) go to 250
 dj=t; kc=j
 250 continue
 if(kc.le.0)go to 400
 ark=0.0;brk=1.0e35;kr=0;i=m
 do 290 l=1,m
 t=a(i,kc);if(t.le.0.)go to 290
 t1=b(i)/t
 if(brk-t1)290,270,280
 270 if(t-ark)290,290,282
 280 brk=t1
 282 ark=t
 kr=i
 290 i=i-1
 if(kr .le.0)go to 292
 if(tol2.le.ark)go to 295
 infs=1
 292 idc(kc)=-idc(kc)
 go to 230
 295 a(kr,kc)=0.0
 do 330 j=1,n
 if(j-kc)300,330,300
 300 if(a(kr,j))310,330,310
 310 zm=a(kr,j)/ark;a(kr,j)=zm
 do 320 i=1,m
 t=a(i,kc)*zm
 t1=a(i,j)-t
 320 a(i,j)=ezero(tol4)
 330 continue
 b(kr)=brk;if(brk.le.0.)go to 350
 do 340 i=1,m
 t=a(i,kc)*brk;t1=b(i)-t
 340 b(i)=ezero(tol4)
 350 a(kr,kc)=-1.0
 do 360 i=1,m
 360 a(i,kc)=-a(i,kc)/ark
 k=idr(kr);l=idc(kc);idr(kr)=l;idc(kc)=k
 t=ci(kr); ci(kr)=cj(kc); cj(kc)=t
 if(infs.le.0)go to 200
 do 380 j=1,n
 if(idc(j)) 370,380,380
 370 idc(j)=-idc(j)
 380 continue
 infs=0
 go to 200
 400 print 910,fmin*zm
 do 440 k=1,m
 nl=3000; i=0
 do 410 l=1,m
 if(nl.le.idr(l))go to 410
 nl=idr(l); i=l
 410 continue
 if(i.le.0)go to 440
 idr(i)=9990;kc=mod(nl,10)+1
 print 920,nl/10,key(kc),fmin*ci(i),b(i)
 440 continue
 print 930
 do 540 k=1,n
 nl=2000; j=0
 do 510 l=1,n
 if(nl.le.iabs(idc(l)))go to 510
 nl=idc(l); j=l
 510 continue
 if(j.le.0)go to 540
 idc(j)=9990;kc=iabs(mod(nl,10))+1
 if(cj(j)+tol5) 530,520,530
 520 cj(j)=0.0
 530 print 920,nl/10,key(kc),fmin*cj(j),zj(m)
 540 continue
  print, " "
 go to 10
 900 format(3i6,1pe15.5)
 910 format(/1x,"  OBJ. FUNCT.    =",1pe15.5//"   BAS VAR",
& 4x,"OBJ. COEFF.",5x,"RESULT")
 920 format(6x,i3,a1,2(1pe15.5))
 930 format(/1x,"N. BAS VAR",4x,"OBJ. COEFF.",5x,"RESULT")
 940 format(4h    ,15a4,/4x,i2," rows x ",i2," cols")
 950 format(v)
999  call detach
  stop
end
c effective zero subroutine
 function ezero(tol4)
 common t,t1
 if(t)100,130,100
 100 if(t1)110,120,110
 110 if(abs(t1/t)-tol4)120,130,130
 120 ezero=0.0;return
 130 ezero=t1;return
 end
c compute (zj-cj) subroutine
cSTART zj
 function zj(m)
 common t(2),j,a(31,51),ci(31),cj(51)
 zj=-cj(j)
 do 100 ne=1,m
 100 zj=zj+a(ne,j)*ci(ne)
 return
 end
 subroutine maxflow
c - maxflow mod. May 1975
 common cap(25,25),flo(25,25)
 dimension lab(25),eps(25),iscan(25)
 logical out;  integer data
 character fname*32
 double precision stat
integer istat
 equivalence (istat,stat)
external com_err_ (descriptors), attach (descriptors), detach
 data=10
 n=25
    print 21
21 format(/"          M u l t i c s   M A X - F L O W / M I N - C U T"/)
 call attach (istat)
 if (istat) 4545,4141,4545
 4545 call com_err_(stat,"maxflow","error in attach")
 stop
 4141 continue
    print 23
 23 format (/" Type 0 for trace if iterations, 1 for optimal solution only")
 read 1000, i;out=i.eq.0
 1000 format (v)
 read(10,3)lab
 3 format(25a2)
 print 4,lab
 4 format(/,20x,25a2,/)
c initialization
 do 101 i=1,n
 lab(i)=0;eps(i)=0;iscan(i)=0
 do 101 j=1,n
 cap(i,j)=-9999.
 101 flo(i,j)=0.
 nn=0
c read in arcs in following format:1=supersource node,
c n=supersink node;on each line
c  line#,start node,end node,capacity,flow
 1 format(v)
 102 read(data,1,end=104)i,j,a1,a2
 cap(i,j)=a1
 flo(i,j)=a2
 if(j.gt.nn)nn=j
 go to 102
 104 lab(1)=9999
 eps(1)=9999.0
 n=nn
 if(out)call moutput(1,nn)
 105 do 113 j=1,nn
c find a labeled, unscanned node.
 if(lab(j).eq.0) go to 114
 107 if(iscan(j).eq.1)goto 114
c scan it.
 do 112 i=1,nn
 if(cap(j,i)+999.) 200,112,200
 200 if(lab(i).ne.0) go to 112
 ps=cap(j,i)-flo(j,i)
 if(ps.gt.0.0) go to 1091
 if(cap(i,j)+9999.) 210,112,210
 210 if(flo(i,j).ge.0.0) go to 112
 lab(i)=-j
 if(flo(j,i).lt.eps(j)) go to 110
 eps(i)=eps(j)
 go to 111
 1091 lab(i)=j
 if(ps.lt.eps(j)) go to 1092
 eps(i)=eps(j)
 go to 111
 1092 eps(i)=ps
 go to 111
 110 eps(i)=flo(j,i)
 111 if(lab(nn).ne.0) go to 1145
 112 continue
 iscan(j)=1
 114 continue
 113 continue
c check to see if all nodes are scanned.
 do 1131 i=1,nn
 if(lab(i).eq.0) go to 1131
 if(iscan(i).eq.1) go to 1131
 go to 105
 1131 continue
 go to 990
c sink node labelled, increase flow
 1145 add=eps(nn)
 115 jk=lab(nn)
 if(jk.gt.0)go to 116
 flo(jk,nn)=flo(jk,nn)-add
 go to 117
 116 flo(jk,nn)=flo(jk,nn)+add
 nn=lab(nn)
 117 if(lab(nn).eq.9999) go to 1171
 go to 115
c back at the beginning, reinitialize
 1171 nn=n
c a feasible flow has been found, tell the world
 if(out)call moutput(3,n)
c reinitialize those things which need reinitialization
 do 118 i=2,nn
 lab(i)=0
 118 iscan(i)=0
 iscan(1)=0
 go to 105
 990 call moutput(2,n)
 end file 10
 call detach
 stop
 end
cSTART moutput
 subroutine moutput(k,n)
 common cap(25,25),flo(25,25)
 go to (15,20,30),k
 15 print 4
 go to 990
 20 print 5
 go to 990
 30 print 8
 990 flomax=0.0
 print 6
 do 992 i=1,n
 do 991 j=1,n
 if(cap(i,j)+9999.) 800,991,800
 800 print 7,i,j,cap(i,j),flo(i,j)
 if(j.eq.n) flomax=flomax+flo(i,j)
 991 continue
 992 continue
 print 10,flomax
 10 format (/20x," Max Flow is ",f10.5,//)
 4 format(/1x,"                   Status of network at start"/)
 5 format(/1x,"              Status of network at optimal solution"/)
 6 format(5x,"Starting Node   Ending Node  Capacity Flow   In Arc  ")
 7 format(1h ,3x,i10,3x,i10,5x,f10.2,3x,f10.2)
 8 format(/1x,"             Status of network at a feasible solution",/)
 return
 end
 subroutine mreg1
c     mreg1     Modified Feb 1975
c
       dimension x(90,10),xsum(10),y(90),xmean(10),xpx(25,11),
&      label(25),xpy(10),xout(10),ll(10)
 character fnam*16
 integer fname,istat
 double precision stat
 equivalence (istat,stat)
external com_err_ (descriptors), attach (descriptors), detach
      data iyes,inos/3hyes,2hno/
     print," "
     print,"                    M R E G 1"
     print," "
    1 print, "Do you want instructions? (Type yes or no)"
       read 100,iy


 fname=10
   100 format(a3)
       if(iyes-iy) 193,130,193
   130 print 132
 print 133
 print 134
 print 135
 print 136
 132 format(///"Multiple regression is an attempt to curve fit observed"/
&"data by the model:",//,"          Y-YMEAN = A(1)(X(1)-XMEAN)+....+A(P)(X(P)-XPMEAN)",//)
 133 format("Where:"//"YMEAN = mean of observed data."/"XIMEAN = mean of each independent variable.")
 134 format("The A(I) are unknown coefficients"/
&"determined by the least squares process."/)
 135 format(" The data is entered in matrix form:"//
& "     Y1    X11    X12...X1M"/
& "     Y2    X21    X22...X2m"/
& "      .     .     .       ."/
& "     Yn    XN1    XN2...XNM
&"//)
 136 format( "Y1 is the observed data while X11,X12,X1M are the"/
& "corresponding independent variables."//"
&For the datafile option build  the file without line numbers"/
&"but with the data in the same sequence as requisted in the terminal option."//
&"Restrictions:"/
& "          M must be .LE. 10 and N must be .LE. 90."/
& "          Input is typed in with each element containing a decimal point"/
& "          and commas seperating each element."//
&"Now you try it.
&"//)
   193 print,"Is data read from a file? (Type yes or no)"
       read 100,iy
       if(iyes-iy) 200,358,200
   200 print,"Number of observations?"
       read,n
       print,"Number of variates?"
       read,np
   141 format(/"Enter the data in the following way,"/5x"Y1,X11,",
&      "X12,...,X1M"/5X"Y2,X21,X22,...,X2M"/6X".  .   ."7X"."/5X,
&      "YN,XN1,XN2,...,XNM"//)
       print 141
     print," "
       do 290 i=1,n
   290 read,y(i),(x(i,j),j=1,np)
   299 print 301
   301 format(/"Are any of the above y(n),x(n,m) elements typed ",
&      "incorrectly?"/" (Type yes or no)")
   306 print,"Any corrections?"
   307 read 100,icy
       if(inos-icy) 315,350,315
   315 if(iyes-icy) 316,320,316
   316  print,"Illegal command, try again."
       go to 306
   320 print 321
   321 format(/"Correct element by typing n (row subscript),",
&      "comma,"/" m (column subscript),comma,y-value,",
&      "comma,"/" x-value,carriage return"/)
   325 read,i,j,ytemp,xtemp
       if(i-90) 335,335,336
   335 if(j-10) 342,342,336
   336 print 337
   337 format(/"Illegal subscript, m greater than 10, or n greater ",
&      "than 90, try again"/)
       go to 325
   342 y(i)=ytemp
       x(i,j)=xtemp
   343 print,"Any corrections?"
       read 100,icy
       if(inos-icy) 345,350,345
   345 if(iyes-icy) 346,325,346
   346 print,"illegal command, try again."
       go to 343
358 call attach (istat)
     if (istat) 3434,4545,3434
      3434 call com_err_(istat,"mreg1","error in attach")
 stop
 4545 continue
308 format (v)
  309 format(v)
       read(fname,309)n
       read(fname,309)np
     print," "
       do 348 i=1,n
       read(fname,308,err=348)y(i),(x(i,j),j=1,np)
   348 print 349,y(i),(x(i,j),j=1,np)
   349 format(6f11.6)
 call detach
       go to 299
   350 print 10
    10 format(//21x"MULTIPLE  REGRESSION  PROGRAM"//25x"YMEAN"9x,
&      "XMEAN")
       do 20 i=1,np
    20 xsum(i)=0.
       ysum=0.
       do 30 i=1,n
       ysum=ysum+y(i)
       do 30 k=1,np
    30 xsum(k)=xsum(k)+x(i,k)
       fn=n
       ymean=ysum/fn
       do 35 i=1,np
    35 xmean(i)=xsum(i)/fn
       print 15,ymean,(xmean(i),i=1,np)
    15 format(20x,1p1e14.6/(34x,1p1e14.6))
c      convert data
       do 45 i=1,n
       do 40 j=1,np
    40 x(i,j)=x(i,j)-xmean(j)
    45 y(i)=y(i)-ymean
     50 print 52
     52 format(/,"Number of variates to consider?")
       read,nvr
      if(nvr.le.0)stop
       print,"Which ones?"
       read,(ll(i),i=1,nvr)
       np1=nvr+1
c      clear xpx array
       do 60 i=1,25
       do 60 j=1,11
    60 xpx(i,j)=0.
c      form cross products and sum
       do 72 ii=1,nvr
       i=ll(ii)
       do 72 jj=ii,nvr
       j=ll(jj)
       do 70 k=1,n
    70 xpx(ii,jj)=xpx(ii,jj)+x(k,i)*x(k,j)
    72 xpx(jj,ii)=xpx(ii,jj)
c      form right side of normal equations
       do 90 ii=1,nvr
       i=ll(ii)
       do 80 k=1,n
    80 xpx(ii,np1)=xpx(ii,np1)+y(k)*x(k,i)
    90 xpy(ii)=xpx(ii,np1)
c      solve normal equations
       call mmtinv(xpx,nvr,np1,label)
       sst=0.
       ssr=0.
       do 101 i=1,n
   101 sst=y(i)*y(i)+sst
       do 102 i=1,nvr
   102 ssr=ssr+xpx(i,np1)*xpy(i)
       sse=sst-ssr
       dfr=nvr
       dfe=n-nvr-1
       dft=n-1
       vsr=ssr/dfr
       vse=sse/dfe
       ftest=vsr/vse
       print,"Want to see predicted values? (Type yes or no)"
  4005 format(a3)
       read 4005,iq
       if(iyes-iq) 400,105,400
c      print read and predicted data
   105 print 2000
       do 131 i=1,n
       out=0.
       do 110 kk=1,nvr
       k=ll(kk)
   110 out=out+x(i,k)*xpx(kk,np1)
       out=out+ymean
       yout=y(i)+ymean
   131 print 3000,out,yout
  2000 format(//10x,36h              CALCULATED    OBSERVED/)
  3000 format(20x,1p2e14.6)
   400 nfe=dfe
       nm1=dft
       print 459
   459 format(///)
       print 460
   460 format(5x,65(1h.))
       print 461
   461 format(5x,1h.11x,25h.DEGREE OF FREE.  SUM OF ,
&      "SQUARES.  VARIANCE ESTIMATE.")
       print 500,nvr,ssr,vsr
       print 501,nfe,sse,vse
   500 format(5x,13h. Regression.i9,6h     .,1p1e14.6,3h  .1p1e16.6,
&      4h   .)
   501 format(5x,13h. Remainder .i9,6h     .,
&      1p1e14.6,3h  .1p1e16.6,4h   .)
       print 502,nm1,sst
   502 format(5x,13h. Total     .i9,6h     .1p1e14.6,3h  .18x,2h .)
       print 460
       print 459
       print 465
   465 format(5x,26hLeast Square Coefficients 4x,12h Coefficient,
&      " t Confidence Band")
       do 510 i=1,nvr
       ser=sqrt(xpx(i,i)*vse)
   510 print 466,xpx(i,np1),ser
   466 format(13x,1p1e14.6,17x,1p1e14.6)
       print 459
       print 468,nvr,nfe,ftest
   468 format(5x,8hf Ratio(i3,1h,i3,22h  degrees of freedom)=,
&     1p1e14.6)
       print 459
       print 4003
 4003 format(10x,"Variance-Covariance-Matrix")
       do 469 i=1,nvr
   469 print 505,(xpx(i,j),j=i,nvr)
   505 format(14x,1p1e14.6)
       go to 50
       end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c      matrix inverse subroutine
c
cSTART gecos.Tue1023.62

       subroutine mmtinv(a,nrarg,ncarg,label)
       dimension a(25,999),label(999)
       nr=nrarg
       nc=ncarg
       do 11 j1=1,nr
   11  label(j1)=j1
       do 21 j1=1,nr
       temp=0.0
       do 13 j2=j1,nr
       if(abs(a(j2,j1))-temp) 13,12,12
   12  temp=abs(a(j2,j1))
       ibig=j2
   13  continue
       if(ibig-j1) 27,16,14
   14  do 15 j2=1,nc
       temp=a(j1,j2)
       a(j1,j2)=a(ibig,j2)
   15  a(ibig,j2)=temp
       i=label(j1)
       label(j1)=label(ibig)
       label(ibig)=i
   16  temp=a(j1,j1)
       a(j1,j1)=1.0
       do 17 j2=1,nc
   17  a(j1,j2)=a(j1,j2)/temp
       do 20 j2=1,nr
       if(j2-j1) 18,20,18
   18  temp=a(j2,j1)
       a(j2,j1)=0.0
       do 19 j3=1,nc
   19  a(j2,j3)=a(j2,j3)-temp*a(j1,j3)
   20  continue
   21  continue
       n1=nr-1
       do 26 j1=1,n1
       do 23 j2=j1,nr
       if(label(j2)-j1) 23,22,23
   22  if(j2-j1) 27,26,24
   23  continue
   24  do 25 j3=1,nr
       temp=a(j3,j1)
       a(j3,j1)=a(j3,j2)
   25  a(j3,j2)=temp
       label(j2)=label(j1)
   26  continue
   27  return
       end
       subroutine orpol1
c     orpol Modified Feb 1975
c ********** orpol .. contains *orpol1* and *orpol2* ************
c *****main program segment orpol1
       common y(100),x(100),rho(100),poly(35),alpha(30),
       & beta(30),coef(30),p0(100),p1(100),p2(100),
        &c(30),ll(30)
 character fname*16
 character yes*4,iy*4
 integer istat,orpol
 double precision stat
 equivalence (istat,stat)
external com_err_ (descriptors), attach (descriptors), detach
 orpol=10
       yes="yes"
       print 2999
  2999 format(///" ORTHOGONAL POLYNOMIAL CURVE-FITTING")
  3000 format(//" Number of points, max degree?")
  3010 format(//" Is data to be read from a file? (Type 'yes','no' or
      & 'stop')")
  3001 format(" Type in dependent data?")
  3002 format(" Type in independent data?")
  3003 format(" Type in weights?")
  3004 format(" Dependent Data Mean =",1pe14.6)
  3005 format(//" DEGREE",9x,"ALPHA",9x,"BETA",9x,"COEFF",
&      11x,"SSR")
  3007 format(v)
  3008 format(i4,5x,4(1pe14.6))
     1 print 3010
       read 3007,iy
       if(iy.eq."stop")stop
       if(iy.eq.yes)go to 150
     3 print 3000
       read 3007,n,max
       mp1=max+1
       print 3001
       read 3007,(y(i),i=1,n)
       print 3002
       read 3007,(x(i),i=1,n)
       print 3003
       read 3007,(rho(i),i=1,n)
     4 fn=0.
       ta=0.
       ba=0.
       ysum=0.
       do 5 i=1,n
       fn=fn+rho(i)
     5 ysum=ysum+y(i)*rho(i)
       ymean=ysum/fn
       do 10 i=1,n
       ta=ta+rho(i)*x(i)
    10 ba=ba+rho(i)
       alpha(1)=ta/ba
       beta(1)=0.
       k=1
       ss=0.
       xpy=0.
       do 20 i=1,n
       p0(i)=1.
       p1(i)=p0(i)*x(i)-alpha(1)*p0(i)
       ss=ss+p1(i)*p1(i)*rho(i)
       y(i)=y(i)-ymean
    20 xpy=xpy+p1(i)*y(i)*rho(i)
       coef(1)=xpy/ss
       ssr=coef(1)*xpy
       print 3004,ymean
       print 3005
       print 3008,k,alpha(1),beta(1),coef(1),ssr
       if(max-1)25,61,25
    25 do 60 k=2,max
       bb=ba
       ba=0.
       ta=0.
       tb=0.
       do 30 i=1,n
       ta=ta+rho(i)*p1(i)*p1(i)*x(i)
       ba=ba+rho(i)*p1(i)*p1(i)
    30 tb=tb+rho(i)*x(i)*p1(i)*p0(i)
       alpha(k)=ta/ba
       beta(k)=tb/bb
       ss=0.
       xpy=0.
       do 40 i=1,n
       p2(i)=x(i)*p1(i)-alpha(k)*p1(i)-beta(k)*p0(i)
       ss=ss+p2(i)*p2(i)*rho(i)
    40 xpy=xpy+p2(i)*y(i)*rho(i)
       coef(k)=xpy/ss
       ssr=coef(k)*xpy
       print 3008,k,alpha(k),beta(k),coef(k),ssr
       do 50 j=1,n
       p0(j)=p1(j)
    50 p1(j)=p2(j)
    60 continue
    61 continue
       call orpol2(max,mp1,yes,ymean)
       go to 1
 150  call attach (istat)
  if (istat) 3535,4545,3535
 3535 call com_err_(istat,"orpol","error in attaching")
 stop
 4545 continue
       read(orpol,3007)n,max
       read(orpol,3007)(y(i),i=1,n)
       read(orpol,3007)(x(i),i=1,n)
       read(orpol,3007)(rho(i),i=1,n)
 call detach
       mp1=max+1
       go to 4
       end
c *****orpol2 segment******
cSTART orpol2
       subroutine orpol2(max,mp1,yes,ymean)
       common y(100),x(100),rho(100),poly(35),alpha(30),
       & beta(30),coef(30),p0(100),p1(100),p2(100),
        &c(30),ll(30)
 character yes*4,iy*4
  1000 format(//" Number and polynomials?")
  1010 format (/" Number, input points?")
  1020 format(//7x,"INPUT",10x,"OR OUTPUT",6x,"REG OUTPUT")
  1030 format (//" Finished? (yes or no)")
  1040 format(//" Regular Polynomial in Decreasing Degree")
  1050 format(//" Number, set of values")
  1060 format(v)
     1 continue
       print 1000
       read 1060,m,(ll(i),i=1,m)
       do 20 i=1,max
    20 c(i)=0.
       do 25 i=1,m
       j=ll(i)
    25 c(j)=coef(j)
       do 30 i=1,mp1
       p0(i)=0.
       p1(i)=0.
       p2(i)=0.
    30 rho(i)=0.
       p0(i)=1.
       p1(i)=1.
       p1(2)=-alpha(1)
       k=2
       rho(1)=coef(1)
        rho(2)=-alpha(1)*coef(1)
       if (max-1)40,120,40
    40 continue
       kp1=k+1
       km1=k-1

    50 do 60 i=1,k
    60 p2(i)=p1(i)
       do 70 i=1,k
    70 p2(i+1)=p2(i+1)-alpha(k)*p1(i)
       do 80 i=1,km1
    80 p2(i+2)=p2(i+2)-beta(k)*p0(i)
       lk =k+1
       do 85 i=1,k
       lk=lk-1
    85 rho(lk+1)=rho(lk)
       rho(1)=0.
       do 90 i=1,kp1
    90 rho(i)=rho(i)+p2(i)*c(k)
       if (max-k) 100,120,100
   100 do 105 i=1,k
   105 p0(i)=p1(i)
       do 106 i=1,kp1
   106 p1(i)=p2(i)
       do 107 i=1,max
   107 p2(i)=0.
       k=k+1
       go to 40
   120 print 1040
       rho(mp1)=rho(mp1)+ymean
       print 1090, (rho(i),i=1,mp1)
   1090 format((e16.8))
       print 1050
       read 1060,mt,(x(i),i=1,mt)
       print 1020
       do 220 i=1,mt
       ans=rho(1)
       do 130 l=2,mp1
   130 ans=ans*x(i)+rho(l)
       poly(1)=1.
       poly(2)=x(i)-alpha(1)
       if(max-1) 180,195,180
   180 do 190 k1=3,mp1
   190 poly(k1)=x(i)*poly(k1-1)-alpha(k1-1)*poly(k1-1)-beta(k1-1)*poly(k1-2)
   195 pred=ymean
       do 200 i1=1,m
       j=ll(i1)
   200 pred=pred+coef(j)*poly(j+1)
   220 print 1022, x(i),pred,ans
  1022 format(e16.8,e16.8,e16.8)
       print 1030
        read 1060, iy
       if(iy.ne.yes) go to 1
   250 return
       end
 subroutine shortest
c     Shortest  -- modified Feb 1975
 common a(25,25),d(25,25),nicid(25,25)
 dimension icon(25)
 character fname*16
 integer istat
 double precision stat
 equivalence (istat,stat)
external com_err_ (descriptors), attach (descriptors), detach
     print 5
  5  format(//,20x,"S H O R T E S T"/)
 iin=10
 print,"1=TRACE ITER., 0=ANSWERS ONLY, -1=MIN SPAN TREE ONLY"
 read,iout
 call attach (istat)
 if (istat) 3535,4545,3535
 3535 call com_err_(istat,"shortest","error in attach")
 stop
 4545 read(iin,3)icon
 3 format(25a2)
 print 4,icon
  4 format(//,2x,25a2/)
 do 1 i=1,25
 do 1 j=1,25
 nicid(i,j)=0
 1 d(i,j)=0.0
 nn=0
 100 read(iin,2,end=1000)i,j,a1
 2 format(v)
 nicid(i,j)=1
 d(i,j)=a1
 if(j.gt.nn)nn=j
 go to 100
 1000 if(-1.eq.iout)go to 1060
 do 101 i=1,25
 101 a(1,i)=0.0
 do 102 i=2,nn
 if(nicid(1,i).ne.0) go to 1012
 a(i,1)=9999.0
 go to 102
 1012 a(i,1)=d(1,i)
 102 continue
 k=1
 if(iout.eq.1)call soutput(1,nn,3,k)
 103 k=k+1
 ifin=1
 do 105 i=2,nn
 xlst=a(i,k-1)
 valu=9999.0
 do 104 j=1,nn
 kk=nicid(i,j)
 if(i.gt.j)kk=nicid(j,i)
 if(kk.eq.0) go to 104
 r=d(j,i)
 if(j.gt.i) r=d(i,j)
 hold=a(j,k-1)+r
 if(hold.ge.valu) go to 104
 valu=hold
 104 continue
 a(i,k)=valu
 if(valu-xlst) 106,105,106
 106  ifin = 0
 105 continue
 if(iout.eq.1)call soutput(1,nn,3,k)
 if(ifin)1051,103,1051
 1051 call soutput(2,nn,3,k)
c minimal spanning tree algorithm

 1060 m=1
 icon(1)=1
 1071 hold=999999.0
 do 108 i=1,m
 k=icon(i)
 do 108 j=1,nn
 do 1072 l=1,m
 if(j.eq.icon(l)) go to 108
 1072 continue
 if(k-j)1073,108,1074
 1073 if(nicid(k,j).eq.0) go to 108
 if(d(k,j).gt.hold) go to 108
 hold=d(k,j)
 go to 1075
 1074 if(nicid(j,k).eq.0) go to 108
 if(d(j,k).gt.hold) go to 108
 hold=d(j,k)
 1075 io=k
 jo=j
 108 continue
 m=m+1
 if(m.gt.nn) go to 1081
 if(io.gt.jo) go to 1082
 nicid(io,jo)=2
 1082 nicid(jo,io)=2
 icon(m)=jo
 go to 1071
 1081 call soutput(0,nn,4,0)
 call detach
 stop
 end
cSTART soutput
c  Modified Feb 1975
 subroutine soutput(k,nn,iwhr,kk)
 common a(25,25),d(25,25),nicid(25,25)
 if(iwhr .eq.4)goto 300
c soutput routine for shortest path algorithm
 go to (1041,1040),k
 1040 print 42,kk
 go to 1042
 1041 print 41,kk
 1042 print 34
 j=1
 do 1043 i=2,nn
 1043 print 7,j,i,a(i,kk)
 7 format(/3x,i10,3x,i10,f10.2,3x,f10.2)
 go to(1045,1044),k
 1044 print 43
 43 format(/4x,"THE SHORTEST PATHS FROM ARC ONE ARE:")
 print 44
 44 format(/3x,13hSTARTING NODE,2x,11hENDING NODE,2x,10hARC LENGTH)
 do 108 i=1,nn
 do 107 j=1,nn
 kx=nicid(i,j)
 if(i.gt.j)kx=nicid(j,i)
 if(kx.eq.0)go to 107
 r=d(j,i)
 if(j.gt.i) r=d(i,j)
 if(abs(a(i,kk)-a(j,kk))-r) 107,105,107
 105 if(i.lt.j)print 7,i,j,d(i,j)
 107 continue
 108 continue
 1045 return
 34 format(/3x,13hSTARTING NODE,2x,11hENDING NODE,2x,8hDISTANCE)
 41 format(/3x,32h STATUS OF NETWORK AT ITERATION ,i2)
 42 format (/"SHORTEST DISTANCE NODE 1 TO ALL OTHERS, ITERATION ",i2)
 431 format(/"THE MINIMAL SPANNING TREE CONSISTS OF THE FOLLOWING ARCS:")
 300 print 431;print 34
 do 302 i=1,nn
 do 302 j=1,nn
 if(i.le.j.and.nicid(i,j).eq.2)print 7,i,j,d(i,j)
 302 continue
 return
 end
 subroutine smlrp
c  smlrp - stepwise multiple linear regression Feb 1975
cSTART gecos.Fri1636.11
 common c(20,20),avg(20),sos(20)
 common vl(20),vh(20),obs(30),vlab(20),ph(20)
 common nr,np,ns,id,pi,nd,sgn,ks3,tol
 common kr,fe,fc,fr,dof,cl,see
 common fname
 external  com_err_(descriptors), detach, smlrpattach (descriptors)
 character fname*16
 character vlab*4,kx*4,nam*4,nal1*4,nal2*4,temp*4
 character kod*3(20)
 double precision stat
 dimension ktr(80),iv1(80),iv2(80),iv3(80),temp(20)
 real na3
 integer ph,vi(29),yes,no,pi,sct,dof,four
  integer istat
  equivalence(istat,stat)
 data bad,yes,no,sct,lbk/-1e12,89,78,65536,"  ["/

data kod/"ADD","SUB","MPY","DVD","RCP","MOV","EXP","TEN","LGN","LOG","SIN","COS","TNH","SQR","CON","END",4*"0"/
lrdata=10; four=11
c  initialize
 3939 continue
      print 120
 120 format(/10x" STEPWISE MULTIPLE LINEAR REGRESSION PROGRAM")
 print 1
 1 format (//)
 ks4=1
 ks6=1
 10 ks1=1
 ks2=1
 ks3=1
 ks5=1
 kr=0
 rewind four
 nal2=" "
c
c  read data file name
c
fname = " "
 call smlrpattach(istat, fname)
 if (istat) 4545,3535,4545
 4545 call com_err_(istat,"smlrp","error attaching file")
c
c  read no. of initial variables
c
 3535 read(lrdata,21)nv
 21 format (v)
 nr=nv
 print 901, nr
 if(nr.gt.20)go to 800
 na3=0
 nt=0
c
c  test for transformation code
c
 idum=0
 227 read(lrdata,21) nam,kdum
 idum=idum+1
 if(nam.eq."LABL" .or. nam.eq."LABE")go to 38
 if(nam.eq."TRNF" .or. nam.eq."TRAN")go to 229
 go to 46
 229 nr=kdum
 print 905, nr
 if ((nr .lt. 2) .or. (nr .gt. 20)) go to 815
 ks5=2
 ker=0
c
c  read and process transformation codes
c
 do 34 l=1,81
 20 na1=0
 na2=0
 na3=0
 na4=0
 read (lrdata,21) nal1,na2,na3,na4
 idum=idum+1
 do 25 i=1,16
 if (nal1 .eq. kod(i)) go to 30
 25 continue
 27 ker=1
 print 907, nal1,na2,na3,na4

 go to 20
 30 if (nal1 .eq. "END") go to 36
 if (nal1 .ne. "CON") go to 32
 obs(na2)=na3
 go to 20
 32 nt=nt+1
 ktr(nt)=i
 if ((na2 .lt. 1) .or. (na2 .gt. 30)) go to 27
 iv1(nt)=na2
 if ((na3 .lt. 1) .or. (na3 .gt. 30)) go to 27
 iv2(nt)=na3
 if ((i .lt. 5) .and. ((na4 .lt. 1) .or. (na4 .gt. 30)))
&  go to 27
 iv3(nt)=na4
 34 continue
 print 2
 2 format (39h *** no "end" following transformations)
 go to 540
 36 if (ker .gt. 0) go to 540
 nal2=" "
 go to 227
c
c  test for variable labels
c
 38 itv=kdum
 do 40 i=1,nr
 vlab(i)="    "
 40 continue
c
c  read and store labels
c
 42 ks7=3

 k=itv+1
 read (lrdata,21) (temp(i), i=1,k)
 idum=idum+1
 do 44 i=1,k
 if (temp(i) .eq. "ELAB") go to 46
 vlab(i)=temp(i)
 44 continue
 print 3
 3 format (33h *** no "endlbl" following labels)
 go to 540
c
c  initialize arrays
c
 46 do 48 i=1,nr
 avg(i)=0.0
 vl(i)=1.0e35
 vh(i)=-1.0e35
 48 continue
 ker=0
c
c  read raw observed data
c
 do 145 k=1,999
 do 50 i=1,nr
 50 obs(i)=1e12
 read (lrdata,21,end=150) (obs(i), i=1,nv)
 go to (100,55), ks5
c
c  transform raw data
c
 55 do 90 j=1,nt
 l=ktr(j)
 na1=iv1(j)
 na2=iv2(j)
 na3=iv3(j)
 go to (61,62,63,64,65,66,67,68,69,70,71,72,73,74),l
 61 obs(na3)=obs(na1)+obs(na2)
    go to 90
 62 obs(na3)=obs(na1)-obs(na2)
    go to 90
 63 obs(na3)=obs(na1)*obs(na2)
    go to 90
 64 obs(na3)=obs(na1)/obs(na2)
    go to 90
 65 obs(na2)=1.0/obs(na1)
    go to 90
 66 obs(na2)=obs(na1)
    go to 90
 67 obs(na2)=exp(obs(na1))
    go to 90
 68 obs(na2)=10.**(obs(na1))
    go to 90
 69 obs(na2)=alog(obs(na1))
    go to 90
 70 obs(na2)=0.4343*alog(obs(na1))
    go to 90
 71 obs(na2)=sin(obs(na1))
    go to 90
 72 obs(na2)=cos(obs(na1))
    go to 90
 73 obs(na2)=tanh(obs(na1))
    go to 90
 74 obs(na2)=sqrt(obs(na1))
 90 continue
 100 write (four,21) (obs(i), i=1,nr)
c
c  accumulate variable sums
c
 do 145 i=1,nr
 x=obs(i)
 if(x-1e12) 140,139,140
 139 ker=1
 print 903, i, k
 go to 145
 140 avg(i)=avg(i)+x
 if (x .lt. vl(i)) vl(i)=x
 if (x .gt. vh(i)) vh(i)=x
 145 continue
 150 nob=k-1
 print 902, nob
 if (ker .gt. 0) go to 540
c
c  read means print switch
c
 nal1="no"
 print 4
 4 format (/" MEANS?")
 read 1001,  nal1
if ((nal1 .eq. "no") .or. (nal1 .eq. "NO")) ks1=2
c
c  compute averages
c
 x=nob
 do 155 i=1,nr
 155 avg(i)=avg(i)/x
 160 rewind four
 do 165 j=1,nr
 do 165 i=1,nr
 165 c(i,j)=0.0
c
c  compute covariance matrix
c
 do 170 k=1,nob
 read (four,21) (obs(i), i=1,nr)
 do 170 i=1,nr
 x=obs(i)-avg(i)
 do 170 j=1,nr
 c(i,j)=c(i,j)+x*(obs(j)-avg(j))
 170 continue
 dof=nob-1
 ker=0
 go to (175,180), ks1
 175 print 904
c
c  compute std-deviations
c
 180 do 190 i=1,nr
 x=c(i,i)
 sos(i)=c(i,i)
 std=sqrt(x/dof)
 if (std .ge. 1.0e-5) go to 5
 ker=1
 ks1=1
 5 go to (185,190), ks1
c
c  print means & std-dev
c
 185 print 906, i,vlab(i),avg(i),std,vl(i),vh(i)
 190 continue
 if (ker .gt. 0) go to 805
c
c  read correlation matrix print switch
c
 go to (192,200), ks2
 192 nal1="no"
 print 6
 6 format(//" CORR. MATRIX")
 read 1001,  nal1
if ((nal1 .eq. "no") .or. (nal1 .eq. "NO")) ks2=2
 go to (195,200), ks2
 195 print 908, (lbk,j,j=1,nr)
c
c  compute correlation matrix
c
 200 do 215 i=1,nr
 x=sos(i)
 do 205 j=1,nr
 205 c(i,j)=c(i,j)/sqrt(x*sos(j))
 go to (210,215), ks2
c
c  print correlation matrix
c
 210 print 910, i,(c(i,j), j=1,i)
 215 continue
 ks1=2
 ks2=2
 220 do 222 i=1,nr
 ph(i)=-1
 vi(i)=0
 222 continue

 print 7
 7 format(/" ENTER NO. of 'X' VARS., THEN INDEX of 'Y' VAR.")
 print 8
 8 format (" FOLLOWED by INDICES of all 'X' VARS.")
c
c  select variables for regression analysis
c
 read 1020,  ni,id,(vi(k), k=1,ni)
 1020 format(v)
 k=ni+1
 if (k .gt. nr) go to 220
 if ((nob-k) .lt. 1) go to 810
 do 228 k=1,ni
 i=vi(k)
 l=iabs(i)
 if (l .gt. nr) go to 220
 if (i) 226,220,224
 224 ph(i)=0
 go to 228
 226 ph(l)=-2
 228 continue
 if ((id .lt. 1) .or. (id .gt. nr)) go to 220
 if (-1 .ne. ph(id)) go to 220
 ph(id)=99
 np=0
 ns=0
 230 fc=0.0
 fe=0.0
 cl=0.0
 x=0.0
 tol=1.0e-4
     print 126
 126 format(/" ENTER F OR CL")
 go to (231,232), ks6
 231 print 9
 9 format(" FOR: CRITICAL F-VALUE or CONFIDENCE LEVEL")
 ks6=2
c
c  read critical-f or confidence level
c
 232 read 1001,  nal1
 if ((nal1 .eq. "CL")  .or. (nal1 .eq. "cl")) go to 233
   print 127
 127 format(/" ENTER CRTICAL-F")
 read 1010,  fe
 go to 234
 233 print 125
 125 format(/" ENTER CONFIDENCE LEVEL")
 read 1010,  cl
 1010 format(v)
 if (cl .ge. 1.0) go to 233
 234 fc=frat(dof)
 ks3=1
 kr=kr+1
 nal1="no"
c
c  read stepwise printout switch
c
 print 128
 128 format (/" STEPWISE RESULTS?")
 read 1001, nal1
 call print2
if ((nal1 .eq. "no") .or. (nal1 .eq. "NO")) ks3 =2
 go to (240,235), ks3
 235 print 924
c
c  pivot on forced variables
c
 240 do 245 i=1,nr
 if (-2 .ne. ph(i)) go to 245
 if (c(i,i) .lt. tol) go to 245
 pi=i
 call pivot(3)
 245 continue
 if (np) 260,260,250
c
c  select variable for deletion
c
 250 call selpvt(2)
 if (pi .eq. 0) go to 260
 call pivot(2)
 go to 250
c
c  select variable for addition
c
 260 call selpvt(1)
 if (pi .eq. 0) go to 300
 call pivot(1)
 go to 250
 300 if (np .gt. 0) go to 310
 print 22
 22 format (29h *** no significant variables)
 go to 230
 310 go to (330,320), ks3
c
c  print regression results
c
 320 call print
 330 call print2
 nal1="no"
 print 11
 11 format (/" RESIDUALS?")
c
c  read residual printout switch
c
 read 1001, nal1
 1001 format (v)
 if ((nal1 .eq. "no").or.(nal1 .eq. "NO")) go to 500
 rewind four
c
c  compute & print residual summary
c
 print 918
 s=0.0
 do 410 k=1,nob
 read (four,21) (obs(i), i=1,nr)
 x=vl(id)
 do 400 j=1,nr

 l=ph(j)
 if ((l .lt. 1) .or. (l .gt. 2)) go to 400
 x=x+vl(j)*obs(j)
 400 continue
 dr=obs(id)-x
 s=s+dr*dr
 pd=100.0*dr/x
 print 920, k,obs(id),x,dr,pd,s
 410 continue
 500 print 12
 12 format (/" ENTER A 1,2,3 or 4")
 go to (510,520), ks4
 510 print 13
 13 format(" FOR: 1 - SELECT NEW DATA FILE",/6x,"2 - SELECT NEW VARIABLES",
&/6x,"3 - SELECT NEW 'F' or 'CL'",/6x,"4 - END OF RUN")
 ks4=2
c
c  read "loop" back switch
c
 520 read 2000,  k
 2000 format (v)
 go to (530,160,230,540), k
 530 rewind lrdata
  continue
 go to 10
 540 continue
rewind lrdata
rewind four
 call detach
 stop
 800 print 14
 14 format (38h *** too many variables, maximum is 20)
 go to 540
 805 print 15
 15 format (42h *** one or more std-dev = 0.0, check data)
 go to 540
 810 k=k+1
 print 950, ni,k
 go to 540
 815 print 16
 16 format (33h *** no. of tranf. vars. in error)
 go to 540
 901 format (/,i4," Initial Variables")
 902 format (i4," Observations")
 903 format(" *** Variable",i3,", Observation",i3,", is missing")
 904 format (/" VAR  LABEL",7x,"MEAN",7x,"STD-DEV",8x,"MIN",10x,"MAX"/)
 905 format (i4," Transformed Variables")
 906 format (i4,a8,2f12.4,1p2e13.4)
 907 format (/" *** TRANF. ERROR",a10,i5,f8.4,i5)
 908 format (/" CORRELATION MATRIX"//(4x,8(a5,i2,"]")))
 910 format (" [",i2,"]"/(4x,8f8.3))
 918 format(/"  OBS",6x,"Y-OBS",7x,"Y-CALC",7x,"ERROR",5x,"-ERR",6x,"C-ER^2"/)
 920 format (i4,3f12.4,f9.2,f12.4)
 924 format (/" STEP   VAR:LABEL    F-CRIT   DOF",4x,"R-SQ",9x,"SEE"/)
 950 format (/" *** MINIMUM NO. OF OBSERVATIONS FOR",i3,"VARIABLES IS",i3)
 end
c  compute f-value routine
c
 function frat(d)
 common c(20,20),avg(20),sos(20)
 common vl(20),vh(20),obs(30),vlab(20),ph(20)
 common nr,np,ns,id,pi,nd,sgn,ks3,tol
 common kr,fe,fc,fr,dof,cl,see
 integer ph,pi,dof,d
c
c  d - degrees of freedom
c  cl - confidence level
c  fe - entered f-value
c
 if(cl) 100,100,200
 100 frat=fe
 return
 200 if(fc) 300,300,400
 300 p=(1.0-cl)/2.0
 t=sqrt(alog(1.0/p**2))
 x2=t*(t*(0.001308*t+0.189269)+1.432788)+1.0
 x=t-(t*(0.010328*t+0.802853)+2.515517)/x2
 x2=x*x
 c1=(x2+1.)/4.
 c2=(x2*(5.*x2+16.)+3.)/96.
 c3=(x2*(x2*(3.*x2+19.)+17.)-15.)/384.
 c4=(x2*(x2*(x2*(79.*x2+776.)+1482.)-1920.)-945.)/92160.
 400 t=x*((((c4/d+c3)/d+c2)/d+c1)/d+1.0)
 frat=t*t
 return
end

c  matrix pivot routine
c
 subroutine pivot(l)
 common c(20,20),avg(20),sos(20)
 common vl(20),vh(20),obs(30),vlab(20),ph(20)
 common nr,np,ns,id,pi,nd,sgn,ks3,tol
 common kr,fe,fc,fr,dof,cl,see
 character vlab*4,kx*4
 integer ph,pi,dof
c
c l=1, forward pivot
c l=2, backward pivot
c l=3, forced pivot
c pi - pivot index
c
 ks=l
 k=pi
 ark=c(k,k)
 c(k,k)=0.0
 do 300 j=1,nr
 if(ph(j)) 300,100,100
 100 if(c(k,j)) 150,300,150
 150 t=c(k,j)/ark
 c(k,j)=c(k,j)/ark
 do 200 i=1,nr
 200 c(i,j)=c(i,j)-c(i,k)*t
 300 continue
 c(k,k)=-1.0
 do 400 i=1,nr
 c(i,k)=-c(i,k)/ark
 400 continue
 go to (410,420,405),ks
 405 kx="   *"
 ph(k)=2
 go to 415
 410 kx="   +"
 ph(k)=1
 415 np=np+1
 dof=dof-1.0
 go to 430
 420 kx="   -"
 ph(k)=0

 np=np-1
 dof=dof+1.0
 430 ns=ns+1
 fc=frat(dof)
 dr=c(id,id)
 r2=1.0-dr
 see=sqrt(dr*sos(id)/dof)
 go to (440,450),ks3
 440 print 900
 450 print 910,ns,kx,k,vlab(k),fc,dof,r2,see
 go to (460,500),ks3
 460 call print
 500 return
 900 format (//" STEP   VAR:LABEL    F-CRIT   DOF",6x,"R-SQ",9x,"SEE"/)
 910 format (i4,a4,i2,":",a6,f9.2,i6,f10.4,f12.4)
 end
c  print regression results routine
c
 subroutine print
 common c(20,20),avg(20),sos(20)
 common vl(20),vh(20),obs(30),vlab(20),ph(20)
 common nr,np,ns,id,pi,nd,sgn,ks3,tol
 common kr,fe,fc,fr,dof,cl,see
 character vlab*4,kx*4
 integer ph,pi,dof
 print 900
 s=0.0
 do 100 i=1,nr
 l=ph(i)
c
 if((l .lt. 1) .or. (l .gt. 2)) go to 100
 x=c(i,id)*sqrt(sos(id)/sos(i))
 vl(i)=c(i,id)*sqrt(sos(id)/sos(i))
 serc=see*sqrt(c(i,i)/sos(i))
 fi=x/serc
 fi=fi*fi
 if((serc .lt. 1.0e-8) .and. (serc .ge. 0.0)) fi=9999.99
 print 910,i,vlab(i),x,serc,fi,c(i,id)
 s=s+x*avg(i)
 100 continue
 vl(id)=avg(id)-s
 print 920, vl(id)
 return
 900 format(//"  VAR  LABEL",6x,"COEFFICIENT",5x,"STD-ERR",5x,"F-RATIO",6x,"BETA-WT"/)
 910  format(i4,a8,f16.6,f13.4,f11.2,f13.4)
 920 format (" CONS",8x,f16.6)
 end
c print summary line routine
c
 subroutine print2
 common c(20,20),avg(20),sos(20)
 common vl(20),vh(20),obs(30),vlab(20),ph(20)
 common nr,np,ns,id,pi,nd,sgn,ks3,tol
 common kr,fe,fc,fr,dof,cl,see
 common fname
 integer dof,ph,pi
 character fname*16
 lrdata=10
 x=sos(id)*c(id,id)
 print 700,fname,kr,id,vlab(id),cl,dof,x
 700 format(//2x,a16,":",i2,"  DEP-VAR=",i2,":",a6,"  CL=",f4.2,"  DOF=",i3,"  RSS=",f13.4)
 return
 end
c  select pivot variable routine
c
 subroutine selpvt(n)
 common c(20,20),avg(20),sos(20)
 common vl(20),vh(20),obs(30),vlab(20),ph(20)
 common nr,np,ns,id,pi,nd,sgn,ks3,tol
 common kr,fe,fc,fr,dof,cl,see
 integer ph,pi,dof
c
c  n=1, select variable for forward addition
c  n=2, select variable for backward deletion
c

 10 pi=0
 go to (20,30),n
 20 l=0
 fr=frat(dof-1)
 go to 100
 30 l=1
 fr=fc
 100 do 350 i=1,nr
 if(ph(i) .ne. l) go to 350
 if(c(i,i) .lt. tol) go to (350,400),n
 dr=c(i,id)*c(i,id)/c(i,i)
 go to (200,250),n
 200 x=c(id,id)-dr
 ft=dr*(dof-1.0)/x
 if((x .lt. 1.0e-8) .and. (x .ge. 0.0)) ft=9999.99
 if(ft .lt. fr) go to 350
 go to 300
 250 ft=dr*dof/c(id,id)
 x=c(id,id)
 if((x .lt. 1.0e-8) .and. (x .ge. 0.0)) ft=9999.99
 if(ft .ge. fr) go to 350
 300 fr=ft
 pi=i
 350 continue
 return
 400 pi=i
 call pivot(2)
 ph(i)=-1
 go to 10
end
      subroutine stat
c     stat
c     Modified Feb 1975
       dimension x(1000),y(1000)
 character key*4
 character fname*16
 double precision xstati
  integer istati,stati
 equivalence(istati,xstati)
external com_err_ (descriptors), attach (descriptors), detach
  1     print 9000
        read 1000, key
     if((key .eq. "yes") .or. (key .eq. "YES")) go to 2
     if((key .eq. "no") .or. (key .eq. "NO")) go to 3
     go to 80
  2    stati=10
 call attach (istati)
 if (istati) 3535,4545,3535
 3535 call com_err_(istati,"stati","error in attach")
stop
 4545 read(stati,1000)n,idelt,(x(k),k=1,n)
  1000 format(v)
 call detach
        go to 4
     3 print 9001
       read 1000, n,idelt
       print 9002
       read 1000, (x(k),k=1,n)
     4 i = 0
       sum=0.
       svr=0.
       sd2=0.
       sd3=0.
       sd4=0.
       if (n-1) 80,80,6
     6 m = n-1
       do 8 k=1,n
       y(k) = x(k)
     8 sum = sum+x(k)
       do 24 k=1,m
       i=k+1
       do 24 j=i,n
       if (y(k)-y(j)) 24,24,12
    12 str = y(k)
       y(k) = y(j)
       y(j) = str
    24 continue
       avg = sum/n
       do 28 k=1,n
       d1 = (x(k)-avg)
       d2 = d1*d1
       d3 = d1*d2
       d4 = d2*d2
       sd2 = sd2+d2
       sd3 = sd3+d3
    28 sd4 = sd4+d4
       uvr = sd2/n
       s=sqrt(uvr)
       avr = uvr*n/(n-1)
       asd=sqrt(avr)
       xn=n
       skw=sqrt(xn)*sd3/sd2**1.5
       xkt = n*sd4/sd2**2
       l = n/2+1
       xmd = y(l)
       rng = y(n)-y(1)
        z1=s/sqrt(float(n-1))
       uclv = n*s**2
    40 print 9004,n,avg,xmd,rng,uvr,s
       print 9005,avr,asd,skw,xkt
       print 9006,avg,z1,avg,z1
       print 9007,uclv,uclv
       z2 = 100./(n+1)
       z3=(-8.)*s+avg
       print 9008
       do 60 k=1,n,idelt
       cpr = k*z2
        cnm=anpf(z3,y(k),avg,s)
       dif = cnm*100.-cpr
    60 print 9012,k,x(k),y(k),cpr,cnm,dif
       print 9016
       go to 1
    80 stop
  9000 format(///" Is data to be read from a file? (Type yes, no or stop)")
  9001 format(/" Type number of data points (n), and the print increment ")
  9002 format(/" Type in n data points"/)
  9004 format(//15x"DESCRIPTIVE STATISTICS FOR A SINGLE"
&        ," SAMPLE"///" Number of observations ",45(1h.),i4,//,
&           " Mean = xbr = sum(x(i))/n ",34(1h.),1pe13.6//" Median ",
&        52(1h.)1pe13.6//" Range ",53(1h.)1pe13.6//" Unadjusted ",
&       "(biased) variance ="/2x,"sum((x(i)-xbr)**2)/n = uvr ",
&           31(1h.),1pe13.6//" Unadjusted standard deviation = ",
&        "sqrt(uvr) = s ",13(1h.),1pe13.6/)
  9005 format(" Adjusted (unbiased) variance = uvr*n/(n-1) = avr ",
&        10(1h.)1pe13.6//" Adjusted standard deviation = sqrt(avr) "
&        19(1h.)1pe13.6//" Skewness = sum((x(i)-xbr)**3)/(n*",
&         "uvr**1.5) ",16(1h.),1pe13.6//" Kurtosis = ",
&        "sum((x(i)-xbr)**4)/(n*uvr**2) ",18(1h.),1pe13.6//)
  9006 format(" Upper confidence limit on mean ="/9x,3hxbr,8x,
&        14h+ t(n-1 d.f.)*2x,11hs/sqrt(n-1)/3x,1pe14.6,17x,1pe14.6//,
&        " Lower confidence limit on mean ="/9x,3hxbr,8x,9h- t(n-1 d,
&        5h.f.)*,2x,11hs/sqrt(n-1)/3x,1pe14.6,17x,1pe14.6//," Upper confidence limit of variance =")
        9007 format(2x,"n*s**2 divided by chi-square",
&        " (right"/2x,"tail critical region) with n-1 d.f.,  ",
&        "n*s**2 ",12(1h.),1pe13.6//" Lower confidence limit on ",
&        "variance =",/2x,"n*s**2 divided by chi-square (left"/2x,
&        "tail critical region) with n-1 d.f.,  n*s**2 ",
&        12(1h.),1pe13.6)
  9008 format(/16x,40(1h*)/,
&        /,23x,7hORDERED,7x,10hCUMULATIVE,3x,10hCUMULATIVE/,
&        10x,4hDATA,10x,4hDATA,9x,10hPERCENTAGE,5x,6hNORMAL,4x,
&        10hDIFFERENCE/)
  9012 format(i4,1p2e14.6,0pf13.3,2pf13.3,0pf14.4)
  9016 format(//)
      end


		    bull_copyright_notice.txt       08/30/05  1008.4r   08/30/05  1007.3    00020025

                                          -----------------------------------------------------------


Historical Background

This edition of the Multics software materials and documentation is provided and donated
to Massachusetts Institute of Technology by Group Bull including Bull HN Information Systems Inc. 
as a contribution to computer science knowledge.  
This donation is made also to give evidence of the common contributions of Massachusetts Institute of Technology,
Bell Laboratories, General Electric, Honeywell Information Systems Inc., Honeywell Bull Inc., Groupe Bull
and Bull HN Information Systems Inc. to the development of this operating system. 
Multics development was initiated by Massachusetts Institute of Technology Project MAC (1963-1970),
renamed the MIT Laboratory for Computer Science and Artificial Intelligence in the mid 1970s, under the leadership
of Professor Fernando Jose Corbato.Users consider that Multics provided the best software architecture for 
managing computer hardware properly and for executing programs. Many subsequent operating systems
incorporated Multics principles.
Multics was distributed in 1975 to 2000 by Group Bull in Europe , and in the U.S. by Bull HN Information Systems Inc., 
as successor in interest by change in name only to Honeywell Bull Inc. and Honeywell Information Systems Inc. .

                                          -----------------------------------------------------------

Permission to use, copy, modify, and distribute these programs and their documentation for any purpose and without
fee is hereby granted,provided that the below copyright notice and historical background appear in all copies
and that both the copyright notice and historical background and this permission notice appear in supporting
documentation, and that the names of MIT, HIS, Bull or Bull HN not be used in advertising or publicity pertaining
to distribution of the programs without specific prior written permission.
    Copyright 1972 by Massachusetts Institute of Technology and Honeywell Information Systems Inc.
    Copyright 2006 by Bull HN Information Systems Inc.
    Copyright 2006 by Bull SAS
    All Rights Reserved