SUBROUTINE lvmrun(nvpar, parval, parnum, parind, ftstat, qercal, & ntrial, crtdel, qask, qansw, statnm, nbins, & tparval, tparnum, qzrnmp, qpgpar, alpha, array, & vmat, beta, wvec, savpar, ainv, ier) INTEGER nvpar DOUBLE PRECISION parval(6,nvpar), tparval(6,nvpar) DOUBLE PRECISION alpha(nvpar,nvpar), array(nvpar,nvpar) DOUBLE PRECISION vmat(nvpar,nvpar), beta(nvpar), wvec(nvpar) DOUBLE PRECISION savpar(nvpar), ainv(nvpar,nvpar) INTEGER parnum(nvpar), parind(nvpar), tparnum(nvpar) LOGICAL qzrnmp(nvpar), qpgpar(nvpar) DOUBLE PRECISION crtdel REAL ftstat INTEGER ntrial, nbins, ier CHARACTER statnm*11 LOGICAL qercal, qask, qansw c lvmrun rashafer 29 Oct 1983 c XSPEC subroutine that initiates the call to the c Marquadt process minimization routine crfitx c nvpar I i: number of parameters in fit c parval D i/r: parameter values (value, low hard limit, c low soft limit, high soft limit, high c hard limit, sigma) c parnum I i: fit parameter numbers c parind I i: parameter index : 0=normalization, -1=not c related to one normalization, N=related to c normalization fit parameter N. c ftstat R r: fit statistic c qercal L4 i: indicates that this is called by the c error calculation routine, therefore the c chattyness is reduced and NO sigmas are c calculated. c ntrial i i: number of trials before prompting c crtdel d i: critical delta fit statistic c qask l i: if true prompt user c qansw l i: if !qask specify whether to continue (T) c statnm c*11 i: name of statistic in use c nbins i4 i: total number of PHA bins - required to test c for an overdetermined fit c tparval D w: workspace array c tparnum I w: workspace array c qzrnmp L w: workspace array c qpgpar L w: workspace array c alpha D w: workspace array c array D w: workspace array c vmat D w: workspace array c beta D w: workspace array c wvec D w: workspace array c savpar D w: workspace array c ainv r D w: workspace array c ier i4 r: CRFITX error flag - Plus if the answer c to the Continue Fitting question is EOF c IER = -1 DOUBLE PRECISION dlamda REAL sttold, curdel INTEGER ivpar, jvpar, ivpart INTEGER kfpar, i, itrial, nvpart INTEGER ier2, ifppeg, iunfrz, logf, iunpeg, lenact LOGICAL qend, qdone, qunpeg CHARACTER wrtstr*255 DOUBLE PRECISION adjval EXTERNAL adjval dlamda = 1.0d-3 sttold = ftstat WRITE(wrtstr, *) 'Number of trials and critical delta : ', & ntrial, crtdel CALL xwrite(wrtstr, 25) DO ivpar = 1, nvpar WRITE(wrtstr,*) ivpar, parnum(ivpar), parind(ivpar), & parval(1,ivpar) CALL xwrite(wrtstr, 25) ENDDO qdone = .FALSE. iunpeg = 0 iunfrz = 0 qunpeg = .FALSE. DO ivpar = 1, nvpar qpgpar(ivpar) = .FALSE. ENDDO c output header kfpar = min(nvpar, 4) WRITE (wrtstr, & '(1x,a11,a20,i2,4x,4(i2,10x))') statnm, & ' Lvl Fit param #', (parnum(i), i=1, kfpar) CALL xwrite(wrtstr, 10) DO WHILE (kfpar.LT.nvpar) WRITE (wrtstr, '(6x,5i12)') (parnum(i), i=kfpar+1, & min(kfpar+5,nvpar)) CALL xwrite(wrtstr, 10) kfpar = kfpar + 5 ENDDO c set initial number of trials and delta itrial = ntrial curdel = 2.*SNGL(crtdel) qdone = .FALSE. c do loop until fit conditions satisfies, or error breakout DO WHILE (.NOT.qdone) itrial = itrial - 1 c freeze any zero-norm parameters and perform the next iteration c reset number of variable parameters for any unpegged or thawed c on previous iteration DO ivpar = 1, nvpar qzrnmp(ivpar) = .FALSE. IF ( parind(ivpar) .GT. 0 ) THEN DO jvpar = 1, nvpar IF ( parnum(jvpar) .EQ. parind(ivpar) ) THEN IF ( parval(1,jvpar) .EQ. 0. ) THEN qzrnmp(ivpar) = .TRUE. WRITE(wrtstr, '(a,i4,a)') & ' Due to zero model norms, fit parameter ', parnum(ivpar), & ' is temporarily frozen' CALL xwrite(wrtstr, 11) ENDIF ENDIF ENDDO ENDIF ENDDO c set up temporary parval and parnum arrays for variables that are not c frozen due to zero norms or pegged due to hitting limits. nvpart = 0 DO ivpar = 1, nvpar IF ( .NOT.qzrnmp(ivpar) .AND. & .NOT.qpgpar(ivpar) ) THEN nvpart = nvpart + 1 DO i = 1, 6 tparval(i,nvpart) = parval(i,ivpar) ENDDO tparnum(nvpart) = parnum(ivpar) ENDIF ENDDO CALL crfitx(nvpart, tparval, tparnum, nbins, & alpha, beta, array, dlamda, & wvec, vmat, savpar, ainv, ftstat, ier, ier2) c set new parameter values in output array ivpart = 0 DO ivpar = 1, nvpar IF ( .NOT.qzrnmp(ivpar) .AND. & .NOT.qpgpar(ivpar) ) THEN ivpart = ivpart + 1 parval(1,ivpar) = tparval(1,ivpart) ENDIF ENDDO c check for any 'zero norm' freezings current and thaw DO ivpar = 1, nvpar IF (qzrnmp(ivpar)) THEN iunfrz = iunfrz + 1 qzrnmp(ivpar) = .FALSE. ENDIF ENDDO c error checking section. On errors 1-3, 5, 7-8 then if there are pegged c parameters unpeg them and try again otherwise exit unless error 1 and c there are no variable parameters left in which case we exit anyhow. On c errors 4 and 6 peg offending parameter and try again. IF (ier.NE.0) CALL wfterr(ier, ier2, nbins) IF (ier .EQ. 4 .OR. ier .EQ. 6) THEN DO ivpar = 1, nvpar IF ( ier2 .EQ. parnum(ivpar) ) THEN qpgpar(ivpar) = .TRUE. ENDIF ENDDO qdone = .FALSE. GOTO 200 ELSEIF (ier .GT. 0) THEN qdone = .TRUE. IF ( ier .EQ. 1 .AND. ier2 .EQ. 0 ) THEN GOTO 300 ELSE DO ivpar = 1, nvpar IF ( qpgpar(ivpar) ) qdone = .FALSE. ENDDO GOTO 100 ENDIF ENDIF c having disposed of all possible error conditions to the end of the loop c continue by check for any pegged parameters IF (.NOT.qunpeg) THEN CALL chkpeg(nvpar, parval, parnum, parind, qpgpar, ifppeg) ENDIF c write out current status in fit logf = log10(dlamda) + 0.5 kfpar = min(nvpar, 4) WRITE (wrtstr, '(1pg14.6,i3,3x,5g12.4)') ftstat, logf, & (adjval(parval(1,i)), i=1, kfpar) CALL xwrite(wrtstr, 10) DO WHILE (kfpar.LT.nvpar) WRITE (wrtstr, '(12x,1p5g12.4)') (adjval(parval(1,i)), & i=kfpar+1, min(kfpar+5,nvpar)) CALL xwrite(wrtstr, 10) kfpar = kfpar + 5 ENDDO curdel = sttold - ftstat sttold = ftstat c check for end of fitting condition IF ( curdel .GE. 0. .AND. curdel.LE.SNGL(crtdel) ) THEN qdone = .TRUE. ELSE qunpeg = .FALSE. GOTO 200 ENDIF c breakout point for error conditions except 4 and 6. unpeg any pegged c variables and reset qdone if necessary (ie fit will continue if c variables are unpegged). If parameters are unpegged then reset dlamda c so effectively the fit is restarted with the current parameters as c the initial guess. 100 CONTINUE iunpeg = 0 DO ivpar = 1, nvpar IF (qpgpar(ivpar)) THEN qpgpar(ivpar) = .FALSE. iunpeg = iunpeg + 1 ENDIF ENDDO IF (qunpeg .AND. iunpeg.NE.0) THEN GOTO 200 ENDIF qunpeg = (iunpeg.NE.0) IF (iunpeg.GT.0) THEN WRITE (wrtstr, '(a,i3,a)') ' Unpegged', iunpeg, & ' parameters' CALL xwrite(wrtstr, 15) qdone = .FALSE. dlamda = 1.0d-3 ENDIF c breakout point for errors 4 and 6 and if curdel is not less than critical c value 200 CONTINUE c if we have exceeded the number of trials then answer the user c whether to carry on. Note the check for the querying flag. If c the user asks not to continue then ensure that the error flag c returned is zero unless they hit EOF in which case it is -1. IF ( itrial .LE. 0 ) THEN IF ( qask ) THEN WRITE (wrtstr, '(a,1pg12.4)') & ' Number of trials exceeded - last iteration delta = ' & , curdel CALL xwrtpr(wrtstr(:lenact(wrtstr))) CALL xquest(' Continue fitting?', 'Y', qdone, qend) ELSE qend = .false. IF ( qansw ) THEN qdone = .true. ELSE qdone = .false. ENDIF ENDIF IF (.NOT.qdone ) ier = 0 IF (qend) ier = -1 qdone = (.NOT.qdone) .OR. qend itrial = ntrial ENDIF c if ier=-1 then absolute override to end loop. IF (ier.EQ.-1) qdone = .TRUE. 300 CONTINUE ENDDO IF (.NOT.qercal) THEN c Calculate the parameter approximate sigmas (under the assumption that c fit statistic space is locally quadratic). IF ( nvpart .GT. 0 ) THEN CALL ivermt(nvpart, tparval, tparnum, array, alpha, wvec, & vmat, ainv) c and write them into the output array DO ivpar = 1, nvpar parval(6,ivpar) = -1.0d0 ENDDO DO ivpar = 1, nvpar DO jvpar = 1, nvpart IF ( parnum(ivpar) .EQ. tparnum(jvpar) ) THEN parval(6,ivpar) = tparval(6,jvpar) ENDIF ENDDO ENDDO ENDIF ENDIF RETURN END