Sources Esope | Cast3M

glsq
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