c
c  compute means and STD of columns of numbers
c
      Parameter (no=25, np=(no*(no+1))/2, nog=np)

      Real*8 x(nog),zm(nog),sd(nog),dx,s,si,dx2,smi
      real*8 eigwrk(nog),eigvec(no,no),eigval(no),c,d,a,p
      integer nob


      open(10,file='ttMi.d',form='unformatted',
     x status='old')
      open(15,file='newRI.d',form='formatted',status='unknown')
      open(17,file='diaRI.d',form='formatted',status='unknown')
      zm=0.d0
      sd=0.d0
      nob=0
c
c  initialize arrays to zero
c
 20   read(10,end=100)ix,nm,x
      if(nm.ne.4)go to 20
c      write(17,1701)ix,(x(k),k=1,5)
 1701 format(1x,i8,5f12.5)
c      if(ix.lt.1000)go to 20
      nob=nob+1
      s=nob
      si=1.d0/s
      smi=1.d0
      if(nob.gt.1)smi=1.d0/(s-1.d0)
      
      do 25 i=1,np
         dx = x(i) - zm(i)
         zm(i) = zm(i) + dx*si
         dx2 = dx*(x(i) - zm(i)) - sd(i)
         sd(i)=sd(i) + dx2*smi
 25   CONTINUE
      GO TO 20
c
c
 100  write(17,4001)nog,nob
 4001 format(1x,i10,' Number of variables'/
     x 1x,i10,' Number of records', i10)
      close(10)
c
      call eigen(no,no,zm,eigval,eigwrk,eigvec)
      nneg=0
      dx=0.0
      do 32 k=1,no
         if(eigval(k).lt.0.d0)then
            nneg=nneg+1
            dx=dx+eigval(k)+eigval(k)
         endif
 32   continue
      dx2 = (dx*dx*100.d0)+1.d0
      if(nneg.gt.0)then
         write(17,4301)nneg
 4301 format(1x,i10,' negative eigenvalues')   
         p=eigval(nneg+1)
         do 124 k=1,nneg
            c=eigval(k)
            eigval(k)=p*(dx-c)*(dx-c)/dx2
 124    continue
        m=0
        do 125 k=1,nt
           do 125 L=k,nt
              m=m+1
              a=0.d0
              do 127 ma=1,nt
                 a=a+eigvec(k,ma)*eigval(ma)*eigvec(L,ma)
 127          continue
          zm(m) = a
 125   continue
      endif
c
c  save results
c
      m = 0
      do 105 i=1,no
         do 106 j=i,no
            m = m + 1
        dx = dsqrt(sd(m))
        write(15,4002)i,j,zm(m),dx
 4002    format(1x,2i5,2f18.10)
 106   continue
       ma=ihmssf(i,i,no)
       dx=dsqrt(sd(ma))
       write(17,4002)i,i,zm(ma),dx
 105  CONTINUE
      close(15)
      close(17)
      stop
      end