      subroutine ainv21(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=21)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      integer function nrand(sseed,nlow,nhigh)
      implicit real*8 (a-h,o-z)

      call ggubs(sseed,1,uni)

      nrand = nlow + idint(uni*dble(real(nhigh-nlow+1)))

      return
      end

      subroutine ols2(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=2)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv2(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine ols3(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=3)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv3(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine locidx(x,xpts,npts,idx,nlow,nhigh,niter)
      implicit real*8 (a-h,o-z)
      parameter (mxiter=50)
      dimension xpts(npts)

      nhigh = 0
      nlow = 0

      if (x .ge. xpts(npts)) then 
         nhigh = 1
         idxtmp = 0
         niter = 0
      elseif (x .lt. xpts(1)) then 
         nlow = 1
         idxtmp = 0
         niter = 0
      else
         
         ilow = 1
         ihigh = npts
         do 10 i = 1,mxiter
            icur = (ilow+ihigh)/2
            if (x .ge. xpts(icur+1)) then
               ilow = icur
            elseif (x .lt. xpts(icur)) then
               ihigh = icur
            else
               idxtmp = icur
               niter = i
               goto 999
            endif
 10      continue

      endif

999   idx = idxtmp

      return
      end

      double precision function dmean2(x,ndim,n)
      implicit real*8 (a-h,o-z)
      dimension x(ndim)
 
      sum = 0.0D+00
      do 10 i = 1,n
         sum = sum + x(i)
 10   continue
 
      dmean2 = sum/dble(real(n))
 
      return
      end
 
      double precision function dstd2(x,ndim,n)
      implicit real*8 (a-h,o-z)
      dimension x(ndim)
 
      xmean = dmean(x,n)

      sum = 0.0D+00
      do 10 i = 1,n
         sum = sum + (x(i)-xmean)*(x(i)-xmean)
 10   continue
 
      dstd2 = dsqrt(sum/dble(real(n)))
 
      return
      end

      double precision function dmax2(x,ndim,n)
      implicit real*8 (a-h,o-z)
      dimension x(ndim)
      
      xmax = x(1)
      do 10 i = 2,n
         if (x(i) .gt. xmax) xmax = x(i)
 10   continue

      dmax2 = xmax

      return
      end

      double precision function dmin2(x,ndim,n)
      implicit real*8 (a-h,o-z)
      dimension x(ndim)
      
      xmin = x(1)
      do 10 i = 2,n
         if (x(i) .lt. xmin) xmin = x(i)
 10   continue

      dmin2 = xmin

      return
      end

      subroutine ols4(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=4)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv4(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine ols5(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=5)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv5(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine ols6(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=6)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv6(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine ols7(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=7)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv7(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine ols10(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=10)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv10(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine ols8(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=8)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv8(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine ols9(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=9)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv9(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine ols11(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=11)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv11(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine ols1(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=1)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv1(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine ols21(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=21)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv21(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine ols36(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=36)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv36(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end

      subroutine ols37(ydata,xdata,ndata,bhat,rsq,stderr,npts)
      implicit real*8 (a-h,o-z)
      parameter (np=37)
      dimension xdata(ndata,np),ydata(ndata),xpx(np,np),
     *          xpxinv(np,np),xpy(np,1),bhat(np,1),temp1(np,1),
     *          temp2(np,1),xt(np,1),xtp(1,np),xtxtp(np,np),
     *          temp(np,np)

      call inimat(xpx,np,np)
      do 51 i = 1,npts
         do 52 j = 1,np
            xt(j,1) = xdata(i,j)
 52      continue
         call tr(xt,xtp,np,1)
         call mult2(xt,xtp,xtxtp,np,1,np)
         call add2(xpx,xtxtp,temp,np,np)
         call copmat(temp,xpx,np,np)
 51   continue
      call inv37(xpx,xpxinv,nsing)

      call inimat(xpy,np,1)
      do 10 i = 1,npts
         do 20 j = 1,np
            temp1(j,1) = xdata(i,j)*ydata(i)
 20      continue
         call add2(xpy,temp1,temp2,np,1)
         call copmat(temp2,xpy,np,1)
 10   continue
    
      call mult2(xpxinv,xpy,bhat,np,np,1)

      sse = 0.0D+00
      sumy = 0.0D+00
      do 12 i = 1,npts
         diff = ydata(i)
         do 14 j = 1,np
            diff = diff - bhat(j,1)*xdata(i,j) 
 14      continue
         sse = sse + diff*diff
         sumy = sumy + ydata(i)
 12   continue
      ymean = sumy/dble(real(npts))
      stderr = dsqrt(sse/dble(real(npts)))

      tss = 0.0D+00
      do 13 i = 1,npts
         tss = tss + (ydata(i)-ymean)*(ydata(i)-ymean)
 13   continue

      rsq = 1.0D+00 - sse/tss

      return
      end