parameter(iii=145)
ccccc********   iii=f1*100=f2*100  ***************
      parameter(n=5)
      parameter(m=13)
cc    dimension xx((n-1)*m),b(n-1),x(n,m),sgm(n)
cc    dimension y(m),kti(n-1),kyi(n-1),a(n,n),xba(n)
      dimension xx(52),b(4),x(n,m),sgm(n)
      dimension y(m),kti(4),kyi(4),a(n,n),xba(n)
      data xx/7.0,1.0,11.0,11.0,7.0,11.0,3.0,1.0,2.0,21.0,
     *       1.0,11.0,10.0,26.0,29.0,
     *        56.0,31.0,52.0,55.0,71.0,31.0,54.0,47.0,
     *       40.0,66.0,68.0,6.0,15.0,8.0,
     *        8.0,6.0,9.0,17.0,22.0,18.0,4.0,23.0,9.0,
     *        8.0,60.0,52.0,20.0,47.0,33.0,
     *        22.0,6.0,44.0,22.0,26.0,34.0,12.0,
     *        12.0/
      data y/78.5,74.3,104.3,87.6,95.9,109.2,
     *        102.7,72.5,93.1,115.9,83.8,113.3,109.4/
      open(6,file='reg2.dat')
      call regres(n,m,xx,y,b,x,kti,kyi,a,sgm,xba,iii)
      close(6)
      stop
      end
c******************************
      subroutine regres(n,m,xx,y,b,x,kti,kyi,a,sgm,xba,iii)
cc    dimension xx((n-1)*m),b(n-1),x(n,m),sgm(n)
cc    dimension y(m),kti(n-1),kyi(n-1),a(n,n),xba(n)
      dimension xx(52),b(4),x(n,m),sgm(n)
      dimension y(m),kti(4),kyi(4),a(n,n),xba(n)
      do  101 j=1,m
      do  100 i=1,n-1
100   x(i,j)=xx((i-1)*m+j)
101   x(n,j)=y(j)
      write(6,95)
95    format(//4x,'original data')
      call prin(x,n,m)
1000  fisn=float(iii)/100.0
      call  pp(n,m,x,xba,sgm,h)
      call rela(n,m,x,xba,sgm,a)
      write(6,96)
96    format(/4x,'original correlation matrix')
      call prin(a,n,n)
      do  200 i=1,n-1
      kti(i)=i
200   kyi(i)=0
      l=0
      c=0.0
888   if (l.eq.n-1) goto 999
      call im(kti,0,v1,k,n,a)
      f1=v1*(m-l-2)/(a(n,n)-v1)
      write(6,7) k,v1,f1
7     format(2x,'kmax=',i2,'   vmax=',f10.4,'   f1=',f10.4)
      if( f1.ge.fisn) then
      write(6,8) k
8     format(1x,'  introduce factor',i4)
      kti(k)=0
      kyi(k)=k
      l=l+1
      call  ma(n,m,l,k,a,r,yn,1,h)
      else
      goto   999
      end if
898    if(l.eq.1)  goto  888
      call  im(kyi,1,v1,k,n,a)
      f1=-v1*(m-l-1)/a(n,n)
      write(6,11) k,v1,f1
 11   format(2x,'kmin=',i2,'  vmin=',f10.4,'   f2=',f10.4)
      if(f1.lt.fisn) then
      write(6,12) k
 12    format(1x,'   cancel factor',i4)
      kti(k)=k
      kyi(k)=0
       l=l-1
      call  ma(n,m,l,k,a,r,yn,1,h)
      goto  898
      else
      goto  888
      end if
999   do  400 i=1,n-1
      if (kyi(i).ne.0) then
      b(i)=a(i,n)*sgm(n)/sgm(i)
      c=b(i)*xba(i)+c
      end if
400    continue
      s1=xba(n)-c
      write(6,777) s1
777   format(2x,'a',1x,f10.4)
666   format(2x,'b',i3,f10.4)
      do  500 i=1,n-1
       if(kyi(i).ne.0) write(6,666) i,b(i)
500   continue
      write(6,333) fisn
333   format(/2x,'fisn=',f10.4)
      write(6,555) r
555   format(/2x,'compound relative coefficient r',f10.4)
      write(6,444) yn
444   format(/2x,'residual error yn',3x,f10.4)
      do  605  j=1,m
      d1=s1
      do  600  i=1,n-1
      if (kyi(i).ne.0) d1=b(i)*x(i,j)+d1
600   continue
605   y(j)=d1
      call prin2(x,y,n,m)
       return
       end
c*******************************************************
c     to translate the matrix a(l) into the matrix a(l+1)
      subroutine ma(n,m,l,k,a,r,yn,iiia,h)
      dimension a(n,n)
      ers=10.0e-5
      if(a(k,k).lt.ers) write(6,5) k,a(k,k)
5     format(2x,'k=',i2,'   a(k,k)=',e10.2,' !!!!!!!!!!!!')
      do  10 i=1,n
      do  10  j=1,n
      if(i.ne.k.and.j.ne.k) a(i,j)=a(i,j)-a(i,k)*a(k,j)/a(k,k)
10    continue
      do  20  j=1,n
       if (j.ne.k) then
       a(k,j)=a(k,j)/a(k,k)
       a(j,k)=-1.0*a(j,k)/a(k,k)
       end if
20    continue
      a(k,k)=1.0/a(k,k)
       r=sqrt(1.0-a(n,n))
       yn=h*sqrt(a(n,n)/(m-l-1))
       if  (iiia.ne.0) then
       write(6,*) '  transformed matrix:'
      do  30  i=1,n
      write(6,14) (a(i,j),j=1,n)
14    format(1x,10f7.2)
30     continue
       end if
       return
       end
c*******************************************************
c     to print out two array
      subroutine prin2(x,y,n,m)
      dimension x(n,m),y(m)
      do 110  ib=1,m,10
      ie=ib+9
       if (ie-m) 2,2,1
  1   ie=m
  2   write(6,33) (x(n,i),i=ib,ie)
      write(6,44) (y(i),i=ib,ie)
110    continue
 33   format(/1x,'observed value y',10f6.1)
 44   format(/1x,'estimate value y',10f6.1)
       return
       end
c*******************************************************
c     to print out one array
      subroutine prin(a,n,m)
      dimension a(n,m)
      do 100  ib=1,m,10
      ie=ib+9
       if (ie-m) 2,2,1
  1   ie=m
  2   write(6,22) (i,i=ib,ie)
      do  150  j=1,n
      write(6,24)j,(a(j,k),k=ib,ie)
150    continue
100    continue
 22   format(/1x,10i7)
 24   format(/1x,i2,10f7.2)
       return
       end
c*******************************************************
c     to calculate the corrlation matrix
      subroutine rela(n,m,x,xba,sgm,a)
      dimension x(n,m),xba(n),sgm(n),a(n,n)
      do  10  i=1,n-1
      do  20  j=i+1,n
       c=0.0
       do  30  l=1,m
      c=(x(i,l)-xba(i))*(x(j,l)-xba(j))+c
30     continue
      a(i,j)=c/(sgm(i)*sgm(j))
      a(j,i)=a(i,j)
20    continue
      a(i,i)=1.0
10     continue
       a(n,n)=1.0
       return
       end
c*******************************************************
c     to calculate average value and variance
      subroutine pp(n,m,x,xba,sgm,h)
      dimension x(n,m),xba(n),sgm(n)
      do  10  i=1,n
       c=0.0
       do  20  l=1,m
      c=x(i,l)+c
20     continue
      xba(i)=c/m
       c=0.0
       do  30  l=1,m
      c=c+(x(i,l)-xba(i))**2
30     continue
      sgm(i)=sqrt(c)
10     continue
        h=sgm(n)
       return
       end
c**************************************************
c     to choose the maximun (or minimun) v(i)
      subroutine im(ks,ii,vi,k,n,a)
ccc   dimension ks(n-1),a(n,n)
      dimension ks(4),a(n,n)
      integer ii
      vi=-10.0
      do  10  i=1,n-1
      if(ks(i).eq.0) then
      else
      v2=a(i,n)*a(i,n)/a(i,i)
      if(ii.eq.1) v2=-v2
       if(v2.ge.vi) then
       vi=v2
       k=i
      end if
      end if
10     continue
       return
       end