C @(#)fttntr.for 17.1.1.1 (ES0-DMD) 01/25/02 17:10:47 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 FTTNTR(NITER,RELAX,PREC,V1,V2,NACT,ACTCH,ISTAT) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 17:25 - 13 JAN 1988 C C.LANGUAGE: F77+ESOext C C.AUTHOR: J.D.PONZ C C.IDENTIFICATION C C FITLIB.FOR VERSION 1.0 27 MAR 1984 C C.PURPOSE C C INTERFACE ROUTINES FOR THE FITTING STRUCTURES C ITERATE ON THE FITTING PROCESS C THE ITERATION IS FINISHED IF C - NITER IS REACHED OR C - CHISQ IS REACHED C - CONVERGENCE IS ACHIEVED C C C.ALGORITHM C C USE MIDAS I/O INTERFACES TO FRAMES AND TABLES C C.KEYWORDS C C NON LINEAR FITTING C C C---------------------------------------------------------------- C C INPUT PARAMETERS C NITER INTG NO. OF ITERATIONS TO BE DONE C <0 USE PREVIOUS COMPUTED VALUES C >0 USE GUESS C RELAX REAL RELAXATION FACTOR C CHISQ REAL CHI SQUARE C V1,V2 REAL OPTIONAL FITTING INTERVAL (IF V1.NE.V2) C C OUTPUT PARAMETERS C NACT INTG ACTUAL NUMBER OF ITERATIONS C ACTCH REAL ACTUAL CHI SQ. C ISTAT INTG STATUS RETURN C INTEGER NITER,NACT,ISTAT,IPRINT,I,J,NITER1,ICOUNT,NN,NP,NNP INTEGER ISTAR,NCOL,NS,KZ6941,N,IDAT,NOFFSET,IFUN,K,IP,NILAST REAL RELAX,PREC,V1,V2,ACTCH,VS,VE,CHIOLD,CHIPERC,AC1,AC2,AC3 REAL DELTA,C1,C2 CHARACTER*80 LINE CHARACTER*6 TYPE DOUBLE PRECISION CHI,PP,Q,Y,Y1,Y2,YOUT,CORRNRM DOUBLE PRECISION AL,B(50),A(50,50) INTEGER ICOL(10),ACOL,AROW,NC,NROW REAL XVAL(10),X(8),W,FSUMSQ LOGICAL NULL(10),VALID,NEXT,AUTO,ISEL,VCMATPR INCLUDE 'MID_INCLUDE:FITI.INC/NOLIST' DOUBLE PRECISION FZDERIV(FZPARMAX) INCLUDE 'MID_INCLUDE:FITC.INC/NOLIST' EQUIVALENCE (XVAL(1),W) C EQUIVALENCE (XVAL(2),Y) EQUIVALENCE (XVAL(3),X(1)) DATA TYPE/'R*4'/ DO 200 I=1,10 NULL(I) = .FALSE. 200 CONTINUE CALL STTPUT(' ',ISTAT) CALL STTPUT('Non-linear Least Square Fitting : '// + 'Newton-Raphson Method (defaulted)',ISTAT) CALL STTPUT('_______'// + '___________________________________________________________' + ,ISTAT) CALL STTPUT(' ',ISTAT) CALL STTPUT(' ',ISTAT) CALL STTPUT( + ' Print Max. Fct. eval. Prec. on Param. Relax. Fact.' + ,ISTAT) WRITE (LINE,9000) NINT(FZMETPAR(1)),NINT(FZMETPAR(2)), + FZMETPAR(3),FZMETPAR(4) CALL STTPUT(LINE,ISTAT) CALL STTPUT(' ',ISTAT) VCMATPR = FZMETPAR(1) .LE. 0.1 IPRINT = NINT(ABS(FZMETPAR(1))) IF (IPRINT.EQ.0) IPRINT = ABS(NITER) AUTO = RELAX .LT. 0. RELAX = ABS(RELAX) DO 10 I = 1,8 NULL(I) = .FALSE. 10 CONTINUE NEXT = .TRUE. IF (NITER.GT.0) THEN DO 20 J = 1,FZNPTOT FZVALUE(J) = FZGUESS(J) 20 CONTINUE FZNITER = 0 NITER1 = NITER ELSE NITER1 = -NITER END IF VS = MIN(V1,V2) VE = MAX(V1,V2) FZCCHIS = 9999. FZCHISQ = 1.D-6 FZRELAX = RELAX CORRNRM = 1.D6 CHI = 0.D0 ICOUNT = 0 IF (FZWEIGHT.EQ.0) THEN ISTAR = 2 W = 1. NCOL = FZNIND + 1 ICOL(1) = 0 ELSE ISTAR = 1 NCOL = FZNIND + 2 ICOL(1) = FZWEIGHT END IF ICOL(2) = FZDVAR NULL(1) = .FALSE. DO 30 I = 1,FZNIND ICOL(I+2) = FZIVAR(I) 30 CONTINUE CALL TBIGET(FZIDEN,NC,NROW,NS,ACOL,AROW,ISTAT) KZ6941 = 0 99 KZ6941 = KZ6941 + 1 IF ( .NOT. (NEXT)) GO TO 180 DO 50 I = 1,50 DO 40 J = 1,50 A(J,I) = 0.D0 40 CONTINUE B(I) = 0.D0 50 CONTINUE CHIOLD = CHI CHI = 0.D0 C C ITERATION FOR TABLES C N = 0 DO 100 IDAT = 1,NROW CALL TBSGET(FZIDEN,IDAT,ISEL,ISTAT) C XVAL(2) = Y ! to avoid problems with the equivalence CALL TBRRDR(FZIDEN,IDAT,NCOL,ICOL(ISTAR), + XVAL(ISTAR),NULL(ISTAR),ISTAT) Y=XVAL(2) VALID = ISEL .AND. ( .NOT. NULL(1)) .AND. + ( .NOT. NULL(2)) .AND. ( .NOT. NULL(3)) C C ... COMPUTE FUNCTION VALUES AND DERIVATIVES (0 IF FIXED PARAM) C IF (V1.NE.V2) VALID = VALID .AND. XVAL(1) .GE. VS .AND. + XVAL(1) .LE. VE IF (VALID) THEN N = N + 1 Y1 = 0. NOFFSET = 1 DO 60 IFUN = 1,FZNFUN CALL FTFUNC(FZFCODE(IFUN),FZNIND,X, + FZACTPAR(IFUN),FZVALUE(NOFFSET), + YOUT,FZDERIV(NOFFSET)) Y1 = Y1 + YOUT NOFFSET = NOFFSET + FZACTPAR(IFUN) 60 CONTINUE DO 70 J = 1,FZNPTOT IP = FZFIXED(J) IF (IP.EQ.0) FZDERIV(J) = 0.D0 IF (IP.GT.0) THEN FZDERIV(IP) = FZDERIV(IP) + + FZDERIV(J)/FZPFAC(J) FZDERIV(J) = 0.D0 END IF 70 CONTINUE C C ... COMPUTE RESIDUALS, WEIGHTED RESIDUALS AND CHI**2 C W = 1./ (W*W) Y2 = Y - Y1 PP = DBLE(W) IF (FZFLAG.EQ.1) PP = 1.D0/AMAX1(1.0,SNGL(Y1)) Q = Y2 CHI = CHI + Q*Q*PP C C ... BUILD NORMAL EQUATION MATRIX C DO 90 J = 1,FZNPTOT AL = FZDERIV(J) DO 80 K = 1,FZNPTOT A(J,K) = A(J,K) + AL*PP*FZDERIV(K) 80 CONTINUE B(J) = B(J) + PP*Q*AL 90 CONTINUE END IF 100 CONTINUE C C ... ADD (PROBABILISTIC) PARAMETER CONSTRAINTS C NN = N NP = 0 DO 110 I = 1,FZNPTOT IF (FZFIXED(I).GE.0) THEN A(I,I) = 1.0D0 ELSE NP = NP + 1 IF (FZUNCER(I).GT.0.) THEN NN = NN + 1 Q = 1.0D0/DBLE(FZUNCER(I)*FZUNCER(I)) A(I,I) = A(I,I) + Q B(I) = B(I) + Q* (FZGUESS(I)-FZVALUE(I)) C ^ C To be modified for soft constraints END IF END IF 110 CONTINUE C C ... INVERT NORMAL EQUATION MATRIX TO GET THE COVARIANCE MATRIX C CALL DMATIN(A,FZNPTOT,50,FZFIXED,ISTAT) NNP = FZNPTOT IF (ISTAT.NE.0) THEN CALL STTPUT('*** ERR-NR : Problems inverting matrix ***', + ISTAT) RETURN END IF C C ... COMPUTE REDUCED CHI**2 AND CORRECTIONS FOR THE UNKNOWNS C FSUMSQ = CHI CHI = DMAX1(CHI/DBLE(NN-NP),0.D0) CORRNRM = 0.D0 DO 130 I = 1,FZNPTOT Q = 0.0D0 DO 120 J = 1,FZNPTOT Q = Q + A(I,J)*B(J) 120 CONTINUE C C ... COMPUTE ERRORS OF THE CORRECTIONS (I.E. THE UNKNOWNS) C ... AND APPLY CORRECTIONS (DO THE ADJUSTMENT!) C IF (FZFIXED(I).EQ.0) THEN FZERROR(I) = 0.D0 ELSE Q = RELAX*Q CORRNRM = CORRNRM + Q*Q FZERROR(I) = DSQRT(A(I,I)*CHI) FZVALUE(I) = FZVALUE(I) + Q END IF 130 CONTINUE CORRNRM = DSQRT(CORRNRM) C C ... LINEAR CONSTR. BETWEEN TWO PARAMS C DO 140 I = 1,FZNPTOT IP = FZFIXED(I) IF (IP.GT.0) THEN FZERROR(I) = FZPFAC(I)*FZERROR(IP) FZVALUE(I) = FZPFAC(I)*FZVALUE(IP) END IF 140 CONTINUE ICOUNT = ICOUNT + 1 FZNITER = FZNITER + 1 FZCCHIS = SNGL(CHI) NEXT = ICOUNT .LT. NITER1 .AND. CORRNRM .GT. PREC IF ( .NOT. NEXT .OR. MOD(ICOUNT,IPRINT).EQ.0) THEN CALL STTPUT(' ',ISTAT) CALL STTPUT( + ' Iter Fct. Eval. Sum of Squares Red. Chisq. % Decr.' + ,ISTAT) IF (ICOUNT.GT.1) THEN CHIPERC = 100.* (CHIOLD-CHI)/CHIOLD WRITE (LINE,9010) FZNITER,FZNITER,FSUMSQ, + FZCCHIS,CHIPERC ELSE WRITE (LINE,9010) FZNITER,FZNITER,FSUMSQ, + FZCCHIS END IF CALL STTPUT(LINE,ISTAT) CALL STTPUT(' ',ISTAT) CALL STTPUT(' Parameters ',ISTAT) DO 150 I = 1,FZNPTOT WRITE (LINE,9020) FZVALUE(I) CALL STTPUT(LINE,ISTAT) 150 CONTINUE CALL STTPUT(' ',ISTAT) END IF C C ... STORE THE CHISQ C IF (FZNITER.LE.100) THEN FZCHI(FZNITER) = CHI NILAST = FZNITER ELSE DO 160 I = 1,99 FZCHI(I) = FZCHI(I+1) 160 CONTINUE FZCHI(100) = CHI NILAST = 100 END IF C C ... C IF (AUTO) THEN IF (NILAST.GT.10) THEN C C ... CHECK FOR CONVERGENCE C IF (FZCHI(NILAST-1).EQ.FZCHI(NILAST)) GO TO 190 AC1 = ALOG(AMAX1(FZCHI(NILAST-2),0.0001)) AC2 = ALOG(AMAX1(FZCHI(NILAST-1),0.0001)) AC3 = ALOG(AMAX1(FZCHI(NILAST),0.0001)) C1 = AC1 - AC2 C2 = AC2 - AC3 IF (ABS(C2).LE.1.E-5) GO TO 190 DELTA = ABS((C1-C2)/C2) IF (C2.GT.0.0 .AND. DELTA.LE.0.25) RELAX = (1.0- + DELTA)**4 END IF END IF GOTO 99 180 CONTINUE 190 CALL STTPUT(' ',ISTAT) NACT = FZNITER ACTCH = FZCCHIS FZNDAT = N IF (CORRNRM.LE.PREC) THEN CALL STTPUT(' --> NR : Convergence achieved <--',ISTAT) ELSE CALL STTPUT('*** WARN-NR : No convergence reached ***',ISTAT) END IF CALL STTPUT(' ',ISTAT) C C ... PRINT CORREL. MATRIX C IF (VCMATPR) CALL DCPRIN(A,NNP,50) RETURN 9000 FORMAT (I6,7X,I7,7X,1PE9.1,8X,0PF4.2) 9010 FORMAT (I5,6X,I5,3X,1PE12.4,2X,1PE12.4,4X,0PF5.2) 9020 FORMAT (2X,1PD15.7) END