C @(#)fitunc.for 17.1.1.1 (ESO-DMD) 01/25/02 17:10:43 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 SUBROUTINE FITUNC(ALGOR1,ISTAT) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.MODULE C FIT C C.NAME C FIT_UNC C C.PURPOSE C makes the interface between the MIDAS data structures and the NAG C routines for non linear LSQ problems without constraints C C.KEYWORDS C Non-linear Least Square, Non-linear Fitting, Non-linear Regression C FIT/TABLE,FIT/IMAGE, Newton Method, Corrected Gauss-Newton Method, C Quasi-Newton Method, Modified Gauss-Newton. C C.DESCRIPTION C Recovers qualifiers and parameters from MIDAS keywords; C Give virtual memory for calculations in NAG routines; C Call the NAG routines for non-linear least squares problems; C Display information messages; C Compute errors on parameters. C C.RESTRICTIONS C ... C C.LANGUAGE C FORTRAN C C.CALLING SEQUENCE C C C.INPUT PARAMETERS C ALGOR (*) CHARACTER algorithm to be used. C (CGNND,QN,MGN) C ISTAT INTEGER exit status C C.MODIFIED PARAMETERS C None C C.OUTPUT PARAMETERS C None C C.FILES C FIT_NAG.INC/NOLIST C C.MODULES CALLED C NAG library (E04...) C LSQMSQ (message display) C Standard Interface (....._ST) C C.AUTHOR C Ph. DEFERT, Feb 1986 C C.MODIFICATIONS C C C.BIBLIOGRAPHY C NAG Mark II Vol 3 (E04) C Algorithms for the Solution of non-linear least squares problem C SIAM J. on Num. An.,15,977-992,1978 C C-unsupported IMPLICIT NONE C IMPLICIT NONE C .. C .. Scalar Arguments .. CHARACTER*(*) ALGOR1 INTEGER ISTAT C .. C .. Scalars in Common .. INTEGER NRCOL,ISTAR,NAC,NAR,NBYT,IFAIM2 CHARACTER WGTTYP*1 CHARACTER ALGOR*10 C .. C .. Arrays in Common .. INTEGER MADRID(1),ICOL(10) C .. C .. Local Scalars .. INTEGER I,IDAT,IFAIL1,IFAIL2,J,K,LIW,LW,NC,NROW,NRVM,NS,STATUS, + PTRW,PTRIW,PTRSV,PTRSVD,PTRFJA,PTRFVC,NRFEVAL,JOB INTEGER KUN,KNUL DOUBLE PRECISION FSUMSQ C DOUBLE PRECISION XTOL,ETA,STEPMX LOGICAL VALID,ISEL,NOMPAR,VCMATPR CHARACTER WGTKEY*10,LINE*78 C .. C .. Local Arrays .. REAL XVAL(10) LOGICAL NULL(10) C INCLUDE 'MID_INCLUDE:FITNAGI.INC/NOLIST' C .. C .. Common blocks .. COMMON /LSQFUN/ICOL,NRCOL,ISTAR,WGTTYP COMMON/VMR/MADRID C .. C .. C .. External Files .. INCLUDE 'MID_INCLUDE:FITNAGC.INC/NOLIST' EXTERNAL LSFUN3,LSQMON,E04HEV,LSFUN4 DATA WGTKEY(1:10)/'FITCHAR '/ C C Executable Statements C C C Initialize parameters and compute number of fixed parameters C ALGOR = ALGOR1 CALL FORUPC(ALGOR,ALGOR) VCMATPR = METPAR(1) .LT. 0. METPAR(1) = ABS(METPAR(1)) NRPFIX = 0 DO 10 J = 1,NRPARAM IF (METPAR(2).GE.0) PARAM(J) = PARINI(J) IF (FIXPAR(J).GE.0) NRPFIX = NRPFIX + 1 10 CONTINUE METPAR(2) = ABS(METPAR(2)) C C Get the type of weighting from FITCHAR 13:13 C CALL STKRDC(WGTKEY,1,FZIWGT,1,K,WGTTYP,KUN,KNUL,STATUS) CALL UPCAS(WGTTYP,WGTTYP) IF (WGTTYP(1:1).NE.'C' .AND. WGTTYP(1:1).NE.'W' .AND. . WGTTYP(1:1).NE.'I' .AND. WGTTYP(1:1).NE.'S') . CALL STETER(12,'FIT/TABLE: Weight type mismatch') C IF (FILTYP(1:1).EQ.'T') THEN C C Get columns of indep and dep variables, of the weight if needed C IF (WGTCOL.EQ.0) THEN ISTAR = 2 NRCOL = NRIND + 1 ICOL(1) = 0 ELSE ISTAR = 1 NRCOL = NRIND + 2 ICOL(1) = WGTCOL END IF ICOL(2) = DEPCOL NULL(1) = .FALSE. DO 20 I = 1,NRIND ICOL(I+2) = INDCOL(I) 20 CONTINUE C C Get number of selected rows in the table (nr of residuals) C NRDATA = 0 CALL TBIGET(FZIDEN,NC,NROW,NS,NAC,NAR,STATUS) DO 40 IDAT = 1,NROW CALL TBSGET(FZIDEN,IDAT,ISEL,STATUS) IF ( .NOT. ISEL) GO TO 40 CALL TBRRDR(FZIDEN,IDAT,NRCOL,ICOL(ISTAR), + XVAL(ISTAR),NULL,STATUS) VALID = ISEL DO 30 K = 1,NRCOL VALID = VALID .AND. ( .NOT. NULL(K)) 30 CONTINUE IF (VALID) NRDATA = NRDATA + 1 40 CONTINUE ELSE C C Get the number of residuals for FIT/IMAG C NRDATA = NRPIX(1)*NRPIX(2)*NRPIX(3) NROW = NRDATA END IF C C Get qualifiers of the FIT command : PRINT,MAX. NR. Of F. EVAL, C PRECISION,ETA and STEP. NOMPAR detects if method parameters C have been given or are all to be defaulted. C NOMPAR = METPAR(1) .LT. 0.1 IF (ABS(ABS(METPAR(3))-999999.0).LT.1.E-6) THEN METPAR(3) = 1.E-5 ELSE NOMPAR = .FALSE. END IF IF (ABS(ABS(METPAR(2))-999999.0).LT.1.E-6) THEN METPAR(2) = 400*NRPARAM ELSE NOMPAR = .FALSE. END IF IF (ABS(ABS(METPAR(4))-999999.0).LT.1.E-6) THEN IF (NRPARAM.EQ.1) THEN METPAR(4) = 0. ELSE IF (ALGOR(1:5).EQ.'CGNND') THEN METPAR(4) = 0.5 ELSE IF (ALGOR(1:3).EQ.'MGN' .OR. ALGOR(1:2).EQ.'QN') THEN METPAR(4) = 0.9 END IF END IF ELSE NOMPAR = .FALSE. END IF IF (ABS(ABS(METPAR(5))-999999.0).LT.1.E-6) THEN METPAR(5) = 1.E05 ELSE NOMPAR = .FALSE. END IF C C Call to the NAG routines after having displayed the parameters and C get the necessary V.M. space. C IFAIL1 = 1 XVAL(1) = 1. IF (ALGOR(1:5).EQ.'CGNND') THEN CALL STTPUT(' ',STATUS) CALL STTPUT( + ' Non-linear Least Squares Fitting : Corrected Gauss-Newton ' + //'Method (no der)',STATUS) CALL STTPUT( + ' __________________________________________________________' + //'_______________',STATUS) CALL STTPUT(' ',STATUS) CALL STTPUT(' ',STATUS) CALL STTPUT( + ' Print Max. Fct. eval. Prec. on Param. Lin. Search Domain' + //' Weight',STATUS) WRITE (LINE,9000) NINT(METPAR(1)),NINT(METPAR(2)),METPAR(3), + METPAR(4),METPAR(5),WGTTYP CALL STTPUT(LINE,STATUS) CALL STTPUT(' ',STATUS) IF (NOMPAR) THEN IF (NRPARAM.EQ.1) THEN LW = 9 + 5*NROW ELSE LW = (7+NRPARAM+2*NROW+NRPARAM/2)*NRPARAM + + 3*NROW + 1 END IF LIW = 100 NRVM = 2*LW + LIW NBYT = NRVM*4 CALL TDMGET(NBYT,PTRW,STATUS) PTRIW = PTRW + 2*LW CALL E04FDF(NRDATA,NRPARAM,PARAM,FSUMSQ,MADRID(PTRIW), + LIW,MADRID(PTRW),LW,IFAIL1) ELSE IF (NRPARAM.EQ.1) THEN LW = 7 + 3*NRPARAM ELSE LW = (6+NROW+NRPARAM/2)*NRPARAM + 2*NROW + 1 END IF LIW = 100 NRVM = 2* (LW+ (NROW+NRPARAM)* (NRPARAM+1)) + LIW NBYT = NRVM*4 CALL TDMGET(NBYT,PTRW,STATUS) PTRFVC = PTRW + 2*LW PTRFJA = PTRFVC + 2*NROW PTRSV = PTRFJA + 2*NROW*NRPARAM PTRSVD = PTRSV + 2*NRPARAM PTRIW = PTRSVD + 2*NRPARAM*NRPARAM CALL E04FCF(NRDATA,NRPARAM,LSFUN3,LSQMON,NINT(METPAR(1)), + NINT(METPAR(2)),DBLE(METPAR(4)), + DBLE(METPAR(3)),DBLE(METPAR(5)),PARAM,FSUMSQ, + MADRID(PTRFVC),MADRID(PTRFJA),NRDATA, + MADRID(PTRSV),MADRID(PTRSVD),NRPARAM,NRITER, + NRFEVAL,MADRID(PTRIW),LIW,MADRID(PTRW),LW, + IFAIL1) END IF ELSE IF (ALGOR(1:2).EQ.'QN') THEN CALL STTPUT(' ',STATUS) CALL STTPUT( + ' Non-linear Least Squares Fitting : Quasi-Newton Method' + ,STATUS) CALL STTPUT( + ' ______________________________________________________' + ,STATUS) CALL STTPUT(' ',STATUS) CALL STTPUT(' ',STATUS) CALL STTPUT( + ' Print Max. Fct. eval. Prec. on Param. Lin. Search Domain' + //' Weight',STATUS) WRITE (LINE,9000) NINT(METPAR(1)),NINT(METPAR(2)),METPAR(3), + METPAR(4),METPAR(5),WGTTYP CALL STTPUT(LINE,STATUS) CALL STTPUT(' ',STATUS) IF (NOMPAR) THEN IF (NRPARAM.EQ.1) THEN LW = 11 + 5*NROW + 1 ELSE LW = (8+2*NRPARAM+2*NROW)*NRPARAM + 3*NROW + 1 END IF LIW = 100 NRVM = LW*2 + LIW NBYT = NRVM*4 CALL TDMGET(NBYT,PTRW,STATUS) PTRIW = PTRW + 2*LW CALL E04GCF(NRDATA,NRPARAM,PARAM,FSUMSQ,MADRID(PTRIW), + LIW,MADRID(PTRW),LW,IFAIL1) ELSE IF (NRPARAM.EQ.1) THEN LW = 9 + 3*NRPARAM ELSE LW = (7+NROW+NRPARAM)*NRPARAM + 2*NROW + 1 END IF LIW = 100 NRVM = (LW+ (NROW+NRPARAM)* (NRPARAM+1))*2 + LIW NBYT = NRVM*4 CALL TDMGET(NBYT,PTRW,STATUS) PTRFVC = PTRW + 2*LW PTRFJA = PTRFVC + 2*NROW PTRSV = PTRFJA + 2*NROW*NRPARAM PTRSVD = PTRSV + 2*NRPARAM PTRIW = PTRSVD + 2*NRPARAM*NRPARAM CALL E04GBF(NRDATA,NRPARAM,E04HEV,LSFUN4,LSQMON, + NINT(METPAR(1)),NINT(METPAR(2)), + DBLE(METPAR(3)),DBLE(METPAR(4)), + DBLE(METPAR(5)),PARAM,FSUMSQ,MADRID(PTRFVC), + MADRID(PTRFJA),NRDATA,MADRID(PTRSV), + MADRID(PTRSVD),NRPARAM,NRITER,NRFEVAL, + MADRID(PTRIW),LIW,MADRID(PTRW),LW,IFAIL1) END IF ELSE IF (ALGOR(1:3).EQ.'MGN') THEN CALL STTPUT(' ',STATUS) CALL STTPUT( +' Non-linear Least Squares Fitting : Modified Gauss-Newton Method' + ,STATUS) CALL STTPUT( +' _______________________________________________________________' + ,STATUS) CALL STTPUT(' ',STATUS) CALL STTPUT(' ',STATUS) CALL STTPUT( + ' Print Max. Fct. eval. Prec. on Param. Lin. Search Domain' + //' Weight',STATUS) WRITE (LINE,9000) NINT(METPAR(1)),NINT(METPAR(2)),METPAR(3), + METPAR(4),METPAR(5),WGTTYP CALL STTPUT(LINE,STATUS) CALL STTPUT(' ',STATUS) IF (NOMPAR) THEN IF (NRPARAM.EQ.1) THEN LW = 11 + 5*NROW + 1 ELSE LW = (8+2*NRPARAM+2*NROW)*NRPARAM + 3*NROW + 1 END IF LIW = 100 NRVM = LW*2 + LIW NBYT = NRVM*4 CALL TDMGET(NBYT,PTRW,STATUS) PTRIW = PTRW + 2*LW CALL E04GEF(NRDATA,NRPARAM,PARAM,FSUMSQ,MADRID(PTRIW), + LIW,MADRID(PTRW),LW,IFAIL1) ELSE IF (NRPARAM.EQ.1) THEN LW = 9 + 3*NRPARAM ELSE LW = (7+NROW+NRPARAM)*NRPARAM + 2*NROW + 1 END IF LIW = 100 NRVM = (LW+ (NROW+NRPARAM)* (NRPARAM+1))*2 + LIW NBYT = NRVM*4 CALL TDMGET(NBYT,PTRW,STATUS) PTRFVC = PTRW + 2*LW PTRFJA = PTRFVC + 2*NROW PTRSV = PTRFJA + 2*NROW*NRPARAM PTRSVD = PTRSV + 2*NRPARAM PTRIW = PTRSVD + 2*NRPARAM*NRPARAM CALL E04GDF(NRDATA,NRPARAM,LSFUN4,LSQMON,NINT(METPAR(1)), + NINT(METPAR(2)),DBLE(METPAR(3)), + DBLE(METPAR(4)),DBLE(METPAR(5)),PARAM,FSUMSQ, + MADRID(PTRFVC),MADRID(PTRFJA),NRDATA, + MADRID(PTRSV),MADRID(PTRSVD),NRPARAM,NRITER, + NRFEVAL,MADRID(PTRIW),LIW,MADRID(PTRW),LW, + IFAIL1) END IF END IF C C Display message after NAG call C CALL LSQMSG(ALGOR(1:5),IFAIL1) CALL STTPUT(' ',STATUS) C C Call to the NAG routine to compute the errors on parameters after C having computed the address of the relevant information to pass. C IF (IFAIL1.NE.1 .AND. IFAIL1.NE.4 .AND. IFAIL1.NE.9) THEN IF (NOMPAR) THEN IF (ALGOR(1:5).EQ.'CGNND') THEN PTRSV = PTRW + 2* (6*NRPARAM+2*NRDATA+NRDATA*NRPARAM+ + MAX0(1, (NRPARAM* (NRPARAM-1)/2))) PTRSVD = PTRSV + 2*NRPARAM ELSE PTRSV = PTRW + 2* (7*NRPARAM+2*NRDATA+NRDATA*NRPARAM+ + NRPARAM* (NRPARAM+1)/2+ + MAX0(1, (NRPARAM* (NRPARAM-1)/2))) PTRSVD = PTRSV + 2*NRPARAM END IF END IF IFAIL2 = 1 IF (VCMATPR) THEN JOB = -1 ELSE JOB = 0 END IF CALL E04YCF(JOB,NRDATA,NRPARAM,FSUMSQ,MADRID(PTRSV), + MADRID(PTRSVD),NRPARAM,ERROR,MADRID(PTRW),IFAIL2) IF (IFAIL2.NE.0) THEN CALL STTPUT(' ',STATUS) IFAIM2 = IFAIL2-2 WRITE (LINE,9010) IFAIM2,NRPARAM CALL STTPUT(LINE,STATUS) CALL STTPUT(' ',STATUS) END IF IF (VCMATPR) THEN CALL DCPRIN(MADRID(PTRSVD),NRPARAM,NRPARAM) CALL MATDIA(MADRID(PTRSVD),ERROR,NRPARAM,NRPARAM) END IF DO 50 K = 1,NRPARAM ERROR(K) = SQRT(ERROR(K)) 50 CONTINUE ELSE CALL STTPUT( + ' *** Unable to evaluate the errors on the parameters ***' + ,STATUS) END IF C C Pass the final reduced chi squares to the COMMON block C NRITER = NRFEVAL FINCHI = FSUMSQ/DBLE(NRDATA-NRPARAM+NRPFIX) C C Exit C CALL TDMFRE(NBYT,PTRW,STATUS) RETURN 9000 FORMAT (I6,8X,I7,7X,1PE9.1,8X,0PF4.2,6X,0PF8.0,5X,A1) 9010 FORMAT ('*** Detected only ',I3,' linear independant parameters ', + 'on the total of ',I3,' ***') END