C  This is the matrix sweep algorithm, written by Dr. Newton.
C  This is what I use for linear regression and matrix inversion.
C
      subroutine sweep(a,mdim,m,k1,k2,ier)
c*****************************************************************
c
c   Subroutine to sweep the mxm matrix a on its k1 st thru k2 th
c   diagonals.  mdim is the row dimension of a in the calling
c   routine. ier is output and is 1 if the algorithm fails.
c
c******************************************************************
c
      implicit real*8(a-h,o-z)
      real*8 a(mdim,mdim)
      integer iwk(100)
c
      ier=0
      do 1 k=k1,k2
 1       iwk(k)=k
      do 50 k=k1,k2
c
c   find first nonzero diagonal element and swap index with kth:
c
         do 70 kk=k,k2
            k3=iwk(kk)
 70      if(abs(a(k3,k3)).gt.1.e-25) go to 71
         go to 99
 71      kc=iwk(k)
         iwk(k)=iwk(kk)
         iwk(kk)=kc
         kk=iwk(k)
         d=1./a(kk,kk)
         a(kk,kk)=1.
         do 10 i=1,m
 10         a(kk,i)=d*a(kk,i)
            do 20 j=1,m
               if(j.eq.kk) go to 20
               a(j,kk)=-a(j,kk)*d
 20         continue
            do 40 j=1,m
               if(j.eq.kk) go to 40
               do 30 i=1,m
                  if(i.eq.kk) go to 30
                  a(j,i)=a(j,i)+a(j,kk)*a(kk,i)/d
 30            continue
 40         continue
 50    continue
       return
 99    ier=1
       return
       end