occam一维反演

C----------------------------------------------------------------------- program runocc C----------------------------------------------------------------------- c C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Program Revision 2.01, 13 Jan 1993 c c runocc is the calling program for OCCAM v2.0. C C REFERENCES: CONSTABLE, PARKER & CONSTABLE, 1987: GEOPHYSICS 52, 289-300. C DEGROOT-HEDLIN & CONSTABLE, 1990: GEOPHYSICS 55, 1613-1624. C CONSTABLE, 1991: GEOPHYS. J. INT. 106, 387-388. C CONSTABLE, 1992: OCCAM DISTRIBUTION NOTES. C C REVISION HISTORY: C MARCH 1986, (VERSION 1.2) C OCTOBER 1987, (VERSION 1.3) C AUGUST 1988, (VERSION 1.4) C JANUARY 1989 (WITH THANKS TO C.DEGROOT-HEDLIN), (VERSION 1.5) C SEPTEMBER 1989, (VERSION 1.5.2) C APRIL 1992 (MORE THANKS TO CATHERINE AND ALSO WSE) (VERSION 2.0) C C IF YOU OBTAIN THIS CODE FROM A THIRD PERSON, PLEASE SEND YOUR NAME AND ADDRESS C TO S. CONSTABLE. YOU WILL THEN RECEIVE UPDATES, NEWS ON BUGS, ETC. C C VERSION 2.0 IS A RE-WRITE TO MAKE OCCAM INDEPENDENT OF THE DIMENSIONALITY C OF THE FORWARD PROBLEM, AND TO MAKE MAXIMUM USE OF DYNAMIC MEMORY ALLOCATION, C AS THE ARRAY SIZES FOR THE 2D PROBLEM ARE NOW GETTING TOO LARGE FOR COMFORT. C INCLUDE FILES ARE NOW USED TO STREAMLINE CHANGES OF DIMENSION. C IMPLICIT DOUBLE PRECISION HAS BEEN REMOVED, AGAIN FOR STORAGE REASONS. C CALLING PROGRAM IS NOW AN INTEGRAL PART OF THE PACKAGE AS ALL THE MODEL- C DEPENDENT STUFF HAS BEEN SPUN OFF. C C IF YOU OBTAIN THIS CODE FROM A THIRD PERSON, PLEASE SEND YOUR NAME AND ADDRESS C TO S. CONSTABLE. YOU WILL THEN RECEIVE UPDATES, NEWS ON BUGS, ETC. C C SUBROUTINES WHICH MUST BE SUPPLIED BY THE USER: C C FORMOD(NP,ND,PM,DP,DM), COMPUTES THE FORWARD FUNCTION FOR MODEL PM() AT C THE DATA PARAMETERS DP() AND RETURNS IT IN DM(). C FORDIV(NP,ND,PM,DP,DM,WJ), COMPUTES THE FORWARD FUNCTION AS DOES FORMOD() C BUT ALSO RETURNS THE MATRIX OF PARTIALS, WJ(ND,NP). C OBJMAT(IRUF, DEL, NPEN), ASSEMBLES THE PENALTY MATRIX (THIS WILL BE DEPENDENT ON MODEL C TYPE AND DIMENSION). OBJMAT IS ALSO RESPONSIBLE FOR SETTING ANY C WEIGHTS, /RESULT/PREWTS(), ASSOCIATED WITH THE MODEL PREJUDICE, ALTHOUGH C THESE DEFAULT TO ZERO IF NOT SET. c DATIME(DATETM), RETURNS AN 80 CHARACTER STRING WITH CURRENT DATE AND TIME. CAN C BE A DUMMY ROUTINE IF SYSTEM PROGRAMMING IS TOO MUCH. C INPUTD(FILENAME), TAKES THE NAMED FILE AND USES IT TO GET THE DATA IMAGE INTO C THE COMMON BLOCK /DATA/. SEE OCCAM HEADER. C INPUTM(FILENAME), TAKES THE NAMED FILE AND READS IN ANY FIXED ATTRIBUTES C OF THE MODEL, SUCH AS LAYER THICKNESSES, FINITE ELEMENT MESH, ETC. COMMUNICATION C BETWEEN INPUTM() AND FORMOD() AND FORDIV() MUST BE VIA THE USER'S OWN COMMON C BLOCKS. WITH STATIC COMMONS (CHECK YOU COMPILER!) THESE NEED NOT BE DECLARED C IN THE CALLING PROGRAMS. INPUTM IS ALSO RESPONSIBLE FOR SETTING ANY C MODEL PREJUDICE IN /RESULT/ PREMOD(), ALTHOUGH THIS DEFAULTS TO ZERO C IF NOT SET. c c Files to be supplied by the user: c c STARTUP must be supplied, in the correct format, to provide inversion parameters. c data and model files pointed at by STARTUP c c Files generated by OCCAM: c c LOGFIL provides a blow by blow account of the inversion c ITERxx provides the inversion parameters, including the model, at every iteration. c ITER00 is a copy of the startup file at iteration 0. Any ITER file can be re- c named startup to continue an inversion. c c As the response of the model in dm() is not necessarily the response of the final c model found, but rather the response of the last model tested in the univariate c optimization routines, we leave it to the user to find the response of c the final model, using low cunning or the program 'FindResponse.f'. c c includes: include 'imp.inc' include 'occamdim.inc' c uses npp include 'occam.inc' c uses np, pm(), idebug, iout c c local variables: integer lerr, maxitr, iruf, nit, i, itmax, konv, ifftol c lerr = error flag for file operations c maxitr = maximum number of subsequent iterations by this run c iruf = roughness measure, 0 (size penalty), 1 (1st deriv.), 2 (2nd deriv.) c nit = current iteration number (0=start) c i = loop index c itmax = maximum iteration number c konv = convergence flag from OCCAM c ifftol = convergence info from OCCAM (has the desired misfit been achieved?) real tolreq, tobt, pmu, stepsz, rlast c tolreq = desired misfit, in RMS standard errors c tobt = current misfit c pmu = current Lagrange multiplier c stepsz = current RMS step though model space c rlast = last RMS measure of roughness character*80 itform, descr, modfil, datfil c itform = format number of iteration files c descr = description line in iteration files c modfil = model file name c datfil = data file name c c parameters integer iof1, iof2 parameter (iof1 = 13, iof2 = 15) c iof1 = unit number for the LOGFIL (passed to /result/iout) c iof2 = unit number for reading and writing ITER files c c calls: c occam, filerr, itrout, itrin, inputm, inputd c C----------------------------------------------------------------------- c Open startup file and read stuff. c (Startup and iteration files should be machine written, and so we do c not necessarily test the validity of the input data unless arrays are c involved (np and pm)). c open (iof2, file='startup', status='old', iostat=lerr) if (lerr .ne. 0) then write(*,*) ' Error opening startup file' stop end if call itrin (itform,descr,modfil,datfil,maxitr,tolreq,iruf, * nit,pmu,rlast,tobt,ifftol,iof2) close(iof2) c c if this is zeroth iteration copy the startup file into ITER00 if (nit .eq. 0) then c initialize the inversion parameters pmu = 5.0 rlast = 1.0e+10 tobt = 1000.0 ifftol = 0 call itrout (itform,descr,modfil,datfil,maxitr,tolreq, * iruf,nit,pmu,rlast,tobt,ifftol,iof2) end if c c zero out the model prejudice and weights in case inputm() and objmat() c do not take care of setting these. do 2 i = 1, np prewts(i) = 0.0 2 premod(i) = 0.0 c itmax = maxitr + nit c c input the data if (abs(idebug) .eq. 1) write(*,*) 'reading data..' call inputd(datfil) c input the fixed model characteristics if (abs(idebug) .eq. 1) write(*,*) 'reading model..' call inputm(modfil) c iout = iof1 open (iout, file='LOGFIL', iostat=lerr) if (lerr .ne. 0) then write(*,*) ' Error opening log file' stop end if write(*,*) 'Major I/O accomplished:' write(*,*) ' sit back and relax.......' c c do it 5 call occam(iruf,tolreq,tobt,ifftol,nit,pmu,rlast,stepsz,konv) c c write an iteration file call itrout (itform,descr,modfil,datfil,maxitr,tolreq,iruf, * nit,pmu,rlast,tobt,ifftol,iof2) c c let the console know roughly what is going on write(*,*) ' **ITERATION ', nit write(*,*) ' RMS MISFIT = ', tobt write(*,*) ' STEPSIZE = ', stepsz write(*,*) ' ' c c test for convergence: continue inversion if c 1) iteration number does not exceed maximum, and c 2) stepsize is still significant, and c 3) rms misfit larger than required, and c 4) there are no irrecoverable convergence problems c if (konv .ne. 0) then write (iof1,*) ' Convergence problems of type ', konv close(iof1) stop else if (nit .ge. itmax) then write (iof1,*) ' Stop on maximum iterations' close(iof1) stop else if ((tobt .ge. 1.01*tolreq) .or. (stepsz .ge. 0.0001)) then goto 5 else write (iof1,*) ' Stop on normal convergence' close(iof1) stop end if c c end c C----------------------------------------------------------------------- SUBROUTINE OCCAM(IRUF,TOLREQ,TOBT,IFFTOL,NIT,PMU,RLAST,STEPSZ, * KONV) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.01, 20 Jan 1993 C C OCCAM EXECUTES ONE ITERATION OF A SMOOTH MODEL FINDER C C ON INPUT: C C INCLUDE FILE 'OCCAM.INC': C /DATA/ CONTAINS C ND = NUMBER OF DATA C D(NDD) = VECTOR OF DATA VALUES C SD(NDD) = VECTOR OF ABSOLUTE STANDARD ERRORS IN THE DATA C DP(NDD,4) = MATRIX CONTAINING UP TO FOUR PARAMETERS ASSOCIATED WITH C EACH DATUM (RANGE, FREQUENCY, DATA TYPE ETC) C NPD = NUMBER OF PARAMETERS ASSOCIATED WITH EACH DATUM (E.G. FOR FREQ. ONLY C NPD = 1, FOR FREQ. AND RANGE NPD = 2) C C /MODEL/ CONTAINS C NP = NUMBER OF MODEL PARAMETERS, USUALLY THE NUMBER OF LAYERS C PM(NPP) = VECTOR OF INITIAL MODEL PARAMETERS, USUALLY LOG10(LAYER RES'IES) C C /RESULT/ CONTAINS C IDEBUG = 0 NORMAL, FASTEST AND QUIETEST EXECUTION C = 1 A LOW LEVEL OF VERBOSITY LETS YOU KNOW HOW THINGS ARE GOING C = 2 A TABLE OF RMS MISFIT VERSUS LAGRANGE MULTIPLIER IS PRODUCED C (USEFUL FOR TRACKING DOWN CONVERGENCE PROBLEMS). THIS TABLE IS USED C TO BRACKET THE MINIMUM C = 3 IF ONE ALSO WANTS A DUMP OF THE MODEL AT EACH LAGRANGE MULT IN C THE TABLE (NOT LIKELY TO BE REQUIRED BY END USERS) c if idebug is negative, the above absolute values are operational but c the step-reduction algorithm is disabled. C IOUT = FORTRAN UNIT NUMBER FOR OUTPUT, CAN EITHER BE TERMINAL OR A FILE C C ARGUMENTS: C IRUF = 0 FOR MODEL MINIMIZATION C = 1 FOR FIRST DERIVATIVE MINIMIZATION C = 2 FOR SECOND DERIVATIVE MINIMIZATION C TOLREQ = REQUIRED RMS MISFIT C TOBT = RMS MISFIT FROM LAST ITERATION (LARGE INITIALLY) C IFFTOL = RECORD OF WHETHER FEASIBLE MODEL HAS BEEN ATTAINED (0=NO, C 1=YES). INITIALLY SET 0 AND THEN LEFT ALONE. C NIT = THE NUMBER OF PREVIOUS ITERATIONS IN THIS INVERSION. SET TO 0 ON FIRST C CALL, IT WILL BE UPDATED BY OCCAM. C PMU = LOG10(LAGRANGE) MULTIPLIER TO START SEARCHES WITH. INITIALLY C SET LARGE AND THEN LEFT FROM LAST CALL FOR EFFICIENCY C RLAST = ROUGHNESS OF LAST MODEL. INITIALLY SET LARGE AND THEN LEFT. C C ON OUTPUT: C C /MODEL/ CONTAINS C PM(NPP) = VECTOR OF UPDATED MODEL PARAMETERS C DM(NDD) = VECTOR CONTAINING *APPROXIMATE* RESPONSE OF NEW MODEL C ARGUMENTS: C TOBT = RMS MISFIT OF NEW MODEL RESPONSE. C NIT = NUMBER OF PREVIOUS CALLS TO OCCAM1 DURING THIS INVERSION C STEPSZ = RMS MEASURE OF THE CHANGES IN THE MODEL PARAMETERS C KONV = STATUS FLAG: C 0 = NORMAL EXIT FOR A SUCCESSFUL ITERATION, C 1 = A PERFECTLY SMOOTH MODEL HAS BEEN FOUND FOR THE REQUIRED TOLERANCE C 2 = CONVERGENCE PROBLEMS. UNABLE TO FIND A MODEL WITH AN R.M.S. MISFIT C LESS THAN OR EQUAL TO THAT OF THE LAST ITERATION'S C 3 = CONVERGENCE PROBLEMS. UNABLE TO FIND A MODEL WITH AN R.M.S. MISFIT C OF TOLREQ AND SMOOTHER THAN THE LAST ITERATION'S C C SUBROUTINES REQUIRED AND SUPPLIED: C C MAKPTP, A SUBROUTINE WHICH CREATES THE ROUGHENING MATRIX TIMES ITS TRANS. C TOFMU(AMU), A FUNCTION WHICH RETURNS THE RMS MISFIT OF THE RESPONSE C OF THE MODEL CONSTRUCTED USING THE LAGRANGE MULTIPLIER AMU. C (CALLS CHOLIN,CHOLSL,FORMOD,ANORM) C INITM, TRMULT, MULT, ATAMULT AND ANORM, SUBROUTINES FOR MATRIX OPERATIONS C FMIN, A SUBROUTINE WHICH MINIMISES A UNIVARIATE FUNCTION C FROOT, A SUBROUTINE WHICH FINDS THE ROOT OF A UNIVARIATE FUNCTION C MINBRK, A SUBROUTINE WHICH BRACKETS A UNIVARIATE MINIMUM C SCANMU, A SUBROUTINE WHICH GENERATES THE MISFIT TABLE WHEN IBUG .GT. 0 C C SUBROUTINES WHICH MUST BE SUPPLIED BY THE USER: C C FORMOD(NP,ND,PM,DP,DM), COMPUTES THE FORWARD FUNCTION FOR MODEL PM() AT C THE DATA PARAMETERS DP() AND RETURNS IT IN DM(). C FORDIV(NP,ND,PM,DP,DM,WJ), COMPUTES THE FORWARD FUNCTION AS DOES FORMOD() C BUT ALSO RETURNS THE MATRIX OF PARTIALS, WJ(ND,NP). C OBJMAT(IRUF, DEL, NPEN), ASSEMBLES THE PENALTY MATRIX (THIS WILL BE DEPENDENT ON MODEL C TYPE AND DIMENSION) C C C INCLUDES: INCLUDE 'imp.inc' INCLUDE 'occamdim.inc' INCLUDE 'occam.inc' C C ARGUMENTS: INTEGER IRUF,IFFTOL,NIT,KONV REAL TOLREQ,TOBT,PMU,RLAST,STEPSZ C LOCAL VARIABLES: INTEGER I REAL TOL0, TMIN, TINT REAL AMU1, AMU2, AMU3, TAMU1, TAMU2, TAMU3 REAL RUF, PMU2 C FUNCTIONS REAL TOFMU, FMIN, FROOT, ANORM, FNDRUF EXTERNAL TOFMU C C LOCAL PARAMETERS REAL TOLM, TOLI C TOLM AND TOLI ARE THE TOLERANCES FOR FMIN AND FROOT RESPECTIVELY C DECREASING THEM WILL IMPROVE ACCURACY AT THE COST OF MORE FORWARD C CALCULATIONS PARAMETER (TOLM= 0.1, TOLI = 0.001) C C----------------------------------------------------------------------- IF (abs(idebug) .EQ. 1) WRITE(*,*) ' ENTERING OCCAM..' C SET THINGS UP FOR THIS ITERATION: KONV = 0 C FRAC CONTROLS THE STEP SIZE; NORMALLY WILL REMAIN AT 1.0 UNLESS WE HAVE C CONVERGENCE PROBLEMS FRAC = 1.0 C CONSTRUCT WJTWJ, WJTWD, AND FIND INITIAL MISFIT IF (abs(idebug) .EQ. 1) WRITE(*,*) * ' CONSTRUCTING DERIVATIVE DEPENDENT MATRICES..' CALL MAKJTJ(TOL0) C IF (NIT .EQ. 0) THEN WRITE(IOUT,*) 'STARTING R.M.S. = ',TOL0 if ((TOL0 .le. TOLREQ) .and. (IFFTOL .eq. 0)) then write (*,*) ' This model is feasible' IFFTOL = 1 end if END IF IF (abs(idebug) .EQ. 1) WRITE(*,*) 'STARTING R.M.S. = ',TOL0 WRITE(IOUT,*) ' ' WRITE(IOUT,*) '** ITERATION ',NIT+1,' **' C C CONSTRUCT THE PENALTY MATRIX MULTIPLIED INTO ITS TRANSPOSE IF (abs(idebug) .EQ. 1) WRITE(*,*)' CONSTRUCTING PENALTY MATRIX..' CALL MAKPTP(IRUF) C C THE NEXT BLOCK OF CODE CONTROLS THE SELECTION OF THE LAGRANGE MULTIPLIER, C PMU: IF (abs(idebug) .EQ. 1) * WRITE(*,*) ' SEARCHING LAGRANGE MULTIPLIER:' 120 NFOR = 0 C PRODUCE THE MISFIT FUNCTION IF REQUIRED (THIS WILL BE USED TO BRACKET MIN) IF (abs(idebug) .EQ. 1) WRITE(*,*) ' ..BRACKETING MINIMUM..' IF (abs(idebug) .GE. 2) THEN CALL SCANMU(AMU1,AMU2,AMU3,TAMU2,TOLREQ) ELSE C BRACKET THE MINIMUM USING MINBRK AND TWO GUESSES AMU1 = PMU - 1.0 AMU2 = PMU CALL MINBRK(AMU1,AMU2,AMU3,TAMU1,TAMU2,TAMU3,TOFMU) END IF C FIND THE MINIMUM IF (abs(idebug) .EQ. 1) WRITE(*,*) ' ..FINDING MINIMUM..' IF (TAMU2 .LT. TOLREQ) THEN C WE'VE BEEN LUCKY AND FOUND AN ACCEPTABLE MINIMUM USING MINBRK TMIN = TAMU2 PMU = AMU2 WRITE(IOUT,*) 'MINIMUM TOL FROM MINBRK IS AT MU =',PMU ELSE TMIN = FMIN(AMU1,AMU2,AMU3,TAMU2,TOFMU,TOLM,PMU) WRITE(IOUT,*) 'MINIMUM TOL FROM FMIN IS AT MU =',PMU END IF WRITE(IOUT,*) 'AND IS =',TMIN WRITE(IOUT,*) 'USING ',NFOR,' EVALUATIONS OF FUNCTION' C IF THE NEW MINIMUM TOLERANCE IS GREATER THAN THE TOLERANCE FROM THE C PREVIOUS MODEL, WE ARE HAVING CONVERGENCE PROBLEMS. CUT THE STEP SIZE IF (((TMIN .GT. TOL0) .AND. (IFFTOL .EQ. 0) .OR. * (TMIN .GT. TOLREQ) .AND. (IFFTOL .EQ. 1)) .and. * idebug .ge. 0) THEN WRITE(IOUT,*) 'DIVERGENCE PROBLEMS, CUTTING STEP SIZE' IF (abs(idebug) .EQ. 1) WRITE(*,*) ' ..DIVERGENCE PROBLEMS..' FRAC = FRAC*2.0 IF (FRAC .GT. 1.0E+05) THEN C WE HAVE CUT STEP SIZE A LOT TO NO EFFECT: GIVE UP KONV = 2 RETURN ELSE GOTO 120 END IF END IF C IF (TMIN .LT. TOLREQ) THEN IF (abs(idebug) .EQ. 1) WRITE(*,*) ' ..FINDING INTERCEPT..' C TOLERANCE IS BELOW THAT REQUIRED; FIND INTERCEPT. C THE LOWER VALUE OF MU BRACKETING THE INTERCEPT IS EASILY FOUND: IT IS C JUST THE MINIMUM AMU1 = PMU TAMU1 = TMIN C THE UPPER BOUND IS FOUND BY TESTING SUCCESSIVELY GREATER MU AMU2 = PMU TAMU2 = TMIN NFOR = 0 140 IF (TAMU2 .LT. TOLREQ) THEN C SUCCESSIVELY DOUBLE MU: AMU2 = AMU2 + 0.30103 TAMU2 = TOFMU(AMU2) GOTO 140 END IF PMU2 = FROOT(TOFMU,AMU1,AMU2,TAMU1,TAMU2,TOLREQ,TOLI) TINT = TAMU2 + TOLREQ WRITE(IOUT,*) 'INTERCEPT IS AT MU = ',PMU2 WRITE(IOUT,*) 'AND IS = ',TINT WRITE(IOUT,*) 'USING ',NFOR,' FUNCTION EVALUATIONS' ELSE C TOLERANCE IS NO YET SMALL ENOUGH. WE WILL KEEP THE MINIMUM. C SINCE FMIN RETURNS WITH PWK1() INSTEAD OF PWK2() WE NEED THIS COPY DO 150 I = 1,NP 150 PWK2(I) = PWK1(I) END IF C***END LAGRANGE MULTIPLIER SELECTION C C COMPUTE ROUGHNESS. WE DO THIS BY A FUNCTION CALL TO AVOID HAVING THE PENALTY C MATRIX HANG AROUND TAKING UP MEMORY. IF (abs(idebug) .EQ. 1) WRITE(*,*) 'COMPUTING ROUGHNESS..' RUF = FNDRUF(IRUF) C C IF WE ATTAINED THE INTERCEPT LAST ITERATION BUT THE MODEL IS GETTING C ROUGHER WE HAVE PROBLEMS WITH CONERGENCE. IF ((RUF.GT.1.01*RLAST) .AND. (IFFTOL.EQ.1) .and. (idebug.ge.0)) * THEN WRITE(IOUT,*) 'ROUGHNESS PROBLEMS, CUTTING STEP SIZE' FRAC = FRAC*2.0 C CHECK TO SEE IF ALL IS HOPELESS IF (FRAC .GT. 1.0E+05) THEN KONV = 3 RETURN END IF C OTHERWISE PLOUGH ON GOTO 120 END IF C C SAVE NEW MODEL AND COMPUTE STEP SIZE IF (abs(idebug) .EQ. 1) WRITE(*,*) 'COMPUTING STEPSIZE..' DO 180 I = 1,NP PWK3(I) = PWK2(I) - PM(I) 180 PM(I) = PWK2(I) C THE STEPSIZE IS THE ACTUAL CHANGE IN THE MODEL, NORMALISED BY NP: STEPSZ = SQRT(ANORM(NP,PWK3)/NP) WRITE(IOUT,*) 'STEPSIZE IS = ',STEPSZ WRITE(IOUT,*) 'ROUGHNESS IS = ',RUF C TIDY UP NIT = NIT + 1 RLAST = RUF IF (TMIN .LT. TOLREQ) THEN IFFTOL = 1 TOBT = TINT ELSE TOBT = TMIN END IF C SEE IF WE HAVE A PERFECLY SMOOTH MODEL THAT FITS DATA: IF (RUF .LT. 1.0E-5 .AND. TOBT .LE. 1.01*TOLREQ) KONV = 1 IF (abs(idebug) .EQ. 1) WRITE(*,*) 'LEAVING OCCAM' RETURN END C C----------------------------------------------------------------------- SUBROUTINE SCANMU(AM1,AM2,AM3,T2,T0) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.00, 13 May 1992 C C SCANMU PRODUCES A REPORT OF MISFIT VERSUS LAGRANGE MULTIPLIER. C IT CALLS THE FORWARD ROUTINE A LOT SO IS REALLY ONLY USED FOR DEBUGGING C AND PARTICULARLY NASTY CASES C IT LOOKS FOR THE GLOBAL MINIMUM AS FOLLOWS: C STARTING AT THE LARGEST MU USED FOR THE DEBUG REPORT, WE LOOK FOR THE C FIRST MINIMUM. WE ALSO LOOK FOR THE GLOLBAL MIN. WE WANT THE GLOBAL C MIN UNLESS THE LOCAL MIN AT LARGER MU IS LOWER THAN THE REQUIRED C TOLERANCE. C C ON INPUT: C T0 IS THE REQUIRED TOLERANCE C /RESULT/ IDEBUG IS THE PRINT LEVEL FROM OCCAM1 C /RESULT/ IOUT IS A UNIT NUMBER FOR OUTPUT C ON OUTPUT: C AM1, AM2, AM3 ARE THREE LAGRANGE MULTIPLIERS WHICH BRACKET A MINIMUM C T2 IS THE TOLERANCE AT AM2. C /RESULT/ ..PMW1().. HAS THE MODEL ASSOCIATED WITH AM2,T2. C SUBROUTINES CALLED: C TOFMU C C INCLUDES: INCLUDE 'imp.inc' INCLUDE 'occamdim.inc' INCLUDE 'occam.inc' C C ARGUMENTS: REAL AM1,AM2,AM3,T2,T0 C LOCAL VARIABLES LOGICAL SEARCH INTEGER I,J,K REAL T1,TG2,T3,TL2,AML3,AML2,AML1,AMG3,AMG2,AMG1,AMU,T C FUNCTIONS REAL TOFMU C C----------------------------------------------------------------------- WRITE(IOUT,*) 'MISFIT AS A FUNCTION OF MU:' T1 = 1.0E+09 T2 = 1.1E+09 TG2 = 0.9E+09 AM2 = 16. AM1 = 16. SEARCH = .TRUE. C C START AT THE TOP (SMOOTHEST AND PRESUMABLY WORST FIT), SWEEPING FROM C LOG10(MU) = 16 TO LOG10(MU) = -4 AT TWO VALUES PER DECADE DO 100 K = 16,-4,-1 C KEEP THE LAST TWO VALUES FOR MINIMUM SEARCH T3 = T2 T2 = T1 AM3 = AM2 AM2 = AM1 C KEEP THE MODEL ASSOCIATED WITH THE MIDDLE VALUE IF WE ARE SEARCHING C FOR THE MINIMUM IF (SEARCH) THEN DO 10 I = 1,NP 10 PWK1(I) = PWK2(I) END IF C KEEP THE MODEL ASSOCIATED WITH THE MIDDLE VALUE ANYWAY IN CASE C GLOBAL AND LOCAL MINIMA ARE FOUND DO 20 I = 1,NP 20 PWK3(I) = PWK2(I) C GET THE NEXT MU AND ASSOCIATED TOLERANCE AM1 = FLOAT(K)/2. T1 = TOFMU(AM1) C IF ((T2 .LE. T1) .AND. (T2 .LE. T3) .AND. SEARCH) THEN C WE HAVE FOUND THE FIRST MINIMUM. KEEP IT AND TURN OFF SEARCH C (MODEL IS ALREADY STASHED IN PM1()) SEARCH = .FALSE. TL2 = T2 AML3 = AM3 AML2 = AM2 AML1 = AM1 END IF C KEEP AN EYE OUT FOR THE GLOBAL MIN IN CASE IT IS DIFFERENT IF (T2 .LE. TG2) THEN TG2 = T2 AMG3 = AM3 AMG2 = AM2 AMG1 = AM1 C COPY ACCROSS MODEL ASSOCIATED WITH GLOBAL MIN. DO 30 I = 1,NP 30 PWK1(I) = PWK3(I) END IF C RECORD IN THE LOG FILE WRITE(IOUT,*) AM1,T1 C C OUTPUT THE MODEL AS WELL IF REQUIRED IF (abs(idebug) .GE. 3) THEN DO 60 J = 1,NP 60 WRITE(IOUT,*) ' ',PWK2(J) END IF 100 CONTINUE C WRITE OUT THE GAUSS STEP AT HIGHEST DEBUG LEVEL IF (abs(idebug) .GE. 3) THEN AMU = -99.0 T = TOFMU(AMU) WRITE(IOUT,*) AM1,T1 DO 110 J = 1,NP 110 WRITE(IOUT,*) ' ',PWK2(J) END IF C C NOW CHECK OUT THE MINIMA C IF THE FIRST MINIMUM ACHIEVES THE TOLERANCE REQUIRED WE WANT IT IF (TL2 .LE. T0) THEN T2 = TL2 AM1 = AML1 AM2 = AML2 AM3 = AML3 ELSE C OTHERWISE THE GLOBAL MIN IS JUST FINE T2 = TG2 AM1 = AMG1 AM2 = AMG2 AM3 = AMG3 END IF RETURN END C C----------------------------------------------------------------------- REAL FUNCTION TOFMU(ALOGMU) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.01, 13 Jan 1993 C C FUNCTION TOFMU RETURND THE RMS MISFIT OF THE RESPONSE OF A MODEL CONSTRUCTED C USING THE GIVEN VALUE OF LAGRANGE MULTIPLIER C C ON INPUT: C ALOGMU = LOG10 OF THE LAGRANGE MULTIPLIER C (-99. IF A ZERO LAGRANGE MULTIPLIER NEEDED) C /RESULT/ WJTWJ(), WJTWD(), PTP(), PREWTS(), PREMOD(), NFOR C /DATA/ DP(), D(), SD(), ND C /MODEL/ NP C C ON OUTPUT C TOFMU = R.M.S. MISFIT OR 1.0E+10 IF CHOLESKY DECOMPOSITION FAILED C /RESULT/ PWK2() CONTAINS THE MODEL WHICH PRODUCES THE MISFIT TOFMU C /MODEL/ DM() CONTIANS THE RESPONSE OF THE MODEL C C SUBROUTINES CALLED: C CHOLIN, CHOLSL ARE R.L. PARKER'S SUBROUTINES TO PERFORM CHOLESKY DECOMP. C AND THEN SOLVE A LINEAR SYSTEM BY BACK SUBSTITUTION. C ANORM, FORMOD ARE EXPLAINED IN OCCAM1 C C AT VERSION 1.4 LOG10(MU) USED TO IMPROVE PERFORMANCE OF 1D OPTIMISATION c C INCLUDES: INCLUDE 'imp.inc' INCLUDE 'occamdim.inc' INCLUDE 'occam.inc' C C ARGUMENTS: REAL ALOGMU C LOCAL VARIABLES REAL A(NPP,NPP+1), DWK1(NDD), AMU, PART, PWK4(NPP) INTEGER I, J, ISTAT C FUNCTIONS: REAL ANORM C C----------------------------------------------------------------------- IF (ALOGMU. EQ. -99.) THEN C THE GAUSS STEP IS REQUIRED AMU = 0.0 ELSE AMU = 10.**(ALOGMU) END IF NFOR = NFOR + 1 C ADD AMU.TRANS(PARTIAL).PARTIAL TO TRANS(W.J).W.J C ALSO ADD MODEL PREJUDICE TO WJTWD() DO 20 I = 1,NP DO 20 J = 1,NP PWK4(J) = WJTWD(J) + AMU*PREWTS(J)*PREMOD(J) 20 A(I,J) = AMU*PTP(I,J) + WJTWJ(I,J) C C SOLVE THE LINEAR SYSTEM CALL CHOLIN(NPP,NP,A,1,ISTAT) IF (ISTAT .NE. 0) THEN C IF DECOMP FAILED WE GIVE WARNING BUT MIGHT STILL GET BY IF TOFMU C IS SET TO A LARGE VALUE WRITE(*,*) * 'Warning:- Cholesky decomp. failed at mu = ',AMU TOFMU = 1.0E+10 RETURN END IF CALL CHOLSL(NPP,NP,A,PWK2,PWK4,0) C CUT STEP SIZE IF NECESSARY IF (FRAC .GT. 1.001) THEN PART = 1./FRAC DO 30 I = 1,NP 30 PWK2(I) = (1.-PART)*PM(I) + PART*PWK2(I) END IF C CALCULATE THE MODEL RESPONSE AND MISFIT CALL FORMOD(NP,ND,PWK2,DP,DM) DO 40 I = 1,ND 40 DWK1(I) = (D(I) - DM(I))/SD(I) TOFMU = SQRT(ANORM(ND,DWK1)/ND) IF (abs(idebug) .EQ. 1) * WRITE(IOUT,*) 'TOFMU =', TOFMU,' AMU =', ALOGMU RETURN END C C C----------------------------------------------------------------------- REAL FUNCTION ANORM(N,D) C----------------------------------------------------------------------- C RETURNS THE SQUARE OF THE EUCLIDEAN NORM OF A VECTOR INCLUDE 'imp.inc' C INTEGER I, N REAL D(N) C ANORM = 0.0 DO 10 I = 1,N ANORM = ANORM + D(I)*D(I) 10 CONTINUE RETURN END C C C----------------------------------------------------------------------- SUBROUTINE INITM(MD,M,N,A) C----------------------------------------------------------------------- C ZEROS A 2D MATRIX INCLUDE 'imp.inc' C INTEGER MD, M, N, I, J REAL A(MD,N) C DO 10 I = 1,M DO 10 J = 1,N 10 A(I,J) = 0.0 RETURN END C C C----------------------------------------------------------------------- SUBROUTINE MULT(MAD,MA,NAD,NA,NB,A,B,C) C----------------------------------------------------------------------- C MULTIPLIES TWO 2D MATRICES INCLUDE 'imp.inc' C INTEGER MAD, MA, NAD, NA, NB, I, J, K REAL A(MAD,NA),B(NAD,NB),C(MAD,NB), CIJ C DO 12 I = 1,MA DO 12 J = 1,NB CIJ = 0.0 DO 10 K = 1,NA 10 CIJ = A(I,K)*B(K,J) + CIJ 12 C(I,J) = CIJ RETURN END C C C----------------------------------------------------------------------- SUBROUTINE TRMULT(MAD,MA,NAD,NA,NB,A,B,C) C----------------------------------------------------------------------- C MULTIPLIES THE TRANSPOSE OF A 2D MATRIX BY ANOTHER MATRIX INCLUDE 'imp.inc' C INTEGER MAD, MA, NAD, NA, NB, I, J, K REAL A(MAD,NA),B(MAD,NB),C(NAD,NB), CIJ DO 12 I = 1,NA DO 12 J = 1,NB CIJ = 0.0 DO 10 K = 1,MA 10 CIJ = A(K,I)*B(K,J) + CIJ 12 C(I,J) = CIJ RETURN END C C C----------------------------------------------------------------------- SUBROUTINE MINBRK(AX,BX,CX,FA,FB,FC,FUNC) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.00, 13 May 1992 C C MINBRK BRACKETS A UNIVARIATE MINIMUM OF A FUNCTION. TO BE USED PRIOR TO C A UNIVARIATE MINIMISATION ROUTINE. C BASED ON A ROUTINE IN NUMERICAL RECIPES BY PRESS ET AL C MODIFIED SO THAT THE MODEL ASSOCIATED WITH THE MISFITS IS CARRIED AROUND C FOR USE IN THE MINIMISATION ROUTINES AND POSSIBLY ULTIMATELY KEPT AS THE C RESULT OF THIS ITERATION; C C ON INPUT: C AX, BX = TWO DISTINCT ESTIMATES OF THE MINIMUM'S WHEREABOUTS C FUNC = THE FUNCTION IN QUESTION C /MODEL/ NP = NUMBER OF MODEL PARAMETERS C ON OUTPUT: C AX,BX,CX = THREE NEW POINTS WHICH BRACKET THE MINIMUM C FA,FB,FC = FUNC(AX), FUNC(BX) ETC. C /RESULT/ PWK1() = THE MODEL ASSOCIATED WITH BX AND FB C C INCLUDES: INCLUDE 'imp.inc' INCLUDE 'occamdim.inc' INCLUDE 'occam.inc' C C ARGUMENTS: REAL AX,BX,CX,FA,FB,FC C LOCAL VARIABLES: INTEGER I REAL DUM,R,Q,ULIM,U,FU C LOCAL PARAMETERS REAL GOLD,GLIMIT,TINY PARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.E-32) C FUNCTIONS: REAL FUNC C C----------------------------------------------------------------------- FB = FUNC(BX) C MAKE PWK1 ASSOCIATED WITH B DO 4 I = 1,NP 4 PWK1(I) = PWK2(I) C PM2 WILL BE ASSOCIATED WITH A FA = FUNC(AX) IF (FB.GT.FA) THEN C SWITCH A AND B DUM = AX AX = BX BX = DUM DUM = FB FB = FA FA = DUM C KEEP PWK1 ASSOCIATED WITH B (A WILL BE LOST SO DOESN'T MATTER) DO 5 I = 1,NP 5 PWK1(I) = PWK2(I) ENDIF CX = BX + GOLD*(BX - AX) FC = FUNC(CX) C MAIN TEST. IF IT FAILS WE EXIT WITH B 10 IF (FB.GE.FC) THEN C KEEP PWK1 ASSOCIATED WITH C DO 11 I = 1,NP 11 PWK1(I) = PWK2(I) R = (BX - AX)*(FB - FC) Q = (BX - CX)*(FB - FA) U = BX-((BX-CX)*Q-(BX-AX)*R)/(2.*SIGN(MAX(ABS(Q-R),TINY),Q-R)) ULIM = BX + GLIMIT*(CX - BX) IF ((BX - U)*(U - CX).GT.0.) THEN FU = FUNC(U) IF (FU.LT.FC) THEN AX = BX FA = FB BX = U FB = FU C KEEP PWK1 ASSOCIATED WITH B DO 12 I = 1,NP 12 PWK1(I) = PWK2(I) GO TO 10 ELSE IF (FU.GT.FB) THEN CX = U FC = FU GO TO 10 ENDIF U = CX + GOLD*(CX - BX) FU = FUNC(U) ELSE IF ((CX - U)*(U - ULIM).GT.0.) THEN FU = FUNC(U) IF (FU.LT.FC) THEN BX = CX CX = U U = CX + GOLD*(CX - BX) FB = FC FC = FU C KEEP PWK1 ASSOCIATED WITH B DO 13 I = 1,NP 13 PWK1(I) = PWK2(I) FU = FUNC(U) ENDIF ELSE IF ((U - ULIM)*(ULIM - CX).GE.0.) THEN U = ULIM FU = FUNC(U) ELSE U = CX + GOLD*(CX - BX) FU = FUNC(U) ENDIF AX = BX BX = CX CX = U FA = FB FB = FC FC = FU GO TO 10 ENDIF RETURN END C C C C----------------------------------------------------------------------- REAL FUNCTION FMIN(AX,BX,CX,FBX,F,TOL,XMIN) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.00, 13 May 1992 C C FMIN RETURNS THE MINIMUM VALUE OF A FUNCTION C BASED ON A ROUTINE IN NUMERICAL RECIPES BY PRESS ET AL C MODIFIED SO THAT THE MODEL ASSOCIATED WITH THE MISFITS IS CARRIED AROUND C FOR USE POSSIBLE ULTIMATE USE AS THE C RESULT OF THIS ITERATION; C C ON INPUT: C AX,BX,CX = INDEPENDENT VARIABLE WHICH BRACKET THE MINIMUM C FBX = F(BX) (USUALLY AVAILABLE FROM THE BRACKETING PROCEDURE) C F = FUNCTION IN QUESTION C TOL = FRACTIONAL TOLERANCE REQUIRED IN THE INDEPENDENT VARIABLE C /RESULT/ PWK1 = MODEL ASSOCIATED WITH BX AND FBX C /MODEL/ NP = NUMBER OF MODEL PARAMETERS C ON OUTPUT: C XMIN = ABSCISSA OF MINIMUM C FMIN = F(XMIN) C /RESULT/ PWK1 = MODEL ASSOCIATED WITH XMIN AND FMIN C C INCLUDES: INCLUDE 'imp.inc' INCLUDE 'occamdim.inc' INCLUDE 'occam.inc' C C ARGUMENTS: REAL AX,BX,CX,FBX,TOL,XMIN C LOCAL VARIABLES REAL AA,B,V,W,X,E,FX,FV,FW,TOL1,Q,R,P,ETEMP,DD,FU,U,XM,TOL2 INTEGER I,ITER C LOCAL PARAMETERS INTEGER ITMAX REAL CGOLD,ZEPS PARAMETER (ITMAX=100,CGOLD=.3819660,ZEPS=1.0E-10) C FUNCTION: REAL F C C----------------------------------------------------------------------- AA = MIN(AX,CX) B = MAX(AX,CX) V = BX W = V X = V E = 0. FX = FBX FV = FX FW = FX DO 11 ITER = 1,ITMAX XM = 0.5*(AA + B) TOL1 = TOL*ABS(X) + ZEPS TOL2 = 2.*TOL1 IF (ABS(X - XM).LE.(TOL2 - .5*(B - AA))) GOTO 30 IF (ABS(E).GT.TOL1) THEN R = (X - W)*(FX - FV) Q = (X - V)*(FX - FW) P = (X - V)*Q - (X - W)*R Q = 2.*(Q - R) IF (Q.GT.0.) P = - P Q = ABS(Q) ETEMP = E E = DD IF (ABS(P).GE.ABS(.5*Q*ETEMP).OR.P.LE.Q*(AA - X) .OR. * P.GE.Q*(B - X)) GOTO 1 DD = P/Q U = X + DD IF (U - AA.LT.TOL2 .OR. B - U.LT.TOL2) DD = SIGN(TOL1,XM - X) GOTO 2 ENDIF 1 IF (X.GE.XM) THEN E = AA - X ELSE E = B - X ENDIF DD = CGOLD*E 2 IF (ABS(DD).GE.TOL1) THEN U = X + DD ELSE U = X + SIGN(TOL1,DD) ENDIF FU = F(U) IF (FU.LE.FX) THEN IF (U.GE.X) THEN AA = X ELSE B = X ENDIF V = W FV = FW W = X FW = FX X = U FX = FU DO 5 I = 1,NP 5 PWK1(I) = PWK2(I) ELSE IF(U.LT.X) THEN AA = U ELSE B = U ENDIF IF (FU.LE.FW .OR. W.EQ.X) THEN V = W FV = FW W = U FW = FU ELSE IF (FU.LE.FV .OR. V.EQ.X .OR. V.EQ.W) THEN V = U FV = FU ENDIF ENDIF 11 CONTINUE WRITE(*,*) 'MAXIMUM ITERATIONS EXCEEDED IN FMIN' 30 XMIN = X FMIN = FX RETURN END C C C----------------------------------------------------------------------- REAL FUNCTION FROOT(FUNC,X1,X2,FA,FB,OFF,TOL) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.00, 13 May 1992 C C FROOT FINDS THE POINT AT WHICH A UNIVARIATE FUNCTION ATTAINS A GIVEN C VALUE C BASED ON A ROUTINE IN NUMERICAL RECIPES BY PRESS ET AL C MODIFIED SO THAT THE MODEL ASSOCIATED WITH THE MISFITS IS CARRIED AROUND C FOR USE POSSIBLE ULTIMATE USE AS THE C RESULT OF THIS ITERATION; C C ON INPUT: C FUNC = THE FUNCTION IN QUESTION C X1,X2 = INDEPENDENT VARIABLES BRACKETING THE ROOT C FA,FB = FUNC(X1),FUNC(X2) (USUALLY AVAILABLE FROM THE BRACKETING) C OFF = VALUE OF THE FUNCTIONAL AT THE REQUIRED ABSCISSA C TOL = FRACTIONAL TOLERANCE REQUIRED FOR ABSCISSA C /RESULT/ PWK2 = MODEL ASSOCIATED WITH X2 AND FB C /RESULT/ PWK1 = MODEL ASSOCIATED WITH X1 AND FA C /MODEL/ NP = NUMBER OF MODEL PARAMTERS C C ON OUTPUT: C FB = FUNC(FROOT) - OFF C FROOT = ABSCISSA REQUIRED C /RESULT/ PWK2 = MODEL ASSOCIATED WITH FB, FROOT C C INCLUDES: INCLUDE 'imp.inc' INCLUDE 'occamdim.inc' INCLUDE 'occam.inc' C C ARGUMENTS: REAL X1,X2,FA,FB,OFF,TOL C LOCAL VARIABLES REAL AA,B,FC,C,E,DD,TOL1,P,Q,R,S,XM INTEGER I,ITER C LOCAL PARAMETERS INTEGER ITMAX REAL EPS PARAMETER (ITMAX = 100,EPS = 3.E-8) C FUNCTIONS: REAL FUNC C C----------------------------------------------------------------------- AA = X1 B = X2 FA = FA - OFF FB = FB - OFF IF (FB*FA.GT.0.) WRITE(*,*) 'ROOT NOT BRACKETED IN FROOT' FC = FB C MODELS WILL BE CARRIED AROUND AS: C AA: PWK1() C B: PWK2() C C: PWK3() C SO COPY PM2 (WITH B) INTO WK (WITH C) DO 1 I = 1,NP 1 PWK3(I) = PWK2(I) DO 11 ITER = 1,ITMAX C IF B AND C ARE NOT ON DIFFERENT SIDES OF ROOT THEN MAKE THEM SO IF (FB*FC.GT.0.) THEN C = AA FC = FA DD = B - AA E = DD C COPY PMF (WITH AA) INTO WK (WITH C) DO 2 I = 1,NP 2 PWK3(I) = PWK1(I) ENDIF IF (ABS(FC).LT.ABS(FB)) THEN C IF B IS NOT CLOSER THAN C MAKE IT SO AA = B B = C C = AA FA = FB FB = FC FC = FA C SWITCH PM2 (WITH B) AND WK (WITH C) (AA WILL BE LOST IF WE DON'T EXIT) DO 3 I = 1,NP PWK1(I) = PWK2(I) PWK2(I) = PWK3(I) 3 PWK3(I) = PWK1(I) ENDIF TOL1 = 2.*EPS*ABS(B) + 0.5*TOL XM = .5*(C - B) IF (ABS(XM).LE.TOL1 .OR. FB.EQ.0.) THEN FROOT = B RETURN ENDIF IF (ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN S = FB/FA IF (AA.EQ.C) THEN P = 2.*XM*S Q = 1. - S ELSE Q = FA/FC R = FB/FC P = S*(2.*XM*Q*(Q - R) - (B - AA)*(R - 1.)) Q = (Q - 1.)*(R - 1.)*(S - 1.) ENDIF IF (P.GT.0.) Q = - Q P = ABS(P) IF (2.*P .LT. MIN(3.*XM*Q - ABS(TOL1*Q),ABS(E*Q))) THEN E = DD DD = P/Q ELSE DD = XM E = DD ENDIF ELSE DD = XM E = DD ENDIF AA = B FA = FB C COPY PM2 (WITH B) INTO PMF (WITH AA) DO 4 I = 1,NP 4 PWK1(I) = PWK2(I) IF (ABS(DD) .GT. TOL1) THEN B = B + DD ELSE B = B + SIGN(TOL1,XM) ENDIF C THE NEW FUNCION EVALUATION AUTOMATICALLY USES PM2 FOR B FB = FUNC(B) - OFF 11 CONTINUE WRITE(*,*) 'MAXIMUM ITERATIONS EXCEEDED IN FROOT' FROOT = B RETURN END C C C----------------------------------------------------------------------- SUBROUTINE CHOLSL(LA,N, A, X, B, ISTAGE) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 C WRITTEN BY: R.L. PARKER C C$$$$$ CALLS NO OTHER ROUTINES C ROUTINE FOR BACK-SUBSTITUTION SOLUTION OF A LINEAR SYSTEM THAT HAS C ALREADY BEEN FACTORED INTO CHOLESKY FACTORS BY 'CHOLIN' (Q.V.) C LA LEADING DIMENSION OF A C N THE ORDER OF MATRIX C A ARRAY CONTAINING ON AND ABOVE ITS DIAGONAL THE TRANSPOSE OF C THE CHOLESKY FACTOR. THE NORMAL WAY OF FINDING THE FACTORS C IS TO CALL CHOLIN AS FOLLOWS C CALL CHOLIN(N, A, 1, IFBAD) C CALL CHOLSL(N, A, X, B, 0) C X THE UNKNOWN VECTOR. SOLUTION TO EQUATIONS IS RETURNED HERE C B THE RIGHT SIDE OF SYSTEM C ISTAGE THE MODE OF THE SOLUTION. IF=1 DO NOT COMPLETE BACK- C SUBSTITUTION. INSTEAD RETURN INVERSE(L-TRANS)*B . OTHER VALUES C FORM THE COMPLETE SOLUTION. C R.L. PARKER C C INCLUDES: INCLUDE 'imp.inc' C ARGUMENTS INTEGER LA,N,ISTAGE REAL A(LA,N+1),X(N),B(N) C LOCAL VARIABLES INTEGER I,K,II,KK,I1 REAL Z C C----------------------------------------------------------------------- DO 2500 I=1,N Z=B(I) IF (I.EQ.1) GO TO 2500 DO 2200 KK=2,I K=I-KK+1 2200 Z=Z - A(K,I+1)*X(K) 2500 X(I)=Z/A(I,I+1) IF (ISTAGE.EQ.1) RETURN C COMPLETE BACK-SUBSTITUTION DO 2700 II=1,N I=N-II+1 Z=X(I) I1=I+1 IF (I.EQ.N) GO TO 2700 DO 2600 K=I1,N 2600 Z= Z - A(I,K+1)*X(K) 2700 X(I)=Z/A(I,I+1) RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE CHOLIN(LA,N, A, JSTAGE, IFBAD) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 C WRITTEN BY: R.L. PARKER C C$$$$$ CALLS NO OTHER ROUTINES C INVERSION AND FACTORIZATION OF POSITIVE DEFINITE SYMMETRIC MATRICES C FROM ALGOL ROUTINE OF WILKINSON & REINSCH 'LINEAR ALGEBRA' PP16-17. C EXPLANATION C CHOLIN CALCULATES THE FACTOR L IN THE CHOLESKY FACTORIZATION C A = L L-TRANS, OR THE INVERSE OF L , OR THE INVERSE OF A. C TRANSPOSED STORAGE IS FOR THE CONVENIENCE OF FORTRAN USERS. HOWEVER, C IT MEANS THAT L-TRANS AND INVERSE OF L-TRANS ARE COMPUTED. C ARGUMENTS C A AN ARRAY DIMENSION N,N+1. ON ENTRY THE INPUT MATRIX SHOULD C SHOULD OCCUPY THE 1ST N COLS, OR IF DESIRED, JUST THE LOWER C LOWER LEFT TRIANGLE. ON EXIT THE UPPER TRIANGLE OF THE ARRAY C DEFINED BY COLS 2,3...N+1 CONTAINS THE DESIRED RESULT. C THIS CAN BE CONVENIENTLY ACCESSED AS A(1,2), OR A(N+1) IF C SINGLY DIMENSIONED IN CALLING PROGRAM. C N THE ORDER OF THE MATRIX C LA LEADING DIMENSION OF A C JSTAGE INDICATES DESIRED CALCULATION - IF=1 UPPER TRIANGLE CONTAINS C L-TRANS. IF=2 UPPER TRIANGLE CONTAINS INVERSE OF L-TRANS. C IF=0 OR 3 CONTAINS THE UPPER HALF OF THE INVERSE OF A. NOTE C LOWER TRIANGLE IS NOT DESTROYED. IF=-2,-3, THE CHOLESKY C FACTORS ARE ASSUMED TO BE IN PLACE ON ENTRY. C IFBAD ERROR FLAG, SET TO 1 IF MATRIX NOT POSITIVE DEFINITE. SET TO C 0 IF OK. C R.L. PARKER C C INCLUDES: INCLUDE 'imp.inc' C PASS PARAMETERS INTEGER LA,N,JSTAGE,IFBAD REAL A(LA,N+1) C LOCAL VARIABLES INTEGER I,J,K,KBACK,ISTAGE,I1,J1,J2,N1 REAL X,Y,Z C C----------------------------------------------------------------------- C STAGE 1 - FORM FACTORS (HERE L-TRANS IS STOTED ABOVE DIAGONAL) ISTAGE=IABS(JSTAGE) IF (JSTAGE.LT.0) GO TO 1550 IFBAD=0 DO 1500 I=1,N I1=I + 1 DO 1400 J=I,N J1=J + 1 X=A(J,I) IF (I.EQ.1) GO TO 1200 DO 1100 KBACK=2,I K=I-KBACK+1 1100 X=X - A(K,J1)*A(K,I1) 1200 IF (I.NE.J) GO TO 1300 IF (X.LE.0.0) GO TO 4000 Y=1.0/SQRT(X) A(I,I1)=Y GO TO 1400 1300 A(I,J1)=X*Y 1400 CONTINUE 1500 CONTINUE IF (ISTAGE.NE.1) GO TO 2000 1550 DO 1600 I=1,N 1600 A(I,I+1)=1.0/A(I,I+1) IF (ISTAGE.EQ.1) RETURN C STAGE 2 - INVERSION OF L-TRANS 2000 DO 2600 I=1,N I1=I+1 IF (I1.GT.N) GO TO 2600 DO 2500 J=I1,N Z=0.0 J1=J + 1 J2=J - 1 DO 2200 KBACK=I,J2 K=J2-KBACK+I 2200 Z=Z - A(K,J1)*A(I,K+1) A(I,J1)=Z*A(J,J1) 2500 CONTINUE 2600 CONTINUE IF (ISTAGE.EQ.2) RETURN C STAGE 3 - CONSTRUCTION OF IVERSE OF A ABOVE DIAGONAL. DO 3500 I=1,N DO 3500 J=I,N Z=0.0 N1=N + 1 J1=J + 1 DO 3200 K=J1,N1 3200 Z=Z + A(J,K)*A(I,K) A(I,J+1)=Z 3500 CONTINUE RETURN C ERROR EXIT 4000 IFBAD=1 C ERROR MESSAGE SUPRESSED - S.CONSTABLE C WRITE(*,400) C 400 FORMAT(44H0CHOLIN FAILS - MATRIX NOT POSITIVE DEFINITE ) RETURN END C C C----------------------------------------------------------------------- SUBROUTINE MAKJTJ(TOL0) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.00, 13 May 1992 C C CONSTRUCTS MATRICES THAT DEPEND ON THE JACOBIAN (WJTWJ AND WJTWD). C THE JACOBIAN ONLY NEEDS TO BE AROUND LONG ENOUGH TO DO C THIS, SO BY STUFFING THIS INTO A SUBROUTINE DYNAMIC MEMORY ALLOCATION C WILL MAKE WJ DISSAPEAR AFTER USE. C C ON INPUT: C CURRENT MODEL AND DATA STORED IN COMMON /MODEL/ AND /DATA/ C ON OUTPUT: C /RESULT/ WJTWJ = W.J.TRANS(W.J) C /RESULT/ WJTWD = WEIGHTED, TRANSLATED DATA PREMULTIPLIED BY TRANS(W.J) C TOL0 = RMS MISFIT ASSOCIATED WITH CURRENT MODEL C CALLS: C ATAMUL, MULT, TRMULT, FORDIV, ANORM C C INCLUDES: INCLUDE 'imp.inc' INCLUDE 'occamdim.inc' INCLUDE 'occam.inc' C C ARGUMENTS: REAL TOL0 C LOCAL VARIABLES REAL WJ(NDD,NPP), DWK1(NDD), DHAT(NDD) INTEGER I, J C FUNCTIONS REAL ANORM C C----------------------------------------------------------------------- C CALCULATE THE MATRIX OF PARTIALS AND MODEL RESPONSE CALL FORDIV(NP,ND,PM,DP,DM,WJ) C CALC MISFIT VECTOR AND MISFIT DO 60 I = 1,ND 60 DWK1(I) = (D(I) - DM(I))/SD(I) TOL0 = SQRT(ANORM(ND,DWK1)/ND) C WEIGHT THE JACOBIAN MATRIX DO 70 I = 1,ND DO 70 J = 1,NP 70 WJ(I,J) = WJ(I,J)/SD(I) C FORM W.J.TRANS(W.J) CALL ATAMUL(NDD,ND,NPP,NP,WJ,WJTWJ) C FORM THE WEIGHTED, TRANSLATED DATA AND PREMULTIPLY BY TRANS(W.J) CALL MULT(NDD,ND,NPP,NP,1,WJ,PM,DHAT) DO 80 I = 1,ND 80 DHAT(I) = DWK1(I) + DHAT(I) CALL TRMULT(NDD,ND,NPP,NP,1,WJ,DHAT,WJTWD) RETURN END C C C----------------------------------------------------------------------- SUBROUTINE MAKPTP(IRUF) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.00, 13 May 1992 C C CONSTRUCTS PENALTY MATRIX (DEL) MULTIPLIED BY ITS TRANSPOSE. C THE PENALTY MATRIX ONLY NEEDS TO BE AROUND LONG ENOUGH TO DO C THIS, SO BY STUFFING THIS INTO A SUBROUTINE DYNAMIC MEMORY ALLOCATION C CAN GET RID OF IT AFTERWARDS C C ON INPUT: C IRUF = 0 FOR NO ROUGHNENING (I.E. MIN. NORM MODELS) (ADDED IN V1.5) C = 1 FOR 1ST DERIVATIVE ROUGHNESS C 2 FOR 2ND DERIVATIVE ROUGHNESS C /MODEL/ NP = NUMBER OF MODEL PARAMETERS C ON OUTPUT: C /RESULT/ PTP() = DEL(TRANS)*DEL C CALLS: C OBJMAT, TRMULT C C INCLUDES: INCLUDE 'imp.inc' INCLUDE 'occamdim.inc' INCLUDE 'occam.inc' C C ARGUMENTS: INTEGER IRUF C LOCAL VARIABLES INTEGER NPEN, I REAL DEL(NNP,NPP) C C----------------------------------------------------------------------- C MAKE PARTIALS MATRIX: CALL OBJMAT(IRUF, DEL, NPEN) IF (abs(idebug) .EQ. 1) * WRITE(*,*) ' NUMBER OF PENALTY TERMS =', NPEN C FORM TRANS(DEL).DEL CALL TRMULT(NNP,NPEN,NPP,NP,NP,DEL,DEL,PTP) C ADD THE WEIGHTS FOR THE MODEL PREJUDICE DO 10 I = 1,NP 10 PTP(I,I) = PTP(I,I) + PREWTS(I)*PREWTS(I) RETURN END C C C----------------------------------------------------------------------- REAL FUNCTION FNDRUF(IRUF) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.00, 13 May 1992 C C COMPUTES THE ROUGHNESS OF THE MODEL. FOR MEMORY ALLOCATION CONSIDERATIONS, C THE PENALTY MATRIX HAS NOT BEEN KEPT AROUND, SO WE NEED TO REFORM IT. C C ON INPUT: C IRUF = 0 FOR NO ROUGHNENING (I.E. MIN. NORM MODELS) (ADDED IN V1.5) C = 1 FOR 1ST DERIVATIVE ROUGHNESS C 2 FOR 2ND DERIVATIVE ROUGHNESS C /RESULT/ PWK2() = LOCATION OF THE MODEL WE WANT TO MEASURE C /MODEL/ NP = NUMBER OF MODEL PARAMETERS C ON OUTPUT: C FNDRUF = ROUGHNESS OF MODEL C CALLS: C OBJMAT, TRMULT C C INCLUDES: INCLUDE 'imp.inc' INCLUDE 'occamdim.inc' INCLUDE 'occam.inc' C C ARGUMENTS: INTEGER IRUF C LOCAL VARIABLES INTEGER I, K, NPEN REAL DEL(NNP,NPP), SUM, VECT(NNP) C FUNCTIONS REAL ANORM C C----------------------------------------------------------------------- C MAKE PARTIALS MATRIX: CALL OBJMAT(IRUF, DEL, NPEN) DO 170 I = 1,NPEN SUM = 0.0 DO 165 K = 1,NP 165 SUM = SUM + DEL(I,K)*PWK2(K) 170 VECT(I) = SUM FNDRUF = ANORM(NPEN,VECT) RETURN END C C C----------------------------------------------------------------------- SUBROUTINE ATAMUL(MAD,MA,NAD,NA,A,C) C----------------------------------------------------------------------- INCLUDE 'imp.inc' C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 C WRITTEN BY: C. DEGROOT-HEDLIN C C DESCRIPTION: C C MULTIPLIES THE TRANSPOSE OF A 2D MATRIX BY ITSELF C TAKES ADVANTAGE OF SYMMETRY ABOUT DIAGONAL C C C ARGUMENTS INTEGER MAD,MA,NAD,NA REAL A(MAD,NA),C(NAD,NA) C C LOCAL VARIABLES INTEGER I,J,K REAL CIJ C C----------------------------------------------------------------------- DO 12 I = 1,NA DO 12 J = I,NA CIJ = 0.0 DO 10 K = 1,MA 10 CIJ = A(K,I)*A(K,J) + CIJ 12 C(I,J) = CIJ DO 15 I=2,NA DO 15 J=1,I-1 15 C(I,J) = C(J,I) RETURN END C C C----------------------------------------------------------------------- subroutine itrout (itform,descr,modfil,datfil,maxitr,tolreq, * iruf,nit,pmu,rlast,tobt,ifftol,iof2) C----------------------------------------------------------------------- c C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.00, 13 May 1992 C c writes an iteration file for OCCAM c c includes: include 'imp.inc' include 'occamdim.inc' include 'occam.inc' c uses np, pm(), idebug c arguments: integer maxitr, iruf, nit, ifftol,iof2 real tolreq, pmu, rlast, tobt character*80 modfil, datfil, itform, descr, datetm c local variables integer lerr, i c lerr = error flag on file operation c i = loop index character*2 num c C----------------------------------------------------------------------- c create file name suffix write (num, '(i2)') nit if (nit .le. 9) num(1:1) = '0' c open iteration file open (unit=iof2, file='ITER'//num, iostat=lerr) if (lerr .ne. 0) then c we could be in the middle of an iteration, so don't die on error, but c send a message and forge on regardless write(*,*) ' Error opening iteration file' return end if call datime(datetm) c c write stuff write(unit=iof2,fmt='(a)') 'FORMAT: '//itform write(unit=iof2,fmt='(a)') 'DESCRIPTION: '//descr write(unit=iof2,fmt='(a)') 'MODEL FILE: '//modfil write(unit=iof2,fmt='(a)') 'DATA FILE: '//datfil write(unit=iof2,fmt='(a)') 'DATE/TIME: '//datetm write(unit=iof2,fmt='(a,i3)') 'MAX ITER: ',maxitr write(unit=iof2,fmt='(a,f6.3)') 'REQ TOL: ',tolreq write(unit=iof2,fmt='(a,i3)') 'IRUF: ',iruf write(unit=iof2,fmt='(a,i3)') 'DEBUG LEVEL: ',idebug write(unit=iof2,fmt='(a,i3)') 'ITERATION: ',nit write(unit=iof2,fmt='(a,g15.7)') 'PMU: ',pmu write(unit=iof2,fmt='(a,g15.7)') 'RLAST: ',rlast write(unit=iof2,fmt='(a,g15.7)') 'TLAST: ',tobt write(unit=iof2,fmt='(a,i3)') 'IFFTOL: ',ifftol write(unit=iof2,fmt='(a,i4)') 'NO. PARMS: ',np write(unit=iof2,fmt='(4g15.7)') (pm(i), i=1,np) C close(unit=iof2) C return end c C C----------------------------------------------------------------------- subroutine itrin (itform,descr,modfil,datfil,maxitr,tolreq, * iruf,nit,pmu,rlast,tobt,ifftol,iof2) C----------------------------------------------------------------------- c C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.01, 13 Jan 1993 C c reads a STARTUP or ITER file for OCCAM c c includes: include 'imp.inc' include 'occamdim.inc' c uses npp include 'occam.inc' c uses np, pm(), idebug c arguments: integer maxitr, iruf, nit, ifftol,iof2 real tolreq, pmu, rlast, tobt character*80 modfil, datfil, itform, descr c local variables integer i, idcode c i = loop index c idcode = function to perform the equivalent of a free-format internal c read for an integer real rdcode c rdcode = function to perform the equivalent of a free-format internal c read for a real character*80 string C C----------------------------------------------------------------------- read(unit=iof2,fmt='(18x,a)', end=198, err=199) itform read(unit=iof2,fmt='(18x,a)', end=198, err=199) descr read(unit=iof2,fmt='(18x,a)', end=198, err=199) modfil read(unit=iof2,fmt='(18x,a)', end=198, err=199) datfil read(unit=iof2,fmt='(18x,a)', end=198, err=199) string c read(unit=iof2,fmt='(18x,a)', end=198, err=199) string maxitr = idcode(string,iof2) read(unit=iof2,fmt='(18x,a)', end=198, err=199) string tolreq = rdcode(string,iof2) read(unit=iof2,fmt='(18x,a)', end=198, err=199) string iruf = idcode(string,iof2) read(unit=iof2,fmt='(18x,a)', end=198, err=199) string idebug = idcode(string,iof2) read(unit=iof2,fmt='(18x,a)', end=198, err=199) string nit = idcode(string,iof2) read(unit=iof2,fmt='(18x,a)', end=198, err=199) string pmu = rdcode(string,iof2) read(unit=iof2,fmt='(18x,a)', end=198, err=199) string rlast = rdcode(string,iof2) read(unit=iof2,fmt='(18x,a)', end=198, err=199) string tobt = rdcode(string,iof2) read(unit=iof2,fmt='(18x,a)', end=198, err=199) string ifftol = idcode(string,iof2) read(unit=iof2,fmt='(18x,a)', end=198, err=199) string np = idcode(string,iof2) if (np .gt. npp) call filerr (' Too many parameters', iof2) read (unit=iof2,fmt= *, end=198, err=199) (pm(i), i=1,np) c return c 198 call filerr (' Startup file ended prematurely', iof2) 199 call filerr (' Error reading startup file', iof2) end c C C----------------------------------------------------------------------- real function rdcode(string, iof2) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.00, 13 Jan 1993 C c rdcode provides a patch to get around the lack of an internal c free-format (list-directed) read on some compilers, notably the c cray. This version for reals. Butchered from decode by Bob Parker. c c on input: c string = character string containing number to be read c iof2 = unit number of open file, to be closed on error c on output: c rdcode contains a real number from string c calls: c filerr c c includes: include 'imp.inc' c input arguments: character*(*) string integer iof2 c local variables character*40 local integer k, k1, kn, l c c Terminate line at % sign kn=index(string//'%', '%') - 1 k1=1 do 1100 k=k1, kn if (string(k:k) .ne. ' ') goto 1200 1100 continue call filerr (' No number found by rdcode', iof2) 1200 do 1300 l=k, kn if (string(l:l).eq. ',' .or. string(l:l) .eq. ' ') goto 1500 1300 continue 1500 local=string(k:l-1) read (local, '(f40.0)', err=1900) rdcode return 1900 call filerr (' Error reading string in rdcode', iof2) end C----------------------------------------------------------------------- integer function idcode(string, iof2) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.00, 13 Jan 1993 C c idcode provides a patch to get around the lack of an internal c free-format (list-directed) read on some compilers, notably the c cray. This version for integers. Butchered from decode by Bob Parker. c c on input: c string = character string containing number to be read c iof2 = unit number of open file, to be closed on error c on output: c idcode contains an integer number from string c calls: c filerr c c includes: include 'imp.inc' c input arguments: character*(*) string integer iof2 c local variables character*40 local integer k, k1, kn, l real realno c c Terminate line at % sign kn=index(string//'%', '%') - 1 k1=1 do 1100 k=k1, kn if (string(k:k) .ne. ' ') goto 1200 1100 continue call filerr (' No number found by idcode', iof2) 1200 do 1300 l=k, kn if (string(l:l).eq. ',' .or. string(l:l) .eq. ' ') goto 1500 1300 continue 1500 local=string(k:l-1) read (local, '(f40.0)', err=1900) realno idcode = int(realno) return 1900 call filerr (' Error reading string in idcode', iof2) end C----------------------------------------------------------------------- subroutine filerr(mssg, io1) C----------------------------------------------------------------------- C C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 2.00, 13 May 1992 C c filerr prints an error message, closes file io1, and stops c c on input: c mssg = character string containing error message c i01 = unit number of file to close (0 = no file open) c on output: c outputs nothing c calls: c no other routines c c c includes: include 'imp.inc' c input arguments: character*(*) mssg integer io1 c C----------------------------------------------------------------------- write(*,*) mssg if (io1 .gt. 0) close (io1) stop end C----------------------------------------------------------------------- subroutine datime(datetm) c----------------------------------------------------------------------- c OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225 c Subroutine Revision 1.00, 18 May 1992 c c Dummy date and time utility for OCCAM 2.0. c c on input: c nothing c on output: c datetm = 80 character string containing message saying that date and time c are not available c character*80 datetm c c----------------------------------------------------------------------- datetm = 'No date and time available on this system' return end c-----------------------------------------------------------------------