C @(#)lsqsol.for	17.1.1.1 (ES0-DMD) 01/25/02 17:15:42
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENTIFICATION
C  subroutine LSQSOL         version 1.0       870130
C  A. Kruszewski             Obs. de Geneve
C.PURPOSE
C  Solves a set of normal equations given by array A.
C  Results are stored in array X. The last element of array X
C  contains a sum of squared deviations. Mean errors are stored in
C  array S. The last element of array S contains the mean error
C  of an observation of unit weight. The dimension N1 is equal to
C  number of unknowns plus one.
C.INPUT/OUTPUT
C  input arguments
C  A           real*8 array     symmetrical array
C  N1          integer*4        dimension of array
C  M           integer*4        number of conditional equations
C  output arguments
C  X           real*8 array     unknowns + vv
C  S           real*8 array     mean errors + sigma0
C-----------------------------------------------------------------------
      SUBROUTINE LSQSOL ( N1 , A , M , X , S )
C
      IMPLICIT          NONE
      INTEGER           N1
      DOUBLE PRECISION  A(N1,N1)
      INTEGER           M
      DOUBLE PRECISION  X(N1)
      DOUBLE PRECISION  S(N1)
C
      INTEGER            N, NN
      INTEGER            I, I1
      INTEGER            J, J1
      INTEGER            K
      DOUBLE PRECISION   TEMP
C
C *** N - number of unknowns.
C
      N = N1 - 1
C
C *** Calculate triangular square root of the input array.
C
      IF ( A(1,1) .LE. 0.0 ) THEN
          S(N1) = -1.0
          RETURN
      ENDIF
      A(1,1) = SQRT( A(1,1) )
      DO 10 I = 2 , N1
          I1 = I - 1
          DO 20 J = 1 , I1
              IF ( A(J,J) .LE. 0.0 ) THEN
                  S(N1) = -1.0
                  RETURN
              ENDIF 
              TEMP = A(I,J)
              J1 = J - 1
              DO 30 K = 1 , J1
                  TEMP = TEMP - A(I,K) * A(J,K)
   30         CONTINUE
              A(I,J) = TEMP / A(J,J)
   20     CONTINUE
          TEMP = A(I,I)
          DO 40 K = 1 , I1
              TEMP = TEMP - A(I,K) * A(I,K)
   40     CONTINUE
          IF ( TEMP .LE. 0.0 ) THEN
              S(N1) = -1.0
              RETURN
          ENDIF
          A(I,I) = SQRT( TEMP )
   10 CONTINUE
C
C *** Invert the triangular square root matrix of size N. 
C
      DO 50 I = 1 , N
          A(I,I) = 1.0 / A(I,I)
   50 CONTINUE
      NN = N - 1
      DO 60 I = 1 , NN
          I1 = I + 1
          DO 70 J = I1 , N
              J1 = J - 1
              TEMP = 0.0
              DO 80 K = I , J1
                  TEMP = TEMP + A(I,K) * A(J,K)   
   80         CONTINUE
              A(I,J) = -TEMP * A(J,J)
   70     CONTINUE
   60 CONTINUE
C
C *** Calculate unknowns. 
C
      DO 90 I = 1 , N
          TEMP = 0.0
          DO 100 J = I , N
              TEMP=TEMP + A(I,J) * A(N1,J)
  100     CONTINUE
          X(I) = SNGL( TEMP )
   90 CONTINUE
C
C *** Calculate vv and sigma0.
C
      TEMP = A(N1,N1)
      X(N1) = SNGL( TEMP * TEMP )
      S(N1) = SNGL( TEMP / SQRT( DBLE( M-N ) ) )
C
C *** Calculate mean errors of unknowns.
C
      DO 110 I = 1 , N
          TEMP = 0.0
          DO 120 J = I , N
              TEMP = TEMP + A(I,J) * A(I,J)
  120     CONTINUE
          S(I) = SNGL( SQRT(TEMP) * S(N1) )
  110 CONTINUE
C
      RETURN
C
      END