c****************************************************************       
c                                                                       
c subroutine gauss appears as figure 14.21 in "introduction to          
c computing" by dyck, v.a., lawson, j.d., and smith, j.a.,              
c copyright 1979 by reston publishing company.                          
c                                                                       
C*****************************************************************
C  Slightly Modified for use in UNIX f77 environment, August 15,87
C  The Modifications do not change the nature or programming
C  Philosophy.  (Just makes them compatible with f77)
C*****************************************************************
c subprogram to reduce a system of n equations in n unknowns,           
c t*x = c, to triangular form.                                          
c                                                                       
c*****************************************************************      
c                                                                       
c ndim  - dimensioned(available) size of all arrays                     
c n     - array size(actual) for problem to be solved                   
c t     - import matrix(general) and export matrix(triangular)          
c c     - import and export form of right-hand-side vector              
c flag  - flag = 1 if a diagonal element of (triangular) t is 0         
c max   - largest (in magnitude) subdiagonal element                    
c index - row containing the largest subdiagonal element                
c ratio - ratio of subdiagonal to diagonal entry                        
c                                                                       
c*****************************************************************      
c                                                                       
      subroutine gauss(ndim,n,t,c,flag)                                 
c                                                                       
      integer ndim, n, i, j, k, index, flag                             
      real t(ndim,ndim), c(ndim), max, temp, ratio, abs                 
c                                                                       
c create zeros in subdiagonal positions of columns 1 - n-1.             
c check for zero diagonal elements in the reduction.                    
c                                                                       
      flag = 0                                                          
      i = 1                                                             
1     if(i .le. n-1) then                                              
c                                                                       
c        find the maximum subdiagonal element in column i.              
c                                                                       
         index = i                                                      
         max = abs(t(i,i))                                              
         j = i + 1                                                      
10       if(j .le. n) then                                             
            if(abs(t(j,i)) .gt. max) then                               
                  max = abs(t(j,i))                                     
                  index = j                                             
            end if                                                      
            j = j + 1                                                   
C The Second while loop
         go to 10
	 endif
c                                                                       
c        if index exceeds i, exchange rows index and i.                 
c                                                                       
         if(index .gt. i) then                                          
               k = i                                                    
20             if(k .le. n) then                                       
                  temp = t(i,k)                                         
                  t(i,k) = t(index,k)                                   
                  t(index,k) = temp                                     
                  k = k + 1                                             
C The third while loop
               go to 20
	       endif
               temp = c(i)                                              
               c(i) = c(index)                                          
               c(index) = temp                                          
         end if                                                         
c                                                                       
c        if the maximum pivot element is zero, set flag = 1 and         
c        bypass the row operations at stage i.                          
c                                                                       
         if(max .eq. 0) then                                            
               flag = 1                                                 
               i = i + 1                                                
            else                                                        
c                                                                       
c              for nonzero pivot, create zeros in subdiagonal           
c              positions of column i by row operations.                 
c                                                                       
               j = i + 1                                                
30             if(j .le. n) then                                       
                  k = i + 1                                             
                  ratio = t(j,i)/t(i,i)                                 
300               if(k .le. n) then                                    
C  The innermost while loop, 3rd level nesting
                     t(j,k) = t(j,k) - t(i,k)*ratio                     
                     k = k + 1                                          
                  go to 300
		  endif
                  c(j) = c(j) - c(i)*ratio                              
                  j = j + 1                                             
               go to 30
C  A second level nesting for the while loop
	       endif
               i = i + 1                                                
         end if                                                         
      go to 1
      endif
C     End of the outermost while
      return                                                            
      end