C @(#)tdmregr.for	17.1.1.1 (ES0-DMD) 01/25/02 17:47:15
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  17:54 - 11 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C.IDENTIFICATION        TDMREGR.FOR
C.KEYWORDS           TABLE, APPLICATIONS
C.ENVIRONMENT  MIDAS
C.PURPOSE
C     THE SUBROUTINE DOES A MULTIPLE REGRESION ANALYZES ON THE MATRIX
C     'A(N,N)' WITH THE WEIGHTS 'WM(N)'
C
C.ALGORITHM
C
C M.A.EFROYMSON IN
C A.RALSTON & H.S.WILF,
C MATHEMATICAL METHODS FOR DIGITAL COMPUTERS, VOL 1
C JOHN WILEY & SONS, 1967
C
C
C------------------------------------------------------------------

      SUBROUTINE TDMLRG(A,N,WM,TW,COE,RMS,F1,F2,TOL)
      IMPLICIT NONE
      INTEGER N  !   IN : number of variables (dep + indep.)
      DOUBLE PRECISION A(N,N)  !   IN : matrix with correlation coeffs.
      DOUBLE PRECISION WM(N)  !   IN : weighted sum of variables
      DOUBLE PRECISION TW  !   IN : total weight
      DOUBLE PRECISION COE(N)  !   OUT: regression coeffs.
      REAL RMS(N)  !   OUT : standard errors on the coeffs.
      REAL F1  !   IN :  F value to enter variable
      REAL F2  !   IN :  F value to remove variable
      DOUBLE PRECISION TOL ! IN : the zero level of the matrix
C
      INTEGER NM1,I,J,IP1,NO,NMIN,NMAX,K
      DOUBLE PRECISION DF1,DF2,FI,SY,BO,BX,VMIN,VMAX,VI,AKK
      DOUBLE PRECISION B(10),RO(10),S(10)
C
C     INITIATE COUNTERS AND ARRAYS
C
      IF (N.LT.2 .OR. N.GT.10) RETURN
      NM1    = N - 1
      DO 10 I = 1,N
          S(I)   = 0.0D0
          COE(I) = 0.D0
          RMS(I) = 0.0
   10 CONTINUE
      IF (F1.LT.F2) F1     = F2
      DF1    = DBLE(F1)
      DF2    = DBLE(F2)
C
C     NORMALIZE MATRIX AND WEIGHTED SUMS
C
      DO 30 I = 1,N
          DO 20 J = 1,N
              A(I,J) = A(I,J) - WM(I)*WM(J)/TW
   20     CONTINUE
   30 CONTINUE
      DO 40 I = 1,N
          WM(I)  = WM(I)/TW
          IF (A(I,I).LT.0.0D0) RETURN
          RO(I)  = DSQRT(A(I,I))
   40 CONTINUE
      DO 60 I = 1,NM1
          IP1    = I + 1
          DO 50 J = IP1,N
              A(I,J) = A(I,J)/ (RO(I)*RO(J))
              A(J,I) = A(I,J)
   50     CONTINUE
   60 CONTINUE
      DO 70 I = 1,N
          A(I,I) = 1.0D0
   70 CONTINUE
C
C     DO THE MATRIX MANIPULATIONS
C
      FI     = TW - 1.0D0
      NO     = 0
   80 IF (FI.LE.0.0D0) RETURN
      SY     = RO(N)*DSQRT(A(N,N)/FI)
      NO     = NO + 1
      IF (NO.GT.1000) GO TO 190
      VMIN   = 1.0D30
      VMAX   = 0.0D0
      NMIN   = 0
      NMAX   = 0
C
      DO 110 I = 1,NM1
          B(I)   = 0.0D0
          S(I)   = 0.0D0
          IF (A(I,I).LE.TOL) GO TO 110
          VI     = A(I,N)*A(N,I)/A(I,I)
          IF (VI) 100,110,90
   90     IF (VI.LE.VMAX) GO TO 110
          VMAX   = VI
          NMAX   = I
          GO TO 110

  100     B(I)   = A(I,N)*RO(N)/RO(I)
          S(I)   = SY*DSQRT(A(I,I))/RO(I)
          IF (DABS(VI).GE.DABS(VMIN)) GO TO 110
          VMIN   = VI
          NMIN   = I
  110 CONTINUE
C
      BX     = 0.0D0
      DO 120 I = 1,NM1
          BX     = BX + B(I)*WM(I)
  120 CONTINUE
      BO     = WM(N) - BX
      IF (DABS(VMIN)*FI.GE.DF2*A(N,N)) GO TO 130
      K      = NMIN
      FI     = FI + 1.0D0
      GO TO 140

  130 IF (VMAX* (FI-1.0D0).LE.DF1* (A(N,N)-VMAX)) GO TO 190
      K      = NMAX
      FI     = FI - 1.0D0
C
  140 CONTINUE
      AKK    = A(K,K)
      DO 160 I = 1,N
          IF (I.EQ.K) GO TO 160
          DO 150 J = 1,N
              IF (J.EQ.K) GO TO 150
              A(I,J) = A(I,J) - A(I,K)*A(K,J)/AKK
  150     CONTINUE
  160 CONTINUE
C
      DO 170 I = 1,N
          IF (I.EQ.K) GO TO 170
          A(I,K) = -A(I,K)/AKK
  170 CONTINUE
      DO 180 J = 1,N
          IF (J.EQ.K) GO TO 180
          A(K,J) = A(K,J)/AKK
  180 CONTINUE
      A(K,K) = 1.0D0/AKK
      GO TO 80
C
C     REGRATION ANALYZES IS FINISHED GET RESULTS OUT
C
  190 CONTINUE
      DO 200 I = 1,NM1
          COE(I) = B(I)
          RMS(I) = SNGL(S(I))
  200 CONTINUE
      COE(N) = BO
      RMS(N) = SNGL(SY)
      RETURN

      END