
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**********************************
C**   LSM_EQNSOL                 **
C**********************************
C
C **** LSM SERVICE ROUTINE ONLY ****
C
      SUBROUTINE LSM_EQNSOL (N,M,C,G,WT,X,SIGX,IERR)
C
      INCLUDE 'implicit_none.finc'
      INCLUDE 'limits.finc'
C
C Purpose: Solve equations for multiple deconvolution
C          (both sides also multiplied by weighting factor)
C
C
C Author:  John W. Campbell, March 1994
C
C Arguments:
C
      INTEGER N, M, IERR
      REAL C(MXEQN,*),G(*),WT(*),X(*),SIGX(*)

C PARAMETERS 
C
C C(1,1)*X(1) + C(1,2)*X(2) + ... C(1,N)*X(N) = G(1)
C C(2,1)*X(1) + C(2,2)*X(2) + ... C(2,N)*X(N) = G(2)
C      .            .                  .         .
C      .            .                  .         .
C C(M,1)*X(1) + C(M,2)*X(2) + ... C(M,N)*X(N) = G(M)
C        N (I)   NUMBER OF UNKNOWNS
C        M (I)   NUMBER OF EQUATIONS
C        C M(I)   EQUATIONS COEFFICIENTS MATRIX
C     MXEQN (I)   FIRST DIMENSION OF 'C' ARRAY
C        G (I)   OBSERVED VALUES
C       WT (I)   WEIGHTS OF OBSERVEED VALUES
C        X (O)   RETURNS X VALUES
C     SIGX (O)   RETURNS SIGMA(X) VALUES
C     IERR (O)   ERROR FLAG 0 = OK
C                           1 = SINGULAR MATRIX
C                           -100 = FEWER EQUATIONS THAN TERMS
C                           -200 = NO UNKNOWNS
C
C
C     PARAMETER MXEQN   MAXIMUM NO. OF EQUATIONS ALLOWED FOR A REFLECTION
C     PARAMETER MXTRM   MAXIMUM NO> OF TERMA ALLOWED IN A EQUATION
C

C
C====== Parameter statements
C 
C
C====== GLOBALS:
C
C
C====== FUNCTION TYPES:
C
C
C====== LOCALS:
C
      INTEGER I, J, K
      DOUBLE PRECISION A(MXEQN,MXTRM),AT(MXTRM,MXEQN)
      DOUBLE PRECISION AN(MXTRM,MXTRM),AI(MXTRM,MXTRM),RHS(MXTRM)
      DOUBLE PRECISION DET

C
C INITIALISE
C
      IERR = 0
      IF (N.LE.0) THEN
         IERR = -200
         RETURN
      END IF
      IF (M.LT.N) THEN
         IERR = -100
         RETURN
      END IF
C
C SET UP 'A' MATRIX
C
      DO 20 I = 1, M
         DO 10 J = 1, N
            A(I,J) = C(I,J)
 10      CONTINUE
 20   CONTINUE
C
C GET TRANSPOSE OF 'A' IN 'AT'
C
      DO 40 I = 1, N
         DO 30 J = 1, M
            AT(I,J) = A(J,I)
30       CONTINUE
40    CONTINUE
C
C GET WEIGHTED 'A'
C
      DO 45 I = 1, M
         DO 42 J = 1, N
            A(I,J) = WT(I)*A(I,J)
 42      CONTINUE
 45   CONTINUE
C
C SET NORMAL MATRIX
C
      DO 70 I = 1,N
         DO 60 J = 1,N
            AN(I,J) = 0.0
            DO 50 K = 1, M
               AN(I,J)=AN(I,J)+AT(I,K)*A(K,J)
 50         CONTINUE   
 60      CONTINUE
 70   CONTINUE
C
C GET RHS ELEMENTS
C
      DO 90 I =1,N
         RHS(I) = 0.0
         DO 80 K = 1, M
            RHS(I) = RHS(I) + AT(I,K)*WT(K)*G(K)
 80      CONTINUE
 90   CONTINUE
C
C INVERT NORMAL MATRIX
C
      CALL INVMAT(AN,AI,N,MXTRM,6,DET,IERR)
      IF (IERR.NE.0) RETURN
C
C GET SOLUTIONS
C
      DO 120 I = 1, N
         X(I) = 0.0
         DO 110 K = 1, N
            X(I) = X(I) + AI(I,K)*RHS(K)
 110     CONTINUE
 120  CONTINUE
C
C GET SIGMAS
C
      DO 130 I = 1, N
         SIGX(I) = AI(I,I)
         IF (SIGX(I).GE.0) THEN
            SIGX(I) = SQRT(SIGX(I))
         ELSE
            SIGX(I) = 0.0
         END IF
 130  CONTINUE
      RETURN
      END
