      block data
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/terprm/maxtry,maxpro,ppconv,r2max
      common/pltprm/nxplt,nyplt,plotrm
      common/pursut/mode
      common/inmode/inmode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      data iprint/2/
      data nprint,maxit/1,2/
      data itap1,itap2,itap3,itap4/3*6,0/
      data maxtry,maxpro,ppconv/5,10,0.15/
      data nxplt,nyplt,plotrm/60,40,0./
      data mode,inmode/2,2/
      data irest/0/
      data rbig,reps,ratio/1.e30,1.e-7,1.e-5/
      data korder,maxkno,banfac/3,8,1.5/
      data sambw,varbw,smofac(1),smofac(2),smofac(3)/0.5,1.0, 0.5,1.0,2.
     *0/
      end
      subroutine plosur(direct)
      dimension direct(1)
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      common/iplopr/iplopr
      iplopr=0
10010 do 10011 i=1,nobs
      s=0.
10020 do 10021 j=1,npred
      s=s+direct(j)*predic((i-1)*npred+j)
10021 continue
10022 continue
      stor4(i)=s
      stor4(2*nobs+i)=smodel(predic((i-1)*npred+1))
10011 continue
10012 continue
      call smooth(stor4(1),stor4(2*nobs+1),nobs,bandw,stor4(1),stor4(nob
     *s+1),  stor4(4*nobs+1),itag)
      ntrim=plotrm*nobs
      imin=ntrim+1
      imax=nobs-ntrim
      xmin=stor4(imin)
      xmax=stor4(imax)
      xmax=xmax+1.e-4*(xmax-xmin)
      rmin=rbig
      rmax=-rmin
10030 do 10031 i=imin,imax
      s=stor4(2*nobs+itag(i))
      if(s .ge. rmin)goto 10051
      rmin=s
10051 continue
      if(s .le. rmax)goto 10071
      rmax=s
10071 continue
10031 continue
10032 continue
      rmax=rmax+1.e-4*(rmax-rmin)
      call psetup(min0(nxplt,100),xmin,xmax,  min0(nyplt,50),rmin,rmax)
10080 do 10081 i=1,nobs
      call plot(stor4(i),stor4(2*nobs+itag(i)))
10081 continue
10082 continue
      call pprint
      return
      end
      subroutine plospl(ipro)
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      common/iplopr/iplopr
      iplopr=ipro
      write(itap3,10090)ipro,(axis((ipro-1)*npred+j),j=1,npred)
10100 do 10101 i=1,nobs
      proj=0.
10110 do 10111 j=1,npred
      proj=proj+axis((ipro-1)*npred+j)*predic((i-1)*npred+j)
10111 continue
10112 continue
      stor4(i)=proj
      itag(i)=i
10101 continue
10102 continue
      call sort(stor4(1),itag,1,nobs)
      ntrim=plotrm*nobs
      imin=ntrim+1
      imax=nobs-ntrim
      xmin=stor4(imin)
      xmax=stor4(imax)
      xmax=xmax+1.e-4*(xmax-xmin)
10120 do 10121 i=imin,imax
      r=bvalue(rkno((ipro-1)*maxkno+1),coeff((ipro-1)*maxkno+1),nkno(ipr
     *o)-korder,  korder,stor4(i),0)
      stor4(nobs+i)=r
10121 continue
10122 continue
      rmin=rbig
      rmax=-rmin
10130 do 10131 i=imin,imax
      if(stor4(nobs+i) .ge. rmin)goto 10151
      rmin=stor4(nobs+i)
10151 continue
      if(stor4(nobs+i) .le. rmax)goto 10171
      rmax=stor4(nobs+i)
10171 continue
10131 continue
10132 continue
      rmin=rmin-0.1*(rmax-rmin)
      rmax=rmax+1.e-4*(rmax-rmin)
      call psetup(min0(nxplt,100),xmin,xmax,min0(nyplt,50),rmin,rmax)
10180 i=1
      goto 10183
10181 i=i+1
10183 if((i).gt.(nkno(ipro)))goto 10182
      call plot(rkno((ipro-1)*maxkno+i),rmin)
      goto 10181
10182 continue
      call pprint
      iplopr=0
      write(itap1,10190)ipro,(axis((ipro-1)*npred+j),j=1,npred)
      rmin=rbig
      rmax=-rmin
10200 do 10201 i=imin,imax
      r=stor4(3*nobs+itag(i))
      if(r .ge. rmin)goto 10221
      rmin=r
10221 continue
      if(r .le. rmax)goto 10241
      rmax=r
10241 continue
      stor4(2*nobs+i)=r
10201 continue
10202 continue
      rmin=rmin-0.1*(rmax-rmin)
      rmax=rmax+1.e-4*(rmax-rmin)
      call smooth(stor4(1),stor4(2*nobs+1),nobs,bandw,stor4(1),stor4(nob
     *s+1),stor4(4*nobs+1),itag)
      call psetup(min0(nxplt,100),xmin,xmax,min0(nyplt,50),rmin,rmax)
10250 do 10251 i=1,nobs
      call plot(stor4(i),stor4(2*nobs+i))
10251 continue
10252 continue
10260 i=1
      goto 10263
10261 i=i+1
10263 if((i).gt.(nkno(ipro)))goto 10262
      call plot(rkno((ipro-1)*maxkno+i),rmin)
      goto 10261
10262 continue
      call pprint
      return
10090 format(   46h1spline function and knots along direction nr ,i3/  2
     *0(7(2x,f7.4)/))
10190 format(   40h1residuals and knots along direction nr ,i3/  20(7(2x
     *,f7.4)/))
      end
      subroutine save
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      if(itap4 .ne. 0)goto 10281
      write(itap1,10290)itap4
      return
10281 continue
      no=nobs
      np=npred
      mk=maxkno
      mp=maxpro
      rewind itap4
      write(itap4) npred,nobs,iprint,nprint,maxit,itap1,itap2,itap3,itap
     *4, ntry,npro,maxtry,maxpro,ppconv,korder,maxkno,banfac,bandw, samb
     *w,varbw,(smofac(i),i=1,3),nxplt,nyplt,plotrm 
      write(itap4) ((predic((j-1)*npred+i),i=1,np),j=1,no),(respon(i),i=
     *1,no),((xlim(i,j), i=1,2),j=1,mp),((axis((j-1)*npred+i),i=1,np),j=
     *1,mp),((rkno((j-1)*maxkno+i),i=1,mk), j=1,mp),((coeff((j-1)*maxkno
     *+i),i=1,mk),j=1,mp),(stor4(3*nobs+i),i=1,no), (nkno(i),i=1,mp),asr
     *old,mode,inmode,irest,rbig,reps,ratio 
      write(itap1,10290)itap4
      return
10290 format(   39h0information for restart saved on tape ,i3)
      end
      subroutine restor
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      if(itap4 .eq. 0)goto 10311
      rewind itap4
      read(itap4) npred,nobs,iprint,nprint,maxit,itap1,itap2,itap3,itap4
     *, ntry,npro,maxtry,maxpro,ppconv,korder,maxkno,banfac,bandw, sambw
     *,varbw,(smofac(i),i=1,3),nxplt,nyplt,plotrm 
      no=nobs
      np=npred
      mk=maxkno
      mp=maxpro
      read(itap4) ((predic((j-1)*npred+i),i=1,np),j=1,no),(respon(i),i=1
     *,no),((xlim(i,j), i=1,2),j=1,mp),((axis((j-1)*npred+i),i=1,np),j=1
     *,mp),((rkno((j-1)*maxkno+i),i=1,mk), j=1,mp),((coeff((j-1)*maxkno+
     *i),i=1,mk),j=1,mp),(stor4(3*nobs+i),i=1,no), (nkno(i),i=1,mp),asro
     *ld,mode,inmode,irest,rbig,reps,ratio 
10311 continue
      irest=1
      write(itap1,10290)itap4
      return
10290 format(   40h0information for restart read from tape ,i3)
      end
      subroutine masa
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      external smuuth
      write(itap1,10090)
      write(itap1,10190)nobs,npred
      write(itap1,10320)mode,maxtry,maxpro,ppconv,maxit
      write(itap1,10330)korder,maxkno,banfac
      write(itap1,10340)iprint,nprint,plotrm
      if(irest .ne. 0)goto 10361
      ntry=1
      npro=1
10370 do 10371 i=1,nobs
      stor4(3*nobs+i)=respon(i)
10371 continue
10372 continue
      inmode=mode
      if(mode .ne. 3)goto 10391
      inmode=1
10391 continue
      sum=0.
10400 do 10401 i=1,nobs
      sum=sum+respon(i)
10401 continue
10402 continue
      sum=sum/float(nobs)
      sstot=0.
10410 do 10411 i=1,nobs
      sstot=sstot+(respon(i)-sum)**2
10411 continue
10412 continue
      asrold=sstot/float(nobs)
      write(itap1,10420)asrold
      goto 10431
10361 continue
      npro=npro+1
      ntry=1
      irest=0
10431 continue
10351 continue
      bandw=banfac/(maxkno-korder)
10440 continue
10441 continue
      write(itap1,10450)npro,ntry
      if(inmode .ne. 1)goto 10471
      call select(asr)
      goto 10481
10471 continue
      call gensta
      call trimin(axis((npro-1)*npred+1),npred,asr,smuuth,maxit,nprint, 
     * stor4(1))
10481 continue
10461 continue
      write(itap1,10490)asr,(axis((npro-1)*npred+i),i=1,npred)
      call range(npro)
      if(iprint .lt. 2)goto 10511
      call plores
10511 continue
10520 do 10521 i=1,npro
      call knopla(i)
10521 continue
10522 continue
      if(iprint .ne. 3)goto 10541
10550 do 10551 i=1,npro
      nki=itnkno(i)
      write(itap1,10560)i,nki,(stor2(maxkno*maxpro+(i-1)*maxkno+j),j=1,n
     *ki)
10551 continue
10552 continue
10541 continue
      call fit(asr)
      write(itap1,10570)asr
      call accept(iacc,asr)
      if(iacc .ne. 0)goto 10591
      write(itap1,10600)
      ntry=ntry+1
      goto 10611
10591 continue
      write(itap1,10620)
10630 do 10631 i=1,nobs
      stor4(3*nobs+i)=stor4(4*nobs+i)
10631 continue
10632 continue
10640 do 10641 i=1,maxkno
10650 do 10651 j=1,npro
      coeff((j-1)*maxkno+i)=stor2(2*maxpro*maxkno+(j-1)*maxkno+i)
      rkno((j-1)*maxkno+i)=stor2(maxkno*maxpro+(j-1)*maxkno+i)
10651 continue
10652 continue
10641 continue
10642 continue
10660 do 10661 i=1,npro
      nkno(i)=itnkno(i)
10661 continue
10662 continue
      npro=npro+1
      ntry=1
      asrold=asr
10611 continue
10581 continue
      call decide(iterm)
      if(iterm .ne. 1)goto 10681
      goto 10442
10681 continue
      goto 10441
10442 continue
      npro=npro-1
      if(iprint .lt. 2)goto 10701
10710 do 10711 i=1,npro
      call plospl(i)
10711 continue
10712 continue
10701 continue
      return
10090 format(   47h1multidimensional additive spline approximation    10
     *h (4/19/80)/   24h0parameters for this run)
10190 format(   13h nobs        ,i5/    13h npred       ,i5)
10320 format(   13h mode        ,i5/    13h maxtry      ,i5/    13h maxp
     *ro      ,i5/    13h ppconv      ,g12.6/    13h maxit       ,i5)
10330 format(    13h korder      ,i5/    13h maxkno      ,i5/    13h ban
     *fac      ,g12.6)
10420 format(   43h0average squared residual around the mean  ,g12.6)
10340 format(   13h iprint      ,i5/    13h nprint      ,i5/    13h plot
     *rm      ,g12.6)
10450 format(   16h1projection no. ,i3,   12h  trial no. ,i3)
10490 format(   48h0minimum asr found by search for new direction: ,g12.
     *6/    20h solution direction:/20(7(2x,f7.4)/))
10560 format(   15h0direction no. ,i3,    9h  nknots=,i3/    16h knot po
     *sitions:/20(7(2x,g12.6)/))
10570 format(   43h0average squared residual from spline fit= g12.6)
10600 format(   24h0projection not accepted)
10620 format(   20h0projection accepted)
      end
      function smodel(x)
      dimension x(1)
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      smod=0.
10720 do 10721 i=1,npro
      proj=0.
10730 do 10731 j=1,npred
      proj=proj+x(j)*axis((i-1)*npred+j)
10731 continue
10732 continue
      if(proj .ge. xlim(1,i))goto 10751
      proj=xlim(1,i)
10751 continue
      if(proj .le. xlim(2,i))goto 10771
      proj=xlim(2,i)
10771 continue
      smod=smod+bvalue(rkno((i-1)*maxkno+1),coeff((i-1)*maxkno+1),  nkno
     *(i)-korder,korder,proj,0)
10721 continue
10722 continue
      smodel=smod
      return
      end
      subroutine gensta
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      if(ntry .ne. 1)goto 10791
      fmin=rbig
10800 do 10801 ipred=1,npred
10810 do 10811 i=1,npred
      axis((npro-1)*npred+i)=0.0
10811 continue
10812 continue
      axis((npro-1)*npred+ipred)=1.0
      asr=smuuth(axis((npro-1)*npred+1),npred)
      if(asr .ge. fmin)goto 10831
      ix=ipred
      fmin=asr
10831 continue
      if(npro .ne. 1)goto 10851
      write(itap1,10090)ipred,asr
      if(iprint .lt. 2)goto 10871
      call range(npro)
      call plores
10871 continue
10851 continue
10801 continue
10802 continue
10880 do 10881 i=1,npred
      axis((npro-1)*npred+i)=0.0
10881 continue
10882 continue
      axis((npro-1)*npred+ix)=1.0
      write(itap1,10190)ix,fmin
      goto 10891
10791 continue
      call genpts(1,npred,axis((npro-1)*npred+1))
      sum=0.
10900 do 10901 i=1,npred
      sum=sum+axis((npro-1)*npred+i)**2
10901 continue
10902 continue
      sum=sqrt(sum)
10910 do 10911 i=1,npred
      axis((npro-1)*npred+i)=axis((npro-1)*npred+i)/sum
10911 continue
10912 continue
      write(itap1,10330)
10891 continue
10781 continue
      return
10090 format(   29h0average squared residual fori3,   14hth predictor= ,
     *g12.6)
10190 format(   14h0predictor nr ,i3,   33h selected for starting direct
     *ion;/      6h asr= ,g12.6)
10330 format(   26h0random starting direction)
      end
      subroutine accept(iacc,asr)
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      iacc=1
      if((asrold-asr) .ge. (ppconv*asrold))goto 10931
      iacc=0
10931 continue
      return
      end
      subroutine decide(iterm)
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      iterm=0
      if(mode .ne. 1)goto 10951
      if((npro .le. maxpro) .and. ((npro .le. npred) .and. (ntry .le. 1)
     *))goto 10971
      iterm=1
10971 continue
      return
10951 continue
      if(mode .ne. 2)goto 10991
      if((npro .le. maxpro) .and. (ntry .le. maxtry))goto 11011
      iterm=1
11011 continue
      return
10991 continue
      if(mode .ne. 3)goto 11031
      if(inmode .ne. 1 .or. (npro .le. npred) .and. (ntry .le. 1))goto 1
     *1051
      inmode=2
      ntry=1
11051 continue
      if((npro .le. maxpro) .and. (ntry .le. maxtry))goto 11071
      iterm=1
11071 continue
11031 continue
      return
      end
      subroutine range(ipro)
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      xmin=rbig
      xmax=-xmin
11080 do 11081 i=1,nobs
      proj=0.
11090 do 11091 j=1,npred
      proj=proj+axis((ipro-1)*npred+j)*predic((i-1)*npred+j)
11091 continue
11092 continue
      if(proj .ge. xmin)goto 11111
      xmin=proj
11111 continue
      if(proj .le. xmax)goto 11131
      xmax=proj
11131 continue
11081 continue
11082 continue
      xlim(1,ipro)=xmin
      xlim(2,ipro)=xmax
      return
      end
      subroutine select(asr)
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      asrmin=rbig
      imin=0
11140 do 11141 i=1,npred
      if(npro .le. 1)goto 11161
      iused=0
11170 j=1
      goto 11173
11171 j=j+1
11173 if((j).gt.(npro-1))goto 11172
      if(axis((j-1)*npred+i) .le. 0.1)goto 11191
      iused=1
      goto 11172
11191 continue
      goto 11171
11172 continue
      if(iused .ne. 1)goto 11211
      goto 11141
11211 continue
11161 continue
11220 do 11221 j=1,npred
      axis((npro-1)*npred+j)=0.
11221 continue
11222 continue
      axis((npro-1)*npred+i)=1.
      asr=smuuth(axis((npro-1)*npred+1),npred)
      if(asr .ge. asrmin)goto 11241
      asrmin=asr
      imin=i
11241 continue
      if(npro .ne. 1)goto 11261
      write(itap1,10090)i,asr
      if(iprint .lt. 2)goto 11281
      call range(npro)
      call plores
11281 continue
11261 continue
11141 continue
11142 continue
11290 do 11291 i=1,npred
      axis((npro-1)*npred+i)=0.
11291 continue
11292 continue
      axis((npro-1)*npred+imin)=1.
      asr=smuuth(axis((npro-1)*npred+1),npred)
10090 format(   29h0average squared residual for,i3,   13hth predictor=,
     *g12.6)
      return
      end
      subroutine fit(asr)
      dimension biatx(20)
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      nbspli=0
11300 do 11301 i=1,npro
      nbspli=nbspli+itnkno(i)
11301 continue
11302 continue
      nbspli=nbspli-(korder+1)*npro+1
11310 do 11311 i=1,nbspli
11320 do 11321 j=1,nbspli
      xtx((j-1)*maxpro*maxkno+i)=0.
11321 continue
11322 continue
      stor4(4*nobs+i)=0.
11311 continue
11312 continue
11330 i=1
      goto 11333
11331 i=i+1
11333 if((i).gt.(nobs))goto 11332
      istart=0
11340 do 11341 j=1,nbspli
      stor2(j)=0.
11341 continue
11342 continue
11350 j=1
      goto 11353
11351 j=j+1
11353 if((j).gt.(npro))goto 11352
      proj=0.
11360 do 11361 k=1,npred
      proj=proj+axis((j-1)*npred+k)*predic((i-1)*npred+k)
11361 continue
11362 continue
      call interv(stor2(maxkno*maxpro+(j-1)*maxkno+1),itnkno(j),proj,lef
     *t,mflag)
      call bsplvb(stor2(maxkno*maxpro+(j-1)*maxkno+1),korder,1,proj,left
     *,biatx)
      if(j .ne. 1)goto 11381
      is=istart+left-korder
11390 do 11391 k=1,korder
      stor2(is+k)=biatx(k)
11391 continue
11392 continue
      istart=istart+itnkno(1)-korder
      goto 11401
11381 continue
      is=istart+left-(korder+1)
      if(left .ne. korder)goto 11421
11430 do 11431 k=2,korder
      stor2(is+k)=biatx(k)
11431 continue
11432 continue
      goto 11441
11421 continue
11450 do 11451 k=1,korder
      stor2(is+k)=biatx(k)
11451 continue
11452 continue
11441 continue
11411 continue
      istart=istart+itnkno(j)-(korder+1)
11401 continue
11371 continue
      goto 11351
11352 continue
11460 do 11461 j=1,nbspli
      stor4(4*nobs+j)=stor4(4*nobs+j)+respon(i)*stor2(j)
11470 do 11471 k=1,nbspli
      xtx((k-1)*maxpro*maxkno+j)=xtx((k-1)*maxpro*maxkno+j)+stor2(j)*sto
     *r2(k)
11471 continue
11472 continue
11461 continue
11462 continue
      goto 11331
11332 continue
      call tred2(maxkno*maxpro,nbspli,xtx,stor4(nobs+1),stor4(1),xtx)
      call tql2(maxkno*maxpro,nbspli,stor4(nobs+1),stor4(1),xtx,ierr)
      if(ierr .eq. 0)goto 11491
      write(itap1,10090)npro,ntry,ierr
      stop
11491 continue
      if(iprint .ne. 3)goto 11511
      write(itap1,10190)(stor4(nobs+i),i=1,nbspli)
10190 format(   19h0eigenvalues of xtx/20(5(1x,g12.6)/))
11511 continue
      small=ratio*stor4(nobs+nbspli)
      ibig=1
11520 continue
11521 if(stor4(nobs+ibig) .ge. small) goto 11522
      ibig=ibig+1
      goto 11521
11522 continue
11530 do 11531 i=ibig,nbspli
      stor4(i)=0.
11540 do 11541 j=1,nbspli
      stor4(i)=stor4(i)+xtx((i-1)*maxpro*maxkno+j)*stor4(4*nobs+j)
11541 continue
11542 continue
      stor4(i)=stor4(i)/stor4(nobs+i)
11531 continue
11532 continue
11550 do 11551 i=1,nbspli
      stor4(4*nobs+i)=0.
11560 do 11561 j=ibig,nbspli
      stor4(4*nobs+i)=stor4(4*nobs+i)+xtx((j-1)*maxpro*maxkno+i)*stor4(j
     *)
11561 continue
11562 continue
11551 continue
11552 continue
      is=0
11570 do 11571 i=1,npro
      it4=itnkno(i)-korder
      if(i .ne. 1)goto 11591
11600 do 11601 j=1,it4
      is=is+1
      stor2(2*maxpro*maxkno+(i-1)*maxkno+j)=stor4(4*nobs+is)
11601 continue
11602 continue
      goto 11611
11591 continue
      stor2(2*maxpro*maxkno+(i-1)*maxkno+1)=0.
11620 do 11621 j=2,it4
      is=is+1
      stor2(2*maxpro*maxkno+(i-1)*maxkno+j)=stor4(4*nobs+is)
11621 continue
11622 continue
11611 continue
11581 continue
11571 continue
11572 continue
      rss=0.
11630 do 11631 i=1,nobs
      stor4(4*nobs+i)=0.
11640 do 11641 j=1,npro
      proj=0.
11650 do 11651 k=1,npred
      proj=proj+axis((j-1)*npred+k)*predic((i-1)*npred+k)
11651 continue
11652 continue
      stor4(4*nobs+i)=stor4(4*nobs+i)+bvalue(stor2(maxkno*maxpro+(j-1)*m
     *axkno+1),stor2(2*maxpro*maxkno+(j-1)*maxkno+1),  itnkno(j)-korder,
     *korder,proj,0)
11641 continue
11642 continue
      stor4(4*nobs+i)=respon(i)-stor4(4*nobs+i)
      rss=rss+stor4(4*nobs+i)**2
11631 continue
11632 continue
      asr=rss/float(nobs-nbspli)
      return
10090 format(   15h1projection nr ,i3,   11h  trial nr ,i3/     33h unab
     *le to compute eigenvalue nr ,i3,    7h of xtx)
      end
      subroutine knopla(ipro)
      real*8 sxy,sx,sy,sxx,syy,a,b,davsr
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
11660 do 11661 i=1,nobs
      proj=0.
      itag(i)=i
11670 do 11671 j=1,npred
      proj=proj+axis((ipro-1)*npred+j)*predic((i-1)*npred+j)
11671 continue
11672 continue
      stor4(i)=proj
11661 continue
11662 continue
      call sort(stor4(1),itag,1,nobs)
      i=0
11680 continue
11681 if(i.ge.nobs-1)goto 11682
      i=i+1
      m0=i
11690 continue
11691 if(stor4(i+1).gt.stor4(i))goto 11692
      i=i+1
      if(i.ge.nobs)goto 11692
      goto 11691
11692 continue
      if(i.eq.m0)goto 11681
      call ranums(stor4(2*nobs+1),i-m0)
11700 j=m0
      goto 11703
11701 j=j+1
11703 if((j).gt.(i-1))goto 11702
      k=j+stor4(2*nobs+j-m0+1)*(i-j+1)
      xt=stor4(j)
      it=itag(j)
      stor4(j)=stor4(k)
      itag(j)=itag(k)
      stor4(k)=xt
      itag(k)=it
      goto 11701
11702 continue
      goto 11681
11682 continue
      avvar=0.
      ibw2=sambw*bandw*nobs/2.0
      sxy=0.0
      sx=sxy
      sy=sx
      sxx=sy
      syy=sxx
11710 do 11711 j=1,ibw2
      x=stor4(j)
      y=stor4(4*nobs+itag(j))
      sxy=sxy+x*y
      sx=sx+x
      sy=sy+y
      sxx=sxx+x*x
      syy=syy+y*y
11711 continue
11712 continue
      ibweff=ibw2
11720 do 11721 i=1,nobs
      if(i .gt. nobs-ibw2)goto 11741
      x=stor4(i+ibw2)
      y=stor4(4*nobs+itag(i+ibw2))
      ibweff=ibweff+1
      sxy=sxy+x*y
      sx=sx+x
      sy=sy+y
      sxx=sxx+x*x
      syy=syy+y*y
11741 continue
      if(i .le. ibw2+1)goto 11761
      x=stor4(i-ibw2-1)
      y=stor4(4*nobs+itag(i-ibw2-1))
      ibweff=ibweff-1
      sxy=sxy-x*y
      sx=sx-x
      sy=sy-y
      sxx=sxx-x*x
      syy=syy-y*y
11761 continue
      denom=sxx-sx*sx/ibweff
      a=0.0
      if(denom.gt.0.0) a=(sxy-sx*sy/ibweff)/denom
      b=(sy-a*sx)/ibweff
      var=(syy+a*a*sxx+2.0*(a*(b*sx-sxy)-b*sy))/ibweff +b*b
      stor4(2*nobs+i)=var
11721 continue
11722 continue
      ibw2=varbw*bandw*nobs/2.0
      sx=0.0
11770 do 11771 j=1,ibw2
      sx=sx+stor4(2*nobs+j)
11771 continue
11772 continue
      ibweff=ibw2
11780 do 11781 i=1,nobs
      if(i .gt. nobs-ibw2)goto 11801
      sx=sx+stor4(2*nobs+i+ibw2)
      ibweff=ibweff+1
11801 continue
      if(i .le. ibw2+1)goto 11821
      sx=sx-stor4(2*nobs+i-ibw2-1)
      ibweff=ibweff-1
11821 continue
      stor4(nobs+i)=sx/ibweff
11781 continue
11782 continue
      avvar=0.
11830 do 11831 i=1,nobs
      avvar=avvar+stor4(nobs+i)
11831 continue
11832 continue
      avvar=avvar/float(nobs)
      rhigh=avvar*smofac(3)
      rlow=avvar*smofac(1)
      avvar=0.
11840 do 11841 i=1,nobs
      if(stor4(nobs+i) .le. rhigh)goto 11861
      stor4(nobs+i)=rhigh
      goto 11851
11861 if(stor4(nobs+i) .ge. rlow)goto 11871
      stor4(nobs+i)=rlow
11871 continue
11851 continue
      stor4(nobs+i)=1./stor4(nobs+i)
      avvar=avvar+stor4(nobs+i)
11841 continue
11842 continue
11880 do 11881 i=1,nobs
      stor4(nobs+i)=stor4(nobs+i)/avvar
11881 continue
11882 continue
      nkint=maxkno-2*korder
11890 do 11891 i=1,korder
      stor2(maxkno*maxpro+(ipro-1)*maxkno+i)=stor4(1)
11891 continue
11892 continue
      nk=korder
      nsame=korder
      sum=0.
      ind=0
11900 do 11901 i=1,nkint
11910 continue
11911 continue
      ind=ind+1
      sum=sum+stor4(nobs+ind)
      if(sum.ge.float(i)/float(nkint+1))goto 11912
      goto 11911
11912 continue
      if(stor4(ind) .ne. stor2(maxkno*maxpro+(ipro-1)*maxkno+nk))goto 11
     *931
      nsame=nsame+1
      goto 11941
11931 continue
      nsame=1
11941 continue
11921 continue
      if(nsame .le. korder)goto 11961
      goto 11901
11961 continue
      nk=nk+1
      stor2(maxkno*maxpro+(ipro-1)*maxkno+nk)=stor4(ind)
11901 continue
11902 continue
11970 i=nk
      goto 11973
11971 i=i+(-1)
11973 if((-1)*((i)-(1)).gt.0)goto 11972
      if(stor2(maxkno*maxpro+(ipro-1)*maxkno+i) .eq. stor4(nobs))goto 11
     *991
      goto 11972
11991 continue
      nk=nk-1
      goto 11971
11972 continue
      if(nk .ge. korder)goto 12011
      write(itap1,10090)(axis((ipro-1)*npred+i),i=1,npred)
      stop
12011 continue
      d=(stor4(nobs)-stor4(1))*10.e-6
12020 i=nk+1
      goto 12023
12021 i=i+1
12023 if((i).gt.(nk+korder))goto 12022
      stor2(maxkno*maxpro+(ipro-1)*maxkno+i)=stor4(nobs)+d
      goto 12021
12022 continue
      itnkno(ipro)=nk+korder
      return
10090 format(   43h0no knot placement possible, only one value /    6h a
     *xis:/20(7(2x,f7.4)/))
      end
      subroutine trimin(x,nvar,fmin,fun,maxit,iprint,scrat)
      dimension x(nvar),scrat(nvar,3)
      common/ntapes/itap1,itap2,itap3,itap4
      external fun
      if(iprint .le. 0)goto 12041
      write(itap2,10090)(x(i),i=1,nvar)
12041 continue
      nit=1
      nfun=0
12050 continue
12051 continue
      xnorm=0.
12060 i=1
      goto 12063
12061 i=i+1
12063 if((i).gt.(nvar))goto 12062
      xnorm=xnorm+x(i)**2
      goto 12061
12062 continue
      xnorm=sqrt(xnorm)
      scrat(1,1)=sqrt((1.-x(1)/xnorm)/2.)
      if(scrat(1,1) .ne. 0.)goto 12081
12090 i=2
      goto 12093
12091 i=i+1
12093 if((i).gt.(nvar))goto 12092
      scrat(i,1)=0.
      goto 12091
12092 continue
      goto 12101
12081 continue
12110 i=2
      goto 12113
12111 i=i+1
12113 if((i).gt.(nvar))goto 12112
      scrat(i,1)=-x(i)/(2.*xnorm*scrat(1,1))
      goto 12111
12112 continue
12101 continue
12071 continue
12120 i=2
      goto 12123
12121 i=i+1
12123 if((i).gt.(nvar))goto 12122
12130 j=1
      goto 12133
12131 j=j+1
12133 if((j).gt.(nvar))goto 12132
      scrat(j,2)=-2.*scrat(i,1)*scrat(j,1)
      goto 12131
12132 continue
      scrat(i,2)=scrat(i,2)+1
      call cirmin(x,scrat(1,2),nvar,fun,nfun,scrat(1,3))
      goto 12121
12122 continue
      if(mod(nit,iprint) .ne. 0 .or. iprint .le. 0)goto 12151
      fmin=fun(x,nvar)
      write(itap2,10190)nit,nfun,fmin,(x(i),i=1,nvar)
12151 continue
      nit=nit+1
      if(nit .le. maxit)goto 12171
      goto 12052
12171 continue
      goto 12051
12052 continue
      nit=nit-1
      fmin=fun(x,nvar)
      nfun=nfun+1
      if(mod(nit,iprint) .eq. 0 .or. iprint .le. 0)goto 12191
      write(itap2,10190)nit,nfun,fmin,(x(i),i=1,nvar)
12191 continue
      if(iprint .le. 0)goto 12211
      write(itap2,10320)
12211 continue
      return
10090 format(   14h0starting axis/10(5(2x,g12.6)/))
10190 format(   11h0iteration ,i5,    6h nfun=,i5,    7h fmin= ,g12.6,  
     *  15h current axis =/10(5(2x,g12.6)/))
10320 format(   21h iteration terminated)
      end
      subroutine cirmin(x,rotvec,nvar,fun,nfun,scrat)
      dimension x(nvar),rotvec(nvar),scrat(nvar),delta(6)
      data delta/0.05,0.1,0.2,0.2,0.2,0.3/
      data pi/3.14159/
      idelta=0
      u=0.
      v=0.
      w=0.
      sign=1.
      fu=fun(x,nvar)
      f0=fu
      nfun=nfun+1
      idelta=idelta+1
      v=u+delta(idelta)
      acos=cos(pi*v)
      asin=sin(pi*v)
12220 i=1
      goto 12223
12221 i=i+1
12223 if((i).gt.(nvar))goto 12222
      scrat(i)=acos*x(i)+asin*rotvec(i)
      goto 12221
12222 continue
      fv=fun(scrat,nvar)
      nfun=nfun+1
      if(fv .lt. fu)goto 12241
      sign=-1
      idelta=0
      uu=u
      u=-v
      v=uu
      fuu=fu
      fu=fv
      fv=fuu
12241 continue
12250 continue
12251 continue
      idelta=idelta+1
      w=v+delta(idelta)
      if(w .ge. 1.)goto 12271
      sw=sign*w
      acos=cos(sw*pi)
      asin=sin(sw*pi)
12280 i=1
      goto 12283
12281 i=i+1
12283 if((i).gt.(nvar))goto 12282
      scrat(i)=acos*x(i)+asin*rotvec(i)
      goto 12281
12282 continue
      fw=fun(scrat,nvar)
      nfun=nfun+1
      if(fw .le. fv)goto 12301
      goto 12310
12301 continue
      u=v
      fu=fv
      v=w
      fv=fw
      goto 12321
12271 continue
      w=1.
      fw=f0
      goto 12310
12321 continue
12261 continue
      goto 12251
12252 continue
12310 continue
      a2=(fw-fu)/((w-u)*(w-v))-(fv-fu)/((v-u)*(w-v))
      a1=(fv-fu)/(v-u)-a2*(v+u)
      xmin=-a1/(2*a2)*sign
      acos=cos(pi*xmin)
      asin=sin(pi*xmin)
12330 i=1
      goto 12333
12331 i=i+1
12333 if((i).gt.(nvar))goto 12332
      scrat(i)=acos*x(i)+asin*rotvec(i)
      goto 12331
12332 continue
      fsol=fun(scrat,nvar)
      nfun=nfun+1
      if(fsol .gt. fv)goto 12351
12360 i=1
      goto 12363
12361 i=i+1
12363 if((i).gt.(nvar))goto 12362
      x(i)=scrat(i)
      goto 12361
12362 continue
      return
12351 continue
      acos=cos(sign*pi*v)
      asin=sin(sign*pi*v)
12370 i=1
      goto 12373
12371 i=i+1
12373 if((i).gt.(nvar))goto 12372
      x(i)=acos*x(i)+asin*rotvec(i)
      goto 12371
12372 continue
      return
      end
      subroutine plores
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      common/iplopr/iplopr
      iplopr=0
      write(itap3,10190)(axis((npro-1)*npred+j),j=1,npred)
      ntrim=plotrm*nobs
      imin=ntrim+1
      imax=nobs-ntrim
12380 do 12381 i=1,nobs
      s=0.
12390 do 12391 j=1,npred
      s=s+axis((npro-1)*npred+j)*predic((i-1)*npred+j)
12391 continue
12392 continue
      stor4(i)=s
      itag(i)=i
12381 continue
12382 continue
      call sort(stor4(1),itag,1,nobs)
      xmin=stor4(imin)
      xmax=stor4(imax)
      xmax=xmax+1.e-4*(xmax-xmin)
      rmin=rbig
      rmax=-rmin
12400 do 12401 i=imin,imax
      r=stor4(3*nobs+itag(i))
      if(r .ge. rmin)goto 12421
      rmin=r
12421 continue
      if(r .le. rmax)goto 12441
      rmax=r
12441 continue
      stor4(nobs+i)=r-stor4(4*nobs+itag(i))
12401 continue
12402 continue
      rmax=rmax+1.e-4*(rmax-rmin)
      call psetup(min0(nxplt,100),xmin,xmax,  min0(nyplt,50),rmin,rmax)
12450 do 12451 i=1,nobs
      call plot(stor4(i),stor4(3*nobs+itag(i)))
12451 continue
12452 continue
      call pprint
      write(itap3,10090)
10090 format(1h0/1h0/1h0)
10190 format(   47h1residuals from model of previous iteration and/     
     *38h moving average smooth along direction/  20(7(2x,f7.4)/))
      return
      end
      function ufun1(x,y)
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      common/iplopr/iplopr
      if(iplopr .ne. 0)goto 12471
      if(x .gt. stor4(1))goto 12491
      ufun1=y-stor4(nobs+1)
      return
12491 continue
      if(x .lt. stor4(nobs))goto 12511
      ufun1=y-stor4(nobs+nobs)
      return
12511 continue
      min=1
      max=nobs
12520 continue
12521 if(max-min.le.1)goto 12522
      nt=(max+min)/2
      if(x .ge. stor4(nt))goto 12541
      max=nt
      goto 12551
12541 continue
      min=nt
12551 continue
12531 continue
      goto 12521
12522 continue
      if(stor4(max) .ne. stor4(min))goto 12571
      ufun1=y-stor4(nobs+min)
      return
12571 continue
      ufun1=y-stor4(nobs+min)-(stor4(nobs+max)-stor4(nobs+min))  *(x-sto
     *r4(min))/(stor4(max)-stor4(min))
      goto 12581
12471 continue
      xx=x
      if(xx .ge. xlim(1,iplopr))goto 12601
      xx=xlim(1,iplopr)
12601 continue
      if(xx .le. xlim(2,iplopr))goto 12621
      xx=xlim(2,iplopr)
12621 continue
      r=bvalue(rkno((iplopr-1)*maxkno+1),coeff((iplopr-1)*maxkno+1),  nk
     *no(iplopr)-korder,korder,xx,0)
      ufun1=r-y
12581 continue
12461 continue
      return
      end
      function ufun2(x,y)
      ufun2=1.
      return
      end
      function smuuth(ax,np)
      dimension ax(np)
      dimension ibw(3),ibwt(3)
      real*8 sxy(3),sx(3),sy(3),sxx(3),syy,a,b,davsr
      common/inprm/npred,nobs
      common/iprint/iprint
      common/minpar/nprint,maxit
      common/ntapes/itap1,itap2,itap3,itap4
      common/ntry/ntry
      common/npro/npro
      common/terprm/maxtry,maxpro,ppconv
      common/smothr/korder,maxkno,banfac
      common/smopar/bandw,sambw,varbw,smofac(3)
      common/pltprm/nxplt,nyplt,plotrm
      common/predic/predic(1)
      common/respon/respon(1)
      common/stor1/xlim(2,1)
      common/axis/axis(1)
      common/stor2/stor2(1)
      common/stor3/xtx(1)
      common/stor4/stor4(1)
      common/stor5/itnkno(1)
      common/stor6/itag(1)
      common/asrold/asrold
      common/pursut/mode
      common/irest/irest
      common/consts/rbig,reps,ratio
      common/coeff/coeff(1)
      common/rkno/rkno(1)
      common/nkno/nkno(1)
      common/inmode/inmode
      ic=0
12630 do 12631 i=1,npred
      if(ax(i).ne.0.0) ic=ic+1
12631 continue
12632 continue
      if(ic .ne. 1)goto 12651
12660 id=1
      goto 12663
12661 id=id+1
12663 if((id).gt.(npred))goto 12662
      if(ax(id).ne.0.0)goto 12662
      goto 12661
12662 continue
12651 continue
12670 i=1
      goto 12673
12671 i=i+1
12673 if((i).gt.(nobs))goto 12672
      if(ic .ne. 1)goto 12691
      stor4(nobs+i)=predic((i-1)*npred+id)
      goto 12671
12691 continue
      sum=0.
12700 j=1
      goto 12703
12701 j=j+1
12703 if((j).gt.(npred))goto 12702
      sum=sum+ax(j)*predic((i-1)*npred+j)
      goto 12701
12702 continue
      stor4(nobs+i)=sum
      goto 12671
12672 continue
12710 i=1
      goto 12713
12711 i=i+1
12713 if((i).gt.(nobs))goto 12712
      itag(i)=i
      goto 12711
12712 continue
      call sort(stor4(nobs+1),itag,1,nobs)
      i=0
12720 continue
12721 if(i.ge.nobs-1)goto 12722
      i=i+1
      m0=i
12730 continue
12731 if(stor4(nobs+i+1).gt.stor4(nobs+i))goto 12732
      i=i+1
      if(i.ge.nobs)goto 12732
      goto 12731
12732 continue
      if(i.eq.m0)goto 12721
      call ranums(stor4(2*nobs+1),i-m0)
12740 j=m0
      goto 12743
12741 j=j+1
12743 if((j).gt.(i-1))goto 12742
      k=j+stor4(2*nobs+j-m0+1)*(i-j+1)
      xt=stor4(nobs+j)
      it=itag(j)
      stor4(nobs+j)=stor4(nobs+k)
      itag(j)=itag(k)
      stor4(nobs+k)=xt
      itag(k)=it
      goto 12741
12742 continue
      goto 12721
12722 continue
      avvar=0.
      ibw2=sambw*bandw*nobs/2.0
      sxy(1)=0.0
      sx(1)=sxy(1)
      sy(1)=sx(1)
      sxx(1)=sy(1)
      syy=sxx(1)
12750 do 12751 j=1,ibw2
      x=stor4(nobs+j)
      y=stor4(3*nobs+itag(j))
      sxy(1)=sxy(1)+x*y
      sx(1)=sx(1)+x
      sy(1)=sy(1)+y
      sxx(1)=sxx(1)+x*x
      syy=syy+y*y
12751 continue
12752 continue
      ibweff=ibw2
12760 do 12761 i=1,nobs
      if(i .gt. nobs-ibw2)goto 12781
      x=stor4(nobs+i+ibw2)
      y=stor4(3*nobs+itag(i+ibw2))
      ibweff=ibweff+1
      sxy(1)=sxy(1)+x*y
      sx(1)=sx(1)+x
      sy(1)=sy(1)+y
      sxx(1)=sxx(1)+x*x
      syy=syy+y*y
12781 continue
      if(i .le. ibw2+1)goto 12801
      x=stor4(nobs+i-ibw2-1)
      y=stor4(3*nobs+itag(i-ibw2-1))
      ibweff=ibweff-1
      sxy(1)=sxy(1)-x*y
      sx(1)=sx(1)-x
      sy(1)=sy(1)-y
      sxx(1)=sxx(1)-x*x
      syy=syy-y*y
12801 continue
      denom=sxx(1)-sx(1)*sx(1)/ibweff
      a=0.0
      if(denom.gt.0.0) a=(sxy(1)-sx(1)*sy(1)/ibweff)/denom
      b=(sy(1)-a*sx(1))/ibweff
      var=(syy+a*a*sxx(1)+2.0*(a*(b*sx(1)-sxy(1))-b*sy(1)))/ibweff +b*b
      avvar=avvar+var
      stor4(4*nobs+i)=var
12761 continue
12762 continue
      avvar=avvar/float(nobs)
      ibw2=varbw*bandw*nobs/2.0
      sx(1)=0.0
12810 do 12811 j=1,ibw2
      sx(1)=sx(1)+stor4(4*nobs+j)
12811 continue
12812 continue
      ibweff=ibw2
12820 do 12821 i=1,nobs
      if(i .gt. nobs-ibw2)goto 12841
      sx(1)=sx(1)+stor4(4*nobs+i+ibw2)
      ibweff=ibweff+1
12841 continue
      if(i .le. ibw2+1)goto 12861
      sx(1)=sx(1)-stor4(4*nobs+i-ibw2-1)
      ibweff=ibweff-1
12861 continue
      stor4(2*nobs+i)=sx(1)/ibweff
12821 continue
12822 continue
12870 do 12871 i=1,3
      sxy(i)=0.0
      sx(i)=sxy(i)
      sy(i)=sx(i)
      sxx(i)=sy(i)
      ibw(i)=smofac(i)*bandw*nobs/2.0+0.5
12871 continue
12872 continue
12880 j=2
      goto 12883
12881 j=j+1
12883 if((j).gt.(ibw(3)))goto 12882
      x=stor4(nobs+j)
      y=stor4(3*nobs+itag(j))
      sxy(3)=sxy(3)+x*y
      sx(3)=sx(3)+x
      sy(3)=sy(3)+y
      sxx(3)=sxx(3)+x*x
      if(j.gt.ibw(2))goto 12881
      sx(2)=sx(2)+x
      sy(2)=sy(2)+y
      sxx(2)=sxx(2)+x*x
      sxy(2)=sxy(2)+x*y
      if(j.gt.ibw(1))goto 12881
      sx(1)=sx(1)+x
      sy(1)=sy(1)+y
      sxx(1)=sxx(1)+x*x
      sxy(1)=sxy(1)+x*y
      goto 12881
12882 continue
12890 do 12891 i=1,3
      ibwt(i)=ibw(i)-1
12891 continue
12892 continue
12900 do 12901 i=1,nobs
12910 do 12911 j=1,3
      if(i .le. 1)goto 12931
      x=stor4(nobs+i-1)
      y=stor4(3*nobs+itag(i-1))
      sxy(j)=sxy(j)+x*y
      sx(j)=sx(j)+x
      sy(j)=sy(j)+y
      sxx(j)=sxx(j)+x*x
      x=stor4(nobs+i)
      y=stor4(3*nobs+itag(i))
      sxy(j)=sxy(j)-x*y
      sx(j)=sx(j)-x
      sy(j)=sy(j)-y
      sxx(j)=sxx(j)-x*x
12931 continue
      if(i .gt. nobs-ibw(j))goto 12951
      x=stor4(nobs+i+ibw(j))
      y=stor4(3*nobs+itag(i+ibw(j)))
      ibwt(j)=ibwt(j)+1
      sxy(j)=sxy(j)+x*y
      sx(j)=sx(j)+x
      sy(j)=sy(j)+y
      sxx(j)=sxx(j)+x*x
12951 continue
      if(i .le. ibw(j)+1)goto 12971
      x=stor4(nobs+i-ibw(j)-1)
      y=stor4(3*nobs+itag(i-ibw(j)-1))
      ibwt(j)=ibwt(j)-1
      sxy(j)=sxy(j)-x*y
      sx(j)=sx(j)-x
      sy(j)=sy(j)-y
      sxx(j)=sxx(j)-x*x
12971 continue
12911 continue
12912 continue
      thisfc=stor4(2*nobs+i)/avvar
      if(thisfc .gt. smofac(2))goto 12991
      denom=sxx(1)-sx(1)*sx(1)/ibwt(1)
      a=0.0
      if(denom.gt.0.0) a=(sxy(1)-sx(1)*sy(1)/ibwt(1))/denom
      s1=a*stor4(nobs+i)+(sy(1)-a*sx(1))/ibwt(1)
      goto 13001
12991 continue
      denom=sxx(3)-sx(3)*sx(3)/ibwt(3)
      a=0.0
      if(denom.gt.0.0) a=(sxy(3)-sx(3)*sy(3)/ibwt(3))/denom
      s3=a*stor4(nobs+i)+(sy(3)-a*sx(3))/ibwt(3)
13001 continue
12981 continue
      if(thisfc .gt. smofac(1))goto 13021
      stor4(4*nobs+i)=s1
      goto 12901
13021 continue
      if(thisfc .lt. smofac(3))goto 13041
      stor4(4*nobs+i)=s3
      goto 12901
13041 continue
      denom=sxx(2)-sx(2)*sx(2)/ibwt(2)
      a=0.0
      if(denom.gt.0.0) a=(sxy(2)-sx(2)*sy(2)/ibwt(2))/denom
      s2=a*stor4(nobs+i)+(sy(2)-a*sx(2))/ibwt(2)
      if(thisfc .gt. smofac(2))goto 13061
      wt=(smofac(2)-thisfc)/(smofac(2)-smofac(1))
      stor4(4*nobs+i)=(1.0-wt)*s2+wt*s1
      goto 13071
13061 continue
      wt=(thisfc-smofac(2))/(smofac(3)-smofac(2))
      stor4(4*nobs+i)=(1.0-wt)*s2+wt*s3
13071 continue
13051 continue
12901 continue
12902 continue
13080 do 13081 i=1,nobs
      stor4(nobs+itag(i))=stor4(4*nobs+i)
13081 continue
13082 continue
      avsr=0.
13090 do 13091 i=1,nobs
      r=stor4(3*nobs+i)-stor4(nobs+i)
      avsr=avsr+r**2
      stor4(4*nobs+i)=r
13091 continue
13092 continue
      avsr=avsr/float(nobs)
      smuuth=avsr
      return
      end
      SUBROUTINE PSETUP(NNXIN,XLIN,XHIN,NNYIN,YLIN,YHIN)
C       INITIALIZATION FOR PLOTTING ROUTINES.
C      NNXIN,NNYIN = NUMBER OF X OR Y BINS
C      XLIN = LEFT LIMIT FOR X
C      XHIN = RIGHT LIMIT FOR X
C      YLIN = BOTTOM LIMIT FOR Y
C      YHIN = TOP LIMIT FOR Y
C         (THERE IS NO REQUIREMENT FOR THE SIGN OF THE BINWIDTH,
C          BUT IT MUST BE NONZERO, I.E. THE LIMITS MUST NOT BE EQUAL)
C      MAXIMUM NUMBER OF BINS IS GIVEN BY MAXX, MAXY, AND MAXXY, BELOW
C                            SET UP THE PLOTS
      COMMON/PLTCM1/ STATS(3,3),XPROJ(100),YPROJ(100),ARRAY(5000)
      COMMON /PLOTCM/ NNX,XL,XH,DX,NNY,YL,YH,DY
      EQUIVALENCE (IVAL,XVAL)
C
C                            GET THE RIGHT TYPE FOR INPUT VALUES
C                            NNX IS INTEGER
      COMMON/NTAPES/ITAP1,ITAP2,ITAP3,ITAP4
      DATA LX/1HX/,LY/1HY/
C                  MAXX IS RESTRICTED BY FORMAT STATEMENTS AND
C                  TO WIDTH OF PAPER.  MAXY AND MAXXY ARE LIMITED
C                  ONLY BY THE DIMENSION STATEMENTS IN COMMON.
      DATA MAXX/100/,MAXY/100/,MAXXY/5000/
      IVAL = NNXIN
*  This is not portable Fortran, and apparently is needed by masa anyway.
*      IF (XVAL.EQ.0.0 .OR. IABS(IVAL).GT.250 000 000) IVAL = XVAL
      NNX = IVAL
C                            XL IS REAL
      XVAL = XLIN
*      IF (XVAL.NE.0.0 .AND. IABS(IVAL).LE.250 000 000) XVAL = IVAL
      XL = XVAL
C                            XH IS REAL
      XVAL = XHIN
*      IF (XVAL.NE.0.0 .AND. IABS(IVAL).LE.250 000 000) XVAL = IVAL
      XH = XVAL
C                            NNY IS INTEGER
      IVAL = NNYIN
*      IF (XVAL.EQ.0.0 .OR. IABS(IVAL).GT.250 000 000) IVAL = XVAL
      NNY = IVAL
C                            YL IS REAL
      XVAL = YLIN
*      IF (XVAL.NE.0.0 .AND. IABS(IVAL).LE.250 000 000) XVAL = IVAL
      YL = XVAL
C                            YH IS REAL
      XVAL = YHIN
*      IF (XVAL.NE.0.0 .AND. IABS(IVAL).LE.250 000 000) XVAL = IVAL
      YH = XVAL
C                            CHECK VALUES FOR OBVIOUS ERROR
      IF (XL.EQ.XH) GO TO 55
      IF (YL.EQ.YH) GO TO 60
      IF (NNX.LE.0 .OR. NNX.GT.MAXX) GO TO 40
      IF (NNY.LE.0 .OR. NNY.GT.MAXY) GO TO 40
      NNXY = NNX*NNY
      IF (NNXY.GT.MAXXY) GO TO 40
C                             THEY LOOK OK
      DX = (XH-XL)/FLOAT(NNX)
      DY = (YH-YL)/FLOAT(NNY)
      DO 10 I=1,NNXY
 10   ARRAY(I)=0.
      DO 21 J=1,3
         DO 20 I=1,3
 20      STATS(I,J)=0.
 21   CONTINUE
      DO 22 I=1,NNX
 22   XPROJ(I)=0.
      DO 24 I=1,NNY
 24   YPROJ(I)=0.
      RETURN
C           ---------- ERRORS FROM HERE ON -------------
 40   WRITE(ITAP3,50) NNX,NNY,MAXX,MAXY,MAXXY
 50   FORMAT (' THE PLOTTING ROUTINE WAS GIVEN NNX=',I4,/
     X   '   AND NNY=',I4,'.  IT REQUIRES NNX FROM 1 TO',I4/
     X   '   NNY FROM 1 TO',I4,' AND NNX*NNY.LE.',I5)
      STOP
 55   LY = LX
      YL = XL
 60   WRITE(ITAP3,70) LY,YL
 70   FORMAT (' THE PLOTTING ROUTINE WAS GIVEN EQUAL UPPER ',
     X   ' AND LOWER LIMITS FOR ',A1,/
     X   ' OF',G12.5,'.  ZERO BINWIDTH IS ILLEGAL.')
      STOP
      END
      SUBROUTINE PLOT(XVAL,YVAL)
      COMMON/PLTCM1/ STATS(3,3),XPROJ(100),YPROJ(100),ARRAY(5000)
      COMMON /PLOTCM/ NNX,XL,XH,DX,NNY,YL,YH,DY
      EQUIVALENCE (X,IXX),(Y,IYY)
      DATA WEIGHT /1.0/
C
C                            CHECK FOR INTEGER INPUT
         X=XVAL
         IF (X.NE.0. .AND. IABS(IXX).LE.250 000 000)
     X      X=FLOAT(IXX)
         Y=YVAL
         IF (Y.NE.0. .AND. IABS(IYY).LE.250 000 000)
     X      Y=FLOAT(IYY)
      IX = (X-XL)/DX
      IY = (Y-YL)/DY
C                            FIND AREA (1-2-3 / 4-5-6 / 7-8-9)
      IXI = 2
      IYI = 2
      IF (IY .GE. NNY) IYI = 1
      IF (IY .LT.   0) IYI = 3
      IF (IX .GE. NNX) IXI = 3
      IF (IX .LT.   0) IXI = 1
C                            INCREMENT PLOT STATISTICS
      STATS(IXI,IYI) = STATS(IXI,IYI) + WEIGHT
C                            QUIT IF NOT IN BOUNDS BOTH WAYS
      IF (IXI.NE.2 .OR. IYI.NE.2) R E T U R N
C                            BUMP PLOT CONTENTS
      NN = IX + NNX*IY
      ARRAY(NN+1) = ARRAY(NN+1) + WEIGHT
      XPROJ(IX+1) = XPROJ(IX+1) + WEIGHT
      YPROJ(IY+1) = YPROJ(IY+1) + WEIGHT
      RETURN
      END
      SUBROUTINE PPRINT
      COMMON/PLTCM1/ STATS(3,3),XPROJ(100),YPROJ(100),ARRAY(5000)
      COMMON /PLOTCM/ NNX,XL,XH,DX,NNY,YL,YH,DY
      COMMON/NTAPES/ITAP1,ITAP2,ITAP3,ITAP4
C
      INTEGER LINE1(116),LINE2(116),LINE3(116)
      REAL WD(120)
      EQUIVALENCE (WD(1),IWD(1))
C                            CHARACTERS, FORMAT A4
      DIMENSION BINED(4),IPROJX(5)
C                    CHARACTERS
      INTEGER CHAR(35),ICH1(15),ICH2(15),ICH3(15),ICHRAY(15)
C                            SIZE OF ICH ARRAYS IS ALSO IN NCHAR, BELOW
      INTEGER BLANK,SPLUS,SMINUS,ASTER,EQUALS,EYE,BINED,SLASH
      INTEGER CRVCHR,CRVCH2,IWD(116),IPROJY(4,3)
C
      DATA IOVER /0/
      DATA BLANK,SPLUS,SMINUS,ASTER /1H ,1H+,1H-,1H*/
      DATA EQUALS,EYE,SLASH /1H=,1HI,1H//
C                             CHARACTERS FOR SHADED PLOTS
      DATA NCHAR /15/
      DATA ICH1 /
     D 1H=,1H+,1HH,1H*,1HW,1HO,1HH,1HO,1HO,1HO,1HH,1HO,1HO,1HO,1HO/
      DATA ICH2 /
     D 1H ,1H ,1H ,1H ,1H ,1H-,1H=,1H=,1H*,1HI,1HI,1HX,1HX,1HX,1HX/
      DATA ICH3 /
     D 1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H-,1H=,1H*,1HI/
C                             CHARACTERS FOR NORMAL PLOTS
      DATA CHAR /1H+,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,
     X  1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,
     X  1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ/
C
C                            MISC CHARACTERS
      DATA IPROJY /1HP,1HR,1HO,1HJ,1HO,1HN,1H ,1HY,1HA,1HX,1HI,1HS/
C                            FORMAT A4
      DATA IPROJX /1H ,1H ,4HPROJ,4HON X,4HAXIS/
      DATA BINED /4HLEFT,4HBIN ,4HEDGE,1H /
      DATA HNONE /4HNONE/
C
      CRVCHR = ASTER
      CRVCH2 = EQUALS
C                            CHECK FLAG FOR SHADED PLOTTING
      IF (IOVER.EQ.0) GO TO 15
C                            SHADED PLOTS.  FIND MAX BIN CONTENTS
         CRVCHR = SMINUS
         NXY = NNX*NNY
         VALMAX = NCHAR
         DO 5 IXY=1,NXY
            VALMAX = AMAX1(VALMAX,ARRAY(IXY))
 5       CONTINUE
C                            SET UP ARRAY
         N = NCHAR
         LAST = VALMAX + 1.0
 8       INC = (LAST-2)/N + 1
         LAST = LAST - INC
         ICHRAY(N) = LAST
         N = N - 1
         IF (N.NE.0) GO TO 8
C           --------- DO TITLES ---------
C                            FIRST LINE
 15   WRITE(ITAP3,20)
 20   FORMAT (1H )
C                            PLOT STATISTICS
      WRITE(ITAP3,30) STATS
 30   FORMAT( 18H PLOT STATISTICS =,3F5.0,3X,1H/,3F5.0,3X,1H/,3F5.0)
 40   CONTINUE
      IF (STATS(2,2).EQ.0.) GO TO 200
C               -------------- PRINT CONTENTS OF X CHANNELS  ------
      DO 50 I=1,116
      IWD(I) = BLANK
 50   LINE1(I) = BLANK
      L = NNX + 2
      DO 100 J=1,5
         CALL VPRINT(XPROJ,LINE1,NNX,2,6,0,J,JPRT)
         IF (J.LT.3) GO TO 95
         DO 90 I=1,4
 90      LINE1(L+I) = IPROJY(I,J-2)
         GO TO 98
 95      IF (JPRT.EQ.0) GO TO 100
 98      WRITE(ITAP3,310) IPROJX(J),LINE1
 100  CONTINUE
C                            PUT ON PLOT BOUNDARY
      DO 105 I=NNX,116
      LINE3(I) = BLANK
      LINE2(I) = BLANK
 105  LINE1(I) = BLANK
      DO 110 I=1,NNX
 110  IWD(I) = SMINUS
      DO 111 I=1,NNX,10
 111  IWD(I) = SPLUS
C                            IWD CONTAINS CHARACTERS.  IT WILL BE USED
C                            AGAIN FOR LOWER AXIS(SO DONT USE IWD OR WD)
      WRITE(ITAP3,330) YH,BLANK,IWD
C         --------- MAKE THE PLOT ----------
      UFLAG = UFUN1(XL,YL)
      VFLAG = UFUN2(XL,YL)
C                            J COUNTS LINES, IY COUNTS FROM NNY TO 1
      IY = NNY + 1
      YSTEPS = FLOAT(NNY)
      Y1 = YH
      DO 140 J=1,NNY
         IY = IY - 1
         LL = NNX*(NNY-J)
         Y2 = Y1
         YSTEPS = YSTEPS - 1.
         Y1 = YL + YSTEPS*DY
         LOVER2 = 0
         LOVER3 = 0
C                            LOOP OVER X-BINS
         XSTEPS = 0.
         X2 = XL
C                            SET BIN COUNTER FOR USER CURVE
         IUX = 0
         IVX = 0
         DO 130 I=1,NNX
            XSTEPS = XSTEPS + 1.
            X1 = X2
            X2 = XL + XSTEPS*DX
            LL = LL + 1
            NMBR=ARRAY(LL)
            IF (NMBR.EQ.0) GO TO 120
            IF (IOVER.NE.0) GO TO 114
               NMBR = MIN0(NMBR,35)
               LINE1(I)=CHAR(NMBR)
               GO TO 122
 114        DO 116 K=1,NCHAR
               IF (NMBR.LE.ICHRAY(K)) GO TO 118
 116        CONTINUE
            K = NCHAR
 118        LINE1(I) = ICH1(K)
            LINE2(I) = ICH2(K)
            IF (LINE2(I).NE.BLANK) LOVER2 = 1
            LINE3(I) = ICH3(K)
            IF (LINE3(I).NE.BLANK) LOVER3 = 1
            GO TO 122
 120        LINE1(I) = BLANK
            LINE2(I) = BLANK
            LINE3(I) = BLANK
 122        IF (UFLAG.EQ.HNONE) GO TO 127
               IF (IUX.EQ.I) GO TO 125
C                            LAST BIN WASNT TREATED. GET U-VALUES
                  U1 = SIGN(1.,UFUN1(X1,Y1))
                  U2 = SIGN(1.,UFUN1(X1,Y2))
 125           U1N = SIGN(1.,UFUN1(X2,Y1))
               U2N = SIGN(1.,UFUN1(X2,Y2))
               IF (U1N.EQ.U2 .AND. U2N.EQ.U1) GO TO 126
                  LINE1(I) = CRVCHR
                  LINE2(I) = BLANK
                  LINE3(I) = BLANK
C                            PLACE RIGHT-SIDE VALUES TO BE
C                            LEFT-SIDE ON NEXT TIME THRU LOOP
 126           U1 = U1N
               U2 = U2N
C                            AND SET BIN COUNTER FOR NEXT BIN
               IUX = I+1
 127        IF (VFLAG.EQ.HNONE) GO TO 130
               IF (IVX.EQ.I) GO TO 128
C                            LAST BIN WASNT TREATED. GET U-VALUES
                  V1 = SIGN(1.,UFUN2(X1,Y1))
                  V2 = SIGN(1.,UFUN2(X1,Y2))
 128           V1N = SIGN(1.,UFUN2(X2,Y1))
               V2N = SIGN(1.,UFUN2(X2,Y2))
               IF (V1N.EQ.V2 .AND. V2N.EQ.V1) GO TO 129
                  LINE1(I) = CRVCH2
                  LINE2(I) = BLANK
                  LINE3(I) = BLANK
C                            PLACE RIGHT-SIDE VALUES TO BE
C                            LEFT-SIDE ON NEXT TIME THRU LOOP
 129           V1 = V1N
               V2 = V2N
C                            AND SET BIN COUNTER FOR NEXT BIN
               IVX = I+1
 130     CONTINUE
C                            CHARACTER FOR RIGHT+LEFT BORDER
         IAXIS=EYE
         IF(MOD(IY,10).EQ.1) IAXIS=SPLUS
         LINE1(NNX+1) = IAXIS
C                            ENCODE VALUE FOR Y-AXIS PROJECTION
         L = NNX + 2
         DO 132 I=1,5
         CALL VPRINT(YPROJ(IY),LINE1(L),1,2,6,0,I,JPRT)
 132     L=L+1
         IF (IOVER.EQ.0) GO TO 138
C                            PUT ON SYMBOL DICTIONARY
            IF (J - (NCHAR+1)) 136,133,138
C                            (NO MORE LINES. MAKE SURE SPACES ARE BLANK)
 133        DO 134 I=L,116
               LINE1(I) = BLANK
               LINE2(I) = BLANK
 134           LINE3(I) = BLANK
            GO TO 138
C                            PUT IN VALUE AND SYMBOL FOR DICTIONARY
 136        L = L + 3
            DO 137 I=1,3
            CALL VPRINT(ICHRAY(J),LINE1(L),1,1,3,0,I,JPRT)
 137        L = L + 1
            LINE1(L) = SLASH
            LINE1(L+1) = ICH1(J)
            LINE2(L+1) = ICH2(J)
            IF (ICH2(J).NE.BLANK) LOVER2 = 1
            LINE3(L+1) = ICH3(J)
            IF (ICH3(J).NE.BLANK) LOVER3 = 1
C                            PUT OUT A LINE
 138     WRITE(ITAP3,330) Y1,IAXIS,LINE1
C                            AND ANY OVERPRINTING
         IF (LOVER2.NE.0) WRITE(ITAP3,335) LINE2
         IF (LOVER3.NE.0) WRITE(ITAP3,335) LINE3
C                            LOOP FOR NEXT LINE
 140  CONTINUE
C                            PUT ON LOWER AXIS
      WRITE(ITAP3,310) BLANK,IWD
C                            PUT ON X-BIN EDGE VALUES
      DO 150 I=NNX,116
 150  LINE1(I) = BLANK
      X2 = 0.
      DO 160 I=1,NNX
      WD(I)=XL + X2*DX
 160  X2 = X2 + 1.
      I = 0
      DO 170 J=1,10
         CALL VPRINT(WD,LINE1,NNX,2,10,4,J,JPRT)
         IF (JPRT.NE.0) GO TO 165
C                     NO CHARACTERS THIS LINE.  LEAVE DO-LOOP IF THRU
         IF (I.NE.0) GO TO 200
         GO TO 170
 165     IF (I.LT.4) I=I+1
         WRITE(ITAP3,310) BINED(I),LINE1
 170  CONTINUE
 200  R E T U R N
C
 310  FORMAT (3X,A4,6X,116A1)
 330  FORMAT(1X,F10.4,   1X,A1,116A1)
 335  FORMAT (1H+,12X,116A1)
      END
      SUBROUTINE VPRINT(VIN,VOUT,NB,MODE,JW,JD,JENT,JPRINT)
C                            VIN = ARRAY OF NB INPUT NUMBERS
C                            VOUT = ARRAY OF NB OUTPUT CHARACTERS
C                            NB = NUMBER OF COLUMNS (BINS) TO DECODE.
C                            MODE   1 = INTEGER INPUT AND OUTPUT
C                                   2 = REAL IN, F OUT
C                                   3 = REAL IN, E OUT
C
C                            JW = PLACES
C                            JD = DECIMAL PLACES (F AND E FORMAT)
C                            JENT = THIS PLACE (VPRINT SHOULD BE CALLED
C                                   JW TIMES, WITH JENT=1,JW
C                            JPRINT = 0 FOR ALL BLANK LINE, 1 OTHERWISE
      DIMENSION IADD(3),TENPWR(21)
C            CHARACTERS
      INTEGER VOUT(NB),DIGIT(10),BLANK,SMINUS,PERIOD,SPLUS
C            INPUT VALUES
      REAL VIN(NB),INPXX
      EQUIVALENCE(XX,IXX,INPXX)
      DATA TENPWR /
     X    1.E0,1.E1,1.E2,1.E3,1.E4,1.E5,1.E6,1.E7,1.E8,1.E9,1.E10
     X   ,1.E11,1.E12,1.E13,1.E14,1.E15,1.E16,1.E17,1.E18,1.E19,1.E20/
      DATA DIGIT /1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/
      DATA BLANK,SPLUS,SMINUS,PERIOD / 1H ,1H+,1H-,1H./
      DATA IADD/1,0,-3/
C                            ACCUR1=MACHINE ACCURACY. ACCUR2=0.1/ACCUR1
      DATA ACCUR2/3.0E3/,ACCUR1/3.3E-5/
C                            SMALL IS LIMIT FOR CHECK FOR ZERO
      DATA SMALL /1.E-40/
C
      JPRINT=0
      KD=JD
      IF(MODE.EQ.1) KD=0
C                            FAC5 IS LOWEST SIG FIG
      FAC5 = 1./TENPWR(KD+1)
C                            KPOW IS CURRENT PLACE (1 UNITS, 2 TENS,...)
      KPOW=JW-KD-JENT+IADD(MODE)
      IF (KPOW) 10,200,20
   10 KPOW = KPOW + 1
C                            RND IS HALF OF LOWEST SIG FIG
   20 RND = 0.4999*FAC5
C                            FAC2 IS 10. IN CURRENT DIGIT.
C     FAC2 = 10.**KPOW
      I=IABS(KPOW)
      FAC2=1.
 22   IF (I.LT.20) GO TO 25
      I=I-20
      FAC2=FAC2*TENPWR(21)
      GO TO 22
 25   FAC2=FAC2*TENPWR(I+1)
      IF (KPOW.LT.0) FAC2=1./FAC2
C                            FAC6 IS 0.1 IN CURRENT DIGIT
C                            (FOR CHECKING IF NEXT IS FIRST DIGIT,
C                            SO THIS IS THE PLACE FOR THE MINUS SIGN)
      FAC6 = 0.01*FAC2
      FAC8 = FAC6 * 5.
C                            FAC1 BRINGS CURRENT DIGIT TO UNITS PLACE
      FAC1=10.0001/FAC2
C                            FAC7 IS LARGEST NUMBER FOR WHICH THE
C                            CURRENT DIGIT IS SIGNIFICANT.
      FAC7 = FAC2*ACCUR2
C                            SIG FIG IS IN UNITS PLACE
C                            NOW LOOP OVER ALL NB NUMBERS
      DO 110 IB=1,NB
      INPXX=VIN(IB)
      IF (XX.NE.0.) GO TO 28
         IF (KPOW.EQ.1) GO TO 27
            VOUT(IB)=BLANK
            GO TO 110
   27    VOUT(IB)=DIGIT(1)
         GO TO 100
   28 GO TO (50,60,30),MODE
C                              FIGURE EXPONENTIAL FORMAT
   30 IF(ABS(XX).GT.SMALL) GO TO 31
      XX=0.0
      ITMP=0
      GO TO 32
   31 TMP=ALOG10(ABS(XX))
      IF(TMP.LT.0.0) TMP=TMP-1.0
      ITMP = IFIX(TMP+0.000001)
   32 IF(JENT-JW+2) 33,37,40
C                              PUT IN VALUE TIMES EXPONENTIAL BIAS
C                              SO I WANT XX=XX*10.0**(-ITMP)
   33 XXFAC=1.
      I=IABS(ITMP)
   34 IF (I.LT.20) GO TO 35
      I=I-20
      XXFAC=XXFAC*TENPWR(21)
      GO TO 34
   35 XXFAC=XXFAC*TENPWR(I+1)
      IF (ITMP.GT.0) XXFAC=1./XXFAC
      XX=XX*XXFAC
      GO TO 60
C                            PUT IN SIGN OF EXP
   37 VOUT(IB)=SPLUS
      IF(ITMP.LT.0) VOUT(IB)=SMINUS
      GO TO 100
C                            PUT IN EXPONENTIAL
 40   ITMP = IABS(ITMP)
      IF (JENT.EQ.JW) GO TO 42
C                            TENS PLACE
      JX = ITMP/10
      GO TO 44
C                            UNITS PLACE
 42   JX = MOD (ITMP,10)
 44   VOUT(IB) = DIGIT(JX+1)
      GO TO 100
C                            INTEGER FORMAT.  INSTALL REAL
   50 XX = FLOAT(IXX)
C                            VALUE NOW IN XX FOR INTERPRETATION
   60 X=ABS(XX)
C                            CHECK SIGNIFICANCE
      IF (X.GT.FAC7) GO TO 68
      X = X + RND
C                            KNOCK OFF HIGH DIGITS
      XM = AMOD(X,FAC2)
      IF (XM.LT.FAC8) GO TO 63
C                            BRING TO UNITS PLACE AND TRUNCATE
      JX = XM*FAC1
C                              ACCEPT IF NOT A ZERO
      IF(JX.NE.0) GO TO 95
C                            ITS A ZERO (OR BLANK)
   63 VOUT(IB) = DIGIT(1)
C                                 JUMP IF PRINTING TENS PLACE OR GREATER
C                                 (TO CHECK FOR - AND DELETE LEADING 0)
C                            ACCEPT IF PRINTING UNITS
      IF (KPOW-1) 65,100,70
C                            (ACCEPTING KPOW=0 WOULD ALLOW TENTHS PLACE)
C                            CHECK IF THERE ARE FURTHER NON-ZERO
C                            SIGNIFICANT DIGITS.
   65 IF (XM.LT.FAC5) GO TO 110
      IF (XM.LT.X*ACCUR1) GO TO 110
      GO TO 100
C                            DIGIT BELOW MACHINE ACCURACY.
   68 VOUT(IB) = BLANK
      GO TO 110
C                              HUNDREDS PLACE OR GREATER
C                            SKIP IF HIGHER PLACE EXISTS
   70 IF(X.GE.FAC2) GO TO 100
      VOUT(IB)=BLANK
C                                 CHECK FOR - SIGN
      IF (XX.GE.0.0) GO TO 110
      IF (KPOW.EQ.2) GO TO 90
      IF (X.LT.FAC6) GO TO 110
   90 VOUT(IB)=SMINUS
      GO TO 100
C                            GET CORRESPONDING CHARACTER
   95 VOUT(IB)=DIGIT(JX+1)
  100 JPRINT=1
  110 CONTINUE
      R E T U R N
C
C                   KPOW=0.  THIS PLACE FOR DECIMAL POINT.
 200  DO 210 I=1,NB
 210  VOUT(I) = PERIOD
C                            SET PRINT FLAG IF EXPONENTIAL FORMAT
      IF (MODE.NE.2) GO TO 230
C                            SET PRINT FLAG IF ANY SUCCEEDING
C                            CHARACTERS WILL BE NON-BLANK
      DO 220 I=1,NB
         INPXX=VIN(I)
         X = ABS(XX)
         XM = AMOD(X,1.)
         IF (XM.LT.FAC5) GO TO 220
         IF (XM.GE.X*ACCUR1) GO TO 230
 220  CONTINUE
C                            NO PRINT FLAG
      RETURN
C                            SET PRINT FLAG
 230  JPRINT = 1
      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
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.GE.11) 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
      SUBROUTINE GENPTS(NVECT,NDIM,Y)
C
C
      DIMENSION Y(1),U(2)
C
      N=NVECT*NDIM
      DO 1000 J=1,N,2
  200 CONTINUE
      CALL RANUMS(U,2)
      V1=2.*U(1)-1.
      V2=2.*U(2)-1.
      S=V1*V1+V2*V2
      IF(S.GE.1.) GO TO 200
      S=SQRT(-2.*ALOG(S)/S)
      Y(J)=V1*S
      IF((J+1).LE.N)Y(J+1)=V2*S
1000  CONTINUE
      RETURN
      END
C
C
      SUBROUTINE RANUMS(X,N)
      REAL X(N)
C     GENERATE N UNIFORM RANDOM NUMBERS IN THE INTERVAL (0.0, 1.0).
C     URAND -- A PORTABLE RANDOM NUMBER GENERATOR, BY MICHAEL A.
C     MALCOLM AND CLEVE B. MOLER, STANFORD COMPUTER SCIENCE REPORT
C     STAN-CS-73-334, JANUARY 1974.
C     NOTE THAT THE STATISTICAL PROPERTIES OF THIS ROUTINE HAVE NOT
C     BEEN EXTENSIVELY VERIFIED ON A LARGE NUMBER OF ARCHITECTURES.
C     THIS IS A UNIFORM RANDOM NUMBER GENERATOR BASED ON THEORY AND
C     SUGGESTIONS GIVEN IN D.E. KNUTH (1969), VOL 2.
      INTEGER IA, IC, ITWO, IY, M2, M
      DOUBLE PRECISION HALFM, DATAN, DSQRT
      DATA M2 / 0 /, ITWO / 2 /, IY /123456789/
      IF (M2 .NE. 0) GO TO 20
C IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORK LENGTH
      M = 1
   10 M2 = M
      M = ITWO * M2
      IF (M .GT. M2) GO TO 10
      HALFM = M2
C COMPUTE MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD
      IA = 8 * IDINT(HALFM*DATAN(1.0D0)/8.0D0) + 5
      IC = 2 * IDINT(HALFM*(0.5D0-DSQRT(3.0D0)/6.0D0)) + 1
C S IS THE SCALE FACTOR FOR CONVERTING TO FLOATING POINT
      S = 0.5 / HALFM
C NOW COMPUTE NEXT RANDOM NUMBER
   20 DO 30 I=1,N
      IY = IY * IA + IC
C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE THE WORD LENGTH FOR
C ADDITION IS GREATER THAN THAT FOR MULITPLICATION
      IF (IY/2 .GT. M2) IY = (IY - M2) - M2
C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER OVERFLOW
C AFFECTS THE SIGN BIT
      IF (IY .LT. 0) IY = (IY + M2) + M2
   30 X(I) = FLOAT(IY) * S
      RETURN
      END
C
C
C
      SUBROUTINE SMOOTH(X,Y,N,BANDW,SORTX,SMO,R3,ITAG)
      DIMENSION X(N),Y(N),SORTX(N),SMO(N),R3(N),ITAG(N)
C
C     SORT OBS IN INCREASING VALUES OF X
C
      DO 10 I=1,N
         ITAG(I)=I
         SORTX(I)=X(I)
   10 CONTINUE
      CALL SORT(SORTX,ITAG,1,N)
      DO 20 I=1,N
         SMO(I)=Y(ITAG(I))
   20 CONTINUE
C
C     COMPUTE RUNNING MEDIANS OF 3 AND FILL IN R3
C
      N1=N-1
      DO 30 I=2,N1
         Y1=SMO(I-1)
         Y2=SMO(I)
         Y3=SMO(I+1)
         IF(Y1.LE.Y2)GOTO 40
         TEMP=Y1
         Y1=Y2
         Y2=TEMP
   40    CONTINUE
         RMED=Y2
         IF(Y2.LE.Y3)GOTO 50
         RMED=Y3
         IF(Y3.GT.Y1)GOTO 50
         RMED=Y1
   50    CONTINUE
         R3(I)=RMED
   30 CONTINUE
      R3(1)=SMO(1)
      R3(N)=SMO(N)
C
C     FIT RUNNING STRAIGHT LINES
C
      IBW=N*BANDW
      K=(IBW-1)/2
      IF(K.LT.1)K=1
C
C     INITIALIZE UPDATE
      SX=0.
      SY=0.
      SXX=0.
      SXY=0.
C
      DO 60 I=1,K
         SX=SX+SORTX(I)
         SY=SY+R3(I)
         SXX=SXX+SORTX(I)**2
         SXY=SXY+SORTX(I)*R3(I)
   60 CONTINUE
      IBWEFF=K
C
C     NOW LOOP OVER POINTS TO SMOOTH
C
      DO 100 I=1,N
         IF(I.GT.N-K)GOTO 70
C
C        ADD OBSERVATION TO THE RIGHT
C
         SX=SX+SORTX(I+K)
         SY=SY+R3(I+K)
         SXX=SXX+SORTX(I+K)**2
         SXY=SXY+SORTX(I+K)*R3(I+K)
         IBWEFF=IBWEFF+1
C
   70    CONTINUE
         IF(I.LE.K+1)GOTO 80
C
C        DELETE OBSERVATION ON THE LEFT
C
         SX=SX-SORTX(I-K-1)
         SY=SY-R3(I-K-1)
         SXX=SXX-SORTX(I-K-1)**2
         SXY=SXY-SORTX(I-K-1)*R3(I-K-1)
         IBWEFF=IBWEFF-1
C
   80    CONTINUE
C
C        NOW COMPUTE SLOPE AND INTERCEPT
C
         B=0.
         DENOM=SXX-SX**2/FLOAT(IBWEFF)
         IF(DENOM.NE.0.)B=(SXY-SX*SY/FLOAT(IBWEFF))/DENOM
         A=(SY-B*SX)/FLOAT(IBWEFF)
         SMO(I)=A+B*SORTX(I)
  100 CONTINUE
      RETURN
      END
