← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to prim/table/libsrc/tdregr.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/tdregr.for
 
1
C @(#)tdregr.for        19.1 (ES0-DMD) 02/25/03 14:11:19
 
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  18:24 - 11 DEC 1987
 
34
C
 
35
C.LANGUAGE: F77+ESOext
 
36
C
 
37
C.AUTHOR: J.D.PONZ
 
38
C
 
39
C.IDENTIFICATION        TDREGR.FOR
 
40
C.KEYWORDS           TABLE, APPLICATIONS
 
41
C.ENVIRONMENT  MIDAS
 
42
C.PURPOSE
 
43
C LINEAR REGRESSION
 
44
C.ALGORITHM
 
45
C
 
46
C HOUSEHOLDER TRANSFORMATION
 
47
C REF.: CH.L.LAWSON AND R.J.HANSON
 
48
C       SOLVING LEAST SQUARES PROBLEMS
 
49
C       PRENTICE HALL INC. 1974
 
50
C
 
51
C
 
52
C------------------------------------------------------------------
 
53
 
 
54
      SUBROUTINE TDSET1(LL,X1,Y,NIND,G,X,N,M)
 
55
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
56
C
 
57
C.PURPOSE
 
58
C
 
59
C PREPARE DATA IN WORK SPACE FOR MULTIPLE LINEAL REGRESSION
 
60
C
 
61
C------------------------------------------------------------------
 
62
      IMPLICIT NONE
 
63
      INTEGER LL
 
64
      INTEGER NIND
 
65
      DOUBLE PRECISION X1(NIND)
 
66
      DOUBLE PRECISION Y
 
67
      INTEGER N
 
68
      INTEGER M
 
69
      DOUBLE PRECISION G(M,M)
 
70
      DOUBLE PRECISION X(M)
 
71
C
 
72
      INTEGER IP,L1
 
73
 
 
74
      IP     = 1
 
75
      G(LL,IP) = 1.D0
 
76
      DO 10 L1 = 1,NIND
 
77
          IP     = IP + 1
 
78
          G(LL,IP) = G(LL,IP-1)*X1(L1)
 
79
   10 CONTINUE
 
80
      G(LL,N+1) = Y
 
81
      RETURN
 
82
 
 
83
      END
 
84
 
 
85
      SUBROUTINE TDSET2(LL,X1,X2,Y,K,L,G,X,N,M)
 
86
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
87
C
 
88
C.PURPOSE
 
89
C
 
90
C PREPARE DATA IN WORK SPACE FOR 1 OR 2D POLYNOMIAL
 
91
C
 
92
C------------------------------------------------------------------
 
93
      IMPLICIT NONE
 
94
      INTEGER LL
 
95
      DOUBLE PRECISION X1
 
96
      DOUBLE PRECISION X2
 
97
      DOUBLE PRECISION Y
 
98
      INTEGER K
 
99
      INTEGER L
 
100
      INTEGER N
 
101
      INTEGER M
 
102
      DOUBLE PRECISION G(M,M)
 
103
      DOUBLE PRECISION X(M)
 
104
C
 
105
      INTEGER IP,L1,K1
 
106
      DOUBLE PRECISION DC
 
107
 
 
108
      IP     = 0
 
109
      DC     = 1.D0
 
110
      DO 20 L1 = 0,L
 
111
          IP     = IP + 1
 
112
          G(LL,IP) = DC
 
113
          DO 10 K1 = 1,K
 
114
              IP     = IP + 1
 
115
              G(LL,IP) = G(LL,IP-1)*X1
 
116
   10     CONTINUE
 
117
          DC     = DC*X2
 
118
   20 CONTINUE
 
119
      G(LL,N+1) = Y
 
120
      RETURN
 
121
 
 
122
      END
 
123
 
 
124
      SUBROUTINE TDHHTR(K,L1,G,X,N,M)
 
125
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
126
C
 
127
C.PURPOSE
 
128
C
 
129
C     Apply a Householder transformation to the two rows K and L1 of G.
 
130
C
 
131
C------------------------------------------------------------------
 
132
      IMPLICIT NONE
 
133
      INTEGER K
 
134
      INTEGER L1  !   IN : RANK OF THE MATRIX
 
135
      INTEGER N
 
136
      INTEGER M
 
137
      DOUBLE PRECISION G(M,M)
 
138
      DOUBLE PRECISION X(M)
 
139
C
 
140
      INTEGER KH,J
 
141
      DOUBLE PRECISION S,H,B
 
142
C
 
143
C ... CONSTRUCT THE TRANSFORMATION
 
144
C
 
145
      S      = DSQRT(G(K,K)**2+G(L1,K)**2)
 
146
      IF (G(K,K).GE.0.D0) S      = -S
 
147
      H      = G(K,K) - S
 
148
      G(K,K) = S
 
149
      B      = G(K,K)*H
 
150
C
 
151
C ... APPLY THE TRANSFORMATION TO THE ROWS K AND L1 OF G.
 
152
C
 
153
      KH     = N + 1 - K
 
154
      IF ((KH.NE.0) .AND. (B.NE.0.D0)) THEN
 
155
          DO 10 J = 1,KH
 
156
              S      = (G(K,K+J)*H+G(L1,K+J)*G(L1,K))/B
 
157
              G(K,K+J) = G(K,K+J) + S*H
 
158
              G(L1,K+J) = G(L1,K+J) + S*G(L1,K)
 
159
   10     CONTINUE
 
160
      END IF
 
161
 
 
162
      RETURN
 
163
 
 
164
      END
 
165
 
 
166
      SUBROUTINE TDSOLV(G,X,N,M)
 
167
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
168
C
 
169
C.PURPOSE
 
170
C
 
171
C
 
172
C                Solve the set of linear equations which are
 
173
C                represented by the upper triangular part
 
174
C                of the matrix G(1:n,1:n+1).
 
175
C
 
176
C------------------------------------------------------------------
 
177
      IMPLICIT NONE
 
178
      INTEGER N
 
179
      INTEGER M
 
180
      DOUBLE PRECISION G(M,M)  
 
181
      DOUBLE PRECISION X(M)
 
182
C
 
183
      INTEGER I,II,J
 
184
      DOUBLE PRECISION S
 
185
C
 
186
      DO 10 I = 1,N
 
187
          X(I)   = G(I,N+1)
 
188
   10 CONTINUE
 
189
      DO 30 II = 1,N
 
190
          I      = N - II + 1
 
191
          S      = 0.D0
 
192
          IF (I.NE.N) THEN
 
193
              DO 20 J = 2,N + 1 - I
 
194
                  S      = S + G(I,I+J-1)*X(I+J-1)
 
195
   20         CONTINUE
 
196
          END IF
 
197
 
 
198
          X(I)   = (X(I)-S)/G(I,I)
 
199
   30 CONTINUE
 
200
      RETURN
 
201
 
 
202
      END