GAIM is an Interactive Program for fitting Generalized Additive
Models.

GAIM is constantly under development by Trevor Hastie and Rob
Tibshirani, so if you have any problems with the program or suggested
improvements, please contact us.  Trevor Hastie- AT&T Bell Labs,
2C-261, 600 Mountain Ave, Murray Hill, NJ07974.  (201) 582-5647.  Rob
Tibshirani- Dept.  of Preventive Medicine and Biostatistics, U. of
Toronto, Toronto Canada M5S1A8, (416)-978-4642.

The program is written in MORTRAN, a FORTRAN preprocessor.  The
 distribution tape has both a FORTRAN and a MORTRAN version.  However,
the FORTRAN produced by MORTRAN is rather messy.

ref:  Hastie, T and Tibshirani, R (1986), Generalized Additive Models,
Statistical Science, 1, 297-318
Hastie, T. and Tibshirani, R. (1987), Generalized Additive Models;
  some Applications, JASA, June 1987.

pause


LIMITATIONS
___________ 
The current limits are set at 1000 observations
and 20 variables. This can easily be changed (by you).

SUMMARY of CAPABILITIES 
_______________________ 
GAIM can fit additive models to GAUSSIAN, BINOMIAL, and MULTINOMIAL
(proportional odds) data, as well as CENSORED survival type data and
MATCHED CASE-CONTROL data.  For the binomial and multinomial, the data
can be grouped or ungrouped.  Users can specifiy their own models by
writing their own versions of the link, deviance, inverse link, weight
and derivative routines.  The models can range from all terms linear to
all terms "smooth", although typically one fits some of each.  One can
generate dummy variables quite easily for the levels of a categorical
variable.  Using the SPECIAL option, you can get confidence bands (+-2
sd) for the fitted functions, as well as variances and covariances for
the linear terms (the current implementation does not allow this for
the MULTINOMIAL, COX or MATCHED CASE-CONTROL Model) 
The program can deal with MISSING data explicitly, and
implicitly one can fit interactions between discrete and smooth terms.
The program allows one to build models in a natural GLIM type way by
adding to or deleteing from a current model.  After each fit, a page of
statistics and details are printed giving slopes, deviance etc.  The
program allows one to print the fitted functions to a PLOT file, which
can then be independently plotted.  At present, the subroutine DUMPS
creates a file for each plot, which can then be manipulated.


pause 
DATA ENTRY
__________ 
The program expects 2 data files in a UNIX
directory, and will prompt for the name of this
directory.  
Unit number 1 reads from a  file called "name/data" which
has a N x P block of  data: N observations on P variables in free
format.  One observation, (p numbers) per line.  
Unit number 11 reads
from a "name/datad" file (Data description). The first line
contains N and P. P lines follow, each line containing a variable
name.

In addition, this file (helpunix) should be in the gaim directory.
Unit 2 is the PLOT output, and UNIT 3 gives fitted values when these
 are requested.


pause 
some details...


ADD, DELETE, MODIFY or RUN 
__________________________ 
These options are
rather obvious, and allow for model building.

ADD--displays a list of the variables NOT currently in the model, and
allows the user to select one. Upon selection, you are requested to
give a span which can be any real number. The following actions are
taken

    SPAN          ACTION
    __________________________________________________________ 
    s < 0   | A cubic Spline Smoother is used with lamda =-span                       s = 0   | The span is chosen by cross-validation 
    0<s<1   | 100s% of the data is used in each neighborhood
	    |   except at the ends, which drop down to 100s/2% 
    s=1     | This fits a linear fit for this term.
    s=2,3,..| Categorical variable with 2,3... categories. If there
	    |   are less then s categories, it allows you to exit or
	    |   continue. If there are more than s, it allows you to exit
            |   or combine the last additional categories.

The categorical option creates k-1 dummy variables if there are k
categories.  The first category gets no dummy, and corresponds to the
constant term in the model. All smooth fits are standardized to have
mean zero---linear fits are not.

pause

DELETE-- This allows you to delete terms in the model.

MODIFY-- This allows you to change the thresh-hold for the BACKFITTING
algorithm (set at !RSSOLD-RSSNEW!/RSSOLD < THBAK=0.0001 ) or the LOCAL
SCORING algorithm (THGAM = 0.0001).

It also allow you to alter the spans for variables. Do not try to
change a continuous variable to a discrete one though.

RUN-- Fits the current model.


pause

MISSING DATA 
____________ 
GAIM treats any observation having the value
99999 (5 nines) as missing (SORRY, blame Fortran).  The procedure is
naturally geared to handle missing data; any fit (smooth, or linear) is
obtained via the backfitting algorithm, and is thus the result of a
simple smooth or linear fit. As such, the missing data is only excluded
when each individual function is estimated. i.e. Missing cases are NOT
removed.  The contribution to the fit for a missing variable is the
average contribution for that variable. Cases with missing RESPONSES
are obviously removed.

MISSING DATA AND INTERACTIONS 
_____________________________ 
For linear
terms, interactions between Linear and Dummy variables can be obtained
in the usual way--by multiplication. For non-parametric fits, this does
not work. Instead, one can use the missing value option to fudge it.
i.e. 
 ________________________________________________________________
   v1!       1,     2, 99999,     7,     6, 99999,999999,99999. !
   v2!   99999, 99999,     9, 99999, 99999,    11,     8,   12  !
group!       1      1      2      1      1      2      2     2  !
_____!__________________________________________________________!

sets up two variables v1 and v2 according to wether group is 1 or 2.

pause

SPECIAL VARIABLES 
_________________ 
Allowance is made for a WEIGHT
variable, and for the BINOMIAL model, one can specify N (r out of N).
For the COX model, one has to specify the censoring variable (0,1).
For the MATCHED CASE-CONTROL model, one codes the (SINGLE) case as 1,
the controls as 0, and then is asked for the matched set indicator. This is 
an integer array with the numbers between 1 and K, where K is the number of sets.

SPECIAL OPTION 
______________ 
This is used to get +-2 st deviation
curves for non-parametric fits and  a covariance matrices for the
parameters of the linear fits.  It uses a lot of CPU, since it
basically fits the current model N times. We are busy working on more
efficient approximate methods for calculating these quantities.

pause

SMOOTHER 
________ 
The current version uses either a LOCAL LINEAR smoother,
			     or a CUBIC SPLINE Smoother.

The LOCAL LINEAR smoother has the following features:

Each neighbourhood contains 100s% of the total WEIGHT, with endpoint
modifications. Weight is defined initially, and also in the non-linear
models is defined implicitly.  A preliminary scan thru the data
replaces all Y's tied in X by their average at that X. A post scan
ensures that all those fitted values are the same as well. See the
references for more details.

The CUBIC SPLINE smoother was kindly supplied by Finbarr O'Sullivan, UC
Berkeley.  It uses the deBoor routines with some fancy Choleski work to
cross-validate in O(n) time. The present implementation requires a
fixed parameter, but this is easily changed.  pause


USER SUPPLIED MODELS 
____________________ 
The following are examples of
user supplied functions which operate in a similar way to the GLIM OWN
macros.  They must be compiled and linked with the rest of the
program.

      REAL FUNCTION DEVU(N,FITS,Y,W,NI) REAL FITS(1),Y(1),W(1),NI(1)
      RSS=0.0 DO 1 I=1,N RSS=RSS+W(I)*(Y(I)-FITS(I))*2.0 1     CONTINUE
      DEVU=RSS RETURN END

      REAL FUNCTION LINKU(ETAHAT) LINKU=ETAHAT RETURN END

      REAL FUNCTION LINVU(MUHAT) REAL MUHAT LINVU=MUHAT RETURN END

      REAL FUNCTION DIRIVU(MUHAT) REAL MUHAT DIRIVU=1.0 RETURN END

      REAL FUNCTION WEIGTU(W,MUHAT,NI) REAL W,MUHAT,NI WEIGTU=1.0
      RETURN END



pause

 PARTIAL RESIDUALS

Note: when partial residuals are requested, the GAIM prints out

x  partial residual     smooth(x)   rank(x)    fitted value    y
pause
*************************************************************************
This Help file was created in a hurry, and will evolve in a natural
way.

*************************************************************************


stop




C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine addtm
      real dumlab(30),modlab(20)
      integer indat(20),indmod(20)
      common/extra/dumlab,modlab,indat,indmod
      integer neffob(20)
      common/neff/neffob
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
      integer sltag(1000)
      logical vlink,vflag
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
      logical ok
      ok=.false.
10011 continue
      write(6,10020)
10020 format(   45h the following variables are not in the model)
      icount=0
      j=1
      goto 10033
10031 j=j+1
10033 if((j).gt.(maxp))goto 10032
      if(indat(j) .ne. 0)goto 10051
      icount=icount+1
      write(6,10060)j,(labels(k,j),k=1,10)
10060 format(i3,    5h ....,10a4)
10051 continue
      goto 10031
10032 continue
      if(icount .le. 0)goto 10081
10091 continue
      write(6,10100)
10100 format(   33h index of new covariate (0= quit))
      read(5,*)ir
      if(ir .gt. 0)goto 10121
      return
10121 continue
      if(indat(ir) .ne. 0)goto 10141
      goto 10092
10141 continue
      write(6,10150)ir
10150 format(i3,   13h not allowed!)
      goto 10091
10092 continue
      p=p+1
      indat(ir)=1
      index(p)=ir
      indmod(p)=1
      modlab(p)=dumlab(1)
      do 10161 i=1,n
      x(i,p)=datam(i,ir)
10161 continue
10162 continue
      cross(p)=.false.
      write(6,10170)ir
10170 format(   19h span for covariate,i3,   18h in interval (0,1),/,
     *  23h 0 --- cross-validated ,/,     13h 1 --- linear,/,     41h -
     *--- cubic spline with lambda = - span ,/,     46h -999   cubic spl
     *ine with auto span selection ,/,     32h k --- categorical with k
     *levels)
      read(5,*)spans(p)
      if((spans(p) .ne. 0) .and. (spans(p) .ne. -999))goto 10191
      cross(p)=.true.
      indmod(p)=3
      modlab(p)=dumlab(21)
10191 continue
      if(spans(p) .le. 1.0)goto 10211
      call getdum(ir)
      goto 10221
10211 continue
      if(spans(p) .ne. 1.0)goto 10241
      indmod(p)=2
      modlab(p)=dumlab(22)
10241 continue
      call mysort(x(1,p),tag(1,p),1,n)
      nef=n
10251 if(x(tag(nef,p),p) .ne. 99999) goto 10252
      nef=nef-1
      goto 10251
10252 continue
      neffob(p)=nef
      miss=n-nef
      if(miss .le. 0)goto 10271
      write(6,10280)ir,miss
10280 format (//   10h variable ,i2,    5h has ,i3,   21h missing observ
     *ations,/,     56h they will receive a fitted value of 0 for this v
     *ariable,/)
10271 continue
10221 continue
10201 continue
      goto 10291
10081 continue
      write(6,10300)
10300 format(   24h oops....they are all in)
      return
10291 continue
10071 continue
      goto 10011
10012 continue
      return
      end
      subroutine getdum(ir)
      real dumlab(30),modlab(20)
      integer indat(20),indmod(20)
      common/extra/dumlab,modlab,indat,indmod
      integer neffob(20)
      common/neff/neffob
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real catval(20),misans
      call mysort(x(1,p),tag(1,p),1,n)
      indcat=1
      catval(1)=x(tag(1,p),p)
      ii=2
      goto 10313
10311 ii=ii+1
10313 if((ii).gt.(n))goto 10312
      i=tag(ii,p)
      if(w(i) .le. 0.0)goto 10331
      if(x(i,p) .le. catval(indcat))goto 10351
      indcat=indcat+1
      if(indcat .le. 20)goto 10371
      write(6,10380)
10380 format(   24h more than 20 categories)
      p=p-1
      indat(ir)=0
      return
10371 continue
      catval(indcat)=x(i,p)
10351 continue
10331 continue
      goto 10311
10312 continue
      icat=aint(spans(p))
      ifmiss=0
      if(catval(indcat) .ne. 99999)goto 10401
      icat=icat+1
      nef=n
      do 10411 i=1,n
      if(x(i,p) .ne. 99999)goto 10431
      nef=nef-1
10431 continue
10411 continue
10412 continue
      miss=n-nef
      write(6,10440)ir,miss
10440 format (//   10h variable ,i2,    5h has ,i3,   21h missing observ
     *ations,/,     45h do you want a separate category for them (y),/,
     *    56h or do you want them to receive the average constant (n))
      read(5,10450)misans
      if(misans .ne. no)goto 10471
      ifmiss=1
10471 continue
10401 continue
      ndum=indcat-1
      if(indcat .eq. icat)goto 10491
      i1=indcat-ifmiss
      i2=icat-ifmiss
      write(6,10500)i1,i2
10500 format(   11h there are ,i3,   22h categories; you said ,i3,/,
     *    44h do you wish to continue (yes) or go back to,
     *	14h the menu (no))
      read(5,10450)answer
10450 format(a1)
      if(answer .ne. no)goto 10521
      p=p-1
      indat(ir)=0
      return
      goto 10531
10521 continue
      if(indcat .le. icat)goto 10551
      ndum=icat-1
      write(6,10560)
10560 format(   39h there are more categories than dummies,/,     36h in
     *itial categories will be combined)
10551 continue
10531 continue
10511 continue
10491 continue
      i1=ndum-ifmiss
      write(6,10570)i1,i1
10570 format(   10h creating ,i2,   17h  dummy variables,/,     14h for
     *the last ,i2,   11h categories)
      indat(ir)=indcat+100
      indc=indcat-ndum
      j=p
      goto 10583
10581 j=j+1
10583 if((j).gt.(p+ndum-1))goto 10582
      if(j .le. 20)goto 10601
      write(6,10610)
10610 format (   35h oops, you`ve exceeded $p variables)
      stop
10601 continue
      indc=indc+1
      cross(j)=.false.
      index(j)=ir
      indmod(j)=100+j
      modlab(j)=dumlab(indc)
      do 10621 i=1,n
      x(i,j)=0.0
10621 continue
10622 continue
      do 10631 i=1,n
      if(datam(i,ir) .ne. catval(indc))goto 10651
      x(i,j)=1.0
10651 continue
      tag(i,j)=i
10631 continue
10632 continue
      spans(j)=1.0
      neffob(j)=n
      goto 10581
10582 continue
      p=p+ndum-1
      if(catval(indcat) .ne. 99999 .or. ifmiss .ne. 0)goto 10671
      modlab(p)=dumlab(23)
10671 continue
      if(ifmiss .ne. 1)goto 10691
      j=p-1
      goto 10703
10701 j=j+(-1)
10703 if((-1)*((j)-(p+1-ndum)).gt.0)goto 10702
      do 10711 i=1,n
      x(i,j)=x(i,j)+99999*x(i,p)
10711 continue
10712 continue
      call mysort(x(1,j),tag(1,j),1,n)
      neffob(j)=nef
      goto 10701
10702 continue
      p=p-1
10691 continue
      return
      end
      subroutine readd
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real dumlab(30),modlab(20)
      integer indat(20),indmod(20)
      common/extra/dumlab,modlab,indat,indmod
      character*60 dsname
      common/dnm/dsname,ilen
      logical ds
      character*80 datan,datadn
      character lett
      write(6,10720)
10720 format(   29h Give the data set directory )
      read(5,10730)dsname
10730 format(a60)
      lett=dsname(1:1)
      i=1
10741 if((lett.eq.    1h ) 	) goto 10742
      	i=i+1
      	lett=dsname(i:i)
      	
      goto 10741
10742 continue
      ilen=i
      datan=dsname
      datadn=dsname
      datadn(i:i+6)=    6h/datad
      datan(i:i+5)=    5h/data
      inquire(file=datan, exist=ds)
      if(.not.(.not.ds))goto 10761
      write(6,10770)datan
10770 format(a80,   27h does not exist! exiting...)
      stop
10761 continue
      inquire(file=datadn, exist=ds)
      if(.not.(.not.ds))goto 10791
      write(6,10800)datadn
10800 format(a80,   27h does not exist! exiting...)
      stop
10791 continue
      open(1,file=datan)
      rewind 1
      open(11,file=datadn)
      rewind 11
      do 10811 i=1,20
      indat(i)=-1
10811 continue
10812 continue
      read(11,*)n,maxp
      if(maxp .le. 20)goto 10831
      write(6,10840)
10840 format(   19h too many variables)
      stop
10831 continue
      do 10851 i=1,maxp
      read(11,10860)(labels(j,i),j=1,10)
      indat(i)=0
10851 continue
10852 continue
10860 format(10a4)
      i=1
      goto 10873
10871 i=i+1
10873 if((i).gt.(n))goto 10872
      read(1,*,end=10880)(datam(i,j),j=1,maxp)
      goto 10871
10872 continue
      return
10880 continue
      write(6,10890)
10890 format(   23h you ve run out of data)
      stop
      end
      subroutine deltm
      real dumlab(30),modlab(20)
      integer indat(20),indmod(20)
      common/extra/dumlab,modlab,indat,indmod
      integer neffob(20)
      common/neff/neffob
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
      integer sltag(1000)
      logical vlink,vflag
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
      logical ok
      ok=.false.
10901 continue
      if(p .ne. 0)goto 10921
      write(6,10930)
10930 format(   11h null model)
      return
10921 continue
      write(6,10940)
10940 format(   41h the following variables are in the model)
      j=1
      goto 10953
10951 j=j+1
10953 if((j).gt.(p))goto 10952
      write(6,10960)index(j),modlab(j),(labels(k,index(j)),k=1,10)
10960 format (i3,    5h ....,11a4)
      goto 10951
10952 continue
      write(6,10970)
10970 format(   43h index of covariate to be deleted (0= quit))
      read(5,*)ir
      if(ir .ne. 0)goto 10991
      return
10991 continue
      if((indat(ir) .ne. 1) .and. (indat(ir) .le. 100))goto 11011
      indat(ir)=0
      j=1
      goto 11023
11021 j=j+1
11023 if((j).gt.(p))goto 11022
      if(ir .ne. index(j))goto 11041
      goto 11022
11041 continue
      goto 11021
11022 continue
      if(p .le. 1)goto 11061
      ij=1
      jj=j
11071 if(index(jj+1) .ne. ir .or. (jj+1) .gt. p) goto 11072
      jj=jj+1
      ij=ij+1
      goto 11071
11072 continue
      k=j
      goto 11083
11081 k=k+1
11083 if((k).gt.(p-ij))goto 11082
      index(k)=index(k+ij)
      cross(k)=cross(k+ij)
      spans(k)=spans(k+ij)
      indmod(k)=indmod(k+ij)
      modlab(k)=modlab(k+ij)
      neffob(k)=neffob(k+ij)
      i=1
      goto 11093
11091 i=i+1
11093 if((i).gt.(n))goto 11092
      tag(i,k)=tag(i,k+ij)
      x(i,k)=x(i,k+ij)
      goto 11091
11092 continue
      goto 11081
11082 continue
      p=p-ij
      goto 11101
11061 continue
      p=0
      return
11101 continue
11051 continue
      goto 11111
11011 continue
      write(6,11120)ir
11120 format(   10h variable ,i3,   20h is not in the model)
11111 continue
11001 continue
      goto 10901
10902 continue
      end
C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine bakfit(n,p,y,x,w,tag,spans,dofs,slopes,cross,  etahat,s
     *,s0,rss,thresh,nit,verbos,quiet,itmax)
      integer n,p,tag(1000,1),nit,itmax
      logical verbos,quiet,cross(1)
      integer neffob(20)
      common/neff/neffob
      real dumlab(30),modlab(20)
      integer indat(20),indmod(20)
      common/extra/dumlab,modlab,indat,indmod
      integer bacnit
      common /bacn/bacnit
      real y(1),x(1000,1),w(1),spans(1),dofs(1),slopes(1),  etahat(1),s(
     *1000,1),s0,rss,thresh
      real z(1000),scrat(1000)
      if(p .ne. 1)goto 10021
      itmax=1
10021 continue
      rss=devg(n,etahat,y,w,junk)
      rssold=10*rss+10
      s1=0.0
      sn=0.0
      se=0.0
      i=1
      goto 10033
10031 i=i+1
10033 if((i).gt.(n))goto 10032
      se=se+w(i)*etahat(i)
      s1=s1+w(i)*y(i)
      sn=sn+w(i)
      goto 10031
10032 continue
      s1=s1/sn
      se=se/sn
      do 10041 i=1,n
      etahat(i)=etahat(i)-se+s1
10041 continue
10042 continue
      s0=s0-se+s1
      nit=0
10051 if(abs((rssold-rss)/rssold) .le. thresh .or. nit .ge. itmax) goto
     *10052
      rssold=rss
      j=1
      goto 10063
10061 j=j+1
10063 if((j).gt.(p))goto 10062
      i=1
      goto 10073
10071 i=i+1
10073 if((i).gt.(n))goto 10072
      etahat(i)=etahat(i)-s(i,j)
      z(i)=y(i)-etahat(i)
      goto 10071
10072 continue
      call smooth(x(1,j),z,w,tag(1,j),spans(j),dofs(j),n,neffob(j),  cro
     *ss(j),s(1,j),s1,rss,quiet,scrat)
      slopes(j)=scrat(1)
      do 10081 i=1,n
      etahat(i)=etahat(i)+s(i,j)+s1
10081 continue
10082 continue
      s0=s0+s1
      if(.not.(verbos))goto 10101
      write(6,10110)j,rss
10110 format(10x,   26hinner loop (backfit) var= ,  i3,    5h rss=,f16.5
     *)
10101 continue
      goto 10061
10062 continue
      nit=nit+1
      if(.not.(.not.quiet))goto 10131
      write(6,10140)nit,rss
10140 format(/10x,   26houter loop (backfit) nit= ,  i3,    5h rss=,f16.
     *5,/)
10131 continue
      goto 10051
10052 continue
      bacnit=nit
      return
      end
      subroutine bsplvd ( t, k, x, left, a, dbiatx, nderiv )
calls bsplvb
calculates value and deriv.s of all b-splines which do not vanish at x
c
c******  i n p u t  ******
c  t     the knot array, of length left+k (at least)
c  k     the order of the b-splines to be evaluated
c  x     the point at which these values are sought
c  left  an integer indicating the left endpoint of the interval of
c        interest. the  k  b-splines whose support contains the interval
c               (t(left), t(left+1))
c        are to be considered.
c  a s s u m p t i o n  - - -  it is assumed that
c               t(left) .lt. t(left+1)
c        division by zero will result otherwise (in  b s p l v b ).
c        also, the output is as advertised only if
c               t(left) .le. x .le. t(left+1) .
c  nderiv   an integer indicating that values of b-splines and their
c        derivatives up to but not including the  nderiv-th  are asked
c        for. ( nderiv  is replaced internally by the integer in (1,k)
c        closest to it.)
c
c******  w o r k   a r e a  ******
c  a     an array of order (k,k), to contain b-coeff.s of the derivat-
c        ives of a certain order of the  k  b-splines of interest.
c
c******  o u t p u t  ******
c  dbiatx   an array of order (k,nderiv). its entry  (i,m)  contains
c        value of  (m-1)st  derivative of  (left-k+i)-th  b-spline of
c        order  k  for knot sequence  t , i=m,...,k; m=1,...,nderiv.
c
c******  m e t h o d  ******
c  values at  x  of all the relevant b-splines of order k,k-1,...,
c  k+1-nderiv  are generated via  bsplvb  and stored temporarily
c  in  dbiatx .  then, the b-coeffs of the required derivatives of the
c  b-splines of interest are generated by differencing, each from the
c  preceding one of lower order, and combined with the values of b-
c  splines of corresponding order in  dbiatx  to produce the desired
c  values.
c
      integer k,left,nderiv,   i,ideriv,il,j,jlow,jp1mid,kp1,kp1mm
     *      ,ldummy,m,mhigh
      real a(k,k),dbiatx(k,nderiv),t(1),x,   factor,fkp1mm,sum
      mhigh = max0(min0(nderiv,k),1)
c     mhigh is usually equal to nderiv.
      kp1 = k+1
      call bsplvb(t,kp1-mhigh,1,x,left,dbiatx)
      if (mhigh .eq. 1)                 go to 99
c     the first column of  dbiatx  always contains the b-spline values
c     for the current order. these are stored in column k+1-current
c     order  before  bsplvb  is called to put values for the next
c     higher order on top of it.
      ideriv = mhigh
      do 15 m=2,mhigh
         jp1mid = 1
         do 11 j=ideriv,k
            dbiatx(j,ideriv) = dbiatx(jp1mid,1)
   11       jp1mid = jp1mid + 1
         ideriv = ideriv - 1
         call bsplvb(t,kp1-ideriv,2,x,left,dbiatx)
   15    continue
c
c     at this point,  b(left-k+i, k+1-j)(x) is in  dbiatx(i,j) for
c     i=j,...,k and j=1,...,mhigh ('=' nderiv). in particular, the
c     first column of  dbiatx  is already in final form. to obtain cor-
c     responding derivatives of b-splines in subsequent columns, gene-
c     rate their b-repr. by differencing, then evaluate at  x.
c
      jlow = 1
      do 20 i=1,k
         do 19 j=jlow,k
   19       a(j,i) = 0.
         jlow = i
   20    a(i,i) = 1.
c     at this point, a(.,j) contains the b-coeffs for the j-th of the
c     k  b-splines of interest here.
c
      do 40 m=2,mhigh
         kp1mm = kp1 - m
         fkp1mm = float(kp1mm)
         il = left
         i = k
c
c        for j=1,...,k, construct b-coeffs of  (m-1)st  derivative of
c        b-splines from those for preceding derivative by differencing
c        and store again in  a(.,j) . the fact that  a(i,j) = 0  for
c        i .lt. j  is used.sed.
         do 25 ldummy=1,kp1mm
            factor = fkp1mm/(t(il+kp1mm) - t(il))
c           the assumption that t(left).lt.t(left+1) makes denominator
c           in  factor  nonzero.
            do 24 j=1,i
   24          a(i,j) = (a(i,j) - a(i-1,j))*factor
            il = il - 1
   25       i = i - 1
c
c        for i=1,...,k, combine b-coeffs a(.,i) with b-spline values
c        stored in dbiatx(.,m) to get value of  (m-1)st  derivative of
c        i-th b-spline (of interest here) at  x , and store in
c        dbiatx(i,m). storage of this value over the value of a b-spline
c        of order m there is safe since the remaining b-spline derivat-
c        ive of the same order do not use this value due to the fact
c        that  a(j,i) = 0  for j .lt. i .
   30    do 40 i=1,k
            sum = 0.
            jlow = max0(i,m)
            do 35 j=jlow,k
   35          sum = a(j,i)*dbiatx(j,m) + sum
   40       dbiatx(i,m) = sum
   99                                   return
      end
      subroutine bsplvb ( t, jhigh, index, x, left, biatx )
calculates the value of all possibly nonzero b-splines at  x  of order
c
c               jout  =  max( jhigh , (j+1)*(index-1) )
c
c  with knot sequence  t .
c
c******  i n p u t  ******
c  t.....knot sequence, of length  left + jout  , assumed to be nonde-
c        creasing.  a s s u m p t i o n . . . .
c                       t(left)  .lt.  t(left + 1)   .
c   d i v i s i o n  b y  z e r o  will result if  t(left) = t(left+1)
c  jhigh,
c  index.....integers which determine the order  jout = max(jhigh,
c        (j+1)*(index-1))  of the b-splines whose values at  x  are to
c        be returned.  index  is used to avoid recalculations when seve-
c        ral columns of the triangular array of b-spline values are nee-
c        ded (e.g., in  bvalue  or in  bsplvd ). precisely,
c                     if  index = 1 ,
c        the calculation starts from scratch and the entire triangular
c        array of b-spline values of orders 1,2,...,jhigh  is generated
c        order by order , i.e., column by column .
c                     if  index = 2 ,
c        only the b-spline values of order  j+1, j+2, ..., jout  are ge-
c        nerated, the assumption being that  biatx , j , deltal , deltar
c        are, on entry, as they were on exit at the previous call.
c           in particular, if  jhigh = 0, then  jout = j+1, i.e., just
c        the next column of b-spline values is generated.
c
c  w a r n i n g . . .  the restriction   jout .le. jmax (= 20)  is im-
c        posed arbitrarily by the dimension statement for  deltal  and
c        deltar  below, but is  n o w h e r e  c h e c k e d  for .
c
c  x.....the point at which the b-splines are to be evaluated.
c  left.....an integer chosen (usually) so that
c                  t(left) .le. x .le. t(left+1)  .
c
c******  o u t p u t  ******
c  biatx.....array of length  jout , with  biatx(i)  containing the val-
c        ue at  x  of the polynomial of order  jout  which agrees with
c        the b-spline  b(left-jout+i,jout,t)  on the interval (t(left),
c        t(left+1)) .
c
c******  m e t h o d  ******
c  the recurrence relation
c
c                       x - t(i)              t(i+j+1) - x
c     b(i,j+1)(x)  =  -----------b(i,j)(x) + ---------------b(i+1,j)(x)
c                     t(i+j)-t(i)            t(i+j+1)-t(i+1)
c
c  is used (repeatedly) to generate the (j+1)-vector  b(left-j,j+1)(x),
c  ...,b(left,j+1)(x)  from the j-vector  b(left-j+1,j)(x),...,
c  b(left,j)(x), storing the new values in  biatx  over the old. the
c  facts that
c            b(i,1) = 1  if  t(i) .le. x .lt. t(i+1)
c  and that
c            b(i,j)(x) = 0  unless  t(i) .le. x .lt. t(i+j)
c  are used. the particular organization of the calculations follows al-
c  gorithm  (8)  in chapter x of the text.
c
      parameter(jmax = 20)
      integer index,jhigh,left,   i,j,jp1
      real biatx(jhigh),t(1),x,   deltal(jmax),deltar(jmax),saved,term
c     dimension biatx(jout), t(left+jout)
current fortran standard makes it impossible to specify the length of
c  t  and of  biatx  precisely without the introduction of otherwise
c  superfluous additional arguments.
      data j/1/
c     save j,deltal,deltar (valid in fortran 77)
c
                                        go to (10,20), index
   10 j = 1
      biatx(1) = 1.
      if (j .ge. jhigh)                 go to 99
c
   20    jp1 = j + 1
         deltar(j) = t(left+j) - x
         deltal(j) = x - t(left+1-j)
         saved = 0.
         do 26 i=1,j
            term = biatx(i)/(deltar(i) + deltal(jp1-i))
            biatx(i) = saved + deltar(i)*term
   26       saved = deltal(jp1-i)*term
         biatx(jp1) = saved
         j = jp1
         if (j .lt. jhigh)              go to 20
c
   99                                   return
      end

      function bvalue ( t, bcoef, n, k, x, jderiv )
calls  interv
c
calculates value at  x  of  jderiv-th derivative of spline from b-repr.
c  the spline is taken to be continuous from the right.
c
c******  i n p u t ******
c  t, bcoef, n, k......forms the b-representation of the spline  f  to
c        be evaluated. specifically,
c  t.....knot sequence, of length  n+k, assumed nondecreasing.
c  bcoef.....b-coefficient sequence, of length  n .
c  n.....length of  bcoef  and dimension of s(k,t),
c        a s s u m e d  positive .
c  k.....order of the spline .
c
c  w a r n i n g . . .   the restriction  k .le. kmax (=20)  is imposed
c        arbitrarily by the dimension statement for  aj, dm, dm  below,
c        but is  n o w h e r e  c h e c k e d  for.
c
c  x.....the point at which to evaluate .
c  jderiv.....integer giving the order of the derivative to be evaluated
c        a s s u m e d  to be zero or positive.
c
c******  o u t p u t  ******
c  bvalue.....the value of the (jderiv)-th derivative of  f  at  x .
c
c******  m e t h o d  ******
c     the nontrivial knot interval  (t(i),t(i+1))  containing  x  is lo-
c  cated with the aid of  interv . the  k  b-coeffs of  f  relevant for
c  this interval are then obtained from  bcoef (or taken to be zero if
c  not explicitly available) and are then differenced  jderiv  times to
c  obtain the b-coeffs of  (d**jderiv)f  relevant for that interval.
c  precisely, with  j = jderiv, we have from x.(12) of the text that
c
c     (d**j)f  =  sum ( bcoef(.,j)*b(.,k-j,t) )
c
c  where
c                   / bcoef(.),                     ,  j .eq. 0
c                   /
c    bcoef(.,j)  =  / bcoef(.,j-1) - bcoef(.-1,j-1)
c                   / ----------------------------- ,  j .gt. 0
c                   /    (t(.+k-j) - t(.))/(k-j)
c
c     then, we use repeatedly the fact that
c
c    sum ( a(.)*b(.,m,t)(x) )  =  sum ( a(.,x)*b(.,m-1,t)(x) )
c  with
c                 (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1)
c    a(.,x)  =    ---------------------------------------
c                 (x - t(.))      + (t(.+m-1) - x)
c
c  to write  (d**j)f(x)  eventually as a linear combination of b-splines
c  of order  1 , and the coefficient for  b(i,1,t)(x)  must then
c  be the desired number  (d**j)f(x). (see x.(17)-(19) of text).
c
      parameter(kmax = 20)
      integer jderiv,k,n,   i,ilo,imk,j,jc,jcmin,jcmax,jj,km1,mflag,nmi
      real bcoef(n),t(1),x,   aj(kmax),dm(kmax),dp(kmax),fkmj
c     dimension t(n+k)
current fortran standard makes it impossible to specify the length of
c  t  precisely without the introduction of otherwise superfluous
c  additional arguments.
      bvalue = 0.
      if (jderiv .ge. k)                go to 99
c
c  *** find  i  s.t.  1 .le. i .lt. n+k  and  t(i) .lt. t(i+1) and
c      t(i) .le. x .lt. t(i+1) . if no such i can be found,  x  lies
c      outside the support of  the spline  f  and  bvalue = 0.
c      (the asymmetry in this choice of  i  makes  f  rightcontinuous)
      if( (x.ne.t(n+1)) .or. (t(n+1).ne.t(n+k)) )  go to 700
      i = n
                                        go to 701
  700 call interv ( t, n+k, x, i, mflag )
      if (mflag .ne. 0)                 go to 99
  701 continue
c  *** if k = 1 (and jderiv = 0), bvalue = bcoef(i).
      km1 = k - 1
      if (km1 .gt. 0)                   go to 1
      bvalue = bcoef(i)
                                        go to 99
c
c  *** store the k b-spline coefficients relevant for the knot interval
c     (t(i),t(i+1)) in aj(1),...,aj(k) and compute dm(j) = x - t(i+1-j),
c     dp(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable
c     from input to zero. set any t.s not obtainable equal to t(1) or
c     to t(n+k) appropriately.
    1 jcmin = 1
      imk = i - k
      if (imk .ge. 0)                   go to 8
      jcmin = 1 - imk
      do 5 j=1,i
    5    dm(j) = x - t(i+1-j)
      do 6 j=i,km1
         aj(k-j) = 0.
    6    dm(j) = dm(i)
                                        go to 10
    8 do 9 j=1,km1
    9    dm(j) = x - t(i+1-j)
c
   10 jcmax = k
      nmi = n - i
      if (nmi .ge. 0)                   go to 18
      jcmax = k + nmi
      do 15 j=1,jcmax
   15    dp(j) = t(i+j) - x
      do 16 j=jcmax,km1
         aj(j+1) = 0.
   16    dp(j) = dp(jcmax)
                                        go to 20
   18 do 19 j=1,km1
   19    dp(j) = t(i+j) - x
c
   20 do 21 jc=jcmin,jcmax
   21    aj(jc) = bcoef(imk + jc)
c
c               *** difference the coefficients  jderiv  times.
      if (jderiv .eq. 0)                go to 30
      do 23 j=1,jderiv
         kmj = k-j
         fkmj = float(kmj)
         ilo = kmj
         do 23 jj=1,kmj
            aj(jj) = ((aj(jj+1) - aj(jj))/(dm(ilo) + dp(jj)))*fkmj
   23       ilo = ilo - 1
c
c  *** compute value at  x  in (t(i),t(i+1)) of jderiv-th derivative,
c     given its relevant b-spline coeffs in aj(1),...,aj(k-jderiv).
   30 if (jderiv .eq. km1)              go to 39
      jdrvp1 = jderiv + 1
      do 33 j=jdrvp1,km1
         kmj = k-j
         ilo = kmj
         do 33 jj=1,kmj
            aj(jj) = (aj(jj+1)*dm(ilo) + aj(jj)*dp(jj))/(dm(ilo)+dp(jj))
   33       ilo = ilo - 1
   39 bvalue = aj(1)
c
   99                                   return
      end
      subroutine interv ( xt, lxt, x, left, mflag )
computes  left = max( i ; 1 .le. i .le. lxt  .and.  xt(i) .le. x )  .
c
c******  i n p u t  ******
c  xt.....a real sequence, of length  lxt , assumed to be nondecreasing
c  lxt.....number of terms in the sequence  xt .
c  x.....the point whose location with respect to the sequence  xt  is
c        to be determined.
c
c******  o u t p u t  ******
c  left, mflag.....both integers, whose value is
c
c   1     -1      if               x .lt.  xt(1)
c   i      0      if   xt(i)  .le. x .lt. xt(i+1)
c  lxt     1      if  xt(lxt) .le. x
c
c        in particular,  mflag = 0 is the 'usual' case.  mflag .ne. 0
c        indicates that  x  lies outside the halfopen interval
c        xt(1) .le. y .lt. xt(lxt) . the asymmetric treatment of the
c        interval is due to the decision to make all pp functions cont-
c        inuous from the right.
c
c******  m e t h o d  ******
c  the program is designed to be efficient in the common situation that
c  it is called repeatedly, with  x  taken from an increasing or decrea-
c  sing sequence. this will happen, e.g., when a pp function is to be
c  graphed. the first guess for  left  is therefore taken to be the val-
c  ue returned at the previous call and stored in the  l o c a l  varia-
c  ble  ilo . a first check ascertains that  ilo .lt. lxt (this is nec-
c  essary since the present call may have nothing to do with the previ-
c  ous call). then, if  xt(ilo) .le. x .lt. xt(ilo+1), we set  left =
c  ilo  and are done after just three comparisons.
c     otherwise, we repeatedly double the difference  istep = ihi - ilo
c  while also moving  ilo  and  ihi  in the direction of  x , until
c                      xt(ilo) .le. x .lt. xt(ihi) ,
c  after which we use bisection to get, in addition, ilo+1 = ihi .
c  left = ilo  is then returned.
c
      integer left,lxt,mflag,   ihi,ilo,istep,middle
      real x,xt(lxt)
      data ilo /1/
c     save ilo  (a valid fortran statement in the new 1977 standard)
      ihi = ilo + 1
      if (ihi .lt. lxt)                 go to 20
         if (x .ge. xt(lxt))            go to 110
         if (lxt .le. 1)                go to 90
         ilo = lxt - 1
         ihi = lxt
c
   20 if (x .ge. xt(ihi))               go to 40
      if (x .ge. xt(ilo))               go to 100
c
c              **** now x .lt. xt(ilo) . decrease  ilo  to capture  x .
   30 istep = 1
   31    ihi = ilo
         ilo = ihi - istep
         if (ilo .le. 1)                go to 35
         if (x .ge. xt(ilo))            go to 50
         istep = istep*2
                                        go to 31
   35 ilo = 1
      if (x .lt. xt(1))                 go to 90
                                        go to 50
c              **** now x .ge. xt(ihi) . increase  ihi  to capture  x .
   40 istep = 1
   41    ilo = ihi
         ihi = ilo + istep
         if (ihi .ge. lxt)              go to 45
         if (x .lt. xt(ihi))            go to 50
         istep = istep*2
                                        go to 41
   45 if (x .ge. xt(lxt))               go to 110
      ihi = lxt
c
c           **** now xt(ilo) .le. x .lt. xt(ihi) . narrow the interval.
   50 middle = (ilo + ihi)/2
      if (middle .eq. ilo)              go to 100
c     note. it is assumed that middle = ilo in case ihi = ilo+1 .
      if (x .lt. xt(middle))            go to 53
         ilo = middle
                                        go to 50
   53    ihi = middle
                                        go to 50
c**** set output and return.
   90 mflag = -1
      left = 1
                                        return
  100 mflag = 0
      left = ilo
                                        return
  110 mflag = 1
      left = lxt
                                        return
      end

C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine dumpm
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real mult,group,ym(1000,15)
      real thetak(15)
      integer mcat,indrm(15)
      common/mltcmn/mult,group,thetak,ym,mcat,indrm
      real pim(15),pm(15)
      real ytot(15)
      write(6,10010)(thetak(j),j=1,mcat)
10010 format(5f10.4)
      i=1
      goto 10023
10021 i=i+1
10023 if((i).gt.(n))goto 10022
      call linkm(mcat,thetak,etahat(i),pim,pm)
      ytot(1)=ym(i,1)
      j=2
      goto 10033
10031 j=j+1
10033 if((j).gt.(mcat-1))goto 10032
      ytot(j)=ytot(j-1)+ym(i,j)
      goto 10031
10032 continue
      mcatm1=mcat-1
      write(4,10040)i,(pim(j),j=1,mcatm1)
      write(4,10040)i,(ytot(j),j=1,mcatm1)
      goto 10021
10022 continue
10040 format(i5,10f7.4)
      return
      end
C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine dumpp(n,p,x,lambda)
      integer n,p
      real x(506,p),lambda(506,2)
      real y(6)
      real sgrid(6,15,15),lam1(15),lam2(15)
      integer table(15,15),lrlam1(2,15), lrlam2(2,15)
      common /cont/ sgrid,lam1,lam2,table,lrlam1,lrlam2
      igrids=15
      i=p+2
      write(14,*)n,i,igrids
      do 10011 i=1,n
      write(14,10020)(lambda(i,j),j=1,2),(x(i,j),j=1,p)
10011 continue
10012 continue
10020 format(8f10.4)
      i=1
      goto 10033
10031 i=i+1
10033 if((i).gt.(igrids))goto 10032
      jgrid=lrlam1(2,i)-lrlam1(1,i)+1
      write(14,*)jgrid
      if(lrlam1(1,i) .le. lrlam1(2,i))goto 10051
      goto 10031
10051 continue
      j=lrlam1(1,i)
      goto 10063
10061 j=j+1
10063 if((j).gt.(lrlam1(2,i)))goto 10062
      k=1
      goto 10073
10071 k=k+1
10073 if((k).gt.(p))goto 10072
      y(k)=sgrid(k,i,j)
      goto 10071
10072 continue
      write(14,10020)lam1(i),lam2(j),(y(ij),ij=1,p)
      goto 10061
10062 continue
      goto 10031
10032 continue
      j=1
      goto 10083
10081 j=j+1
10083 if((j).gt.(igrids))goto 10082
      jgrid=lrlam2(2,j)-lrlam2(1,j)+1
      write(14,*)jgrid
      if(lrlam2(1,j) .le. lrlam2(2,j))goto 10101
      goto 10081
10101 continue
      i=lrlam2(1,j)
      goto 10113
10111 i=i+1
10113 if((i).gt.(lrlam2(2,j)))goto 10112
      k=1
      goto 10123
10121 k=k+1
10123 if((k).gt.(p))goto 10122
      y(k)=sgrid(k,i,j)
      goto 10121
10122 continue
      write(14,10020)lam1(i),lam2(j),(y(ij),ij=1,p)
      goto 10111
10112 continue
      goto 10081
10082 continue
      return
      end
C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine dumps(diriv)
      real mult,group,ym(1000,15)
      real thetak(15)
      integer mcat,indrm(15)
      common/mltcmn/mult,group,thetak,ym,mcat,indrm
      integer neffob(20)
      common/neff/neffob
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
      integer sltag(1000)
      logical vlink,vflag
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
      character*60 dsname
      common/dnm/dsname,ilen
      real diriv
      external diriv
      character*90 plotn
      real sep,part
      write(6,10010)
10010 format(   36h want to dump a plot of this model? )
      read(5,10020)ans
      if(ans .ne. no)goto 10041
      return
10041 continue
10020 format(a1)
      nitav=numsum/nitout
      r2=(devnll-devian)*100.0/devnll
      if(p .ne. 1 .or. mod .eq. mult)goto 10061
      call nplot(plotn)
      oldx=-9999
      oldy=-9999
      i=1
      goto 10073
10071 i=i+1
10073 if((i).gt.(neffob(1)))goto 10072
      ii=tag(i,1)
      if((x(ii,1) .eq. oldx) .and. (y(ii) .eq. oldy))goto 10091
      write(2,10100)x(ii,1),y(ii),muhat(ii)
      oldx=x(ii,1)
      oldy=y(ii)
10091 continue
      goto 10071
10072 continue
      close(2)
      goto 10111
10061 continue
      write(6,10120)
10120 format(   32h want to dump plots separately? )
      read(5,10020)sep
      write(6,10130)
10130 format(   25h want partial residuals? )
      read(5,10020)part
      j=1
      goto 10143
10141 j=j+1
10143 if((j).gt.(p))goto 10142
      write(6,10150)j,(labels(k,index(j)),k=1,10)
10150 format(   10h Variable ,i3,    2h: ,10a4)
      if(sep .ne. yes)goto 10171
      write(6,10180)
10180 format(   18h do you want it...)
      read(5,10020)ans
      if(ans .ne. no)goto 10201
      goto 10141
10201 continue
10171 continue
      call nplot(plotn)
      totsm=0.0
      i=1
      goto 10213
10211 i=i+1
10213 if((i).gt.(neffob(j)))goto 10212
      totsm=totsm+s(tag(i,j),j)
      goto 10211
10212 continue
      totsm=totsm/neffob(j)
      if(part .ne. yes)goto 10231
      i=1
      goto 10243
10241 i=i+1
10243 if((i).gt.(neffob(j)))goto 10242
      ii=tag(i,j)
      sput=s(ii,j)
      z=s(ii,j)+(y(ii)-ni(ii)*muhat(ii))*diriv(muhat(ii))
      if(mod .ne. mult)goto 10261
      z=-z
      sput=-sput
10261 continue
      write(2,10100)x(ii,j),z,sput,ii,muhat(ii),y(ii)
      goto 10241
10242 continue
      close(2)
      goto 10271
10231 continue
      oldx=-9999
      i=1
      goto 10283
10281 i=i+1
10283 if((i).gt.(neffob(j)))goto 10282
      ii=tag(i,j)
      if(x(ii,j) .eq. oldx)goto 10301
      sput=s(ii,j)
      if(mod .ne. mult)goto 10321
      sput=-sput
10321 continue
      write(2,10100)x(ii,j),sput
      oldx=x(ii,j)
10301 continue
      goto 10281
10282 continue
      close(2)
10271 continue
10221 continue
      write(6,10330)totsm
10330 format(   24h average of function is ,f10.3)
      goto 10141
10142 continue
10111 continue
10051 continue
      if(.not.(vlink))goto 10351
      write(6,10360)
10360 format(   14h variable link)
      call nplot(plotn)
      iii=0
      i=1
      goto 10373
10371 i=i+1
10373 if((i).gt.(n))goto 10372
      write(2,10100)argf(i),slink(i)
      goto 10371
10372 continue
      close(2)
10351 continue
10380 format(   24h data a; input x y yhat;)
10390 format(    7h data a,i1,   15h; input x yhat;)
10400 format(   39h title .c=blue estimated mean function;  /,   37hproc
     * gplot; plot y*x yhat*x /overlay;,/,     54hsymbol1 v=star c=red;s
     *ymbol2 i=join v=diamond c=green;)
10410 format(   33h label x = eta yhat=f(eta);cards;)
10420 format(   30h title .c=blue estimated link;  /,   26hproc gplot; p
     *lot  yhat*x ;,/,     32hsymbol i=join c=green v=diamond;)
10430 format(   34h title .c=blue estimated function;  /,   26hproc gplo
     *t; plot  yhat*x ;,/,     32hsymbol i=join c=green v=diamond;)
10440 format(   56h title .c=blue estimated function and partial residua
     *ls;  /,   37hproc gplot; plot y*x yhat*x /overlay;,/,     45hsymbo
     *l1 v=star c=red;symbol2 i=join  c=green;)
10450 format(    8hlabel x=,10a4,   17hyhat=muhat;cards;)
10460 format(   23hfootnote2 .c=blue span=,f4.2,    5h dof=,f6.2,     8h
     *n-nmiss=,i3,    13h #iterations=, i3,    1h(,i3,    2h);,/,     27
     *hfootnote3 .c=blue deviance=, f8.2,    4h r2=,f8.2,    2h%;)
10470 format(    8hlabel x=,10a4,   13hyhat=s;cards;)
10480 format(3f10.4)
10100 format(3f15.4,i5,2f10.4)
      return
      end
      subroutine nplot(plotn)
      character*60 dsname
      common/dnm/dsname,ilen
      character*90 plotn
      character*30 plf
      do 10491 i=1,90
      plotn(i:i)=    1h
10491 continue
10492 continue
      do 10501 i=1,30
      plf(i:i)=    1h
10501 continue
10502 continue
      close(2)
      write(6,10510)
10510 format(   30h Give a name for the plot file)
      read(5,10520)plf
10520 format(a30)
      plotn=dsname
      plotn(ilen:ilen)=    1h/
      plotn(ilen+1:ilen+30)=plf(1:30)
      open(2,file=plotn)
      rewind 2
      return
      end
C     mortran 2.79 (reserved keyword macros of 09/28/81)
      block data
      real mult,group,ym(1000,15)
      real thetak(15)
      integer mcat,indrm(15)
      common/mltcmn/mult,group,thetak,ym,mcat,indrm
      real dumlab(30),modlab(20)
      integer indat(20),indmod(20)
      common/extra/dumlab,modlab,indat,indmod
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
      integer sltag(1000)
      logical vlink,vflag
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
     *),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
      integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
      common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
     *icensq,irindq,iriskq,nq,kq
      real match
      integer indms,kms(1000)
      common/mtchco/match,indms,kms
      data yes/    1hy/, no/    1hn/, binom/    1hb/,gauss/    1hg/,cox/
     *    1hc/
      data mult/    1hm/
      data match/    1hr/
      data user/    1hu/
      data n/1000/, p/20/
      data verbos/.false./, quiet/.false./
      data vlink/.false./
      data score/    1hs/,local/    1hl/
      data thresh/0.001/,thbak/0.001/
      data dumlab/    4h    ,     4hd2: ,    4hd3: ,    4hd4: ,    4hd5:
     * ,    4hd6: ,    4hd7: ,    4hd8: ,    4hd9: ,    4hd10:,    4hd11
     *:,     4hd12:,    4hd13:,    4hd14:,    4hd15:,    4hd16:,     4hd
     *17:,    4hd18:,    4hd19:,    4hd20:,    4hcv: ,    4hlin:,    4hm
     *i: ,7*    4h    /
      end
      external devb,devg,linkg,linkb,dirivg,dirivb,weigtg,weigtb, devu,d
     *evc,linku,dirivu,weigtu,linvg,linvb,linvu,linkmc,linvmc,devmc
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real mult,group,ym(1000,15)
      real thetak(15)
      integer mcat,indrm(15)
      common/mltcmn/mult,group,thetak,ym,mcat,indrm
      integer neffob(20)
      common/neff/neffob
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
      integer sltag(1000)
      logical vlink,vflag
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
     *),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
      integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
      common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
     *icensq,irindq,iriskq,nq,kq
      real z(1000),wz(1000)
      common /wvc/z,wz
      character*60 dsname
      common/dnm/dsname,ilen
      real match
      integer indms,kms(1000)
      common/mtchco/match,indms,kms
      real adterm,determ,moterm,runmod,specia,quit
      data adterm/    1ha/, determ/    1hd/, moterm/    1hm/,specia/
     *1hs/, quit/    1hq/
      data help/    1hh/
      data runmod/    1hr/
      call title
      call readd
      call name7
      answer=yes
10011 if(answer .ne. yes) goto 10012
      call model
10021 if(answer .eq. quit) goto 10022
10031 continue
      prevan=answer
      write(6,10040)
10040 format(   50h add, delete, modify, run, specials, help  or quit)
      read(5,10050)answer
      if(answer .ne. adterm)goto 10071
      call addtm
      goto10061
10071 if(answer .ne. determ)goto 10081
      call deltm
      goto10061
10081 if(answer .ne. help)goto 10091
      call helpg
      goto10061
10091 if(answer .ne. moterm)goto 10101
      call modtm
      goto10061
10101 if(answer .ne. specia)goto 10111
      if(prevan .eq. runmod)goto 10131
      write(6,10140)
10140 format(   20h run the model first)
      goto 10151
10131 continue
      if(mod .ne. gauss)goto 10171
      call hat(weigtg)
      goto10161
10171 if(mod .ne. binom)goto 10181
      call hat(weigtb)
      goto10161
10181 if(mod .ne. cox)goto 10191
      call hat(weigtg)
      goto 10201
10191 continue
      write(6,10210)mod
10210 format(   27h sorry, not yet implemented,     16hfor  model type ,
     *a1)
10201 continue
10161 continue
10151 continue
10121 continue
      goto10061
10111 if((answer .ne. quit) .and. (answer .ne. runmod))goto 10221
      goto 10032
10221 continue
10061 continue
      goto 10031
10032 continue
      if(answer .ne. quit)goto 10241
      goto 10022
10241 continue
      if(mod .ne. gauss)goto 10261
      call gam(linkg,dirivg,weigtg,devg,linvg)
      call info
      call dumps(dirivg)
      goto10251
10261 if(mod .ne. binom)goto 10271
      call gam(linkb,dirivb,weigtb,devb,linvb)
      call info
      call dumps(dirivb)
      goto10251
10271 if(mod .ne. cox)goto 10281
      nq=n
      call gam(linkg,dirivg,weigtg,devc,linvg)
      call info
      call dumps(dirivg)
      goto10251
10281 if(mod .ne. user)goto 10291
      call gam(linku,dirivu,weigtu,devu,linvu)
      call info
      call dumps(dirivu)
      goto10251
10291 if(mod .ne. mult)goto 10301
      call multam
      call minfo
      call dumps(dirivb)
      write(6,10310)
10310 format(   38h want fitted cumulative probabilities )
      read(5,10050)ans
      if(ans .ne. yes)goto 10331
      call dumpm
10331 continue
      goto10251
10301 if(mod .ne. match)goto 10341
      call gam(linkmc,dirivb,weigtb,devmc,linvmc)
      call info
      call dumps(dirivb)
      goto 10351
10341 continue
      write(6,10360)mod
10360 format(   32h sorry, we dont have model type ,a1)
10351 continue
10251 continue
      goto 10021
10022 continue
      write(6,10370)
10370 format(   15h another model?)
      read(5,10050)answer
10050 format(a1)
      goto 10011
10012 continue
      stop
      end
      subroutine model
      real dumlab(30),modlab(20)
      integer indat(20),indmod(20)
      common/extra/dumlab,modlab,indat,indmod
      real mult,group,ym(1000,15)
      real thetak(15)
      integer mcat,indrm(15)
      common/mltcmn/mult,group,thetak,ym,mcat,indrm
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
     *),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
      integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
      common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
     *icensq,irindq,iriskq,nq,kq
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
      integer sltag(1000)
      logical vlink,vflag
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
      real match
      integer indms,kms(1000)
      common/mtchco/match,indms,kms
      real v(20),more,quit
      logical ok
      data quit/    1hq/
      ok=.false.
      vlink=.false.
10381 if(ok) goto 10382
      write(6,10390)
10390 format(   48h which model...binomial(b), gaussian(g), cox(c),,/,
     *   17h multinomial(m), ,     44h matched case control(r), user(u)
     *or quit(q))
      read(5,10400)mod
10400 format(a1)
      if((mod .ne. gauss) .and. ((mod .ne. binom) .and. ((mod .ne. cox)
     *.and. ((mod .ne. user) .and. ((mod .ne. mult) .and. (mod .ne. matc
     *h))))))goto 10421
      ok=.true.
      goto10411
10421 if(mod .ne. quit)goto 10431
      stop
      goto 10441
10431 continue
      write(6,10450)mod
10450 format(   32h sorry, we dont have model type ,a1)
10441 continue
10411 continue
      goto 10381
10382 continue
      esttyp=score
      if(mod .ne. mult)goto 10471
      call multin
      return
10471 continue
      j=1
      goto 10483
10481 j=j+1
10483 if((j).gt.(maxp))goto 10482
      indat(j)=0
      write(6,10490)j,(labels(k,j),k=1,10)
10490 format(i3,    5h ....,10a4)
      goto 10481
10482 continue
      ir=maxp+1
10501 if((ir .le. maxp) .and. (ir .gt. 0)) goto 10502
      write(6,10510)
10510 format(   19h response variable?)
      read(5,*)ir
      goto 10501
10502 continue
      indat(ir)=2
      indy=ir
      do 10521 i=1,n
      y(i)=datam(i,ir)
10521 continue
10522 continue
      ir=maxp+1
10531 if((ir .le. maxp) .and. (ir .ne. indy)) goto 10532
      write(6,10540)
10540 format(   55h give index of weight variable  0 is no weight variab
     *le)
      read(5,*)ir
      goto 10531
10532 continue
      indw=ir
      if(ir .gt. 0)goto 10561
      do 10571 i=1,n
      w(i)=1.0
10571 continue
10572 continue
      goto 10581
10561 continue
      totw=0.0
      do 10591 i=1,n
      w(i)=abs(datam(i,ir))
      totw=totw+w(i)
10591 continue
10592 continue
      do 10601 i=1,n
      w(i)=w(i)*n/totw
10601 continue
10602 continue
      indat(ir)=3
10581 continue
10551 continue
      nmiss=0
      do 10611 i=1,n
      if(y(i) .ne. 99999)goto 10631
      w(i)=0.0
      y(i)=0.0
      nmiss=nmiss+1
10631 continue
10611 continue
10612 continue
      if(nmiss .le. 0)goto 10651
      ratmis=(100*nmiss)/n
      write(6,10660)nmiss,ratmis
10660 format (i4,   34h responses missing (=> 99999), or ,f4.1,   14h% o
     *f the data;,/,     31h their weights will be set to 0)
10651 continue
      p=0
      do 10671 i=1,n
      ni(i)=1.0
10671 continue
10672 continue
      if(mod .ne. binom)goto 10691
      write(6,10700)
10700 format(   44h y out of n ... give index of n  0 makes n=1)
      ir=maxp+1
10711 if((ir .le. maxp) .and. (ir .ne. indy)) goto 10712
      read(5,*)ir
      goto 10711
10712 continue
      if(ir .gt. 0)goto 10731
      indni=0
      goto 10741
10731 continue
      indni=ir
      do 10751 i=1,n
      ni(i)=datam(i,ir)
      y(i)=y(i)/ni(i)
10751 continue
10752 continue
      indat(ir)=4
10741 continue
10721 continue
10691 continue
      if(mod .ne. match)goto 10771
      write(6,10780)
10780 format(   36h give index of matched set indicator)
      ir=maxp+1
10791 if((ir .le. maxp) .and. (ir .ne. indy)) goto 10792
      read(5,*)ir
      goto 10791
10792 continue
      if(ir .le. 0)goto 10811
      indms=ir
      do 10821 i=1,n
      kms(i)=datam(i,ir)
10821 continue
10822 continue
      indat(ir)=6
10811 continue
10771 continue
      if(mod .ne. cox)goto 10841
      write(6,10850)
10850 format(   33h give index of censoring variable)
      read(5,*) ic
      indat(ic)=5
      do 10861 i=1,n
      icensq(i)=datam(i,ic)
10861 continue
10862 continue
      write(6,10870)
10870 format(   55h there should be no censored obs before the 1st failu
     *re)
      do 10881 i=1,n
      tq(i)=y(i)
10881 continue
10882 continue
      call mysort(tq,tagy,1,n)
      do 10891 i=1,n
      tq(i)=y(tagy(i))
      aq(i)=-1.0*icensq(tagy(i))
10891 continue
10892 continue
      j=1
10900 continue
      jj=j
10911 if(tq(jj+1) .ne. tq(j)) goto 10912
      jj=jj+1
      goto 10911
10912 continue
      if(jj.gt.j) call psort(aq,tagy,j,jj)
      j=jj+1
      if(j.lt.n) go to 10900
      do 10921 i=1,n
      tq(i)=y(i)
      aq(i)=icensq(i)
      bq(i)=w(i)
      do 10931 j=1,maxp
      xt(i,j)=datam(i,j)
10931 continue
10932 continue
10921 continue
10922 continue
      do 10941 i=1,n
      y(i)=tq(tagy(i))
      icensq(i)=aq(tagy(i))
      w(i)=bq(tagy(i))
      do 10951 j=1,maxp
      datam(i,j)=xt(tagy(i),j)
10951 continue
10952 continue
10941 continue
10942 continue
10841 continue
      return
      end
      subroutine modtm
      real dumlab(30),modlab(20)
      integer indat(20),indmod(20)
      common/extra/dumlab,modlab,indat,indmod
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
      integer sltag(1000)
      logical vlink,vflag
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
10400 format(a1)
      logical ok
      vlink=.false.
      write(6,10960)
10960 format(   45h do you want to change the thresh-holds (y/n))
      read(5,10400)answer
      if(answer .ne. yes)goto 10981
      write(6,10990)thresh
10990 format(   23h change gam threshold (,f10.8,    6h y/n ))
      read(5,10400)answer
      if(answer .ne. yes)goto 11011
      write(6,11020)
11020 format(   26h give new gam thresh-hold )
      read(5,*)thresh
11011 continue
      write(6,11030)thbak
11030 format(   28h change back-fit threshold (,f10.8,    6h y/n ))
      read(5,10400)answer
      if(answer .ne. yes)goto 11051
      write(6,11060)
11060 format(   31h give new back-fit thresh-hold )
      read(5,*)thbak
11051 continue
10981 continue
      ok=.false.
11071 continue
      write(6,11080)
11080 format(   40h the following variables may be modified)
      j=1
      goto 11093
11091 j=j+1
11093 if((j).gt.(p))goto 11092
      if(.not.(cross(j)))goto 11111
      spans(j)=0.0
11111 continue
      if(indmod(j) .gt. 3)goto 11131
      write(6,11140)index(j),spans(j),modlab(j),(labels(k,index(j)),k=1,
     *10)
11140 format (i3,   10h ....span=,f5.2,    6h .... ,11a4)
11131 continue
      goto 11091
11092 continue
      write(6,11150)
11150 format(   44h index of covariate to be modified (0= quit))
      read(5,*)ir
      if(ir .ne. 0)goto 11171
      return
11171 continue
      if(indat(ir) .eq. 1)goto 11191
      write(6,11200)
11200 format(   21h naughty, naughty....)
      goto 11211
11191 continue
      do 11221 j=1,p
      if(index(j) .ne. ir)goto 11241
      goto 11222
11241 continue
11221 continue
11222 continue
      write(6,11250)ir,spans(j)
11250 format (   10h variable i3,   10h has span=,f5.2,/,     62h give n
     *ew span in interval (-?,1)   (1 => linear, - => spline))
      modlab(j)=dumlab(1)
      indmod(j)=1
      read(5,*)spans(j)
      cross(j)=.false.
      if(spans(j) .ne. 1.0)goto 11271
      indmod(j)=2
      modlab(j)=dumlab(22)
11271 continue
      if((spans(j) .ne. 0.0) .and. (spans(j) .ne. -999))goto 11291
      cross(j)=.true.
      indmod(j)=3
      modlab(j)=dumlab(21)
11291 continue
11211 continue
11181 continue
      goto 11071
11072 continue
      return
      end
      subroutine name7
      character*60 dsname
      common/dnm/dsname,ilen
      character*90 devn
      character*6 devn2
      close(7)
      devn=dsname
      devn(ilen:ilen)=    1h/
      devn(ilen+1:ilen+6)=    6hanodev
      open(7,file=devn)
      rewind 7
      return
      end
C     mortran 2.79 (reserved keyword macros of 09/28/81)              
      subroutine gam(link,diriv,weight,dev,linv)                                
      external link,diriv,weight,dev,linv                                       
      real linv,diriv,weight,dev,scrat1(1000),scrat2(1000)                      
      real scrat3(1000),scrat4(1000,3)                                          
      integer q                                                                 
      real datam(1000,20),labels(10,20)                                         
      integer maxp,n                                                            
      common /datac/datam,labels,maxp,n                                         
      real muhat(1000),etahat(1000)                                             
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak                   
      logical verbos,quiet                                                      
      real s0,devnll                                                            
      integer maxit,nitout,numsm                                                
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d        
     *evnll, maxit,verbos,quiet,nitout,numsm
      real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000        
     *),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
      integer icensq(1000),irindq(1000),iriskq(1000),nq,kq                      
      common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,         
     *icensq,irindq,iriskq,nq,kq
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)            
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score            
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc         
      logical cross(20)                                                         
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod        
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real z(1000),wz(1000)                                                     
      common /wvc/z,wz                                                          
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc           
      integer sltag(1000)                                                       
      logical vlink,vflag                                                       
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc          
      itss=30                                                                   
      call init(linv)                                                           
      if(mod .ne. cox)goto 10021                                                
      call initcx                                                               
10021 continue                                                                  
      devnew=dev(n,muhat,y,w,ni)                                                
      devnll=devnew                                                             
      devold=10*devnew+10                                                       
      if(mod .ne. cox)goto 10041                                                
      call stpar                                                                
10041 continue                                                                  
      nit=0                                                                     
      numsm=0                                                                   
      if(.not.(.not.vlink))goto 10061                                           
10071 if(abs((devold-devnew)/devold) .le. thresh .or. nit .ge. maxit) go        
     *to 10072
      nit=nit+1                                                                 
      if(mod .eq. cox)goto 10091                                                
      i=1                                                                       
      goto 10103                                                                
10101 i=i+1                                                                     
10103 if((i).gt.(n))goto 10102                                                  
      z(i)=etahat(i)+(y(i)-muhat(i))*diriv(muhat(i))                            
      wz(i)=weight(w(i),muhat(i),ni(i))                                         
      goto 10101                                                                
10102 continue                                                                  
      goto 10111                                                                
10091 continue                                                                  
      call calcz                                                                
10111 continue                                                                  
10081 continue                                                                  
      call bakfit(n,p,z,x,wz,tag,spans,dofs,slopes,cross,etahat,  s,s0,r        
     *ss,thbak,nitb,verbos,quiet,itss)
      numsm=numsm+nitb                                                          
      call link(n,etahat,muhat)                                                 
      devold=devnew                                                             
      devnew=dev(n,muhat,y,w,ni)                                                
      if(.not.(.not.quiet))goto 10131                                           
      write(6,10140)nit,devnew                                                  
10140 format(/10x,   22houter loop (gam) nit= ,  i3,    5h dev=,f16.5,//        
     *)
10131 continue                                                                  
      goto 10071                                                                
10072 continue                                                                  
10061 continue                                                                  
      if(.not.(vlink))goto 10161                                                
10171 if(abs((devold-devnew)/devold) .le. thresh .or. nit .ge. 5) goto 1        
     *0172
      devold=devnew                                                             
      devnei=devold                                                             
      devoli=devnei*10                                                          
      nitin=0                                                                   
10181 if(abs((devoli-devnei)/devoli) .le. thresh .or. nitin .ge. maxit)         
     *goto 10182
      nitin=nitin+1                                                             
      i=1                                                                       
      goto 10193                                                                
10191 i=i+1                                                                     
10193 if((i).gt.(n))goto 10192                                                  
      fdirin=fpinv(etahat(i))                                                   
      z(i)=etahat(i)+  (y(i)-muhat(i))*diriv(muhat(i))*fdirin                   
      wz(i)=weight(w(i),muhat(i),ni(i))/(fdirin*fdirin)                         
      goto 10191                                                                
10192 continue                                                                  
      call bakfit(n,p,z,x,wz,tag,spans,dofs,slopes,cross,etahat,  s,s0,r        
     *ss,thbak,nitb,verbos,quiet,itss)
      numsm=numsm+nitb                                                          
      i=1                                                                       
      goto 10203                                                                
10201 i=i+1                                                                     
10203 if((i).gt.(n))goto 10202                                                  
      sleta(i)=flink(etahat(i))                                                 
      goto 10201                                                                
10202 continue                                                                  
      call link(n,etahat,muhat)                                                 
      devoli=devnei                                                             
      devnei=dev(n,muhat,y,w,ni)                                                
      if(.not.(.not.quiet))goto 10221                                           
      write(6,10230)nitin,devnei                                                
10230 format(/10x,   27hinner loop  eta (gam) nit= ,  i3,    5h dev=,f16        
     *.5,//)
10221 continue                                                                  
      goto 10181                                                                
10182 continue                                                                  
      if(.not.(.not.vflag))goto 10251                                           
      slinc=devnei                                                              
10251 continue                                                                  
      devnes=devnei                                                             
      devols=devnes*100                                                         
      nitsin=0                                                                  
      vflag=.true.                                                              
      call mysort(etahat,sltag,1,n)                                             
      do 10261 i=1,n                                                            
      argf(i)=etahat(sltag(i))                                                  
10261 continue                                                                  
10262 continue                                                                  
10271 if(abs((devols-devnes)/devols) .le. thresh .or. nitsin .ge. 6) got        
     *o 10272
      nitsin=nitsin+1                                                           
      i=1                                                                       
      goto 10283                                                                
10281 i=i+1                                                                     
10283 if((i).gt.(n))goto 10282                                                  
      z(i)=sleta(i) +(y(i)-muhat(i))*diriv(muhat(i))                            
      wz(i)=weight(w(i),muhat(i),ni(i))                                         
      goto 10281                                                                
10282 continue                                                                  
      call smooth(etahat,z,wz,sltag,0.5,sldof,n,n,  .false.,sleta,sv,rss        
     *,quiet,scrat1)
      do 10291 i=1,n                                                            
      slink(i)=sleta(sltag(i))+sv                                               
10291 continue                                                                  
10292 continue                                                                  
      call montne(slink,n)                                                      
      call link(n,slink,scrat1)                                                 
      i=1                                                                       
      goto 10303                                                                
10301 i=i+1                                                                     
10303 if((i).gt.(n))goto 10302                                                  
      sleta(sltag(i))=slink(i)                                                  
      muhat(sltag(i))=scrat1(i)                                                 
      goto 10301                                                                
10302 continue                                                                  
      devols=devnes                                                             
      devnes=dev(n,muhat,y,w,ni)                                                
      if(.not.(.not.quiet))goto 10321                                           
      write(6,10330)nitsin,devnes                                               
10330 format(/10x,   30hinner loop  s(eta) (gam) nit= ,  i3,    5h dev=,        
     *f16.5,//)
10321 continue                                                                  
      goto 10271                                                                
10272 continue                                                                  
      call der(n,argf,slink,w,0.001,sldir,scrat4)                               
      dif=slink(n)-slink(1)                                                     
      idif=1                                                                    
      if(dif .ge. 0)goto 10351                                                  
      idif=-1                                                                   
10351 continue                                                                  
      do 10361 i=1,n                                                            
      if(sldir(i)*idif .ge. 0.01)goto 10381                                     
      sldir(i)=idif*0.01                                                        
10381 continue                                                                  
10361 continue                                                                  
10362 continue                                                                  
      devold=devnew                                                             
      devnew=devnes                                                             
      nit=nit+1                                                                 
      if(.not.(.not.quiet))goto 10401                                           
      write(6,10410)nit,devnew                                                  
10410 format(10x,30(    1h_),/10x,   22houter loop (gam) nit= ,  i3,            
     *5h dev=,f16.5,/10x,30(    1h-)/)
10401 continue                                                                  
      goto 10171                                                                
10172 continue                                                                  
10161 continue                                                                  
      devian=devnew                                                             
      nitout=nit                                                                
      return                                                                    
      end                                                                       
C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine hat(weight)
      external weight
      real weight
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      integer neffob(20)
      common/neff/neffob
      real z(1000),wz(1000)
      common /wvc/z,wz
      integer bacnit
      common /bacn/bacnit
      character*60 dsname
      common/dnm/dsname,ilen
      logical parm1,parm2
      real zero
      integer its
      real*8 vf(1000,20),vyhat(1000),trudof,temp1,temp2
      real vs(1000,20),veta(1000),vs0,vrss,vscale(1000)
      real varco(20,20)
      integer indco(20),numco
      character*10 plotn
      ans=0
10011 if(ans .eq. yes .or. ans .eq. no) goto 10012
      write(6,10020)
10020 format(   48h this program computes standard deviation curves,/,
     *   53h for the estimated functions, as well as a covariance,/,
     * 49h matrix for the estimated constants in the model.,/,/,     55h
     * the procedure takes n x (time for final gam iteration),/,     47h
     * ( as you ll find out), so dont use carelessly.,/,     26h do you
     *wish to continue ?)
      read(5,10030)ans
      if(ans .ne. yes)goto 10051
      goto 10012
10051 continue
      if(ans .ne. no)goto 10071
      return
10071 continue
      goto 10011
10012 continue
      numco=1
      j=1
      goto 10083
10081 j=j+1
10083 if((j).gt.(p))goto 10082
      if(spans(j) .ne. 1)goto 10101
      numco=numco+1
      indco(numco)=j
10101 continue
      goto 10081
10082 continue
      do 10111 j=1,numco
      do 10121 jj=1,j
      varco(j,jj)=0.0
10121 continue
10122 continue
10111 continue
10112 continue
      parm1=.true.
      parm2=.false.
      its=bacnit
      if(bacnit .le. 10)goto 10141
      its=10
10141 continue
      zero=0.0
      trudof=0d0
      vtotw=0.0
      i=1
      goto 10153
10151 i=i+1
10153 if((i).gt.(n))goto 10152
      vyhat(i)=0d0
      wz(i)=weight(w(i),muhat(i),ni(i))
      vtotw=vtotw+wz(i)
      j=1
      goto 10163
10161 j=j+1
10163 if((j).gt.(p))goto 10162
      vf(i,j)=0d0
      goto 10161
10162 continue
      goto 10151
10152 continue
      ii=1
      goto 10173
10171 ii=ii+1
10173 if((ii).gt.(n))goto 10172
      if(wz(ii) .ne. 0.0)goto 10191
      goto 10171
10191 continue
      vs0=wz(ii)/vtotw
      i=1
      goto 10203
10201 i=i+1
10203 if((i).gt.(n))goto 10202
      z(i)=0.0
      veta(i)=vs0
      j=1
      goto 10213
10211 j=j+1
10213 if((j).gt.(p))goto 10212
      vs(i,j)=0.0
      goto 10211
10212 continue
      goto 10201
10202 continue
      z(ii)=1
      call bakfit(n,p,z,x,wz,tag,spans,dofs,slopes,cross,veta,  vs,vs0,v
     *rss,zero,nitb,parm2,parm1,its)
      j=1
      goto 10223
10221 j=j+1
10223 if((j).gt.(numco))goto 10222
      if(j .ne. 1)goto 10241
      sloup1=vs0
      goto 10251
10241 continue
      sloup1=slopes(indco(j))
10251 continue
10231 continue
      jj=1
      goto 10263
10261 jj=jj+1
10263 if((jj).gt.(j))goto 10262
      if(jj .ne. 1)goto 10281
      sloup2=vs0
      goto 10291
10281 continue
      sloup2=slopes(indco(jj))
10291 continue
10271 continue
      varco(j,jj)=varco(j,jj)+sloup1*sloup2/wz(ii)
      goto 10261
10262 continue
      goto 10221
10222 continue
      i=1
      goto 10303
10301 i=i+1
10303 if((i).gt.(n))goto 10302
      j=1
      goto 10313
10311 j=j+1
10313 if((j).gt.(p))goto 10312
      temp1=vs(i,j)
      vf(i,j)=vf(i,j)+temp1*temp1/wz(ii)
      goto 10311
10312 continue
      temp2=veta(i)
      vyhat(i)=vyhat(i)+temp2*temp2/wz(ii)
      goto 10301
10302 continue
      write(6,10320)ii,n,nitb
10320 format(   11h completed ,i3,    4h of ,i3,    5h with,i3,   11h it
     *erations)
      goto 10171
10172 continue
      i=1
      goto 10333
10331 i=i+1
10333 if((i).gt.(n))goto 10332
      trudof=trudof+vyhat(i)*wz(i)
      goto 10331
10332 continue
      scale=1.0
      write(6,10340)trudof
10340 format(   13h true dof is ,f10.2)
      if(mod .ne. gauss)goto 10361
      scale=sqrt(devian/(n-trudof))
10361 continue
      j=1
      goto 10373
10371 j=j+1
10373 if((j).gt.(numco))goto 10372
      jj=1
      goto 10383
10381 jj=jj+1
10383 if((jj).gt.(j))goto 10382
      varco(j,jj)=varco(j,jj)*scale*scale
      varco(jj,j)=varco(j,jj)
      goto 10381
10382 continue
      goto 10371
10372 continue
      write(6,10390)
10390 format(/   31h covariance matrix of constants,/)
      write(6,10400)(varco(1,j),j=1,numco)
10400 format(   25h intercept s0 -----------,10(e9.2,1x))
      jj=2
      goto 10413
10411 jj=jj+1
10413 if((jj).gt.(numco))goto 10412
      write(6,10420)(labels(k,index(indco(jj))),k=1,6),(varco(jj,j),j=1,
     *numco)
10420 format (1x,6a4,10(e9.2,1x))
      goto 10411
10412 continue
      write(6,10430)
10430 format(   33h want to dump  plots with +-2sd? )
      read(5,10030)ans
      if(ans .ne. no)goto 10451
      return
10451 continue
10030 format(a1)
      r2=(devnll-devian)*100.0/devnll
      j=1
      goto 10463
10461 j=j+1
10463 if((j).gt.(p))goto 10462
      write(6,10470)(labels(k,index(j)),k=1,10)
10470 format (   16h do you want ...,10a4)
      read(5,10030)ans
      if(ans .ne. no)goto 10491
      goto 10461
10491 continue
      call nplot(plotn)
      totsm=0.0
      do 10501 i=1,neffob(j)
      totsm=totsm+s(tag(i,j),j)
10501 continue
10502 continue
      totsm=totsm/neffob(j)
      oldx=-9999
      i=1
      goto 10513
10511 i=i+1
10513 if((i).gt.(neffob(j)))goto 10512
      ii=tag(i,j)
      if(x(ii,j) .eq. oldx)goto 10531
      sput=s(ii,j)-totsm
      if(vf(ii,j) .le. 0.0)goto 10551
      sputd=sqrt(vf(ii,j))
      goto 10561
10551 continue
      sputd=0.0
10561 continue
10541 continue
      sputd=2*sputd*scale
      sputl=sput-sputd
      sputu=sput+sputd
      write(2,*)x(ii,j),sput,sputl,sputu
      oldx=x(ii,j)
10531 continue
      goto 10511
10512 continue
      write(6,10570)spans(j),dofs(j),trudof,devian,r2
      goto 10461
10462 continue
10580 format(    7h data a,i1,   27h; input x yhat yhatl yhatu;)
10590 format(   34h title .c=blue estimated function;  /,   49hproc gplo
     *t; plot  yhat*x yhatl*x yhatu*x/overlay;,/,     33hsymbol1 i=join
     *c=green v=diamond;,/,     28hsymbol2 i=join c=red v=star;,/,     2
     *8hsymbol3 i=join c=red v=star;)
10570 format(    5hspan=,f5.2,    5h dof=,f6.2,    19h true dof (overall
     *),f6.2,      9hdeviance=, f8.2,    4h r2=,f8.2,    2h%;)
10600 format(    8hlabel x=,10a4,   13hyhat=s;cards;)
      return
      end
C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine helpg
      real a(20),paus,halt
      data paus/    4hpaus/, halt/    4hstop/
      data quit/    1hq/
      open(12,file=8hhelpunix)
      rewind 12
10011 continue
      read(12,10020)a
10020 format(20a4)
      if(a(1) .ne. halt)goto 10041
      return
      goto10031
10041 if(a(1) .ne. paus)goto 10051
      write(6,10060)
10060 format(    40h type q for quit, anything else for more)
      read(5,10070)ans
10070 format(a1)
      if(ans .ne. quit)goto 10091
      return
10091 continue
      goto 10101
10051 continue
      write(6,10110)a
10101 continue
10031 continue
10110 format(    1h ,20a4)
      goto 10011
10012 continue
      return
      end
C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine info
      real dumlab(30),modlab(20)
      integer indat(20),indmod(20)
      common/extra/dumlab,modlab,indat,indmod
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
      integer sltag(1000)
      logical vlink,vflag
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real match
      integer indms,kms(1000)
      common/mtchco/match,indms,kms
      totdof=1.0
      do 10011 j=1,p
      totdof=totdof+dofs(j)
10011 continue
10012 continue
      if(mod .ne. cox)goto 10031
      totdof=totdof-1.0
10031 continue
      errn=n-totdof
      scale=1.0
      if(mod .ne. gauss)goto 10051
      scale=devian/errn
10051 continue
      if(mod .ne. binom)goto 10071
      write(6,10080)
10080 format(//   35h logistic additive regression model/)
      write(7,10090)
10090 format(//   35h logistic additive regression model/)
      goto10061
10071 if(mod .ne. gauss)goto 10101
      write(6,10110)
10110 format(//   35h gaussian additive regression model/)
      write(7,10120)
10120 format(//   35h gaussian additive regression model/)
      goto10061
10101 if(mod .ne. match)goto 10131
      write(6,10140)
10140 format(//   47h matched case-control additive regression model/)
      write(7,10150)
10150 format(//   47h matched case-control additive regression model/)
      goto10061
10131 if(mod .ne. cox)goto 10161
      write(6,10170)
10170 format(//   50h additive cox proportional hazard regression model/
     *)
      write(7,10180)
10180 format(//   50h additive cox proportional hazard regression model/
     *)
      goto 10191
10161 continue
10191 continue
10061 continue
      if(esttyp .ne. score)goto 10211
      write(6,10220)
10220 format(//   28h estimation by local scoring/)
      write(7,10230)
10230 format(//   28h estimation by local scoring/)
      goto 10241
10211 continue
      write(6,10250)
10250 format(//   31h estimation by local likelihood/)
      write(7,10260)
10260 format(//   31h estimation by local likelihood/)
10241 continue
10201 continue
      write(6,10270)(labels(k,indy),k=1,10)
10270 format (/   23h response variable ....,10a4)
      write(7,10280)(labels(k,indy),k=1,10)
10280 format (/   23h response variable ....,10a4)
      if(mod .ne. binom)goto 10301
      if(indni .eq. 0)goto 10321
      write(6,10330)(labels(k,indni),k=1,10)
10330 format (   23h binomial denominator..,10a4)
      write(7,10340)(labels(k,indni),k=1,10)
10340 format (   23h binomial denominator..,10a4)
      goto 10351
10321 continue
      write(6,10360)
10360 format(   20h bernoulli response )
      write(7,10370)
10370 format(   20h bernoulli response )
10351 continue
10311 continue
10301 continue
      write(6,10380)devian,nitout,numsm
10380 format (   12h deviance = ,f14.3,   14h # iterations=,i3,   20h #
     *smooths/variable=,i4)
      write(7,10390)devian,nitout,numsm
10390 format (   12h deviance = ,f14.3,   14h # iterations=,i3,   20h #
     *smooths/variable=,i4)
      a=devian/n
      write(6,10400)a
10400 format(   20h average deviance = , f12.3)
      write(7,10410)a
10410 format(   20h average deviance = , f12.3)
      write(6,10420)errn,scale
10420 format(   18h dof of deviance  ,f6.2,     18h scale estimate   ,f1
     *0.3)
      write(7,10430)errn,scale
10430 format(   18h dof of deviance  ,f6.2,     18h scale estimate   ,f1
     *0.3)
      r2=100*(devnll-devian)/devnll
      write(6,10440)r2,devnll
10440 format (   12h r square = ,f5.2,   24h% of a null deviance of ,f14
     *.3)
      write(7,10450)r2,devnll
10450 format (   12h r square = ,f5.2,   24h% of a null deviance of ,f14
     *.3)
      write(6,10460)
10460 format(/,    65h variable    span   dof     slope                n
     *ame            ,/,    65h --------    ----   ---     ------    ---
     *------------------------)
      write(7,10470)
10470 format(/,    65h variable    span   dof     slope                n
     *ame            ,/,    65h --------    ----   ---     ------    ---
     *------------------------)
      if(mod .eq. cox)goto 10491
      write(6,10500)s0
10500 format (   23h   0          --     1 ,f10.4,   26h   s0---the inte
     *rcept term)
      write(7,10510)s0
10510 format (   23h   0          --     1 ,f10.4,   26h   s0---the inte
     *rcept term)
10491 continue
      j=1
      goto 10523
10521 j=j+1
10523 if((j).gt.(p))goto 10522
      if(spans(j) .lt. 1.0)goto 10541
      write(6,10550)index(j),spans(j), dofs(j),slopes(j),modlab(j),(labe
     *ls(k,index(j)),k=1,10)
10550 format (    1h ,i3,7x,f6.2,2x, f5.3, f9.4,3x,11a4)
      write(7,10560)index(j),spans(j), dofs(j),slopes(j),modlab(j),(labe
     *ls(k,index(j)),k=1,10)
10560 format (    1h ,i3,7x,f6.2,2x, f5.3, f9.4,3x,11a4)
      goto 10571
10541 continue
      write(6,10580)index(j),spans(j),dofs(j),modlab(j),(labels(k,index(
     *j)),k=1,10)
10580 format (    1h ,i3,7x,f6.2,2x, f6.3,    8h  smooth,3x,11a4)
      write(7,10590)index(j),spans(j),dofs(j),modlab(j),(labels(k,index(
     *j)),k=1,10)
10590 format (    1h ,i3,7x,f6.2,2x, f6.3,    8h  smooth,3x,11a4)
10571 continue
10531 continue
      goto 10521
10522 continue
      write(6,10600)totdof
10600 format (   29h                    -----    ,/,     19h
     *       ,f6.3)
      write(7,10610)totdof
10610 format (   29h                    -----    ,/,     19h
     *       ,f6.3)
      if(.not.(vlink))goto 10631
      slinc=slinc-devian
      write(6,10640)sldof,slinc
10640 format(   21h variable link (dof= ,f5.2,    8h) caused,  /   21h d
     *eviance to drop by ,f10.3)
      write(7,10650)sldof,slinc
10650 format(   21h variable link (dof= ,f5.2,    8h) caused,  /   21h d
     *eviance to drop by ,f10.3)
10631 continue
      return
      end
      real function devu(n,fits,y,w,ni)
      real fits(1),y(1),w(1),ni(1)
      rss=0.0
      do 10661 i=1,n
      rat=y(i)/abs(fits(i))
      if(rat .gt. 0.0)goto 10681
      write(6,10690)
10690 format(     6hshriek)
      stop
10681 continue
      ratl=alog(rat)
      rss=rss+w(i)*(rat-ratl -1)*2.0
10661 continue
10662 continue
      devu=rss
      return
      end
      real function linku(n,etahat,muhat)
      linku=etahat
      return
      end
      real function linvu(muhat)
      real muhat
      linvu=muhat
      return
      end
      real function dirivu(muhat)
      real muhat
      dirivu=1.0
      return
      end
      real function weigtu(w,muhat,ni)
      real w,muhat,ni
      t=muhat/700.00
      t=t*t
      if(t .gt. 0.0001)goto 10711
      t=0.0001
10711 continue
      weigtu=1.0/t
      return
      end
      real function fpinv(arg)
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
      integer sltag(1000)
      logical vlink,vflag
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      if(.not.(vflag))goto 10731
      fpinv=1.0/yinter(argf,sldir,n,arg)
      goto 10741
10731 continue
      fpinv=1.0
10741 continue
10721 continue
      return
      end
      real function flink(arg)
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
      integer sltag(1000)
      logical vlink,vflag
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      if(.not.(vflag))goto 10761
      flink=yinter(argf,slink,n,arg)
      goto 10771
10761 continue
      flink=arg
10771 continue
10751 continue
      return
      end
      subroutine montne(x,n)
      real x(n)
      integer bb,eb,br,er,bl,el
      bb=0
      eb=bb
10780 continue
10781 if(eb.ge.n) goto 10782
      bb=eb+1
      eb=bb
10791 if(eb .ge. n) goto 10792
      if(x(bb) .ne. x(eb+1))goto 10811
      eb=eb+1
      goto 10821
10811 continue
      goto 10792
10821 continue
10801 continue
      goto 10791
10792 continue
10831 continue
      if(eb .ge. n)goto 10851
      if(x(eb) .le. x(eb+1))goto 10871
      br=eb+1
      er=br
10881 if(er .ge. n) goto 10882
      if(x(er+1) .ne. x(br))goto 10901
      er=er+1
      goto 10911
10901 continue
      goto 10882
10911 continue
10891 continue
      goto 10881
10882 continue
      pmn=(x(bb)*(eb-bb+1)+x(br)*(er-br+1))/(er-bb+1)
      eb=er
      do 10921 i=bb,eb
      x(i)=pmn
10921 continue
10922 continue
10871 continue
10851 continue
      if(bb.le.1) goto 10781
      if(x(bb-1).le.x(bb)) goto 10781
      bl=bb-1
      el=bl
10931 if(bl .le. 1) goto 10932
      if(x(bl-1) .ne. x(el))goto 10951
      bl=bl-1
      goto 10961
10951 continue
      goto 10932
10961 continue
10941 continue
      goto 10931
10932 continue
      pmn=(x(bb)*(eb-bb+1)+x(bl)*(el-bl+1))/(eb-bl+1)
      bb=bl
      do 10971 i=bb,eb
      x(i)=pmn
10971 continue
10972 continue
      goto 10831
10832 continue
      goto 10781
10782 continue
      return
      end
      subroutine der (n,x,s,w,fdel,d,sc)
      real x(n),s(n),w(n),d(n),sc(n,3)
      integer bl,el,bc,ec,br,er
      if(x(n) .gt. x(1))goto 10991
      do 11001 j=1,n
      d(j)=0.0
11001 continue
11002 continue
      return
10991 continue
      i=n/4
      j=3*i
      scale=x(j)-x(i)
11011 if(scale.gt.0.0) goto 11012
      if(j.lt.n) j=j+1
      if(i.gt.1) i=i-1
      scale=x(j)-x(i)
      goto 11011
11012 continue
      del=fdel*scale*2.0
      do 11021 j=1,n
      sc(j,1)=x(j)
      sc(j,2)=s(j)
      sc(j,3)=w(j)
11021 continue
11022 continue
      call pool (n,sc,sc(1,2),sc(1,3),del)
      bc=0
      br=bc
      er=br
11031 continue
      br=er+1
      er=br
11041 if(er .ge. n) goto 11042
      if(sc(br,1) .ne. sc(er+1,1))goto 11061
      er=er+1
      goto 11071
11061 continue
      goto 11042
11071 continue
11051 continue
      goto 11041
11042 continue
      if(br .ne. 1)goto 11091
      bl=br
      el=er
      goto 11031
11091 continue
      if(bc .ne. 0)goto 11111
      bc=br
      ec=er
      do 11121 j=bl,el
      d(j)=(sc(bc,2)-sc(bl,2))/(sc(bc,1)-sc(bl,1))
11121 continue
11122 continue
      goto 11031
11111 continue
      do 11131 j=bc,ec
      d(j)=(sc(br,2)-sc(bl,2))/(sc(br,1)-sc(bl,1))
11131 continue
11132 continue
      if(er .ne. n)goto 11151
      do 11161 j=br,er
      d(j)=(sc(br,2)-sc(bc,2))/(sc(br,1)-sc(bc,1))
11161 continue
11162 continue
      goto 11032
11151 continue
      bl=bc
      el=ec
      bc=br
      ec=er
      goto 11031
11032 continue
      return
      end
      subroutine pool (n,x,y,w,del)
      real x(n),y(n),w(n)
      integer bb,eb,br,er,bl,el
      bb=0
      eb=bb
11171 if(eb.ge.n) goto 11172
      bb=eb+1
      eb=bb
11181 if(eb .ge. n) goto 11182
      if(x(bb) .ne. x(eb+1))goto 11201
      eb=eb+1
      goto 11211
11201 continue
      goto 11182
11211 continue
11191 continue
      goto 11181
11182 continue
      if(eb .ge. n)goto 11231
      if(x(eb+1)-x(eb) .ge. del)goto 11251
      br=eb+1
      er=br
11261 if(er .ge. n) goto 11262
      if(x(er+1) .ne. x(br))goto 11281
      er=er+1
      goto 11291
11281 continue
      goto 11262
11291 continue
11271 continue
      goto 11261
11262 continue
      if(x(er+1)-x(er).lt.x(eb+1)-x(eb))goto 11171
      eb=er
      pw=w(bb)+w(eb)
      px=(x(bb)*w(bb)+x(eb)*w(eb))/pw
      py=(y(bb)*w(bb)+y(eb)*w(eb))/pw
      do 11301 i=bb,eb
      x(i)=px
      y(i)=py
      w(i)=pw
11301 continue
11302 continue
11251 continue
11231 continue
11311 continue
      if(bb.le.1)goto 11312
      if(x(bb)-x(bb-1).ge.del)goto 11312
      bl=bb-1
      el=bl
11321 if(bl .le. 1) goto 11322
      if(x(bl-1) .ne. x(el))goto 11341
      bl=bl-1
      goto 11351
11341 continue
      goto 11322
11351 continue
11331 continue
      goto 11321
11322 continue
      bb=bl
      pw=w(bb)+w(eb)
      px=(x(bb)*w(bb)+x(eb)*w(eb))/pw
      py=(y(bb)*w(bb)+y(eb)*w(eb))/pw
      do 11361 i=bb,eb
      x(i)=px
      y(i)=py
      w(i)=pw
11361 continue
11362 continue
      goto 11311
11312 continue
      goto 11171
11172 continue
      return
      end
      real function yinter(z,fz,n,x)
      real z(n),fz(n),x
      integer n
      if(x .gt. z(1))goto 11381
      ii=1
      jj=2
11391 if(fz(jj) .ne. fz(ii)) goto 11392
      jj=jj+1
      goto 11391
11392 continue
      goto11371
11381 if(x .lt. z(n))goto 11401
      jj=n
      ii=n-1
11411 if(fz(ii) .ne. fz(jj)) goto 11412
      ii=ii-1
      goto 11411
11412 continue
      goto 11421
11401 continue
      ii=1
11431 if((z(ii).le.x).and.(x.le.z(ii+1))) goto 11432
      ii=ii+1
      goto 11431
11432 continue
      jj=ii+1
11421 continue
11371 continue
      if(z(ii) .ne. z(jj))goto 11451
      yinter=(fz(ii)+fz(jj))/2
      goto 11461
11451 continue
      slope=(fz(jj)-fz(ii))/(z(jj)-z(ii))
      yinter=fz(ii)+slope*(x-z(ii))
11461 continue
11441 continue
      return
      end
      subroutine mysort(input,tag,from,till)
      real input(1)
      integer tag(1),from,till
      real scrat(1000)
      i=from
      goto 11473
11471 i=i+1
11473 if((i).gt.(till))goto 11472
      tag(i)=i
      scrat(i)=input(i)
      goto 11471
11472 continue
      call sort(scrat,tag,from,till)
      return
      end
      subroutine psort(input,tag,from,till)
      real input(1)
      integer tag(1),from,till
      real scrat(1000)
      i=from
      goto 11483
11481 i=i+1
11483 if((i).gt.(till))goto 11482
      scrat(i)=input(i)
      goto 11481
11482 continue
      call sort(scrat,tag,from,till)
      return
      end
      subroutine sort (v,a,ii,jj)
c
c     puts into a the permutation vector which sorts v into
c     increasing order.  only elements from ii to jj are considered.
c     arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements
c     v is returned sorted
c     this is a modification of cacm algorithm #347 by r. c. singleton,
c     which is a modified hoare quicksort.
c
      dimension a(jj),v(1),iu(20),il(20)
      integer t,tt
      integer a
      real v
      m=1
      i=ii
      j=jj
 10   if (i.ge.j) go to 80
 20   k=i
      ij=(j+i)/2
      t=a(ij)
      vt=v(ij)
      if (v(i).le.vt) go to 30
      a(ij)=a(i)
      a(i)=t
      t=a(ij)
      v(ij)=v(i)
      v(i)=vt
      vt=v(ij)
 30   l=j
      if (v(j).ge.vt) go to 50
      a(ij)=a(j)
      a(j)=t
      t=a(ij)
      v(ij)=v(j)
      v(j)=vt
      vt=v(ij)
      if (v(i).le.vt) go to 50
      a(ij)=a(i)
      a(i)=t
      t=a(ij)
      v(ij)=v(i)
      v(i)=vt
      vt=v(ij)
      go to 50
 40   a(l)=a(k)
      a(k)=tt
      v(l)=v(k)
      v(k)=vtt
 50   l=l-1
      if (v(l).gt.vt) go to 50
      tt=a(l)
      vtt=v(l)
 60   k=k+1
      if (v(k).lt.vt) go to 60
      if (k.le.l) go to 40
      if (l-i.le.j-k) go to 70
      il(m)=i
      iu(m)=l
      i=k
      m=m+1
      go to 90
 70   il(m)=k
      iu(m)=j
      j=l
      m=m+1
      go to 90
 80   m=m-1
      if (m.eq.0) return
      i=il(m)
      j=iu(m)
 90   if (j-i.gt.10) go to 20
      if (i.eq.ii) go to 10
      i=i-1
 100  i=i+1
      if (i.eq.j) go to 80
      t=a(i+1)
      vt=v(i+1)
      if (v(i).le.vt) go to 100
      k=i
 110  a(k+1)=a(k)
      v(k+1)=v(k)
      k=k-1
      if (vt.lt.v(k)) go to 110
      a(k+1)=t
      v(k+1)=vt
      go to 100
      end
C     mortran 2.79 (reserved keyword macros of 09/28/81)              
      subroutine init(linv)                                                     
      external linv                                                             
      real linv                                                                 
      real muhat(1000),etahat(1000)                                             
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak                   
      logical verbos,quiet                                                      
      real s0,devnll                                                            
      integer maxit,nitout,numsm                                                
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d        
     *evnll, maxit,verbos,quiet,nitout,numsm
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)            
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score            
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc         
      logical cross(20)                                                         
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod        
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real datam(1000,20),labels(10,20)                                         
      integer maxp,n                                                            
      common /datac/datam,labels,maxp,n                                         
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc           
      integer sltag(1000)                                                       
      logical vlink,vflag                                                       
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc          
      real match                                                                
      integer indms,kms(1000)                                                   
      common/mtchco/match,indms,kms                                             
      theta=1.0                                                                 
      alpha=1.0                                                                 
      maxit=20                                                                  
      if(mod .ne. gauss .or. .not.(.not.vlink))goto 10021                       
      maxit=1                                                                   
10021 continue                                                                  
      s0=0                                                                      
      sn=0                                                                      
      i=1                                                                       
      goto 10033                                                                
10031 i=i+1                                                                     
10033 if((i).gt.(n))goto 10032                                                  
      s0=s0+w(i)*y(i)                                                           
      sn=sn+w(i)                                                                
      goto 10031                                                                
10032 continue                                                                  
      p0=s0/sn                                                                  
      i=1                                                                       
      goto 10043                                                                
10041 i=i+1                                                                     
10043 if((i).gt.(n))goto 10042                                                  
      muhat(i)=p0                                                               
      etahat(i)=linv(muhat(i))                                                  
      do 10051 j=1,p                                                            
      s(i,j)=0.0                                                                
10051 continue                                                                  
10052 continue                                                                  
      goto 10041                                                                
10042 continue                                                                  
      s0=etahat(1)                                                              
      if(.not.(vlink))goto 10071                                                
      do 10081 i=1,n                                                            
      sldir(i)=1.0                                                              
10081 continue                                                                  
10082 continue                                                                  
      vflag=.false.                                                             
10071 continue                                                                  
      return                                                                    
      end                                                                       
      real function devg(n,fits,y,w,ni)                                         
      real fits(1),y(1),w(1),ni(1)                                              
      rss=0.0                                                                   
      do 10091 i=1,n                                                            
      rss=rss+w(i)*(y(i)-fits(i))*(y(i)-fits(i))                                
10091 continue                                                                  
10092 continue                                                                  
      devg=rss                                                                  
      return                                                                    
      end                                                                       
      function devb(n,fits,y,w,ni)                                              
      real fits(n),y(n),ni(n),w(1)                                              
      dev=0.0                                                                   
      ent1=0                                                                    
      ent2=0                                                                    
      do 10101 i=1,n                                                            
      pr=fits(i)                                                                
      if(pr .ge. 0.0001)goto 10121                                              
      pr=0.0001                                                                 
10121 continue                                                                  
      if(pr .le. 0.9999)goto 10141                                              
      pr=0.9999                                                                 
10141 continue                                                                  
      if(((1.0-y(i))*y(i)) .gt. 0.0)goto 10161                                  
      entrop=0.0                                                                
      goto 10171                                                                
10161 continue                                                                  
      entrop=2.0*w(i)*ni(i)*(y(i)*alog(y(i))+  (1.0-y(i))*alog(1.0-y(i))        
     *)
10171 continue                                                                  
10151 continue                                                                  
      entadd=2.*w(i)*ni(i)*(y(i)*alog(pr)+(1.0-y(i))*alog(1.0-pr))              
      dev=dev+entrop-entadd                                                     
      ent1=ent1+entrop                                                          
      ent2=ent2+entadd                                                          
10101 continue                                                                  
10102 continue                                                                  
      devb=dev                                                                  
      return                                                                    
      end                                                                       
      function devmc(n,fits,y,w,ni)                                             
      real fits(n),y(n),ni(n),w(1)                                              
      dev=0.0                                                                   
      do 10181 i=1,n                                                            
      pr=fits(i)                                                                
      if(pr .ge. 0.0001)goto 10201                                              
      pr=0.0001                                                                 
10201 continue                                                                  
      if(pr .le. 0.9999)goto 10221                                              
      pr=0.9999                                                                 
10221 continue                                                                  
      entrop=2.*w(i)*ni(i)*y(i)*alog(pr)                                        
      dev=dev-entrop                                                            
10181 continue                                                                  
10182 continue                                                                  
      devmc=dev                                                                 
      return                                                                    
      end                                                                       
      subroutine linkg(n,etahat,muhat)                                          
      real muhat(1),etahat(1)                                                   
      do 10231 i=1,n                                                            
      muhat(i)=etahat(i)                                                        
10231 continue                                                                  
10232 continue                                                                  
      return                                                                    
      end                                                                       
      real function linvg(muhat)                                                
      real muhat                                                                
      linvg=muhat                                                               
      return                                                                    
      end                                                                       
      subroutine linkmc(n,etahat,muhat)                                         
      real muhat(1),etahat(1),work(1000)                                        
      real match                                                                
      integer indms,kms(1000)                                                   
      common/mtchco/match,indms,kms                                             
      do 10241 i=1,n                                                            
      muhat(i)=exp(etahat(i))                                                   
10241 continue                                                                  
10242 continue                                                                  
      do 10251 i=1,n                                                            
      work(i)=0.0                                                               
10251 continue                                                                  
10252 continue                                                                  
      do 10261 i=1,n                                                            
      work(kms(i))=work(kms(i))+muhat(i)                                        
10261 continue                                                                  
10262 continue                                                                  
      do 10271 i=1,n                                                            
      muhat(i)=muhat(i)/work(kms(i))                                            
10271 continue                                                                  
10272 continue                                                                  
      return                                                                    
      end                                                                       
      real function linvmc(muhat)                                               
      real muhat                                                                
      linvmc=alog(muhat)                                                        
      return                                                                    
      end                                                                       
      subroutine linkb(n,etahat,muhat)                                          
      real muhat(1),etahat(1)                                                   
      do 10281 i=1,n                                                            
      pr=exp(etahat(i))                                                         
      muhat(i)=pr/(1+pr)                                                        
10281 continue                                                                  
10282 continue                                                                  
      return                                                                    
      end                                                                       
      real function linvb(muhat)                                                
      real muhat                                                                
      real logit                                                                
      logit=muhat                                                               
      d=1.0-logit                                                               
      if(d .ge. 0.0001)goto 10301                                               
      d=0.0001                                                                  
10301 continue                                                                  
      logit=logit/d                                                             
      if(logit .ge. 0.0001)goto 10321                                           
      linvb=alog(0.0001)                                                        
      goto10311                                                                 
10321 if(logit .le. 0 9999.0)goto 10331                                         
      linvb=alog(9999.0)                                                        
      goto 10341                                                                
10331 continue                                                                  
      linvb=alog(logit)                                                         
10341 continue                                                                  
10311 continue                                                                  
      return                                                                    
      end                                                                       
      real function dirivg(muhat)                                               
      real muhat                                                                
      dirivg=1.0                                                                
      return                                                                    
      end                                                                       
      real function dirivb(muhat)                                               
      real muhat                                                                
      pr=muhat                                                                  
      pr=pr*(1.0-pr)                                                            
      if(pr .gt. 0.01)goto 10361                                                
      pr=0.01                                                                   
10361 continue                                                                  
      dirivb=1.0/pr                                                             
      return                                                                    
      end                                                                       
      real function weigtg(w,muhat,ni)                                          
      real muhat,ni                                                             
      weigtg=w                                                                  
      return                                                                    
      end                                                                       
      real function weigtb(w,muhat,ni)                                          
      real muhat,ni                                                             
      pr=dirivb(muhat)                                                          
      weigtb=w*ni/pr                                                            
      return                                                                    
      end                                                                       
      subroutine title                                                          
      write(6,10370)                                                            
10370 format (//    64h     *******         ***        ******     *             
     * **            ,/,    64h    *********       *****       ******   
     *  **     **            ,/,    64h    **     **      **   **       
     * **       ***   ***            ,/,    64h    **            **     
     ***       **       **** ****            ,/,    64h    **           
     * **     **       **       ** *** **            ,/,    64h    **   
     *****     *********       **       **  *  **            ,/,    64h 
     *   **   ****     *********       **       **     **            ,/,
     *    64h    **     **     **     **       **       **     **       
     *     ,/,    64h    *********     **     **     ******     **     *
     **            ,/,    64h     *******      **     **     ******     
     ***     **            ,/,    62h                                   
     *                           ,/,    51h   generalized    additive   
     * interactive    models,//,    53h      coded by     trevor hastie 
     *                    ,/,    61h         and       robert tibshirani
     *  (SLAC      aug 25 1984),/,    61h      last update              
     *       (AT&T labs apr  8 1987),/, //)
      return                                                                    
      end                                                                       
C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine initcx
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
     *),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
      integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
      common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
     *icensq,irindq,iriskq,nq,kq
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      do 10011 i=1,n
      etahat(i)=0
      muhat(i)=0
10011 continue
10012 continue
      i=1
10021 if(icensq(i) .eq. 1) goto 10022
      i=i+1
      goto 10021
10022 continue
      tq(1)=y(i)
      iriskq(1)=i
      jj=1
      i=i+1
10031 if(i .gt. nq) goto 10032
      if(y(i) .eq. y(i-1) .or. icensq(i) .ne. 1)goto 10051
      jj=jj+1
      tq(jj)=y(i)
      iriskq(jj)=i
10051 continue
      i=i+1
      goto 10031
10032 continue
      kq=jj
      call calcdd
      return
      end
      subroutine calcdd
      real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
     *),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
      integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
      common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
     *icensq,irindq,iriskq,nq,kq
      do 10061 i=1,kq
      ddq(i)=0
      ii=iriskq(i)
      if(i .ge. kq)goto 10081
      ij=i+1
      ik=iriskq(ij)-1
      goto 10091
10081 continue
      ik=nq
10091 continue
10071 continue
      do 10101 j=ii,ik
      ddq(i)=ddq(i)+icensq(j)
10101 continue
10102 continue
10061 continue
10062 continue
      return
      end
      real function devc(n,fits,y,w,ni)
      real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
     *),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
      integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
      common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
     *icensq,irindq,iriskq,nq,kq
      real fits(1),y(1),w(1),ni(1)
      call calcsf(fits)
      call calcll(fits,value)
      devc=value
      return
      end
      subroutine calcll(fits,value)
      real fits(1)
      real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
     *),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
      integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
      common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
     *icensq,irindq,iriskq,nq,kq
      j=1
      sum2=0
10111 if(icensq(j) .ne. 0) goto 10112
      j=j+1
      goto 10111
10112 continue
      do 10121 i=j,nq
      temp=fits(i)
      termq(i)=exp(temp)
      sum2=sum2+termq(i)
10121 continue
10122 continue
      sum1=0
      do 10131 i=1,kq
      if(i .eq. 1)goto 10151
      is=i-1
      it=iriskq(is)
      iu=iriskq(i)-1
      do 10161 j=it,iu
      sum2=sum2-termq(j)
10161 continue
10162 continue
10151 continue
      sum1=sum1+sfitsq(i)-ddq(i)*log(sum2)
10131 continue
10132 continue
      value=-2*sum1
      return
      end
      subroutine calcsf(fits)
      real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
     *),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
      integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
      common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
     *icensq,irindq,iriskq,nq,kq
      real fits(1)
      do 10171 i=1,kq
      sfitsq(i)=0
      ii=iriskq(i)
      if(i .ge. kq)goto 10191
      ij=i+1
      ik=iriskq(ij)-1
      goto 10201
10191 continue
      ik=nq
10201 continue
10181 continue
      do 10211 j=ii,ik
      sfitsq(i)=sfitsq(i)+icensq(j)*fits(j)
10211 continue
10212 continue
10171 continue
10172 continue
      return
      end
      subroutine calcz
      real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
     *),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
      integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
      common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
     *icensq,irindq,iriskq,nq,kq
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real z(1000),wz(1000)
      common /wvc/z,wz
      real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
      integer sltag(1000)
      logical vlink,vflag
      common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
      do 10221 i=1,nq
      irindq(i)=0
10221 continue
10222 continue
      do 10231 i=1,kq
      irindq(iriskq(i))=i
10231 continue
10232 continue
      aq(1)=0
      ii=iriskq(1)
      do 10241 i=ii,nq
      aq(1)=aq(1)+exp(etahat(i))
10241 continue
10242 continue
      do 10251 i=2,kq
      ii=i-1
      aq(i)=aq(ii)
      it=i-1
      iv=iriskq(it)
      is=iriskq(i)-1
      do 10261 j=iv,is
      aq(i)=aq(i)-exp(etahat(j))
10261 continue
10262 continue
10251 continue
10252 continue
      bq(nq)=0
      bbq(nq)=0
      do 10271 i=1,kq
      bq(nq)=bq(nq)+ddq(i)*(1.0/aq(i))
      bbq(nq)=bbq(nq)+ddq(i)*(1.0/aq(i))**2
10271 continue
10272 continue
      ii=nq
      do 10281 i=2,nq
      ii=ii-1
      ij=ii+1
      bq(ii)=bq(ij)
      bbq(ii)=bbq(ij)
      if(irindq(ij) .eq. 0)goto 10301
      bq(ii)=bq(ii)-ddq(irindq(ij))*(1.0/aq(irindq(ij)))
      bbq(ii)=bbq(ii)-ddq(irindq(ij))*(1.0/aq(irindq(ij)))**2
10301 continue
10281 continue
10282 continue
      do 10311 i=1,nq
      r1q(i)=icensq(i)-exp(etahat(i))*bq(i)
      r2q(i)=-exp(etahat(i))*bq(i)+exp(2*etahat(i))*bbq(i)
10311 continue
10312 continue
      do 10321 i=1,nq
      z(i)=etahat(i)-r1q(i)/r2q(i)
10321 continue
10322 continue
      do 10331 i=1,nq
      wz(i)=-1.0*r2q(i)
10331 continue
10332 continue
      return
      end
      subroutine stpar
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
     *),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
      integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
      common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
     *icensq,irindq,iriskq,nq,kq
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      s0=0
      do 10341 j=1,p
      slopes(j)=0
10341 continue
10342 continue
      do 10351 i=1,n
      do 10361 j=1,p
      s(i,j)=slopes(j)*x(i,j)
10361 continue
10362 continue
10351 continue
10352 continue
      do 10371 j=1,p
      temp=0
      do 10381 i=1,n
      temp=temp+s(i,j)/n
10381 continue
10382 continue
      do 10391 i=1,n
      s(i,j)=s(i,j)-temp
10391 continue
10392 continue
10371 continue
10372 continue
      do 10401 i=1,n
      muhat(i)=0
      do 10411 j=1,p
      muhat(i)=muhat(i)+s(i,j)
10411 continue
10412 continue
10401 continue
10402 continue
      do 10421 i=1,n
      etahat(i)=muhat(i)
10421 continue
10422 continue
      do 10431 i=1,n
      etaq(i)=0
10431 continue
10432 continue
      s0q=0
      return
      end
C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine multam
      real dumlab(30),modlab(20)
      integer indat(20),indmod(20)
      common/extra/dumlab,modlab,indat,indmod
      real mult,group,ym(1000,15)
      real thetak(15)
      integer mcat,indrm(15)
      common/mltcmn/mult,group,thetak,ym,mcat,indrm
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real z(1000),wz(1000)
      common /wvc/z,wz
      real wm(15,15),wm1(1000)
      real scrat1(1000),scrat2(1000)
      itss=30
      call initm
      devnew=devm(0)
      devnll=devnew
      devold=10*devnew+10
      nit=0
      numsm=0
10011 if(abs((devold-devnew)/devold) .le. thresh .or. nit .ge. maxit) go
     *to 10012
      nit=nit+1
      call getz
      i=1
      goto 10023
10021 i=i+1
10023 if((i).gt.(n))goto 10022
      z(i)=z(i)+etahat(i)
      goto 10021
10022 continue
      call bakfit(n,p,z,x,wz,tag,spans,dofs,slopes,cross,etahat,  s,s0,r
     *ss,thbak,nitb,verbos,quiet,itss)
      do 10031 i=1,n
      etahat(i)=etahat(i)-s0
10031 continue
10032 continue
      j=1
      goto 10043
10041 j=j+1
10043 if((j).gt.(mcat-1))goto 10042
      thetak(j)=thetak(j)+s0
      goto 10041
10042 continue
      s0=0
      numsm=numsm+nitb
      call thetaf
      devold=devnew
      devnew=devm(1)
      if(.not.(.not.quiet))goto 10061
      write(6,10070)nit,devnew
10070 format(/10x,   22houter loop (gam) nit= ,  i3,    5h dev=,f16.5,//
     *)
10061 continue
      goto 10011
10012 continue
      devian=devnew
      nitout=nit
      return
      end
      subroutine initm
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real mult,group,ym(1000,15)
      real thetak(15)
      integer mcat,indrm(15)
      common/mltcmn/mult,group,thetak,ym,mcat,indrm
      real tots(15)
      real linvb
      maxit=20
      do 10081 k=1,mcat
      tots(k)=0
10081 continue
10082 continue
      totn=0.0
      i=1
      goto 10093
10091 i=i+1
10093 if((i).gt.(n))goto 10092
      k=1
      goto 10103
10101 k=k+1
10103 if((k).gt.(mcat))goto 10102
      tots(k)=tots(k)+w(i)*ni(i)*ym(i,k)
      goto 10101
10102 continue
      totn=totn+w(i)*ni(i)
      goto 10091
10092 continue
      tots(1)=tots(1)/totn
      do 10111 k=2,mcat
      tots(k)=tots(k)/totn+tots(k-1)
10111 continue
10112 continue
      if(abs(tots(mcat)-1.0) .le. 00001)goto 10131
      write(6,10140)
10140 format(    18h oops, check initm)
10131 continue
      k=1
      goto 10153
10151 k=k+1
10153 if((k).gt.(mcat-1))goto 10152
      thetak(k)=linvb(tots(k))
      goto 10151
10152 continue
      i=1
      goto 10163
10161 i=i+1
10163 if((i).gt.(n))goto 10162
      etahat(i)=0.0
      do 10171 j=1,p
      s(i,j)=0.0
10171 continue
10172 continue
      goto 10161
10162 continue
      return
      end
      real function devm(inp)
      real mult,group,ym(1000,15)
      real thetak(15)
      integer mcat,indrm(15)
      common/mltcmn/mult,group,thetak,ym,mcat,indrm
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real pm(15),pim(15)
      real pr
      dev=0.0
      ent1=0
      ent2=0
      k=2
      goto 10183
10181 k=k+1
10183 if((k).gt.(mcat-1))goto 10182
      if(thetak(k) .ge. thetak(k-1))goto 10201
      thetak(k)=thetak(k-1)*(1.00001)
      write(6,10210)
10210 format(   27h oops, thetas not ascending)
10201 continue
      goto 10181
10182 continue
      i=1
      goto 10223
10221 i=i+1
10223 if((i).gt.(n))goto 10222
      call linkm(mcat,thetak,etahat(i),pim,pm)
      k=1
      goto 10233
10231 k=k+1
10233 if((k).gt.(mcat))goto 10232
      pr=pm(k)
      yk=ym(i,k)
      if(pr .ge. 0.0001)goto 10251
      pr=0.0001
10251 continue
      if(yk .gt. 0.000001)goto 10271
      entrop=0.0
      goto 10281
10271 continue
      entrop=2.0*w(i)*ni(i)*yk*alog(yk)
10281 continue
10261 continue
      entadd=2.*w(i)*ni(i)*yk*alog(pr)
      dev=dev+entrop-entadd
      goto 10231
10232 continue
      goto 10221
10222 continue
      devm=dev
      return
      end
      subroutine linkm(mcat,thetak,etahat,pim,pm)
      integer mcat
      real thetak(mcat),etahat,pim(mcat),pm(mcat)
      eta=thetak(1)+etahat
      call linkb(1,eta,pm(1))
      pim(1)=pm(1)
      k=2
      goto 10293
10291 k=k+1
10293 if((k).gt.(mcat-1))goto 10292
      eta=thetak(k)+etahat
      call linkb(1,eta,pim(k))
      pm(k)=pim(k)-pim(k-1)
      goto 10291
10292 continue
      pm(mcat)=1.0-pim(mcat-1)
      return
      end
      subroutine getz
      real mult,group,ym(1000,15)
      real thetak(15)
      integer mcat,indrm(15)
      common/mltcmn/mult,group,thetak,ym,mcat,indrm
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real z(1000),wz(1000)
      common /wvc/z,wz
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real wm1(15,15),ystar(15)
      real wm2(15),wystar(15)
      real pim(15),pm(15)
      do 10301 k=1,mcat
      wystar(k)=0.0
      wm2(k)=0.0
10301 continue
10302 continue
      i=1
      goto 10313
10311 i=i+1
10313 if((i).gt.(n))goto 10312
      call linkm(mcat,thetak,etahat(i),pim,pm)
      call getw(mcat,ni(i),wm1,pm,pim)
      ystar(1)=ym(i,1)-pm(1)
      k=2
      goto 10323
10321 k=k+1
10323 if((k).gt.(mcat-1))goto 10322
      ystar(k)=ystar(k-1)+ym(i,k)-pm(k)
      goto 10321
10322 continue
      k=1
      goto 10333
10331 k=k+1
10333 if((k).gt.(mcat-1))goto 10332
      c=pim(k)*(1.0-pim(k))
      ystar(k)=ystar(k)/c
      goto 10331
10332 continue
      k=1
      goto 10343
10341 k=k+1
10343 if((k).gt.(mcat-1))goto 10342
      wm2(k)=0.0
      kk=1
      goto 10353
10351 kk=kk+1
10353 if((kk).gt.(mcat-1))goto 10352
      wm2(k)=wm2(k)+wm1(k,kk)
      goto 10351
10352 continue
      goto 10341
10342 continue
      z(i)=0.0
      wz(i)=0.0
      k=1
      goto 10363
10361 k=k+1
10363 if((k).gt.(mcat-1))goto 10362
      z(i)=z(i)+ystar(k)*wm2(k)
      wz(i)=wz(i)+wm2(k)
      goto 10361
10362 continue
      z(i)=z(i)/wz(i)
      goto 10311
10312 continue
      return
      end
      subroutine thetaf
      real mult,group,ym(1000,15)
      real thetak(15)
      integer mcat,indrm(15)
      common/mltcmn/mult,group,thetak,ym,mcat,indrm
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real wm1(15,15),ystar(15)
      real wm2(15,15),wystar(15)
      real pim(15),pm(15)
      do 10371 k=1,mcat
      wystar(k)=0.0
      do 10381 kk=1,mcat
      wm2(k,kk)=0.0
10381 continue
10382 continue
10371 continue
10372 continue
      i=1
      goto 10393
10391 i=i+1
10393 if((i).gt.(n))goto 10392
      call linkm(mcat,thetak,etahat(i),pim,pm)
      call getw(mcat,ni(i),wm1,pm,pim)
      ystar(1)=ym(i,1)-pm(1)
      k=2
      goto 10403
10401 k=k+1
10403 if((k).gt.(mcat-1))goto 10402
      ystar(k)=ystar(k-1)+ym(i,k)-pm(k)
      goto 10401
10402 continue
      k=1
      goto 10413
10411 k=k+1
10413 if((k).gt.(mcat-1))goto 10412
      c=pim(k)*(1.0-pim(k))
      ystar(k)=ystar(k)/c + thetak(k)
      goto 10411
10412 continue
      k=1
      goto 10423
10421 k=k+1
10423 if((k).gt.(mcat-1))goto 10422
      kk=1
      goto 10433
10431 kk=kk+1
10433 if((kk).gt.(mcat-1))goto 10432
      wystar(k)=wystar(k)+wm1(k,kk)*ystar(kk)
      wm2(k,kk)=wm2(k,kk)+wm1(k,kk)
      goto 10431
10432 continue
      goto 10421
10422 continue
      goto 10391
10392 continue
      call tholve(mcat,wm2,wystar,thetak)
      return
      end
      subroutine tholve(mcat,wm,wystar,thetak)
      integer mcat
      real wm(15,mcat),wystar(mcat),thetak(mcat)
      k=1
      goto 10443
10441 k=k+1
10443 if((k).gt.(mcat-1))goto 10442
      wm(mcat,k)=wystar(k)
      wm(k,mcat)=wystar(k)
      goto 10441
10442 continue
      wm(mcat,mcat)=1.0
      k=1
      goto 10453
10451 k=k+1
10453 if((k).gt.(mcat-1))goto 10452
      call sweep(15,mcat,k,wm)
      goto 10451
10452 continue
      k=1
      goto 10463
10461 k=k+1
10463 if((k).gt.(mcat-1))goto 10462
      thetak(k)=-wm(mcat,k)
      goto 10461
10462 continue
      return
      end
      subroutine sweep(maxp,p,k,a)
      integer p
      real a(maxp,p)
      d=a(k,k)
      do 10471 j=1,p
      a(k,j)=a(k,j)/d
10471 continue
10472 continue
      i=1
      goto 10483
10481 i=i+1
10483 if((i).gt.(p))goto 10482
      if(i .ne. k)goto 10501
      goto 10481
10501 continue
      b=a(i,k)
      j=1
      goto 10513
10511 j=j+1
10513 if((j).gt.(p))goto 10512
      a(i,j)=a(i,j)-b*a(k,j)
      goto 10511
10512 continue
      a(i,k)=-b/d
      goto 10481
10482 continue
      a(k,k)=1/d
      return
      end
      subroutine minfo
      real dumlab(30),modlab(20)
      integer indat(20),indmod(20)
      common/extra/dumlab,modlab,indat,indmod
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real muhat(1000),etahat(1000)
      real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
      logical verbos,quiet
      real s0,devnll
      integer maxit,nitout,numsm
      common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
     *evnll, maxit,verbos,quiet,nitout,numsm
      real mult,group,ym(1000,15)
      real thetak(15)
      integer mcat,indrm(15)
      common/mltcmn/mult,group,thetak,ym,mcat,indrm
      totdof=mcat-1.0
      do 10521 j=1,p
      totdof=totdof+dofs(j)
10521 continue
10522 continue
      totn=(mcat-1)*n
      errn=totn-totdof
      scale=1.0
      write(6,10530)
10530 format(//   35h logistic additive regression model/,    42h propor
     *tional odds model (mc cullagh 1980)/)
      write(6,10540)
10540 format(//   28h estimation by local scoring/)
      write(6,10550)devian,nitout,numsm
10550 format (   12h deviance = ,f14.3,   14h # iterations=,i3,   20h #
     *smooths/variable=,i4)
      a=devian/totn
      write(6,10560)a
10560 format(   20h average deviance = , f12.3)
      write(6,10570)errn,scale
10570 format(   18h dof of deviance  ,f8.2,     18h scale estimate   ,f1
     *0.3)
      r2=100*(devnll-devian)/devnll
      write(6,10580)r2,devnll
10580 format (   12h r square = ,f5.2,   24h% of a null deviance of ,f14
     *.3)
      write(6,10590)
10590 format(/,    65h variable    span   dof     slope                n
     *ame            ,/,    65h --------    ----   ---     ------    ---
     *------------------------)
      write(6,10600)thetak(1),(labels(k,indrm(1)),k=1,10)
10600 format(    26h theta(1)            1    ,f9.4,4x,10a4)
      if(group .ne. yes)goto 10621
      k=2
      goto 10633
10631 k=k+1
10633 if((k).gt.(mcat-1))goto 10632
      write(6,10640)k,thetak(k),(labels(kk,indrm(k)),kk=1,10)
10640 format(     7h theta(,i2,   17h)           1    ,f9.4,4x,10a4)
      goto 10631
10632 continue
      goto 10651
10621 continue
      k=2
      goto 10663
10661 k=k+1
10663 if((k).gt.(mcat-1))goto 10662
      write(6,10670)k,thetak(k)
10670 format(     7h theta(,i2,   17h)           1    ,f9.4)
      goto 10661
10662 continue
10651 continue
10611 continue
      j=1
      goto 10683
10681 j=j+1
10683 if((j).gt.(p))goto 10682
      if(spans(j) .lt. 1.0)goto 10701
      write(6,10710)index(j),spans(j), dofs(j),slopes(j),modlab(j),(labe
     *ls(k,index(j)),k=1,10)
10710 format (    1h ,i3,7x,f6.2,2x, f5.3,1x, f9.4,5x,11a4)
      goto 10721
10701 continue
      write(6,10730)index(j),spans(j),dofs(j),modlab(j),(labels(k,index(
     *j)),k=1,10)
10730 format (    1h ,i3,7x,f6.2,2x, f6.3,    8h  smooth,5x,11a4)
10721 continue
10691 continue
      goto 10681
10682 continue
      write(6,10740)totdof
10740 format (   29h                    -----    ,/,     19h
     *       ,f6.3)
      return
      end
C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine multin
      real dumlab(30),modlab(20)
      integer indat(20),indmod(20)
      common/extra/dumlab,modlab,indat,indmod
      real mult,group,ym(1000,15)
      real thetak(15)
      integer mcat,indrm(15)
      common/mltcmn/mult,group,thetak,ym,mcat,indrm
      real datam(1000,20),labels(10,20)
      integer maxp,n
      common /datac/datam,labels,maxp,n
      real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
      real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
      integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
      logical cross(20)
      common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
     *,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
     *core
      j=1
      goto 10013
10011 j=j+1
10013 if((j).gt.(maxp))goto 10012
      indat(j)=0
      write(6,10020)j,(labels(k,j),k=1,10)
10020 format(i3,    5h ....,10a4)
      goto 10011
10012 continue
10031 continue
      write(6,10040)
10040 format(   31h how many response categories k)
      read(5,*)mcat
      if(mcat .gt. 15)goto 10061
      goto 10032
10061 continue
      write(6,10070)
10070 format(   19h no more than $mcat)
      goto 10031
10032 continue
      write(6,10080)
10080 format(   47h is the input grouped( ie k response variables),/,
     *  39h or ungrouped (ie y=1 or 2 or ... or k),/,     31h yes = grou
     *ped   no = ungrouped)
      read(5,10090)group
10090 format(a1)
      j=1
      goto 10103
10101 j=j+1
10103 if((j).gt.(maxp))goto 10102
      indat(j)=0
      write(6,10110)j,(labels(k,j),k=1,10)
10110 format(i3,    5h ....,10a4)
      goto 10101
10102 continue
      if(group .ne. yes)goto 10131
      write(6,10140)
10140 format(   33h give k response variable indices)
      read(5,*)(indrm(j),j=1,mcat)
      do 10151 j=1,mcat
      indat(indrm(j))=2
10151 continue
10152 continue
      indni=0
      i=1
      goto 10163
10161 i=i+1
10163 if((i).gt.(n))goto 10162
      ni(i)=0
      j=1
      goto 10173
10171 j=j+1
10173 if((j).gt.(mcat))goto 10172
      ym(i,j)=datam(i,indrm(j))
      ni(i)=ni(i)+ym(i,j)
      goto 10171
10172 continue
      j=1
      goto 10183
10181 j=j+1
10183 if((j).gt.(mcat))goto 10182
      ym(i,j)=ym(i,j)/ni(i)
      goto 10181
10182 continue
      w(i)=1.0
      goto 10161
10162 continue
      goto 10191
10131 continue
      write(6,10200)
10200 format(   35h give index of categorical response)
      read(5,*)ir
      indat(ir)=2
      indni=0
      ymin=100
      ymax=0
      i=1
      goto 10213
10211 i=i+1
10213 if((i).gt.(n))goto 10212
      yy=datam(i,ir)
      if(yy .ge. ymin)goto 10231
      ymin=yy
10231 continue
      if(yy .le. ymax)goto 10251
      ymax=yy
10251 continue
      goto 10211
10212 continue
      icat=int(ymax-ymin+1)
      if(icat .eq. mcat)goto 10271
      write(6,10280)mcat,icat
10280 format(   11h there are ,i3,   13h categories, ,      4hnot ,i3)
10271 continue
      if(icat .ge. mcat)goto 10301
      mcat=icat
10301 continue
      imin=int(ymin)
      ito=imin+mcat-1
      write(6,10310)mcat,imin,ito
10310 format(i3,   22h categories are used: ,i3,    4h to ,i3)
      i=1
      goto 10323
10321 i=i+1
10323 if((i).gt.(n))goto 10322
      ni(i)=1
      w(i)=1
      do 10331 j=1,mcat
      ym(i,j)=0.0
10331 continue
10332 continue
      iy=int(datam(i,ir))-imin+1
      if(iy .le. mcat)goto 10351
      iy=mcat
10351 continue
      ym(i,iy)=1
      goto 10321
10322 continue
10191 continue
10121 continue
      p=0
      return
      end
      subroutine getw(mcat,ni,wm1,pm,pim)
      real wm1(15,mcat),pm(mcat),pim(mcat)
      real temp(15)
      real ni
      k=1
      goto 10363
10361 k=k+1
10363 if((k).gt.(mcat-1))goto 10362
      temp(k)=pim(k)*(1.0-pim(k))
      kk=k
      goto 10373
10371 kk=kk+1
10373 if((kk).gt.(mcat-1))goto 10372
      wm1(k,kk)=0.0
      wm1(kk,k)=0.0
      goto 10371
10372 continue
      goto 10361
10362 continue
      k=1
      goto 10383
10381 k=k+1
10383 if((k).gt.(mcat-1))goto 10382
      wm1(k,k)=1.0/pm(k) + 1.0/pm(k+1)
      wm1(k,k)=ni*wm1(k,k)*temp(k)*temp(k)
      wm1(k,k+1)=ni*(-1.0/pm(k+1))*temp(k)*temp(k+1)
      wm1(k+1,k)=wm1(k,k+1)
      goto 10381
10382 continue
      return
      end
      subroutine saxpy(n,sa,sx,incx,sy,incy)
c
c     constant times a vector plus a vector.
c     uses unrolled loop for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      real sx(1),sy(1),sa
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (sa .eq. 0.0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        sy(iy) = sy(iy) + sa*sx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        sy(i) = sy(i) + sa*sx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        sy(i) = sy(i) + sa*sx(i)
        sy(i + 1) = sy(i + 1) + sa*sx(i + 1)
        sy(i + 2) = sy(i + 2) + sa*sx(i + 2)
        sy(i + 3) = sy(i + 3) + sa*sx(i + 3)
   50 continue
      return
      end

      subroutine sbart(xs,ys,ws,n,knot,nk,coef,sz,lev,crit,icrit,spar,
     &ispar,lspar,uspar,tol,isetup,xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,
     &abd,p1ip,p2ip,ld4,ldnk,ier)
      real xs(n),ys(n),ws(n),knot(nk+4),coef(nk),sz(n),lev(n),crit,spar,
     &lspar,uspar,tol,xwy(nk),hs0(nk),hs1(nk),hs2(nk),hs3(nk), sg0(nk),
     &sg1(nk),sg2(nk),sg3(nk),abd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk)
      integer n,nk,isetup,icrit,ispar,ld4,ldnk,ier
      realt1,t2,ratio, a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w, fu,fv,fw,
     &fx,x,ax,bx
      integer i
      i=1
23000 if(.not.(i.le.n))goto 23002
      if(.not.(ws(i).gt.0))goto 23003
      ws(i)=sqrt(ws(i)) 
23003 continue
      i=i+1
      goto 23000
23002 continue
      if(.not.(isetup.eq.0))goto 23005
      call sgram(sg0,sg1,sg2,sg3,knot,nk)
      call stxwx(xs,ys,ws,n,knot,nk,xwy,hs0,hs1,hs2,hs3)
      t1=0. 
      t2=0.
      do 23007 i=3,nk-3 
      t1 = t1 + hs0(i) 
23007 continue
      do 23009 i=3,nk-3 
      t2 = t2 + sg0(i) 
23009 continue
      ratio = t1/t2
      isetup = 1 
23005 continue
      if(.not.(ispar.eq.1))goto 23011
      call sslvrg(xs,ys,ws,n,knot,nk,coef,sz,lev,crit,icrit,spar,ratio,
     &xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,p2ip,ld4,ldnk,ier)
      return
23011 continue
      ax=lspar 
      bx=uspar
      c = 0.5*(3. - sqrt(5.0))
      eps = 1.00
10    eps = eps/2.00
      tol1 = 1.0 + eps
      if(.not.(tol1 .gt. 1.00))goto 23013
      go to 10
23013 continue
      eps = sqrt(eps)
      a = ax
      b = bx
      v = a + c*(b - a)
      w = v
      x = v
      e = 0.0
      spar = x
      call sslvrg(xs,ys,ws,n,knot,nk,coef,sz,lev,crit,icrit,spar,ratio,
     &xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,p2ip,ld4,ldnk,ier)
      fx = crit
      fv = fx
      fw = fx
20    xm = 0.5*(a + b)
      tol1 = eps*abs(x) + tol/3.0
      tol2 = 2.0*tol1
      if(.not.(abs(x - xm) .le. (tol2 - 0.5*(b - a))))goto 23015
      go to 90
23015 continue
      if(.not.(abs(e) .le. tol1))goto 23017
      go to 40
23017 continue
      r = (x - w)*(fx - fv)
      q = (x - v)*(fx - fw)
      p = (x - v)*q - (x - w)*r
      q = 2.00*(q - r)
      if(.not.(q .gt. 0.0))goto 23019
      p = -p
23019 continue
      q = abs(q)
      r = e
      e = d
30    if(.not.(abs(p) .ge. abs(0.5*q*r)))goto 23021
      go to 40
23021 continue
      if(.not.(p .le. q*(a - x)))goto 23023
      go to 40
23023 continue
      if(.not.(p .ge. q*(b - x)))goto 23025
      go to 40
23025 continue
      d = p/q
      u = x + d
      if(.not.((u - a) .lt. tol2))goto 23027
      d = sign(tol1, xm - x)
23027 continue
      if(.not.((b - u) .lt. tol2))goto 23029
      d = sign(tol1, xm - x)
23029 continue
      go to 50
40    if(.not.(x .ge. xm))goto 23031
      e = a - x
23031 continue
      if(.not.(x .lt. xm))goto 23033
      e = b - x
23033 continue
      d = c*e
50    if(.not.(abs(d) .ge. tol1))goto 23035
      u = x + d
23035 continue
      if(.not.(abs(d) .lt. tol1))goto 23037
      u = x + sign(tol1, d)
23037 continue
      spar = u
      call sslvrg(xs,ys,ws,n,knot,nk,coef,sz,lev,crit,icrit,spar,ratio,
     &xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,p2ip,ld4,ldnk,ier)
      fu = crit
      if(.not.(fu .gt. fx))goto 23039
      go to 60
23039 continue
      if(.not.(u .ge. x))goto 23041
      a = x
23041 continue
      if(.not.(u .lt. x))goto 23043
      b = x
23043 continue
      v = w
      fv = fw
      w = x
      fw = fx
      x = u
      fx = fu
      go to 20
60    if(.not.(u .lt. x))goto 23045
      a = u
23045 continue
      if(.not.(u .ge. x))goto 23047
      b = u
23047 continue
      if(.not.(fu .le. fw))goto 23049
      go to 70
23049 continue
      if(.not.(w .eq. x))goto 23051
      go to 70
23051 continue
      if(.not.(fu .le. fv))goto 23053
      go to 80
23053 continue
      if(.not.(v .eq. x))goto 23055
      go to 80
23055 continue
      if(.not.(v .eq. w))goto 23057
      go to 80
23057 continue
      go to 20
70    v = w
      fv = fw
      w = u
      fw = fu
      go to 20
80    v = u
      fv = fu
      go to 20
90    continue 
      spar = x 
      crit = fx
      return
23012 continue
      return
      end
Caveat receptor.  [Jack] dongarra@anl-mcs, [Eric Grosse] research!ehg
Compliments of netlib   Fri May  2 13:34:23 EDT 1986
      real function sdot(n,sx,incx,sy,incy)
c
c     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      real sx(1),sy(1),stemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      stemp = 0.0e0
      sdot = 0.0e0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        stemp = stemp + sx(ix)*sy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      sdot = stemp
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        stemp = stemp + sx(i)*sy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) +
     *   sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4)
   50 continue
   60 sdot = stemp
      return
      end

      subroutine sgram(sg0,sg1,sg2,sg3,tb,nb)
      real sg0(nb),sg1(nb),sg2(nb),sg3(nb),tb(nb+4),vnikx(4,3),work(16),
     &yw1(4),yw2(4),wpt
      integer nb,ileft,ilo,mflag,i,ii,jj
      do 23000 i=1,nb
      sg0(i)=0.
      sg1(i)=0.
      sg2(i)=0.
      sg3(i)=0.
23000 continue
      ilo = 1
      do 23002 i=1,nb 
      call interv(tb(1),(nb+1),tb(i),ileft,mflag)
      call bsplvd (tb,4,tb(i),ileft,work,vnikx,3)
      do 23004 ii=1,4 
      yw1(ii) = vnikx(ii,3) 
23004 continue
      call bsplvd (tb,4,tb(i+1),ileft,work,vnikx,3)
      do 23006 ii=1,4 
      yw2(ii) = vnikx(ii,3) - yw1(ii) 
23006 continue
      wpt = tb(i+1) - tb(i)
      if(.not.(ileft.ge.4))goto 23008
      do 23010 ii=1,4 
      jj=ii
      sg0(ileft-4+ii) = sg0(ileft-4+ii) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
      jj=ii+1
      if(.not.(jj.le.4))goto 23012
      sg1(ileft+ii-4) = sg1(ileft+ii-4) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23012 continue
      jj=ii+2
      if(.not.(jj.le.4))goto 23014
      sg2(ileft+ii-4) = sg2(ileft+ii-4) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23014 continue
      jj=ii+3
      if(.not.(jj.le.4))goto 23016
      sg3(ileft+ii-4) = sg3(ileft+ii-4) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23016 continue
23010 continue
      goto 23009
23008 continue
      if(.not.(ileft.eq.3))goto 23018
      do 23020 ii=1,3 
      jj=ii
      sg0(ileft-3+ii) = sg0(ileft-3+ii) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
      jj=ii+1
      if(.not.(jj.le.3))goto 23022
      sg1(ileft+ii-3) = sg1(ileft+ii-3) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23022 continue
      jj=ii+2
      if(.not.(jj.le.3))goto 23024
      sg2(ileft+ii-3) = sg2(ileft+ii-3) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23024 continue
23020 continue
      goto 23019
23018 continue
      if(.not.(ileft.eq.2))goto 23026
      do 23028 ii=1,2 
      jj=ii
      sg0(ileft-2+ii) = sg0(ileft-2+ii) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
      jj=ii+1
      if(.not.(jj.le.2))goto 23030
      sg1(ileft+ii-2) = sg1(ileft+ii-2) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23030 continue
23028 continue
      goto 23027
23026 continue
      if(.not.(ileft.eq.1))goto 23032
      do 23034 ii=1,1 
      jj=ii
      sg0(ileft-1+ii) = sg0(ileft-1+ii) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23034 continue
23032 continue
23027 continue
23019 continue
23009 continue
23002 continue
      return
      end
      subroutine sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,flag)
      realabd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk),wjm3(3),wjm2(2),wjm1(1)
     &,c0,c1,c2,c3
      integerflag,ld4,nk,ldnk,i,j,k
      wjm3(1)=0. 
      wjm3(2)=0. 
      wjm3(1)=0.
      wjm2(1)=0. 
      wjm2(2)=0.
      wjm1(1)=0.
      do 23000 i=1,nk 
      j=nk-i+1
      c0 = 1./abd(4,j)
      if(.not.(j.le.nk-3))goto 23002
      c1 = abd(1,j+3)*c0
      c2 = abd(2,j+2)*c0
      c3 = abd(3,j+1)*c0 
      goto 23003
23002 continue
      if(.not.(j.eq.nk-2))goto 23004
      c1 = 0.
      c2 = abd(2,j+2)*c0
      c3 = abd(3,j+1)*c0 
      goto 23005
23004 continue
      if(.not.(j.eq.nk-1))goto 23006
      c1 = 0.
      c2 = 0.
      c3 = abd(3,j+1)*c0 
      goto 23007
23006 continue
      if(.not.(j.eq.nk))goto 23008
      c1 = 0.
      c2 = 0.
      c3 = 0.
23008 continue
23007 continue
23005 continue
23003 continue
      p1ip(1,j) = 0. - (c1*wjm3(1)+c2*wjm3(2)+c3*wjm3(3))
      p1ip(2,j) = 0. - (c1*wjm3(2)+c2*wjm2(1)+c3*wjm2(2))
      p1ip(3,j) = 0. - (c1*wjm3(3)+c2*wjm2(2)+c3*wjm1(1))
      p1ip(4,j) = c0**2 +c1**2*wjm3(1)+2.*c1*c2*wjm3(2)+2.*c1*c3*wjm3(3)
     & +c2**2*wjm2(1)+2.*c2*c3*wjm2(2) +c3**2*wjm1(1)
      wjm3(1)=wjm2(1) 
      wjm3(2)=wjm2(2) 
      wjm3(3)=p1ip(2,j)
      wjm2(1)=wjm1(1) 
      wjm2(2)=p1ip(3,j)
      wjm1(1)=p1ip(4,j)
23000 continue
      if(.not.(flag.eq.0))goto 23010
      return
23010 continue
      do 23012 i=1,nk 
      j=nk-i+1
      k=1
23014 if(.not.(k.le.4.and.j+k-1.le.nk))goto 23016
      p2ip(j,j+k-1) = p1ip(5-k,j) 
      k=k+1
      goto 23014
23016 continue
23012 continue
      do 23017 i=1,nk 
      j=nk-i+1
      k=j-4
23019 if(.not.(k.ge.1))goto 23021
      c0 = 1./abd(4,k) 
      c1 = abd(1,k+3)*c0
      c2 = abd(2,k+2)*c0 
      c3 = abd(3,k+1)*c0
      p2ip(k,j) = 0. - ( c1*p2ip(k+3,j) +c2*p2ip(k+2,j) +c3*p2ip(k+1,j) 
     &) 
      k=k-1
      goto 23019
23021 continue
23017 continue
      return
23011 continue
      end
      subroutine sknotl(x,n,knot,k)
      real x(n),knot(n+6),a1,a2,a3,a4
      integer n,k,ndk,j
      a1 = log(50.)/log(2.) 
      a2 = log(100.)/log(2.)
      a3 = log(140.)/log(2.) 
      a4 = log(200.)/log(2.)
      if(.not.(n.lt.50))goto 23000
      ndk = n 
      goto 23001
23000 continue
      if(.not.(n.ge.50 .and. n.lt.200))goto 23002
      ndk = 2.**(a1+(a2-a1)*(n-50.)/150.) 
      goto 23003
23002 continue
      if(.not.(n.ge.200 .and. n.lt.800))goto 23004
      ndk = 2.**(a2+(a3-a2)*(n-200.)/600.) 
      goto 23005
23004 continue
      if(.not.(n.ge.800 .and. n.lt.3200))goto 23006
      ndk = 2.**(a3+(a4-a3)*(n-800.)/2400.) 
      goto 23007
23006 continue
      if(.not.(n.ge.3200))goto 23008
      ndk = 200. + (n-3200)**.2 
23008 continue
23007 continue
23005 continue
23003 continue
23001 continue
      k = ndk + 6
      do 23010 j=1,3 
      knot(j) = x(1) 
23010 continue
      do 23012 j=1,ndk 
      knot(j+3) = x( 1 + (j-1)*(n-1)/(ndk-1) ) 
23012 continue
      do 23014 j=1,3 
      knot(ndk+3+j) = x(n) 
23014 continue
      return
      end
C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine smooth(x,y,w,itag,span,dof,n,nef, cross,smo,s0,rss,quie
     *t,scrat)
      dimension x(n),y(n),w(n),smo(n),scrat(n),itag(n)
      logical cross,quiet
      real cvrss(6),cvspan(6),cvmin
      integer idmin
      data cvspan/0.3,0.4,0.5,0.6,0.7,1.0/
      if(span .ge. 0.0)goto 10021
      call splsm(x,y,w,itag,span,cross,dof,  nef,n,smo,s0,rss,quiet,scra
     *t)
      goto 10031
10021 continue
      penal=0.01
      cvmin=1e15
      idmin=1
      if(.not.(cross))goto 10051
      k=1
      goto 10063
10061 k=k+1
10063 if((k).gt.(6))goto 10062
      call smth(x,y,w,itag,cvspan(k), dof,n,nef,.true.,smo,s0,cvrss(k),s
     *crat)
      if(cvrss(k) .gt. cvmin)goto 10081
      cvmin=cvrss(k)
      idmin=k
10081 continue
      goto 10061
10062 continue
      span=cvspan(idmin)
      if(penal .le. 0.)goto 10101
      cvmin=(1.+penal)*cvmin
      k=6
      goto 10113
10111 k=k+(-1)
10113 if((-1)*((k)-(1)).gt.0)goto 10112
      if(cvrss(k) .gt. cvmin)goto 10131
      goto 10112
10131 continue
      goto 10111
10112 continue
      span=cvspan(k)
10101 continue
      if(.not.(.not.quiet))goto 10151
      do 10161 k=1,6
      write(6,10170)k,cvspan(k),cvrss(k)
10170 format (i4,    7h span= ,f4.1,    8h cvrss= , f10.3)
10161 continue
10162 continue
      write(6,10180)span
10180 format (   20h span chosen .......,f4.1)
10151 continue
10051 continue
      call smth(x,y,w,itag,span,dof,n,nef,.false.,smo,s0,rss,scrat)
10031 continue
10011 continue
      return
      end
      subroutine smth(x,y,w,itag,span,dof,n,nef,cross,smo,s0,rss,scrat)
      dimension x(n),y(n),w(n),smo(n),scrat(n),itag(n)
      real scrat2(1000)
      real*8 sumw,xbar,ybar,cov,var,wtot,wspan,lspan,rspan
      integer left,right
      logical cross,line
      wtot=0
      do 10191 i=1,nef
      wtot=wtot+w(itag(i))
10191 continue
10192 continue
      line=.true.
      if(span .ge. 1.0)goto 10211
      line=.false.
10211 continue
      istart=1
      cov=0.
      var=0.
10221 if(w(itag(istart)).gt.0.0) goto 10222
      istart=istart+1
      goto 10221
10222 continue
      xbar=x(itag(istart))
      ybar=y(itag(istart))
      sumw=w(itag(istart))
      if(.not.(line))goto 10241
      istart=istart+1
      do 10251 i=istart,nef
      xin=x(itag(i))
      yin=y(itag(i))
      win=w(itag(i))
      xbar=(sumw*xbar+xin*win)/(sumw+win)
      ybar=(sumw*ybar+yin*win)/(sumw+win)
      cov=cov+win*(xin-xbar)*(yin-ybar)*(sumw+win)/sumw
      var=var+win*(xin-xbar)**2*(sumw+win)/sumw
      sumw=sumw+win
10251 continue
10252 continue
      i=1
      goto 10263
10261 i=i+1
10263 if((i).gt.(nef))goto 10262
      if(.not.(cross))goto 10281
      xout=x(itag(i))
      yout=y(itag(i))
      win=w(itag(i))
      cov=cov-win*(xout-xbar)*(yout-ybar)*sumw/(sumw-win)
      var=var-win*(xout-xbar)**2*sumw/(sumw-win)
      xbar=(sumw*xbar-win*xout)/(sumw-win)
      ybar=(sumw*ybar-win*yout)/(sumw-win)
      sumw=sumw-win
10281 continue
      if(var .le. 0.)goto 10301
      smo(itag(i))=cov*x(itag(i))/var
      goto 10311
10301 continue
      smo(itag(i))=0
10311 continue
10291 continue
      if(.not.(cross))goto 10331
      xin=x(itag(i))
      yin=y(itag(i))
      win=w(itag(i))
      xbar=(sumw*xbar+xin*win)/(sumw+win)
      ybar=(sumw*ybar+yin*win)/(sumw+win)
      cov=cov+win*(xin-xbar)*(yin-ybar)*(sumw+win)/sumw
      var=var+win*(xin-xbar)**2*(sumw+win)/sumw
      sumw=sumw+win
10331 continue
      goto 10261
10262 continue
      if(var .le. 0)goto 10351
      s0=ybar-cov*xbar/var
      goto 10361
10351 continue
      s0=ybar
10361 continue
10341 continue
      i=nef+1
      goto 10373
10371 i=i+1
10373 if((i).gt.(n))goto 10372
      if(var .le. 0.)goto 10391
      smo(itag(i))=cov*xbar/var
      goto 10401
10391 continue
      smo(itag(i))=0
10401 continue
10381 continue
      goto 10371
10372 continue
      if(var .le. 0)goto 10421
      scrat(1)=cov/var
      goto 10431
10421 continue
      scrat(1)=0.0
10431 continue
10411 continue
      dof=1.0
      goto 10441
10241 continue
      do 10451 i=1,nef
      scrat(i)=y(i)
      scrat2(i)=w(i)
10451 continue
10452 continue
      if(.not.(.not.cross))goto 10471
      i=0
10481 if(i.ge.nef-1) goto 10482
      i=i+1
      m0=i
10491 if(x(itag(i+1)).gt.x(itag(i))) goto 10492
      i=i+1
      if(i .lt. nef)goto 10491
10492 continue
      if(i.eq.m0)goto 10481
      ntie=i-m0+1
      r=0.
      wt=0.
      do 10501 jj=m0,i
      j=itag(jj)
      r=r+y(j)*w(j)
      wt=wt+w(j)
10501 continue
10502 continue
      if(wt .le. 0.0)goto 10521
      r=r/wt
10521 continue
      wt=wt/(i-m0+1)
      do 10531 j=m0,i
      y(itag(j))=r
      w(itag(j))=wt
10531 continue
10532 continue
      goto 10481
10482 continue
10471 continue
      left=istart
      right=istart
      dof=-1.0
      wspan=span*wtot/2
      lspan=-w(itag(istart))
      rspan=w(itag(istart))
      do 10541 i=istart,nef
      lspan=lspan+w(itag(i))
      rspan=rspan-w(itag(i))
10551 if((rspan+w(itag(right+1))) .ge. wspan .or. right .ge. nef) goto 1
     *0552
      right=right+1
      rspan=rspan+w(itag(right))
      xin=x(itag(right))
      yin=y(itag(right))
      win=w(itag(right))
      xbar=(sumw*xbar+xin*win)/(sumw+win)
      ybar=(sumw*ybar+yin*win)/(sumw+win)
      cov=cov+win*(xin-xbar)*(yin-ybar)*(sumw+win)/sumw
      var=var+win*(xin-xbar)**2*(sumw+win)/sumw
      sumw=sumw+win
      goto 10551
10552 continue
10561 if(lspan .le. wspan .or. left .ge. i) goto 10562
      xout=x(itag(left))
      yout=y(itag(left))
      win=w(itag(left))
      cov=cov-win*(xout-xbar)*(yout-ybar)*sumw/(sumw-win)
      var=var-win*(xout-xbar)**2*sumw/(sumw-win)
      xbar=(sumw*xbar-win*xout)/(sumw-win)
      ybar=(sumw*ybar-win*yout)/(sumw-win)
      sumw=sumw-win
      lspan=lspan-w(itag(left))
      left=left+1
      goto 10561
10562 continue
      if(.not.(cross))goto 10581
      xout=x(itag(i))
      yout=y(itag(i))
      win=w(itag(i))
      cov=cov-win*(xout-xbar)*(yout-ybar)*sumw/(sumw-win)
      var=var-win*(xout-xbar)**2*sumw/(sumw-win)
      xbar=(sumw*xbar-win*xout)/(sumw-win)
      ybar=(sumw*ybar-win*yout)/(sumw-win)
      sumw=sumw-win
10581 continue
      if(var .le. 0.)goto 10601
      smo(itag(i))=ybar+cov*(x(itag(i))-xbar)/var
      dof=dof+w(itag(i))/sumw+  (w(itag(i))*(x(itag(i))-xbar)**2)/var
      goto 10611
10601 continue
      smo(itag(i))=ybar
      dof=dof+w(itag(i))/sumw
10611 continue
10591 continue
      if(.not.(cross))goto 10631
      xin=x(itag(i))
      yin=y(itag(i))
      win=w(itag(i))
      xbar=(sumw*xbar+xin*win)/(sumw+win)
      ybar=(sumw*ybar+yin*win)/(sumw+win)
      cov=cov+win*(xin-xbar)*(yin-ybar)*(sumw+win)/sumw
      var=var+win*(xin-xbar)**2*(sumw+win)/sumw
      sumw=sumw+win
10631 continue
10541 continue
10542 continue
      if(.not.(.not.cross))goto 10651
      i=0
10661 if(i.ge.nef-1) goto 10662
      i=i+1
      m0=i
10671 if(x(itag(i+1)).gt.x(itag(i))) goto 10672
      i=i+1
      if(i .lt. nef)goto 10671
10672 continue
      if(i.eq.m0)goto 10661
      ntie=i-m0+1
      r=0.
      wt=0.
      do 10681 jj=m0,i
      j=itag(jj)
      r=r+smo(j)*w(j)
      wt=wt+w(j)
10681 continue
10682 continue
      if(wt .le. 0.0)goto 10701
      r=r/wt
10701 continue
      do 10711 j=m0,i
      smo(itag(j))=r
10711 continue
10712 continue
      goto 10661
10662 continue
10651 continue
      do 10721 i=1,nef
      y(i)=scrat(i)
      w(i)=scrat2(i)
10721 continue
10722 continue
      totsm=0.0
      ybar=0.0
      sumw=0.0
      do 10731 ii=1,nef
      i=itag(ii)
      totsm=totsm+smo(i)*w(i)
      ybar=ybar+w(i)*y(i)
      sumw=sumw+w(i)
10731 continue
10732 continue
      totsm=totsm/sumw
      do 10741 ii=1,nef
      i=itag(ii)
      smo(i)=smo(i)-totsm
10741 continue
10742 continue
      ii=nef+1
      goto 10753
10751 ii=ii+1
10753 if((ii).gt.(n))goto 10752
      i=itag(ii)
      smo(i)=0.0
      ybar=ybar+w(i)*y(i)
      sumw=sumw+w(i)
      goto 10751
10752 continue
      ybar=ybar/sumw
      s0=ybar
10441 continue
10231 continue
      rss=0.0
      do 10761 i=1,n
      rss=rss+(w(i)/sumw)*(y(i)-s0-smo(i))**2
10761 continue
10762 continue
      return
      end
      subroutine spbfa(abd,lda,n,m,info)
      integer lda,n,m,info
      real abd(lda,1)
c
c     spbfa factors a real symmetric positive definite matrix
c     stored in band form.
c
c     spbfa is usually called by spbco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c
c     on entry
c
c        abd     real(lda, n)
c                the matrix to be factored.  the columns of the upper
c                triangle are stored in the columns of abd and the
c                diagonals of the upper triangle are stored in the
c                rows of abd .  see the comments below for details.
c
c        lda     integer
c                the leading dimension of the array  abd .
c                lda must be .ge. m + 1 .
c
c        n       integer
c                the order of the matrix  a .
c
c        m       integer
c                the number of diagonals above the main diagonal.
c                0 .le. m .lt. n .
c
c     on return
c
c        abd     an upper triangular matrix  r , stored in band
c                form, so that  a = trans(r)*r .
c
c        info    integer
c                = 0  for normal return.
c                = k  if the leading minor of order  k  is not
c                     positive definite.
c
c     band storage
c
c           if  a  is a symmetric positive definite band matrix,
c           the following program segment will set up the input.
c
c                   m = (band width above diagonal)
c                   do 20 j = 1, n
c                      i1 = max0(1, j-m)
c                      do 10 i = i1, j
c                         k = i-j+m+1
c                         abd(k,j) = a(i,j)
c                10    continue
c                20 continue
c
c     linpack.  this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas sdot
c     fortran max0,sqrt
c
c     internal variables
c
      real sdot,t
      real s
      integer ik,j,jk,k,mu
c     begin block with ...exits to 40
c
c
         do 30 j = 1, n
            info = j
            s = 0.0e0
            ik = m + 1
            jk = max0(j-m,1)
            mu = max0(m+2-j,1)
            if (m .lt. mu) go to 20
            do 10 k = mu, m
               t = abd(k,j) - sdot(k-mu,abd(ik,jk),1,abd(mu,j),1)
               t = t/abd(m+1,jk)
               abd(k,j) = t
               s = s + t*t
               ik = ik - 1
               jk = jk + 1
   10       continue
   20       continue
            s = abd(m+1,j) - s
c     ......exit
            if (s .le. 0.0e0) go to 40
            abd(m+1,j) = sqrt(s)
   30    continue
         info = 0
   40 continue
      return
      end

      subroutine spbsl(abd,lda,n,m,b)
      integer lda,n,m
      real abd(lda,1),b(1)
c
c     spbsl solves the real symmetric positive definite band
c     system  a*x = b
c     using the factors computed by spbco or spbfa.
c
c     on entry
c
c        abd     real(lda, n)
c                the output from spbco or spbfa.
c
c        lda     integer
c                the leading dimension of the array  abd .
c
c        n       integer
c                the order of the matrix  a .
c
c        m       integer
c                the number of diagonals above the main diagonal.
c
c        b       real(n)
c                the right hand side vector.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains
c        a zero on the diagonal.  technically this indicates
c        singularity but it is usually caused by improper subroutine
c        arguments.  it will not occur if the subroutines are called
c        correctly and  info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call spbco(abd,lda,n,rcond,z,info)
c           if (rcond is too small .or. info .ne. 0) go to ...
c           do 10 j = 1, p
c              call spbsl(abd,lda,n,c(1,j))
c        10 continue
c
c     linpack.  this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas saxpy,sdot
c     fortran min0
c
c     internal variables
c
      real sdot,t
      integer k,kb,la,lb,lm
c
c     solve trans(r)*y = b
c
      do 10 k = 1, n
         lm = min0(k-1,m)
         la = m + 1 - lm
         lb = k - lm
         t = sdot(lm,abd(la,k),1,b(lb),1)
         b(k) = (b(k) - t)/abd(m+1,k)
   10 continue
c
c     solve r*x = y
c
      do 20 kb = 1, n
         k = n + 1 - kb
         lm = min0(k-1,m)
         la = m + 1 - lm
         lb = k - lm
         b(k) = b(k)/abd(m+1,k)
         t = -b(k)
         if(lm .gt. 0)call saxpy(lm,t,abd(la,k),1,b(lb),1)
   20 continue
      return
      end
C     mortran 2.79 (reserved keyword macros of 09/28/81)
      subroutine splsm(x,y,w,itag,span,cross,dof,nef,n,smo,  s0,rss,quie
     *t,scrat)
      real nx(1000),ny(1000),nw(1000),ns(1000)
      integer point(1000)
      logical quiet,cross
      real scrat(1)
      real x(1),y(1),w(1),smo(1),rss,s0
      real range,min,knot(1020),coef(1020),lev(1000)
      real xwy(1000),hs0(1000),hs1(1000),hs2(1000),hs3(1000),  sg0(1000)
     *,sg1(1000),sg2(1000),sg3(1000)
      real abd(4,1000),p1ip(4,1000),p2ip(1)
      integer itag(1)
      i=1
      j=1
      eps=.1e-9
10011 if(i .gt. nef) goto 10012
      sw=0.0
      sy=0.0
      cx=x(itag(i))
10021 if(x(itag(i))-cx .ge. eps .or. i .gt. nef) goto 10022
      point(itag(i))=j
      it=itag(i)
      sw=sw+w(it)
      sy=sy+w(it)*y(it)
      i=i+1
      goto 10021
10022 continue
      if(sw .le. eps)goto 10041
      nx(j)=cx
      nw(j)=sw
      ny(j)=sy/sw
      j=j+1
10041 continue
      goto 10011
10012 continue
      nspl=j-1
      range=nx(nspl)-nx(1)
      min=nx(1)
      i=1
      goto 10053
10051 i=i+1
10053 if((i).gt.(nspl))goto 10052
      nx(i)=(nx(i)- min)/range
      goto 10051
10052 continue
      call sknotl(nx,nspl,knot,k)
      nk=k-4
      spar=-span
      isetup = 0
      ier = 1
      icrit=1
      if(.not.(cross))goto 10071
      ispar=0
      goto 10081
10071 continue
      ispar=1
10081 continue
10061 continue
      call sbart(nx,ny,nw,nspl,knot,nk, coef,ns,lev,crit,icrit,spar,ispa
     *r,0.,1.,.05,isetup, xwy, hs0,hs1,hs2,hs3, sg0,sg1,sg2,sg3, abd,p1i
     *p,p2ip,4,nk,ier)
      span=-spar
      dof=-1.0
      do 10091 i=1,nspl
      dof=dof+lev(i)
10091 continue
10092 continue
      do 10101 ii=1,nef
      i=itag(ii)
      smo(i)=ns(point(i))
10101 continue
10102 continue
      totsm=0.0
      ybar=0.0
      sumw=0.0
      do 10111 ii=1,nef
      i=itag(ii)
      totsm=totsm+smo(i)*w(i)
      ybar=ybar+w(i)*y(i)
      sumw=sumw+w(i)
10111 continue
10112 continue
      totsm=totsm/sumw
      do 10121 ii=1,nef
      i=itag(ii)
      smo(i)=smo(i)-totsm
10121 continue
10122 continue
      ii=nef+1
      goto 10133
10131 ii=ii+1
10133 if((ii).gt.(n))goto 10132
      i=itag(ii)
      smo(i)=0.0
      ybar=ybar+w(i)*y(i)
      sumw=sumw+w(i)
      goto 10131
10132 continue
      ybar=ybar/sumw
      s0=ybar
      rss=0.0
      do 10141 i=1,n
      rss=rss+(w(i)/sumw)*(y(i)-s0-smo(i))**2
10141 continue
10142 continue
      return
      end
      subroutine sslvrg(x,y,w,n,knot,nk,coef,sz,lev,crit,icrit,spar,
     &ratio,xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,p2ip,ld4,ldnk,
     &info)
      real x(n),y(n),w(n),knot(nk+4),coef(nk),sz(n),lev(n),crit,ratio,
     &spar,xwy(nk),hs0(nk),hs1(nk),hs2(nk),hs3(nk), sg0(nk),sg1(nk),sg2(
     &nk),sg3(nk),abd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk),lambda,b0,b1,
     &b2,b3,eps,vnikx(4,1),work(16),xv,bvalue,rss,df
      integer n,nk,icrit,ld4,ldnk,i,icoef,ileft,ilo,info,j,mflag
      ilo = 1 
      eps = .1e-10
      lambda = ratio*16.**(-2. + spar*(6.))
      do 23000 i=1,nk 
      coef(i) = xwy(i) 
23000 continue
      do 23002 i=1,nk 
      abd(4,i) = hs0(i)+lambda*sg0(i) 
23002 continue
      do 23004 i=1,(nk-1) 
      abd(3,i+1) = hs1(i)+lambda*sg1(i) 
23004 continue
      do 23006 i=1,(nk-2) 
      abd(2,i+2) = hs2(i)+lambda*sg2(i) 
23006 continue
      do 23008 i=1,(nk-3) 
      abd(1,i+3) = hs3(i)+lambda*sg3(i) 
23008 continue
      call spbfa(abd,ld4,nk,3,info)
      if(.not.(info.ne.0))goto 23010
      return
23010 continue
      call spbsl(abd,ld4,nk,3,coef)
      icoef = 1
      do 23012 i=1,n 
      xv = x(i)
      sz(i) = bvalue(knot,coef,nk,4,xv,0)
23012 continue
      if(.not.(icrit.eq.0))goto 23014
      return
23014 continue
      call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,0)
      do 23016 i=1,n 
      xv = x(i)
      call interv(knot(1),(nk+1),xv,ileft,mflag)
      if(.not.(mflag.eq.-1))goto 23018
      ileft = 4 
      xv = knot(4)+eps 
23018 continue
      if(.not.(mflag.eq.1))goto 23020
      ileft = nk 
      xv = knot(nk+1)-eps 
23020 continue
      j=ileft-3
      call bsplvd(knot,4,xv,ileft,work,vnikx,1)
      b0=vnikx(1,1)
      b1=vnikx(2,1)
      b2=vnikx(3,1)
      b3=vnikx(4,1)
      lev(i) = (p1ip(4,j)*b0**2 + 2.*p1ip(3,j)*b0*b1 +2.*p1ip(2,j)*b0*
     &b2 + 2.*p1ip(1,j)*b0*b3 +p1ip(4,j+1)*b1**2 + 2.*p1ip(3,j+1)*b1*b2 
     &+2.*p1ip(2,j+1)*b1*b3 +p1ip(4,j+2)*b2**2 + 2.*p1ip(3,j+2)*b2*b3 +
     &p1ip(4,j+3)*b3**2 )*w(i)**2 
23016 continue
      if(.not.(icrit.eq.1))goto 23022
      rss = 0. 
      df = 0.
      do 23024 i=1,n 
      rss = rss + ((y(i)-sz(i))*w(i))**2
23024 continue
      do 23026 i=1,n 
      df = df + 1.-lev(i)
23026 continue
      crit = (rss/n)/((df/n)**2) 
      goto 23023
23022 continue
      crit = 0.
      do 23028 i=1,n 
      crit = crit +(((y(i)-sz(i))*w(i))/(1-lev(i)))**2 
23028 continue
23023 continue
      return
23015 continue
      end
      subroutine stxwx(x,z,w,k,xknot,n,y,hs0,hs1,hs2,hs3)
      real z(k),w(k),x(k),xknot(n+4),y(n),hs0(n),hs1(n),hs2(n),hs3(n),
     &eps,vnikx(4,1),work(16)
      integer k,n,j,i,ilo,ileft,mflag
      do 23000 i=1,n 
      y(i)=0. 
      hs0(i)=0. 
      hs1(i)=0.
      hs2(i)=0. 
      hs3(i)=0. 
23000 continue
      ilo=1 
      eps = .1e-9
      do 23002 i=1,k 
      call interv(xknot(1),(n+1),x(i),ileft,mflag)
      if(.not.(mflag.eq. 1))goto 23004
      if(.not.(x(i).le.(xknot(ileft)+eps)))goto 23006
      ileft=ileft-1
      goto 23007
23006 continue
      return
23007 continue
23004 continue
      call bsplvd (xknot,4,x(i),ileft,work,vnikx,1)
      j= ileft-4+1
      y(j) = y(j)+w(i)**2*z(i)*vnikx(1,1)
      hs0(j)=hs0(j)+w(i)**2*vnikx(1,1)**2
      hs1(j)=hs1(j)+w(i)**2*vnikx(1,1)*vnikx(2,1)
      hs2(j)=hs2(j)+w(i)**2*vnikx(1,1)*vnikx(3,1)
      hs3(j)=hs3(j)+w(i)**2*vnikx(1,1)*vnikx(4,1)
      j= ileft-4+2
      y(j) = y(j)+w(i)**2*z(i)*vnikx(2,1)
      hs0(j)=hs0(j)+w(i)**2*vnikx(2,1)**2
      hs1(j)=hs1(j)+w(i)**2*vnikx(2,1)*vnikx(3,1)
      hs2(j)=hs2(j)+w(i)**2*vnikx(2,1)*vnikx(4,1)
      j= ileft-4+3
      y(j) = y(j)+w(i)**2*z(i)*vnikx(3,1)
      hs0(j)=hs0(j)+w(i)**2*vnikx(3,1)**2
      hs1(j)=hs1(j)+w(i)**2*vnikx(3,1)*vnikx(4,1)
      j= ileft-4+4
      y(j) = y(j)+w(i)**2*z(i)*vnikx(4,1)
      hs0(j)=hs0(j)+w(i)**2*vnikx(4,1)**2 
23002 continue
      return
      end