← 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/tdregp.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/tdregp.for
 
1
C===========================================================================
 
2
C Copyright (C) 1995-2005 European Southern Observatory (ESO)
 
3
C
 
4
C This program is free software; you can redistribute it and/or 
 
5
C modify it under the terms of the GNU General Public License as 
 
6
C published by the Free Software Foundation; either version 2 of 
 
7
C the License, or (at your option) any later version.
 
8
C
 
9
C This program is distributed in the hope that it will be useful,
 
10
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
C GNU General Public License for more details.
 
13
C
 
14
C You should have received a copy of the GNU General Public 
 
15
C License along with this program; if not, write to the Free 
 
16
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
C MA 02139, USA.
 
18
C
 
19
C Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
C       Internet e-mail: midas@eso.org
 
21
C       Postal address: European Southern Observatory
 
22
C                       Data Management Division 
 
23
C                       Karl-Schwarzschild-Strasse 2
 
24
C                       D 85748 Garching bei Muenchen 
 
25
C                       GERMANY
 
26
C===========================================================================
 
27
C
 
28
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
29
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
30
C                                         all rights reserved
 
31
C
 
32
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  16:02 - 19 NOV 1987
 
33
C
 
34
C.LANGUAGE: F77+ESOext
 
35
C
 
36
C.AUTHOR: J.D.PONZ
 
37
C
 
38
C.IDENTIFICATION
 
39
C
 
40
C  program TDREGP
 
41
C
 
42
C.PURPOSE
 
43
C
 
44
C  Execute the command
 
45
C  REGRESSION/POL table dep_var ind_var1[,ind_var2] deg1[,deg2]  [outkey]
 
46
C  OR
 
47
C  REGRESSION/POL table dep_var ind_var1[,ind_var2,...]   ?      [outkey]
 
48
C
 
49
C.KEYWORDS
 
50
C
 
51
C  Tables, Multiple linear regression, POLYNOMIAL FIT
 
52
C
 
53
C.ALGORITHM
 
54
C
 
55
C  uses table interface routines
 
56
C  fitting algorithm using Householder transformation
 
57
C
 
58
C.VERSION
 
59
C 051019        last modif
 
60
C 
 
61
C-----------------------------------------------------------
 
62
C
 
63
      SUBROUTINE TDREGP
 
64
      IMPLICIT NONE
 
65
C
 
66
      LOGICAL NULL(9), SEL, IPOL, NEXT
 
67
C
 
68
      INTEGER PARVAL
 
69
      INTEGER IVAR(9),IDEG(2),IPAR(13)
 
70
      INTEGER NROW,NCOL, NAC, NAR, I, ISTAT, STATUS
 
71
      INTEGER IAV, TID, NSC1, NN, NVAR, I1, I2, K, L, N, LL, N1
 
72
      INTEGER IROW, KMIN, NPAR, NM1, K1
 
73
      INTEGER INDEX, KUN, KNUL, TYPE, MAXPAR, INDKEY
 
74
      INTEGER NOELM,BYTELM,DEGMAX
 
75
      PARAMETER (MAXPAR=100)
 
76
C
 
77
      REAL ERROR,RPAR1(4),RPAR(13)
 
78
C
 
79
      DOUBLE PRECISION XMIN,XMAX,YMIN,YMAX
 
80
      DOUBLE PRECISION VAL(9)
 
81
      DOUBLE PRECISION G(MAXPAR,MAXPAR),X(MAXPAR)
 
82
C
 
83
      CHARACTER*16 MSG
 
84
      CHARACTER*80 TABLE,LINE,SELE
 
85
      CHARACTER*90 SDUM
 
86
      CHARACTER*17 COLUMN(9)
 
87
      CHARACTER*8  FORM,TYPKEY
 
88
      CHARACTER*20 CPAR, KEYNAM, INPNAM
 
89
C
 
90
      INCLUDE 'MID_INCLUDE:TABLES.INC'
 
91
      INCLUDE 'MID_INCLUDE:TABLED.INC'
 
92
C 
 
93
      DATA MSG/'ERR:TREGRMUL'/
 
94
      DATA PARVAL/5/
 
95
      DATA CPAR/' '/
 
96
      DATA IPAR/13*0/
 
97
      DATA RPAR/13*0./
 
98
C
 
99
C ... GET INPUT PARAMETERS
 
100
C
 
101
      DO 10, I=1,9
 
102
          NULL(I) = .FALSE.
 
103
          IVAR(I) = 0
 
104
          VAL(I) = 0
 
105
   10 CONTINUE
 
106
 
 
107
      CALL TDPGET(PARVAL,NPAR,ISTAT)
 
108
 
 
109
      IF (ISTAT.NE.0) GO TO 70
 
110
      IF (NPAR.LT.4) THEN
 
111
          ISTAT  = ERRPAR
 
112
          GO TO 70
 
113
      END IF
 
114
 
 
115
      TABLE  = TPARBF(1)
 
116
      KEYNAM = TPARBF(5)
 
117
      INDKEY = INDEX(KEYNAM,' ') - 1
 
118
 
 
119
      DO 15 I = 1, MAXPAR
 
120
          X(I) = 0.D0
 
121
   15 CONTINUE
 
122
      INPNAM = KEYNAM(1:INDKEY)//'D'
 
123
      CALL STKFND(INPNAM,TYPKEY,NOELM,BYTELM,ISTAT)
 
124
      CALL STKWRD(INPNAM,X,1,NOELM,KUN,ISTAT)
 
125
 
 
126
      CALL STKRDR('INPUTR',1,2,IAV,RPAR1,KUN,KNUL,ISTAT)
 
127
      IDEG(1) = RPAR1(1)
 
128
      IDEG(2) = RPAR1(2)
 
129
 
 
130
      DEGMAX = 0
 
131
      IF (IDEG(2).GT.0)   DEGMAX = IDEG(1)
 
132
 
 
133
      IF (DEGMAX.GT.8) THEN
 
134
         CALL STETER(9,'REGRES/POLY: degree along X is limited to 8')
 
135
      ENDIF
 
136
 
 
137
      IF (IDEG(2).GT.8) THEN
 
138
         CALL STETER(9,'REGRES/POLY: degree along Y is limited to 8')
 
139
      ENDIF
 
140
C
 
141
C ... READ INPUT TABLE
 
142
C
 
143
      TID = -1
 
144
      CALL TBTOPN(TABLE,F_I_MODE,TID,ISTAT)
 
145
      CALL TBIGET(TID,NCOL,NROW,NSC1,NAC,NAR,ISTAT)
 
146
      CALL TDRSEL(TID,SELE,ISTAT)
 
147
C
 
148
C ... DECODE PARAMETERS
 
149
C
 
150
      COLUMN(1) = TPARBF(2)
 
151
      LINE   = ' '
 
152
      LINE   = TPARBF(3)
 
153
      NN     = INDEX(LINE,' ') - 1
 
154
      NVAR   = 1
 
155
      I1     = 1
 
156
      SDUM   = LINE//','
 
157
   20 CONTINUE
 
158
      I2     = INDEX(SDUM,',')
 
159
      SDUM(I2:I2) = ';'
 
160
      NVAR   = NVAR + 1
 
161
      IF (NVAR.GT.9) THEN
 
162
          ISTAT  = ERRPAR
 
163
          GO TO 70
 
164
      END IF
 
165
      COLUMN(NVAR) = LINE(I1:I2-1)
 
166
      I1     = I2 + 1
 
167
      IF (I1.LT.NN) GO TO 20
 
168
      K      = IDEG(1)
 
169
      L      = IDEG(2)
 
170
      IF (K.EQ.0 .AND. L.EQ.0) THEN
 
171
          IPOL   = .FALSE.
 
172
          N      = NVAR
 
173
          NM1    = N - 1
 
174
      ELSE
 
175
          IPOL   = .TRUE.
 
176
          N      = (K+1)* (L+1)
 
177
      END IF
 
178
      IF (N.GT.NOELM) THEN
 
179
          ISTAT  = ERRPAR
 
180
          GO TO 70
 
181
      END IF
 
182
C
 
183
C ... SEARCH FOR COLUMNS AND CHECK FORMAT
 
184
C
 
185
      DO 30 I = 1,NVAR
 
186
          CALL TBCSER(TID,COLUMN(I),IVAR(I),ISTAT)
 
187
          IF (ISTAT.NE.0) GO TO 70
 
188
          IF (IVAR(I).EQ.-1) THEN
 
189
              ISTAT  = ERRPAR
 
190
              GO TO 70
 
191
          END IF
 
192
          CALL TBFGET(TID,IVAR(I),FORM,LL,TYPE,ISTAT)
 
193
          IF (TYPE.EQ.D_C_FORMAT) THEN
 
194
              ISTAT  = ERRFMT
 
195
              GO TO 70
 
196
          END IF
 
197
   30 CONTINUE
 
198
      CALL STTPUT(TABLE,ISTAT)
 
199
C
 
200
C ... SOLVE LEAST SQUARE PROBLEM
 
201
C
 
202
      XMIN   = 1.D30
 
203
      YMIN   = 1.D30
 
204
      XMAX   = -XMIN
 
205
      YMAX   = -YMIN
 
206
      LL     = 0
 
207
      N1     = 0
 
208
      DO 50 IROW = 1,NROW
 
209
          CALL TBSGET(TID,IROW,SEL,ISTAT)
 
210
          IF (SEL) THEN
 
211
              CALL TBRRDD(TID,IROW,NVAR,IVAR,VAL,NULL,ISTAT)
 
212
              NEXT   = ( .NOT. NULL(1)) .AND. ( .NOT. NULL(2)) .AND.
 
213
     +                 ( .NOT. NULL(3)) .AND. ( .NOT. NULL(4)) .AND.
 
214
     +                 ( .NOT. NULL(5)) .AND. ( .NOT. NULL(6)) .AND.
 
215
     +                 ( .NOT. NULL(7)) .AND. ( .NOT. NULL(8)) .AND.
 
216
     +                 ( .NOT. NULL(9))
 
217
              IF (NEXT) THEN
 
218
                  N1     = N1 + 1
 
219
                  IF (IPOL) THEN
 
220
                      XMIN   = MIN(XMIN,VAL(2))
 
221
                      XMAX   = MAX(XMAX,VAL(2))
 
222
                      YMIN   = MIN(YMIN,VAL(3))
 
223
                      YMAX   = MAX(YMAX,VAL(3))
 
224
                      CALL TDSET2(LL+1,VAL(2),VAL(3),VAL(1),K,L,G,X,N,
 
225
     +                            MAXPAR)
 
226
                  ELSE
 
227
                      CALL TDSET1(LL+1,VAL(2),VAL(1),NM1,G,X,N,MAXPAR)
 
228
                  END IF
 
229
                  IF (LL.NE.0) THEN
 
230
                      KMIN   = MIN(LL,N+1)
 
231
                      DO 40 K1 = 1,KMIN
 
232
                          CALL TDHHTR(K1,LL+1,G,X,N,MAXPAR)
 
233
   40                 CONTINUE
 
234
                  END IF
 
235
                  LL     = MIN(LL+1,N+1)
 
236
              END IF
 
237
          END IF
 
238
   50 CONTINUE
 
239
      CALL TDSOLV(G,X,N,MAXPAR)
 
240
      ERROR  = ABS(G(N+1,N+1))
 
241
C
 
242
C ... COPY RESULT OF REGRESSION
 
243
C
 
244
      CPAR(9:16) = TABLE
 
245
      IPAR(1) = N1
 
246
      IPAR(2) = NVAR - 1
 
247
      IPAR(3) = IVAR(1)
 
248
      IF (IPOL) THEN
 
249
          CPAR(17:20) = 'MULT'
 
250
          IPAR(4) = IVAR(2)
 
251
          IPAR(5) = IVAR(3)
 
252
          IPAR(6) = IDEG(1)
 
253
          IPAR(7) = IDEG(2)
 
254
      ELSE
 
255
          CPAR(17:20) = 'LINE'
 
256
          DO 60 I = 1,NVAR - 1
 
257
              IPAR(I+3) = IVAR(I+1)
 
258
   60     CONTINUE
 
259
      END IF
 
260
      RPAR(1) = XMIN
 
261
      RPAR(2) = XMAX
 
262
      RPAR(3) = YMIN
 
263
      RPAR(4) = YMAX
 
264
      RPAR(5) = SQRT((ERROR*ERROR)/IPAR(1))
 
265
 
 
266
      INPNAM = KEYNAM(1:INDKEY)//'C'
 
267
      CALL STKWRC(INPNAM,1,CPAR,1,20,KUN,ISTAT)
 
268
      INPNAM = KEYNAM(1:INDKEY)//'I'
 
269
      CALL STKWRI(INPNAM,IPAR,1,13,KUN,ISTAT)
 
270
      INPNAM = KEYNAM(1:INDKEY)//'R'
 
271
      CALL STKWRR(INPNAM,RPAR,1,5,KUN,ISTAT)
 
272
      INPNAM = KEYNAM(1:INDKEY)//'D'
 
273
      CALL STKWRD(INPNAM,X,1,N,KUN,ISTAT)
 
274
C
 
275
C ... PRINT RESULT
 
276
C
 
277
      IF (SELE(1:1).NE.'-') THEN
 
278
          LINE   = ' SELECT  '//SELE
 
279
          CALL STTPUT(LINE,ISTAT)
 
280
      END IF
 
281
      RPAR(5) = ERROR
 
282
      CALL TDRDIS(CPAR,IPAR,RPAR,X,ISTAT)
 
283
      IF (ISTAT.NE.0) GO TO 70
 
284
C
 
285
C ... END
 
286
C
 
287
      CALL DSCUPT(TID,TID,' ',ISTAT)
 
288
      CALL TBTCLO(TID,ISTAT)
 
289
   70 IF (ISTAT.NE.0) THEN
 
290
          WRITE (MSG(13:16),9000) ISTAT
 
291
          CALL TDERRR(ISTAT,MSG,STATUS)
 
292
      END IF
 
293
      RETURN
 
294
 9000 FORMAT (I4)
 
295
      END