C @(#)tdregr.for	17.1.1.1 (ES0-DMD) 01/25/02 17:47:16
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.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  18:24 - 11 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.IDENTIFICATION        TDREGR.FOR
C.KEYWORDS           TABLE, APPLICATIONS
C.ENVIRONMENT  MIDAS
C.PURPOSE
C LINEAR REGRESSION
C.ALGORITHM
C
C HOUSEHOLDER TRANSFORMATION
C REF.: CH.L.LAWSON AND R.J.HANSON
C       SOLVING LEAST SQUARES PROBLEMS
C       PRENTICE HALL INC. 1974
C
C
C------------------------------------------------------------------

      SUBROUTINE TDSET1(LL,X1,Y,NIND,G,X,N,M)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.PURPOSE
C
C PREPARE DATA IN WORK SPACE FOR MULTIPLE LINEAL REGRESSION
C
C------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER LL
      INTEGER NIND
      DOUBLE PRECISION X1(NIND)
      DOUBLE PRECISION Y
      INTEGER N
      INTEGER M
      DOUBLE PRECISION G(M,M)
      DOUBLE PRECISION X(M)
C
      INTEGER IP,L1

      IP     = 1
      G(LL,IP) = 1.D0
      DO 10 L1 = 1,NIND
          IP     = IP + 1
          G(LL,IP) = G(LL,IP-1)*X1(L1)
   10 CONTINUE
      G(LL,N+1) = Y
      RETURN

      END

      SUBROUTINE TDSET2(LL,X1,X2,Y,K,L,G,X,N,M)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.PURPOSE
C
C PREPARE DATA IN WORK SPACE FOR 1 OR 2D POLYNOMIAL
C
C------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER LL
      DOUBLE PRECISION X1
      DOUBLE PRECISION X2
      DOUBLE PRECISION Y
      INTEGER K
      INTEGER L
      INTEGER N
      INTEGER M
      DOUBLE PRECISION G(M,M)
      DOUBLE PRECISION X(M)
C
      INTEGER IP,L1,K1
      DOUBLE PRECISION DC

      IP     = 0
      DC     = 1.D0
      DO 20 L1 = 0,L
          IP     = IP + 1
          G(LL,IP) = DC
          DO 10 K1 = 1,K
              IP     = IP + 1
              G(LL,IP) = G(LL,IP-1)*X1
   10     CONTINUE
          DC     = DC*X2
   20 CONTINUE
      G(LL,N+1) = Y
      RETURN

      END

      SUBROUTINE TDHHTR(K,L1,G,X,N,M)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.PURPOSE
C
C     Apply a Householder transformation to the two rows K and L1 of G.
C
C------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER K
      INTEGER L1  !   IN : RANK OF THE MATRIX
      INTEGER N
      INTEGER M
      DOUBLE PRECISION G(M,M)
      DOUBLE PRECISION X(M)
C
      INTEGER KH,J
      DOUBLE PRECISION S,H,B
C
C ... CONSTRUCT THE TRANSFORMATION
C
      S      = DSQRT(G(K,K)**2+G(L1,K)**2)
      IF (G(K,K).GE.0.D0) S      = -S
      H      = G(K,K) - S
      G(K,K) = S
      B      = G(K,K)*H
C
C ... APPLY THE TRANSFORMATION TO THE ROWS K AND L1 OF G.
C
      KH     = N + 1 - K
      IF ((KH.NE.0) .AND. (B.NE.0.D0)) THEN
          DO 10 J = 1,KH
              S      = (G(K,K+J)*H+G(L1,K+J)*G(L1,K))/B
              G(K,K+J) = G(K,K+J) + S*H
              G(L1,K+J) = G(L1,K+J) + S*G(L1,K)
   10     CONTINUE
      END IF

      RETURN

      END

      SUBROUTINE TDSOLV(G,X,N,M)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.PURPOSE
C
C
C                Solve the set of linear equations which are
C                represented by the upper triangular part
C                of the matrix G(1:n,1:n+1).
C
C------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER N
      INTEGER M
      DOUBLE PRECISION G(M,M)  
      DOUBLE PRECISION X(M)
C
      INTEGER I,II,J
      DOUBLE PRECISION S
C
      DO 10 I = 1,N
          X(I)   = G(I,N+1)
   10 CONTINUE
      DO 30 II = 1,N
          I      = N - II + 1
          S      = 0.D0
          IF (I.NE.N) THEN
              DO 20 J = 2,N + 1 - I
                  S      = S + G(I,I+J-1)*X(I+J-1)
   20         CONTINUE
          END IF

          X(I)   = (X(I)-S)/G(I,I)
   30 CONTINUE
      RETURN

      END