text:         subroutine gspiv(b,n,m,eps,rnorms,ipivots,ncols)
          implicit double precision (a-h,o-z)
          save
          dimension rnorms(1),ipivots(1)
          double precision b(n,m),cd
  c      
  c       Set rank to zero just in case.
  c     
          ncols=0
  c
  c        . . . initialize the array of pivots
  c 
          do 1100 i=1,m
          ipivots(i)=i
   1100 continue
  c 
  c       . . . prepare the array of values of norms
  c             of columns
  c 
          done=1
          dtot=0
          do 1400 i=1,m
  c
          d=0
          do 1200 j=1,n
          d=d+b(j,i)*(b(j,i))
   1200 continue
          rnorms(i)=sqrt(d)
          dtot=dtot+d
   1400 continue
  c
          thresh = eps
  
          if (dtot .lt. eps**2) then
          ncols = 0
          return
          endif
  
  c$$$        thresh=(dtot)*eps**2
  c$$$        thresh=sqrt(thresh)
  
  
  cccccccccccccccccccccccccccccccccccccccccccc
  
  c
  c
  
  c 
  c       . . . conduct gram-schmidt iterations
  c 
          do 4000 i=1,min(m,n)
  c 
  c       find the pivot
  c 
          ipivot=i
          rn=rnorms(i)
  c 
          do 2200 j=i+1,m
          if(rnorms(j) .le. rn) goto 2200
          rn=rnorms(j)
          ipivot=j
   2200 continue
   2400 continue
  c 
  c       put the column number ipivot in the i-th place
  c 
          do 2600 j=1,n
          cd=b(j,i)
          b(j,i)=b(j,ipivot)
          b(j,ipivot)=cd
   2600 continue
  c 
          iijj=ipivots(i)
          ipivots(i)=ipivots(ipivot)
          ipivots(ipivot)=iijj
  c 
          d=rnorms(ipivot)
          rnorms(ipivot)=rnorms(i)
          rnorms(i)=d
  c 
  c       orthogonalize the i-th column to all preceding ones
  c 
          if(i .eq. 1) goto 2790
          do 2780 j=1,i-1
  c 
          call gsrleascap(b(1,i),b(1,j),n,cd)
  c 
          do 2770 l=1,n
          b(l,i)=b(l,i)-b(l,j)*cd
   2770 continue
   2780 continue
   2790 continue
  c 
  c       normalize the i-th column
  c 
          call gsrleascap(b(1,i),b(1,i),n,cd)
  c 
          d=cd
          if(d .lt. thresh**2 ) return
  c 
          ncols=i
  c 
          d=done/sqrt(d)
          do 2800 j=1,n
          b(j,i)=b(j,i)*d
   2800 continue
  c 
          if(i .eq. m) goto 3400
  c 
  c        orthogonalize everything else to it
  c 
          do 3200 j=i+1,m
  c
          if(rnorms(j) .lt. thresh) goto 3200
  c 
          call gsrleascap(b(1,i),b(1,j),n,cd)
  c 
          cd=(cd)
  c 
          rrn=0
          do 3000 l=1,n
          b(l,j)=b(l,j)-b(l,i)*cd
          rrn=rrn+b(l,j)*(b(l,j))
   3000 continue
          rnorms(j)=sqrt(rrn)
   3200 continue
   3400 continue
  c 
   4000 continue
  c 
          return
          end
  c
  c
  c
          subroutine gsrleascap(x,y,n,prod)
          implicit double precision (a-h,o-z)
          save
          dimension x(1),y(1)
  c 
          prod=0
          do 1200 i=1,n
          prod=prod+x(i)*(y(i))
   1200 continue
          return
          end
  c
  c
  c
          subroutine gspiv_move(n,a,b)
          implicit double precision (a-h,o-z)
          dimension a(1),b(1)
          do 1000 j=1,n
          b(j)=a(j)
   1000 continue
          end       
  
  
  
  