C GLSQ      SOURCE    CHAT      05/01/13    00:18:31     5004           
*********************************************************************           
      SUBROUTINE GLSQ(A,X,IL,N,M,ALPHA,E1,E2,M50)                               
C                                                                               
C     LEAST SQUARES SOLUTION OF A LINEAR SYSTEM                                 
C-----------------------------------------------------------------------        
C                                                                               
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)                                                  
      DIMENSION A(M50,M50), X(*)                                                
      INTEGER IL(*)                                                             
C                                                                               
C-----------------------------------------------------------------------        
C                                                                               
      MPLUS1 = M + 1                                                            
      I = 1                                                                     
C                                                                               
      DO 50 K = 1, M                                                            
      IL(K) = 0                                                                 
      IF(I - N) 10, 40, 50                                                      
C                                                                               
   10 IPLUS1 = I + 1                                                            
      DO 30 J = IPLUS1, N                                                       
      IF(ABS(A(J,K)) .LE. E1) GO TO 30                                          
      T1 = SQRT(A(J,K)**2 + A(I,K)**2)                                          
      S = A(J,K) / T1                                                           
      C = A(I,K) / T1                                                           
C                                                                               
      DO 20 L = K, MPLUS1                                                       
      T2     = + C * A(I,L) + S * A(J,L)                                        
      A(J,L) = - S * A(I,L) + C * A(J,L)                                        
      A(I,L) = T2                                                               
   20 CONTINUE                                                                  
C                                                                               
   30 CONTINUE                                                                  
C                                                                               
   40 IF(ABS(A(I,K)) .LE. E2) GO TO 50                                          
      IL(K) = I                                                                 
      I = I + 1                                                                 
   50 CONTINUE                                                                  
C                                                                               
      SUM = 0.0                                                                 
      IF(I .GT. N) GO TO 70                                                     
C                                                                               
      DO 60 J = I, N                                                            
      SUM = SUM + A(J,MPLUS1)**2                                                
   60 CONTINUE                                                                  
C                                                                               
   70 ALPHA = SQRT(SUM)                                                         
      X(MPLUS1) = - 1.0                                                         
C                                                                               
      DO 90 J = 1, M                                                            
      I = MPLUS1 - J                                                            
      X(I) = 0.0                                                                
      L = IL(I)                                                                 
      IF(L .EQ. 0) GO TO 90                                                     
      IPLUS1 = I + 1                                                            
      SUM = 0.0                                                                 
C                                                                               
      DO 80 K = IPLUS1, MPLUS1                                                  
      SUM = SUM + A(L,K) * X(K)                                                 
   80 CONTINUE                                                                  
C                                                                               
      X(I) = - SUM / A(L,I)                                                     
   90 CONTINUE                                                                  
C                                                                               
      RETURN                                                                    
      END