← Back to branch summary

~ubuntu-branches/ubuntu/wily/eso-midas/wily-proposed

« back to all changes in this revision

Viewing changes to prim/table/libsrc/tdmregr.for

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
prim/table/libsrc/tdmregr.for
 
1
C @(#)tdmregr.for       19.1 (ES0-DMD) 02/25/03 14:11:18
 
2
C===========================================================================
 
3
C Copyright (C) 1995 European Southern Observatory (ESO)
 
4
C
 
5
C This program is free software; you can redistribute it and/or 
 
6
C modify it under the terms of the GNU General Public License as 
 
7
C published by the Free Software Foundation; either version 2 of 
 
8
C the License, or (at your option) any later version.
 
9
C
 
10
C This program is distributed in the hope that it will be useful,
 
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
C GNU General Public License for more details.
 
14
C
 
15
C You should have received a copy of the GNU General Public 
 
16
C License along with this program; if not, write to the Free 
 
17
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Corresponding concerning ESO-MIDAS should be addressed as follows:
 
21
C       Internet e-mail: midas@eso.org
 
22
C       Postal address: European Southern Observatory
 
23
C                       Data Management Division 
 
24
C                       Karl-Schwarzschild-Strasse 2
 
25
C                       D 85748 Garching bei Muenchen 
 
26
C                       GERMANY
 
27
C===========================================================================
 
28
C
 
29
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
31
C                                         all rights reserved
 
32
C
 
33
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  17:54 - 11 DEC 1987
 
34
C
 
35
C.LANGUAGE: F77+ESOext
 
36
C
 
37
C.AUTHOR: J.D.PONZ
 
38
C.IDENTIFICATION        TDMREGR.FOR
 
39
C.KEYWORDS           TABLE, APPLICATIONS
 
40
C.ENVIRONMENT  MIDAS
 
41
C.PURPOSE
 
42
C     THE SUBROUTINE DOES A MULTIPLE REGRESION ANALYZES ON THE MATRIX
 
43
C     'A(N,N)' WITH THE WEIGHTS 'WM(N)'
 
44
C
 
45
C.ALGORITHM
 
46
C
 
47
C M.A.EFROYMSON IN
 
48
C A.RALSTON & H.S.WILF,
 
49
C MATHEMATICAL METHODS FOR DIGITAL COMPUTERS, VOL 1
 
50
C JOHN WILEY & SONS, 1967
 
51
C
 
52
C
 
53
C------------------------------------------------------------------
 
54
 
 
55
      SUBROUTINE TDMLRG(A,N,WM,TW,COE,RMS,F1,F2,TOL)
 
56
      IMPLICIT NONE
 
57
      INTEGER N  !   IN : number of variables (dep + indep.)
 
58
      DOUBLE PRECISION A(N,N)  !   IN : matrix with correlation coeffs.
 
59
      DOUBLE PRECISION WM(N)  !   IN : weighted sum of variables
 
60
      DOUBLE PRECISION TW  !   IN : total weight
 
61
      DOUBLE PRECISION COE(N)  !   OUT: regression coeffs.
 
62
      REAL RMS(N)  !   OUT : standard errors on the coeffs.
 
63
      REAL F1  !   IN :  F value to enter variable
 
64
      REAL F2  !   IN :  F value to remove variable
 
65
      DOUBLE PRECISION TOL ! IN : the zero level of the matrix
 
66
C
 
67
      INTEGER NM1,I,J,IP1,NO,NMIN,NMAX,K
 
68
      DOUBLE PRECISION DF1,DF2,FI,SY,BO,BX,VMIN,VMAX,VI,AKK
 
69
      DOUBLE PRECISION B(10),RO(10),S(10)
 
70
C
 
71
C     INITIATE COUNTERS AND ARRAYS
 
72
C
 
73
      IF (N.LT.2 .OR. N.GT.10) RETURN
 
74
      NM1    = N - 1
 
75
      DO 10 I = 1,N
 
76
          S(I)   = 0.0D0
 
77
          COE(I) = 0.D0
 
78
          RMS(I) = 0.0
 
79
   10 CONTINUE
 
80
      IF (F1.LT.F2) F1     = F2
 
81
      DF1    = DBLE(F1)
 
82
      DF2    = DBLE(F2)
 
83
C
 
84
C     NORMALIZE MATRIX AND WEIGHTED SUMS
 
85
C
 
86
      DO 30 I = 1,N
 
87
          DO 20 J = 1,N
 
88
              A(I,J) = A(I,J) - WM(I)*WM(J)/TW
 
89
   20     CONTINUE
 
90
   30 CONTINUE
 
91
      DO 40 I = 1,N
 
92
          WM(I)  = WM(I)/TW
 
93
          IF (A(I,I).LT.0.0D0) RETURN
 
94
          RO(I)  = DSQRT(A(I,I))
 
95
   40 CONTINUE
 
96
      DO 60 I = 1,NM1
 
97
          IP1    = I + 1
 
98
          DO 50 J = IP1,N
 
99
              A(I,J) = A(I,J)/ (RO(I)*RO(J))
 
100
              A(J,I) = A(I,J)
 
101
   50     CONTINUE
 
102
   60 CONTINUE
 
103
      DO 70 I = 1,N
 
104
          A(I,I) = 1.0D0
 
105
   70 CONTINUE
 
106
C
 
107
C     DO THE MATRIX MANIPULATIONS
 
108
C
 
109
      FI     = TW - 1.0D0
 
110
      NO     = 0
 
111
   80 IF (FI.LE.0.0D0) RETURN
 
112
      SY     = RO(N)*DSQRT(A(N,N)/FI)
 
113
      NO     = NO + 1
 
114
      IF (NO.GT.1000) GO TO 190
 
115
      VMIN   = 1.0D30
 
116
      VMAX   = 0.0D0
 
117
      NMIN   = 0
 
118
      NMAX   = 0
 
119
C
 
120
      DO 110 I = 1,NM1
 
121
          B(I)   = 0.0D0
 
122
          S(I)   = 0.0D0
 
123
          IF (A(I,I).LE.TOL) GO TO 110
 
124
          VI     = A(I,N)*A(N,I)/A(I,I)
 
125
          IF (VI) 100,110,90
 
126
   90     IF (VI.LE.VMAX) GO TO 110
 
127
          VMAX   = VI
 
128
          NMAX   = I
 
129
          GO TO 110
 
130
 
 
131
  100     B(I)   = A(I,N)*RO(N)/RO(I)
 
132
          S(I)   = SY*DSQRT(A(I,I))/RO(I)
 
133
          IF (DABS(VI).GE.DABS(VMIN)) GO TO 110
 
134
          VMIN   = VI
 
135
          NMIN   = I
 
136
  110 CONTINUE
 
137
C
 
138
      BX     = 0.0D0
 
139
      DO 120 I = 1,NM1
 
140
          BX     = BX + B(I)*WM(I)
 
141
  120 CONTINUE
 
142
      BO     = WM(N) - BX
 
143
      IF (DABS(VMIN)*FI.GE.DF2*A(N,N)) GO TO 130
 
144
      K      = NMIN
 
145
      FI     = FI + 1.0D0
 
146
      GO TO 140
 
147
 
 
148
  130 IF (VMAX* (FI-1.0D0).LE.DF1* (A(N,N)-VMAX)) GO TO 190
 
149
      K      = NMAX
 
150
      FI     = FI - 1.0D0
 
151
C
 
152
  140 CONTINUE
 
153
      AKK    = A(K,K)
 
154
      DO 160 I = 1,N
 
155
          IF (I.EQ.K) GO TO 160
 
156
          DO 150 J = 1,N
 
157
              IF (J.EQ.K) GO TO 150
 
158
              A(I,J) = A(I,J) - A(I,K)*A(K,J)/AKK
 
159
  150     CONTINUE
 
160
  160 CONTINUE
 
161
C
 
162
      DO 170 I = 1,N
 
163
          IF (I.EQ.K) GO TO 170
 
164
          A(I,K) = -A(I,K)/AKK
 
165
  170 CONTINUE
 
166
      DO 180 J = 1,N
 
167
          IF (J.EQ.K) GO TO 180
 
168
          A(K,J) = A(K,J)/AKK
 
169
  180 CONTINUE
 
170
      A(K,K) = 1.0D0/AKK
 
171
      GO TO 80
 
172
C
 
173
C     REGRATION ANALYZES IS FINISHED GET RESULTS OUT
 
174
C
 
175
  190 CONTINUE
 
176
      DO 200 I = 1,NM1
 
177
          COE(I) = B(I)
 
178
          RMS(I) = SNGL(S(I))
 
179
  200 CONTINUE
 
180
      COE(N) = BO
 
181
      RMS(N) = SNGL(SY)
 
182
      RETURN
 
183
 
 
184
      END