Esope source files

   1 : C GLSQ      SOURCE    CHAT      05/01/13    00:18:31     5004           
   2 : *********************************************************************           
   3 :       SUBROUTINE GLSQ(A,X,IL,N,M,ALPHA,E1,E2,M50)                               
   4 : C                                                                               
   5 : C     LEAST SQUARES SOLUTION OF A LINEAR SYSTEM                                 
   6 : C-----------------------------------------------------------------------        
   7 : C                                                                               
   8 :       IMPLICIT INTEGER(I-N)
   9 :       IMPLICIT REAL*8(A-H,O-Z)                                                  
  10 :       DIMENSION A(M50,M50), X(*)                                                
  11 :       INTEGER IL(*)                                                             
  12 : C                                                                               
  13 : C-----------------------------------------------------------------------        
  14 : C                                                                               
  15 :       MPLUS1 = M + 1                                                            
  16 :       I = 1                                                                     
  17 : C                                                                               
  18 :       DO 50 K = 1, M                                                            
  19 :       IL(K) = 0                                                                 
  20 :       IF(I - N) 10, 40, 50                                                      
  21 : C                                                                               
  22 :    10 IPLUS1 = I + 1                                                            
  23 :       DO 30 J = IPLUS1, N                                                       
  24 :       IF(ABS(A(J,K)) .LE. E1) GO TO 30                                          
  25 :       T1 = SQRT(A(J,K)**2 + A(I,K)**2)                                          
  26 :       S = A(J,K) / T1                                                           
  27 :       C = A(I,K) / T1                                                           
  28 : C                                                                               
  29 :       DO 20 L = K, MPLUS1                                                       
  30 :       T2     = + C * A(I,L) + S * A(J,L)                                        
  31 :       A(J,L) = - S * A(I,L) + C * A(J,L)                                        
  32 :       A(I,L) = T2                                                               
  33 :    20 CONTINUE                                                                  
  34 : C                                                                               
  35 :    30 CONTINUE                                                                  
  36 : C                                                                               
  37 :    40 IF(ABS(A(I,K)) .LE. E2) GO TO 50                                          
  38 :       IL(K) = I                                                                 
  39 :       I = I + 1                                                                 
  40 :    50 CONTINUE                                                                  
  41 : C                                                                               
  42 :       SUM = 0.0                                                                 
  43 :       IF(I .GT. N) GO TO 70                                                     
  44 : C                                                                               
  45 :       DO 60 J = I, N                                                            
  46 :       SUM = SUM + A(J,MPLUS1)**2                                                
  47 :    60 CONTINUE                                                                  
  48 : C                                                                               
  49 :    70 ALPHA = SQRT(SUM)                                                         
  50 :       X(MPLUS1) = - 1.0                                                         
  51 : C                                                                               
  52 :       DO 90 J = 1, M                                                            
  53 :       I = MPLUS1 - J                                                            
  54 :       X(I) = 0.0                                                                
  55 :       L = IL(I)                                                                 
  56 :       IF(L .EQ. 0) GO TO 90                                                     
  57 :       IPLUS1 = I + 1                                                            
  58 :       SUM = 0.0                                                                 
  59 : C                                                                               
  60 :       DO 80 K = IPLUS1, MPLUS1                                                  
  61 :       SUM = SUM + A(L,K) * X(K)                                                 
  62 :    80 CONTINUE                                                                  
  63 : C                                                                               
  64 :       X(I) = - SUM / A(L,I)                                                     
  65 :    90 CONTINUE                                                                  
  66 : C                                                                               
  67 :       RETURN                                                                    
  68 :       END                                                                       
  69 :