← Back to branch summary

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

« back to all changes in this revision

Viewing changes to .pc/writting.patch/applic/fit/libsrc/fitcon.for

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-sl34juxohmn4aty4
Tags: 13.09pl1.2+dfsg-1
http://bugs.debian.org/740702
Initial release. (Closes: #740702)

Show diffs side-by-side

added added

removed removed

Lines of Context:
.pc/writting.patch/applic/fit/libsrc/fitcon.for
 
1
C  @(#)fitcon.for       19.1 (ESO-DMD) 02/25/03 13:17:46 
 
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 Massachusetts Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Correspondence 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
      SUBROUTINE FITCON(ALGOR,BNDTYP,ISTAT)                                    
 
30
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++        
 
31
C                                                                               
 
32
C.MODULE                                                                        
 
33
C       FIT                                                                     
 
34
C                                                                               
 
35
C.NAME                                                                          
 
36
C       FIT_CON                                                                 
 
37
C                                                                               
 
38
C.PURPOSE                                                                       
 
39
C   makes the interface between the MIDAS data structures and the NAG           
 
40
C   routines for non linear LSQ problems with constraints                       
 
41
C                                                                               
 
42
C.KEYWORDS                                                                      
 
43
C  Constrained Non-linear Least Square, Constrained Non-linear Fitting,         
 
44
C  Constrained Non-linear Regression, FIT/TABLE, FIT/IMAGE,                     
 
45
C  Corrected Gauss-Newton Method, Quasi-Newton Method, Modified                 
 
46
C  Gauss-Newton.                                                                
 
47
C                                                                               
 
48
C.DESCRIPTION                                                                   
 
49
C       Recovers qualifiers and parameters from MIDAS keywords;                 
 
50
C       Recovers bounds;                                                        
 
51
C       Give virtual memory for calculations in NAG routines;                   
 
52
C       Call the NAG routines for non-linear least squares problems;            
 
53
C       Display information messages;                                           
 
54
C       Compute errors on parameters.                                           
 
55
C ****
 
56
C CHECK CALLING SEQUENCE OF E04KAF,E04KBF,E04KCF,E04KDF
 
57
C     LAST TWO ARGUMENTS...
 
58
C ****
 
59
C                                                                               
 
60
C.RESTRICTIONS                                                                  
 
61
C       ...                                                                     
 
62
C                                                                               
 
63
C.LANGUAGE                                                                      
 
64
C       FORTRAN                                                                 
 
65
C                                                                               
 
66
C.CALLING SEQUENCE                                                              
 
67
C       CALL FITCON(ALGOR,BNDTYP,ISTAT)                                        
 
68
C                                                                               
 
69
C.INPUT PARAMETERS                                                              
 
70
C       ALGOR    (*)        CHARACTER   algorithm to be used.                   
 
71
C                           (CGNND,QN,MGN)                                      
 
72
C       BNDTYP   (*)        CHARACTER   Bounds types                            
 
73
C                           (P,G,I)                                             
 
74
C       ISTAT               INTEGER     exit status                             
 
75
C                                                                               
 
76
C.MODIFIED PARAMETERS                                                           
 
77
C       None                                                                    
 
78
C                                                                               
 
79
C.OUTPUT PARAMETERS                                                             
 
80
C       None                                                                    
 
81
C                                                                               
 
82
C.FILES                                                                         
 
83
C       FIT_NAG.INC/NOLIST                                                      
 
84
C                                                                               
 
85
C.MODULES CALLED                                                                
 
86
C       NAG library     (E04___)                                                
 
87
C       MINMSQ          (message display)                                       
 
88
C                                                                               
 
89
C.AUTHOR                                                                        
 
90
C       Ph. DEFERT,      Feb 1986                                               
 
91
C        M. Peron        Jul 1992                                             
 
92
C.MODIFICATIONS                                                                 
 
93
C                                                                               
 
94
C                                                                               
 
95
C.BIBLIOGRAPHY                                                                  
 
96
C       NAG Mark II Vol 3 (E04)                                                 
 
97
C       Algorithms for the Solution of non-linear least squares problem         
 
98
C          SIAM J. on Num. An.,15,977-992,1978                                  
 
99
C.VERSION
 
100
C 021122        last modif
 
101
C                                                                               
 
102
C-----------------------------------------------------------------------        
 
103
C 
 
104
C      IMPLICIT NONE                                                            
 
105
C     .. Scalar Arguments ..                                                    
 
106
      CHARACTER ALGOR* (*),BNDTYP* (*)                                          
 
107
      INTEGER ISTAT ,KUN,KNUL,TID,TID1
 
108
C     ..                                                                        
 
109
C     .. Scalars in Common ..                                                   
 
110
      INTEGER NRCOL,ISTAR                                                       
 
111
      CHARACTER WGTTYP*1                                                        
 
112
C     ..                                                                        
 
113
C     .. Arrays in Common ..                                                    
 
114
      INTEGER ICOL(10)                 
 
115
      INTEGER MADRID(1)                 
 
116
C     ..                                                                        
 
117
C     .. Local Scalars ..                                                       
 
118
      INTEGER I,IDAT,IFAIL1,IFAIL3,J,K,LIW,LW,NC,NROW,NRFEVA,           
 
119
     +        NRVM,NS,STATUS
 
120
      INTEGER IBOUND,LHSL,IMIN,IMAX,NRBND,ROUTIN                     
 
121
      INTEGER NAC,NAR,STAT,NAC1,NAR1,NBYT
 
122
      INTEGER*8 PTRW,PTRIW,PTRDLT,PTRHSL,PTRHSD,PTRSTA,PTRGRD
 
123
 
 
124
      DOUBLE PRECISION FSUMSQ
 
125
      DOUBLE PRECISION HMIN,HMAX
 
126
 
 
127
      LOGICAL VALID,ISEL,NOMPAR,VCMATPR                                         
 
128
      CHARACTER WGTKEY*10,LINE*78,BNDTAB*60                           
 
129
 
 
130
C     ..                                                                        
 
131
C     .. Local Arrays ..                                                        
 
132
      REAL XVAL(10)                                                             
 
133
      LOGICAL NULL(10)                                                          
 
134
C     ..                                                                        
 
135
C     ..                                                                        
 
136
C     .. External Files ..                                                      
 
137
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
138
      INCLUDE 'MID_INCLUDE:FITNAGI.INC'                                 
 
139
      COMMON/VMR/MADRID
 
140
      COMMON /LSQFUN/ICOL,NRCOL,ISTAR,WGTTYP 
 
141
      INCLUDE 'MID_INCLUDE:FITNAGC.INC'                                 
 
142
      EXTERNAL E04LBS,E04JBQ,FUNCT3,FUNCT4,MINMON
 
143
 
 
144
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
145
      DATA WGTKEY(1:10)/'FITCHAR   '/                                           
 
146
 
 
147
C                                                                               
 
148
C  Executable Statements                                                        
 
149
C                                                                               
 
150
C  Initialize parameters and compute number of fixed parameters                 
 
151
C                                                                               
 
152
      VCMATPR = METPAR(1) .LT. 0.                                               
 
153
      METPAR(1) = ABS(METPAR(1))                                                
 
154
      WRITE(6,*)METPAR(1)
 
155
      NRPFIX = 0                                                                
 
156
      DO 10 J = 1,NRPARAM                                                       
 
157
          IF (METPAR(2).GE.0) PARAM(J) = PARINI(J)                              
 
158
          IF (FIXPAR(J).GE.0) NRPFIX = NRPFIX + 1                               
 
159
   10 CONTINUE                                                                  
 
160
      METPAR(2) = ABS(METPAR(2))                                                
 
161
C                                                                               
 
162
C  Get the type of weighting from FITCHAR 13:13                                 
 
163
C                                                                               
 
164
      CALL STKRDC(WGTKEY,1,FZIWGT,1,K,WGTTYP,
 
165
     +    KUN,KNUL,ISTAT)
 
166
      CALL UPCAS(WGTTYP,WGTTYP)
 
167
      IF (WGTTYP(1:1).NE.'C' .AND. WGTTYP(1:1).NE.'W' .AND.                     
 
168
     +    WGTTYP(1:1).NE.'I' .AND. WGTTYP(1:1).NE.'S') CALL STETER(12,          
 
169
     +    'FIT/TABLE: Weight type mismatch')                                    
 
170
C                                                                               
 
171
      IF (FILTYP(1:1).EQ.'T') THEN                                              
 
172
C                                                                               
 
173
C  Get the selected columns for indep and dep variables                         
 
174
C                                                                               
 
175
          IF (WGTCOL.EQ.0) THEN                                                 
 
176
              ISTAR  = 2                                                        
 
177
              NRCOL  = NRIND + 1                                                
 
178
              ICOL(1) = 0                                                       
 
179
                                                                                
 
180
          ELSE                                                                  
 
181
              ISTAR  = 1                                                        
 
182
              NRCOL  = NRIND + 2                                                
 
183
              ICOL(1) = WGTCOL                                                  
 
184
          END IF                                                                
 
185
                                                                                
 
186
          ICOL(2) = DEPCOL                                                      
 
187
          NULL(1) = .FALSE.                                                     
 
188
          DO 20 I = 1,NRIND                                                     
 
189
              ICOL(I+2) = INDCOL(I)                                             
 
190
   20     CONTINUE                                                              
 
191
C                                                                               
 
192
C  Get number of selected rows in the table (nr of residuals)                   
 
193
C                                                                               
 
194
          NRDATA = 0                                                            
 
195
          CALL TBIGET(FZIDEN,NC,NROW,NS,NAC,NAR,STAT)                                 
 
196
          DO 40 IDAT = 1,NROW                                                   
 
197
              CALL TBSGET(FZIDEN,IDAT,ISEL,ISTAT)                              
 
198
              IF ( .NOT. ISEL) GO TO 40                                         
 
199
              CALL TBRRDR(FZIDEN,IDAT,NRCOL,ICOL(ISTAR),                  
 
200
     +                    XVAL(ISTAR),NULL,ISTAT)                              
 
201
              VALID  = ISEL                                                     
 
202
              DO 30 K = 1,NRCOL                                                 
 
203
                  VALID  = VALID .AND. ( .NOT. NULL(K))                         
 
204
   30         CONTINUE                                                          
 
205
              IF (VALID) NRDATA = NRDATA + 1                                    
 
206
   40     CONTINUE                                                              
 
207
                                                                                
 
208
      ELSE                                                                      
 
209
C                                                                               
 
210
C  Get the number of residuals for FIT/IMAG                                     
 
211
C                                                                               
 
212
          NRDATA = NRPIX(1)*NRPIX(2)*NRPIX(3)                                   
 
213
          NROW   = NRDATA                                                       
 
214
      END IF                                                                    
 
215
C                                                                               
 
216
C     Get qualifiers of the FIT command : PRINT,MAX. NR. Of F. EVAL,            
 
217
C     PRECISION,ETA and STEP. NOMPAR detects if method parameters               
 
218
C     have been given or are all to be defaulted.                               
 
219
C                                                                               
 
220
      NOMPAR = METPAR(1) .LT. 0.1                                               
 
221
      IF (ABS(ABS(METPAR(3))-999999.).LT.1.E-6) THEN                            
 
222
          METPAR(3) = 0.0                                                       
 
223
                                                                                
 
224
      ELSE                                                                      
 
225
          NOMPAR = .FALSE.                                                      
 
226
      END IF                                                                    
 
227
                                                                                
 
228
      IF (ABS(ABS(METPAR(2))-999999.).LT.1.E-6) THEN                            
 
229
          IF (ALGOR(1:5).EQ.'CGNND') THEN                                       
 
230
              METPAR(2) = 40*NRPARAM* (NRPARAM+5)                               
 
231
                                                                                
 
232
          ELSE IF (ALGOR(1:2).EQ.'QN') THEN                                     
 
233
              METPAR(2) = 100*NRPARAM                                           
 
234
                                                                                
 
235
          ELSE IF (ALGOR(1:3).EQ.'MGN') THEN                                    
 
236
              METPAR(2) = 50*NRPARAM                                            
 
237
          END IF                                                                
 
238
                                                                                
 
239
      ELSE                                                                      
 
240
          NOMPAR = .FALSE.                                                      
 
241
      END IF                                                                    
 
242
                                                                                
 
243
      IF (ABS(ABS(METPAR(4))-999999.).LT.1.E-6) THEN                            
 
244
          IF (NRPARAM.EQ.1) THEN                                                
 
245
              METPAR(4) = 0.                                                    
 
246
                                                                                
 
247
          ELSE                                                                  
 
248
              IF (ALGOR(1:5).EQ.'CGNND') THEN                                   
 
249
                  METPAR(4) = 0.5                                               
 
250
                                                                                
 
251
              ELSE IF (ALGOR(1:2).EQ.'QN') THEN                                 
 
252
                  METPAR(4) = 0.9                                               
 
253
                                                                                
 
254
              ELSE IF (ALGOR(1:3).EQ.'MGN') THEN                                
 
255
                  IF (NRPARAM.LT.10) THEN                                       
 
256
                      METPAR(4) = 0.5                                           
 
257
                                                                                
 
258
                  ELSE IF (NRPARAM.LE.20) THEN                                  
 
259
                      METPAR(4) = 0.1                                           
 
260
                                                                                
 
261
                  ELSE                                                          
 
262
                      METPAR(4) = 0.01                                          
 
263
                  END IF                                                        
 
264
                                                                                
 
265
              END IF                                                            
 
266
                                                                                
 
267
          END IF                                                                
 
268
                                                                                
 
269
      ELSE                                                                      
 
270
          NOMPAR = .FALSE.                                                      
 
271
      END IF                                                                    
 
272
                                                                                
 
273
      IF (ABS(ABS(METPAR(5))-999999.).LT.1.E-6) THEN                            
 
274
          METPAR(5) = 1.E05                                                     
 
275
                                                                                
 
276
      ELSE                                                                      
 
277
          NOMPAR = .FALSE.                                                      
 
278
      END IF                                                                    
 
279
C                                                                               
 
280
C   Get the bounds and store them in PARLOW and PARUPP                          
 
281
C                                                                               
 
282
      IF (BNDTYP(1:1).EQ.'G') THEN                                              
 
283
          IBOUND = 3                                                            
 
284
          CALL STKRDC('BNDTAB',1,1,60,K,BNDTAB,KUN,KNUL,ISTAT)                  
 
285
          IF (BNDTAB(1:1).EQ.'+') THEN
 
286
            CALL STTPUT('Bounds table not defined',ISTAT)
 
287
            CALL STSEPI
 
288
          ENDIF 
 
289
          CALL TBTOPN(BNDTAB,F_IO_MODE,TID,ISTAT)                             
 
290
          CALL TBIGET(TID,NC,NRBND,NS,NAC1,NAR1,ISTAT)                           
 
291
          IF (NC.LT.2) THEN                                                     
 
292
              CALL STETER(12,'FIT/... : Error in bounds table')                 
 
293
                                                                                
 
294
          ELSE IF (NC.GT.2) THEN                                                
 
295
              CALL STTPUT('*** WARN-FIT : Bounds table more than'//             
 
296
     +                    ' 2 columns, take first ones ***',ISTAT)
 
297
          END IF                                                                
 
298
                                                                                
 
299
          IF (NRBND.NE.1) THEN                                                  
 
300
              CALL STTPUT('*** WARN-FIT : Bounds table more than'//             
 
301
     +                    ' 1 row, global bounds taken as first one ***'        
 
302
     +                    ,ISTAT)                                              
 
303
          END IF                                                                
 
304
                                                                                
 
305
          CALL TBSGET(TID,1,ISEL,ISTAT)                                     
 
306
          ICOL(1) = 1                                                           
 
307
          ICOL(2) = 2                                                           
 
308
          CALL TBRRDR(TID,1,2,ICOL,                  
 
309
     +                    XVAL,NULL,ISTAT)                        
 
310
          VALID  = ISEL .AND. ((.NOT.NULL(1)) .OR. (.NOT.NULL(2)))              
 
311
          IF ( .NOT. VALID) CALL STETER(8,                                      
 
312
     +                           ' FIT/TABLE : Global bounds not founds'        
 
313
     +                                  )                                       
 
314
          IF ( .NOT. NULL(1)) THEN                                              
 
315
              PARLOW(1) = XVAL(1)                                               
 
316
                                                                                
 
317
          ELSE                                                                  
 
318
              PARLOW(1) = -1.D+6                                                
 
319
          END IF                                                                
 
320
                                                                                
 
321
          IF ( .NOT. NULL(2)) THEN                                              
 
322
              PARUPP(1) = XVAL(2)                                               
 
323
                                                                                
 
324
          ELSE                                                                  
 
325
              PARUPP(1) = 1.D+6                                                 
 
326
          END IF                                                                
 
327
                                                                                
 
328
          PARLOW(1) = XVAL(1)                                                   
 
329
          PARUPP(1) = XVAL(2)                                                   
 
330
          CALL STTPUT('   Lower Bound     Upper Bound',ISTAT)
 
331
          WRITE (LINE,9000) PARLOW(1),PARUPP(1)                                 
 
332
          CALL STTPUT(LINE,ISTAT)                                              
 
333
          CALL STTPUT(' ',ISTAT)                                               
 
334
                                                                                
 
335
      ELSE IF (BNDTYP(1:1).EQ.'I') THEN                                         
 
336
          IBOUND = 0                                                            
 
337
          CALL STKRDC('BNDTAB',1,1,60,K,BNDTAB,KUN,KNUL,ISTAT)   
 
338
          IF (BNDTAB(1:1).EQ.'+') THEN
 
339
            CALL STTPUT('Bounds table not defined',ISTAT)
 
340
            CALL STSEPI
 
341
          ENDIF 
 
342
          CALL TBTOPN(BNDTAB,F_IO_MODE,TID1,ISTAT)                          
 
343
          CALL TBIGET(TID1,NC,NRBND,NS,NAC1,NAR1,ISTAT)
 
344
          IF (NC.LT.2) THEN                                                     
 
345
              CALL STETER(12,'FIT/... : Error in bounds table')                 
 
346
                                                                                
 
347
          ELSE IF (NC.GT.2) THEN                                                
 
348
              CALL STTPUT('*** WARN-FIT : Bounds table more than'//             
 
349
     +                    ' 2 columns, take first ones ***',STATUS)             
 
350
          END IF                                                                
 
351
                                                                                
 
352
          IF (NRBND.GT.NRPARAM) THEN                                            
 
353
              CALL STTPUT('*** WARN-FIT : Bounds table has more '//             
 
354
     +                    'rows than parameters, take first ones ***',          
 
355
     +                    STATUS)                                               
 
356
          END IF                                                                
 
357
                                                                                
 
358
          DO 50 IDAT = 1,MIN0(NRBND,NRPARAM)                                    
 
359
              CALL TBSGET(TID1,IDAT,ISEL,STATUS)                              
 
360
              IF (ISEL) THEN                                                    
 
361
                  ICOL(1) = 1                                                   
 
362
                  ICOL(2) = 2                                                   
 
363
                  CALL TBRRDR(TID1,IDAT,2,ICOL,                  
 
364
     +                    XVAL,NULL,ISTAT)                        
 
365
               IF ( .NOT. NULL(1)) THEN                                      
 
366
                      PARLOW(IDAT) = XVAL(1)                                    
 
367
                                                                                
 
368
                  ELSE                                                          
 
369
                      PARLOW(IDAT) = -1.D+6                                     
 
370
                  END IF                                                        
 
371
                                                                                
 
372
                  IF ( .NOT. NULL(2)) THEN                                      
 
373
                      PARUPP(IDAT) = XVAL(2)                                    
 
374
                                                                                
 
375
                  ELSE                                                          
 
376
                      PARUPP(IDAT) = 1.D+6                                      
 
377
                  END IF                                                        
 
378
                                                                                
 
379
              ELSE                                                              
 
380
                  PARLOW(IDAT) = -1.D+6                                         
 
381
                  PARUPP(IDAT) = 1.D+6                                          
 
382
              END IF                                                            
 
383
                                                                                
 
384
   50     CONTINUE                                                              
 
385
          IF (NRBND.LT.NRPARAM) THEN                                            
 
386
              DO 60 IDAT = NRBND + 1,NRPARAM                                    
 
387
                  PARLOW(IDAT) = -1.D+6                                         
 
388
                  PARUPP(IDAT) = 1.D+6                                          
 
389
   60         CONTINUE                                                          
 
390
          END IF                                                                
 
391
                                                                                
 
392
          CALL STTPUT(                                                          
 
393
     +             '   Initial values     Lower Bounds     Upper Bounds'        
 
394
     +                ,STATUS)                                                  
 
395
          DO 70 K = 1,NRPARAM                                                   
 
396
              WRITE (LINE,9010) PARAM(K),PARLOW(K),PARUPP(K)                    
 
397
              CALL STTPUT(LINE,STATUS)                                          
 
398
   70     CONTINUE                                                              
 
399
          CALL STTPUT(' ',STATUS)                                               
 
400
                                                                                
 
401
      ELSE IF (BNDTYP(1:1).EQ.'P') THEN                                         
 
402
          CALL STTPUT(' Parameters must be positive',STATUS)                    
 
403
          IBOUND = 2                                                            
 
404
      END IF                                                                    
 
405
C                                                                               
 
406
C  Call to the NAG routines after having displayed the parameters and           
 
407
C  get the necessary V.M. space.                                                
 
408
C                                                                               
 
409
      IFAIL1 = 1                                                                
 
410
      XVAL(1) = 1.                                                              
 
411
      IF (ALGOR(1:5).EQ.'CGNND') THEN                                           
 
412
          CALL STTPUT(' ',STATUS)                                               
 
413
          CALL STTPUT(                                                          
 
414
     +          ' Constrained Non-linear Least Squares Fitting : Corr. '        
 
415
     +                //'Gauss-Newton M. (no der)',STATUS)                      
 
416
          CALL STTPUT(                                                          
 
417
     +          ' _____________________________________________________'        
 
418
     +                //'________________________',STATUS)                      
 
419
          CALL STTPUT(' ',STATUS)                                               
 
420
          CALL STTPUT(' ',STATUS)                                               
 
421
          CALL STTPUT(                                                          
 
422
     +  ' Print   Max. Fct. eval.  Prec. on Param.  Lin. Search  Domain'        
 
423
     +                //'  Weight',STATUS)                                      
 
424
          WRITE (LINE,9020) NINT(METPAR(1)),NINT(METPAR(2)),METPAR(3),          
 
425
     +      METPAR(4),METPAR(5),WGTTYP                                          
 
426
          CALL STTPUT(LINE,STATUS)                                              
 
427
          CALL STTPUT(' ',STATUS)                                               
 
428
                                                                                
 
429
          IF (NOMPAR) THEN                                                      
 
430
              ROUTIN = 1                                                        
 
431
              IF (NRPARAM.EQ.1) THEN                                            
 
432
                  LW     = 13                                                   
 
433
                                                                                
 
434
              ELSE                                                              
 
435
                  LW     = 12*NRPARAM + NRPARAM* (NRPARAM-1)/2                  
 
436
              END IF                                                            
 
437
                                                                                
 
438
              LIW    = NRPARAM + 2                                              
 
439
                                                                                
 
440
              NRVM   = 2*LW + LIW                                               
 
441
              NBYT = 4*NRVM
 
442
              CALL TDMGET(NBYT,PTRW,STATUS)
 
443
              PTRIW  = PTRW + 2*LW                                              
 
444
                                                                                
 
445
              CALL E04JAF(NRPARAM,IBOUND,PARLOW,PARUPP,PARAM,FSUMSQ,            
 
446
     +                    MADRID(PTRIW),LIW,MADRID(PTRW),LW,IFAIL1)             
 
447
                                                                                
 
448
          ELSE                                                                  
 
449
              ROUTIN = 2                                                        
 
450
              LW     = 9*NRPARAM                                                
 
451
              LHSL   = MAX0(NRPARAM* (NRPARAM-1)/2,1)                           
 
452
              LIW    = 10                                                       
 
453
                                                                                
 
454
              NRVM   = 2* (LW+LHSL+4*NRPARAM) + LIW                             
 
455
 
 
456
C             CALL STIPUT('WORK ',10,9,1,1,NRVM,1.,XVAL,' ',' ',PTRW,            
 
457
C    +                    STATUS)                                               
 
458
              NBYT = 4*NRVM
 
459
              CALL TDMGET(NBYT,PTRW,STATUS)
 
460
              PTRDLT = PTRW + 2*LW                                              
 
461
              PTRHSL = PTRDLT + 2*NRPARAM                                       
 
462
              PTRHSD = PTRHSL + 2*LHSL                                          
 
463
              PTRSTA = PTRHSD + 2*NRPARAM                                       
 
464
              PTRGRD = PTRSTA + NRPARAM                                       
 
465
              PTRIW  = PTRGRD + 2*NRPARAM                                       
 
466
                                                                                
 
467
              IFAIL3 = 1                                                        
 
468
              CALL E04HBF(NRPARAM,FUNCT3,PARAM,NRFEVA,MADRID(PTRDLT),          
 
469
     +                    MADRID(PTRHSL),LHSL,MADRID(PTRHSD),FSUMSQ,            
 
470
     +                    MADRID(PTRGRD),MADRID(PTRIW),LIW,MADRID(PTRW),        
 
471
     +                    LW,IFAIL3)                                            
 
472
                                                                                
 
473
                                                                                
 
474
              CALL MNMX(MADRID(PTRHSD),NRPARAM,HMIN,HMAX,IMIN,IMAX)             
 
475
              IF (HMAX.LT.1.D+4*HMIN) THEN                                      
 
476
                  WRITE (LINE,9030) HMIN,HMAX                                   
 
477
                  CALL STETER(21,LINE)                                          
 
478
              END IF                                                            
 
479
                                                                                
 
480
              CALL E04JBF(NRPARAM,FUNCT3,MINMON,NINT(METPAR(1)),.TRUE.,         
 
481
     +                    1,E04JBQ,NINT(METPAR(2)),DBLE(METPAR(4)),             
 
482
     +                    DBLE(METPAR(3)),DBLE(METPAR(5)),0.D0,                 
 
483
     +                    MADRID(PTRDLT),IBOUND,PARLOW,PARUPP,PARAM,            
 
484
     +                    MADRID(PTRHSL),LHSL,MADRID(PTRHSD),                   
 
485
     +                    MADRID(PTRSTA),FSUMSQ,MADRID(PTRGRD),                 
 
486
     +                    MADRID(PTRIW),LIW,MADRID(PTRW),LW,IFAIL1)             
 
487
                                                                                
 
488
          END IF                                                                
 
489
                                                                                
 
490
      ELSE IF (ALGOR(1:2).EQ.'QN') THEN                                         
 
491
          CALL STTPUT(' ',STATUS)                                               
 
492
          CALL STTPUT(                                                          
 
493
     +   ' Constrained Non-linear Least Squares Fitting : Quasi-Newton '        
 
494
     +                //'Method',STATUS)                                        
 
495
          CALL STTPUT(                                                          
 
496
     +   ' ____________________________________________________________'        
 
497
     +                //'______',STATUS)                                        
 
498
          CALL STTPUT(' ',STATUS)                                               
 
499
          CALL STTPUT(' ',STATUS)                                               
 
500
          CALL STTPUT(                                                          
 
501
     +  ' Print   Max. Fct. eval.  Prec. on Param.  Lin. Search  Domain'        
 
502
     +                //'  Weight',STATUS)                                      
 
503
          WRITE (LINE,9020) NINT(METPAR(1)),NINT(METPAR(2)),METPAR(3),          
 
504
     +      METPAR(4),METPAR(5),WGTTYP                                          
 
505
          CALL STTPUT(LINE,STATUS)                                              
 
506
          CALL STTPUT(' ',STATUS)                                               
 
507
                                                                                
 
508
          IF (NOMPAR) THEN                                                      
 
509
              ROUTIN = 3                                                        
 
510
              IF (NRPARAM.EQ.1) THEN                                            
 
511
                  LW     = 11                                                   
 
512
                                                                                
 
513
              ELSE                                                              
 
514
                  LW     = 10*NRPARAM + NRPARAM* (NRPARAM-1)/2                  
 
515
              END IF                                                            
 
516
                                                                                
 
517
              LIW    = NRPARAM                                              
 
518
              NRVM   = LW*2 + LIW+ 2*NRPARAM                                               
 
519
C              CALL STIPUT('WORK ',10,9,1,1,NRVM,1.,XVAL,' ',' ',PTRW,            
 
520
C     +                    STATUS)                                               
 
521
              NBYT=NRVM*4
 
522
              CALL TDMGET(NBYT,PTRW,STATUS)
 
523
              PTRIW  = PTRW + 2*LW                                              
 
524
              PTRGRD = PTRIW + 2*NRPARAM
 
525
              CALL E04KAF(NRPARAM,IBOUND,PARLOW,PARUPP,PARAM,FSUMSQ,            
 
526
     +                    MADRID(PTRGRD),MADRID(PTRIW),LIW,MADRID(PTRW),        
 
527
     +                    LW,IFAIL1)                                            
 
528
                                                                                
 
529
          ELSE                                                                  
 
530
              ROUTIN = 4                                                        
 
531
              LW     = 9*NRPARAM                                                
 
532
              LIW    = 10                                                       
 
533
              LHSL   = MAX(NRPARAM* (NRPARAM-1)/2,1)                            
 
534
                                                                                
 
535
              NRVM   = (LW+3*NRPARAM+LHSL+1)*2 + LIW + NRPARAM                  
 
536
              NBYT = NRVM*4
 
537
              CALL TDMGET(NBYT,PTRW,STATUS)
 
538
C              CALL STIPUT('WORK ',10,9,1,1,NRVM,1.,XVAL,' ',' ',PTRW,            
 
539
C     +                    STATUS)                                               
 
540
              PTRHSL = PTRW + 3*NRPARAM                                         
 
541
              PTRHSD = PTRHSL + 2*LHSL                                          
 
542
              PTRSTA = PTRHSD + 2*NRPARAM                                       
 
543
              PTRGRD = PTRSTA + NRPARAM                                       
 
544
              PTRIW  = PTRGRD + 2*NRPARAM                                       
 
545
                                                                                
 
546
                                                                                
 
547
              IFAIL3 = 1                                                        
 
548
              CALL E04HCF(NRPARAM,FUNCT4,PARINI,FSUMSQ,MADRID(PTRGRD),          
 
549
     +                    MADRID(PTRIW),LIW,MADRID(PTRW),LW,IFAIL3)             
 
550
                                                                                
 
551
              IF (IFAIL3.EQ.2) THEN                                             
 
552
                  CALL MINMSG(ALGOR(1:5),ROUTIN,10)                             
 
553
                  GO TO 100                                                     
 
554
                                                                                
 
555
              END IF                                                            
 
556
                                                                                
 
557
              CALL SETZER(MADRID(PTRHSL),LHSL)                                 
 
558
              CALL SETONE(MADRID(PTRHSD),NRPARAM)                               
 
559
                                                                                
 
560
              CALL E04KBF(NRPARAM,FUNCT4,MINMON,NINT(METPAR(1)),.TRUE.,         
 
561
     +                    1,E04LBS,NINT(METPAR(2)),DBLE(METPAR(4)),             
 
562
     +                    DBLE(METPAR(3)),DBLE(METPAR(5)),0.D0,IBOUND,          
 
563
     +                    PARLOW,PARUPP,PARAM,MADRID(PTRHSL),LHSL,              
 
564
     +                    MADRID(PTRHSD),MADRID(PTRSTA),FSUMSQ,                 
 
565
     +                    MADRID(PTRGRD),MADRID(PTRIW),LIW,MADRID(PTRW),        
 
566
     +                    LW,IFAIL1)                                            
 
567
                                                                                
 
568
          END IF                                                                
 
569
                                                                                
 
570
      ELSE IF (ALGOR(1:3).EQ.'MGN') THEN                                        
 
571
          CALL STTPUT(' ',STATUS)                                               
 
572
          CALL STTPUT(                                                          
 
573
     +       ' Constrained Non-linear Least Squares Fitting : Modified '        
 
574
     +                //'Gauss-Newton Method',STATUS)                           
 
575
          CALL STTPUT(                                                          
 
576
     +       ' ________________________________________________________'        
 
577
     +                //'___________________',STATUS)                           
 
578
          CALL STTPUT(' ',STATUS)                                               
 
579
          CALL STTPUT(' ',STATUS)                                               
 
580
          CALL STTPUT(                                                          
 
581
     +  ' Print   Max. Fct. eval.  Prec. on Param.  Lin. Search  Domain'        
 
582
     +                //'  Weight',STATUS)                                      
 
583
          WRITE (LINE,9020) NINT(METPAR(1)),NINT(METPAR(2)),METPAR(3),          
 
584
     +      METPAR(4),METPAR(5),WGTTYP                                          
 
585
          CALL STTPUT(LINE,STATUS)                                              
 
586
          CALL STTPUT(' ',STATUS)                                               
 
587
                                                                                
 
588
          IF (NOMPAR) THEN                                                      
 
589
              ROUTIN = 5                                                        
 
590
              IF (NRPARAM.EQ.1) THEN                                            
 
591
                  LW     = 11                                                   
 
592
                                                                                
 
593
              ELSE                                                              
 
594
                  LW     = 10*NRPARAM + NRPARAM* (NRPARAM-1)/2                  
 
595
              END IF                                                            
 
596
                                                                                
 
597
              LIW    = NRPARAM + 4                                              
 
598
              NRVM   = LW*2 + LIW                                               
 
599
C             CALL STIPUT('WORK ',10,9,1,1,NRVM,1.,XVAL,' ',' ',PTRW,            
 
600
C     +                    STATUS)                  
 
601
              NBYT = NRVM*4
 
602
              CALL TDMGET(NBYT,PTRW,STATUS)                             
 
603
              PTRIW  = PTRW + 2*LW                                              
 
604
              PTRGRD = PTRIW + 2 + NRPARAM
 
605
              CALL E04KCF(NRPARAM,IBOUND,PARLOW,PARUPP,PARAM,FSUMSQ,            
 
606
     +                    MADRID(PTRGRD),MADRID(PTRIW),LIW,MADRID(PTRW),        
 
607
     +                    LW,IFAIL1)                                            
 
608
                                                                                
 
609
          ELSE                                                                  
 
610
              ROUTIN = 6                                                        
 
611
              IF (NRPARAM.EQ.1) THEN                                            
 
612
                  LW     = 8                                                    
 
613
                                                                                
 
614
              ELSE                                                              
 
615
                  LW     = 7*NRPARAM + NRPARAM* (NRPARAM-1)/2                   
 
616
              END IF                                                            
 
617
                                                                                
 
618
              LIW    = 10                                                       
 
619
              LHSL   = MAX(NRPARAM* (NRPARAM-1)/2,1)                            
 
620
                                                                                
 
621
              NRVM   = (LW+3*NRPARAM+LHSL+1)*2 + LIW + NRPARAM                  
 
622
              NBYT = NRVM*4
 
623
              CALL TDMGET(NBYT,PTRW,STATUS)
 
624
C              CALL STIPUT('WORK ',10,9,1,1,NRVM,1.,XVAL,' ',' ',PTRW,            
 
625
C     +                    STATUS)                                               
 
626
              PTRHSL = PTRW + 2*LW                                              
 
627
              PTRHSD = PTRHSL + 2*LHSL                                          
 
628
              PTRSTA = PTRHSD + 2*NRPARAM                                       
 
629
              PTRGRD = PTRSTA + NRPARAM                                       
 
630
              PTRIW  = PTRGRD + 2*NRPARAM                                       
 
631
                                                                                
 
632
                                                                                
 
633
              IFAIL3 = 1                                                        
 
634
              CALL E04HCF(NRPARAM,FUNCT4,PARINI,FSUMSQ,MADRID(PTRGRD),          
 
635
     +                    MADRID(PTRIW),LIW,MADRID(PTRW),LW,IFAIL3)             
 
636
                                                                                
 
637
              IF (IFAIL3.EQ.2) THEN                                             
 
638
                  CALL MINMSG(ALGOR(1:5),ROUTIN,10)                             
 
639
                  GO TO 100                                                     
 
640
                                                                                
 
641
              END IF                                                            
 
642
                                                                                
 
643
              CALL E04KDF(NRPARAM,FUNCT4,MINMON,NINT(METPAR(1)),                
 
644
     +                    NINT(METPAR(2)),DBLE(METPAR(4)),                      
 
645
     +                    DBLE(METPAR(3)),0.D0,DBLE(METPAR(5)),IBOUND,          
 
646
     +                    PARLOW,PARUPP,PARAM,MADRID(PTRHSL),LHSL,              
 
647
     +                    MADRID(PTRHSD),MADRID(PTRSTA),FSUMSQ,                 
 
648
     +                    MADRID(PTRGRD),MADRID(PTRIW),LIW,MADRID(PTRW),        
 
649
     +                    LW,IFAIL1)                                            
 
650
                                                                                
 
651
          END IF                                                                
 
652
                                                                                
 
653
      END IF                                                                    
 
654
                                                                                
 
655
      CALL MINMSG(ALGOR(1:5),ROUTIN,IFAIL1)                                     
 
656
      CALL STTPUT(' ',STATUS)                                                   
 
657
                                                                                
 
658
      IF (IFAIL1.NE.1 .AND. IFAIL1.NE.4 .AND. IFAIL1.NE.9) THEN                 
 
659
          IF (NOMPAR) THEN                                                      
 
660
              IF (ALGOR(1:5).EQ.'CGNND') THEN                                   
 
661
                  PTRHSD = PTRW + 8* (6*NRPARAM+2*NRDATA+NRDATA*NRPARAM+        
 
662
     +                     MAX0(1, (NRPARAM* (NRPARAM-1)/2)))                   
 
663
                  PTRHSD = PTRHSD + 8*NRPARAM                                   
 
664
                                                                                
 
665
              ELSE                                                              
 
666
                  PTRHSD = PTRW + 8* (7*NRPARAM+2*NRDATA+NRDATA*NRPARAM+        
 
667
     +                     NRPARAM* (NRPARAM+1)/2+                              
 
668
     +                     MAX0(1, (NRPARAM* (NRPARAM-1)/2)))                   
 
669
                  PTRHSD = PTRHSD + 8*NRPARAM                                   
 
670
                                                                                
 
671
              END IF                                                            
 
672
                                                                                
 
673
          END IF                                                                
 
674
                                                                                
 
675
                                                                                
 
676
      ELSE                                                                      
 
677
          CALL STTPUT(                                                          
 
678
     +        ' *** Unable to evaluate the errors on the parameters ***'        
 
679
     +                ,STATUS)                                                  
 
680
      END IF                                                                    
 
681
                                                                                
 
682
  100 CONTINUE                                                                  
 
683
C                                                                               
 
684
C  Pass the final reduced chi squares to the COMMON block                       
 
685
C                                                                               
 
686
      NRITER = NRFEVA                                                          
 
687
      FINCHI = FSUMSQ/DBLE(NRDATA-NRPARAM+NRPFIX)                             
 
688
C                                                                               
 
689
C  Exit                                                                         
 
690
C                                                                               
 
691
      RETURN                                                                    
 
692
                                                                                
 
693
 9000 FORMAT (5X,G12.5,10X,G12.5)                                               
 
694
 9010 FORMAT (5X,G12.5,10X,G12.5,10X,G12.5)                                     
 
695
 9020 FORMAT (I6,8X,I7,7X,1PE9.1,8X,0PF4.2,6X,0PF8.0,5X,A1)                     
 
696
 9030 FORMAT ('FIT/TABLE : Bad Scaling of the problem ','HSD from ',            
 
697
     +       1PD9.2,' to ',1PD9.2)                                              
 
698
 9040 FORMAT ('*** Detected only ',I3,' linear independant parameters ',        
 
699
     +       'on the total of ',I3,' ***')                                      
 
700
      END