C++++++++++++++++ DOUBLE PRECISION VERSION 3DP (JUN 1988) ++++++++++++++MAIN0001 C SPLMOD. MAIN SUBPROGRAM. MAIN0002 C FOR THE ANALYSIS OF SUMS OF ONE-PARAMETER FUNCTIONS. MAIN0003 C APPROXIMATES OPTIONALLY THE MODEL WITH CUBIC B-SPLINES MAIN0004 C (SPLINE=.TRUE.) OR USES EXACT MODEL (SPLINE=.FALSE.). MAIN0005 C MAIN0006 C Author: Robert H. Vogel C C REFERENCES - R.H.VOGEL (1988) SPLMOD USERS MANUAL, EMBL TECHNICAL MAIN0007 C REPORT DA09, (EUROPEAN MOLECULAR BIOLOGY MAIN0008 C LABORATORY, HEIDELBERG, F.R. OF GERMANY) MAIN0009 C S.W.PROVENCHER, R.H.VOGEL (1983) IN PROGRESS IN MAIN0010 C SCIENTIFIC COMPUTING, VOL. 2, PP. 304-319, MAIN0011 C P.DEUFLHARD AND E.HAIRER EDS. (BIRKHAEUSER, BOSTON) MAIN0012 C MAIN0013 C DISTRIBUTION AND USE OF THIS PROGRAM IS UNRESTRICTED AS LONG AS MAIN0014 C NO CHARGE IS MADE AND PROPER REFERENCE GIVEN. HOWEVER, AS A NEW MAIN0015 C USER YOU SHOULD NOTIFY Stephen W. Provencher at: C sp@S-provencher.COM C http://S-provencher.COM C C SO THAT YOU CAN OBTAIN THE USERS MANUAL (WHICH IS ESSENTIAL) C AND HAVE YOUR NAME PUT ON A MAILING LIST FOR UPDATES, ETC. C-----------------------------------------------------------------------MAIN0026 C SUBPROGRAMS CALLED - ANALYZ, CHOTRD, ERRMES, INICAP, INIT, MAIN0027 C INPUT, SETWT, USERSI, USERTR, WRITYT MAIN0028 C WHICH IN TURN CALL - BCOEFF, BETAIN, BSPLIN, CHOSOL, DIFF, MAIN0029 C ERRMES, FISHNI, FULLSM, GAMLN, GETPRU, MAIN0030 C HESEXA, HESSEN, HESSPL, INTERP, INTERV, MAIN0031 C LINEQS, LINPAR, LSTSQR, NOKSUB, NPASCL, MAIN0032 C OUTALL, OUTCOR, OUTITR, OUTPAR, PGAUSS, MAIN0033 C PIVOT, PIVOT1, PLPRIN, PLRES, RANDOM, MAIN0034 C READYT, RGAUSS, STORIN, USEREX, USERFL, MAIN0035 C USERFN, USERIN, USEROU, USERST, USERTR, MAIN0036 C USERWT, VARIAN, WRITIN, WRITYT MAIN0037 C COMMON USED - DBLOCK, SBLOCK, IBLOCK, LBLOCK MAIN0038 C-----------------------------------------------------------------------MAIN0039 DOUBLE PRECISION PRECIS, RANGE, MAIN0040 1 DZ, DZINV, VARMIN, VAROLD, VARSAV, ZSTART, MAIN0041 2 ZERO, ONE, TWO, THREE, FOUR, SIX, TEN, MAIN0042 3 TWOTHR, FORTHR, SIXTEL MAIN0043 DOUBLE PRECISION EK, EPK, ETE, R, DELTAL, DELTAN, GTG, H, MAIN0044 1 DELTAH, EKEK, ASAVE, ATOT, DTOLER, DX, BSPL0, MAIN0045 2 BSPL1, BSPL2, CCAP, ECAP, YCAP, SGGCAP, MAIN0046 3 SGYCAP, YW, BETA, ZGRID, D, U, WORK, YKN, MAIN0047 4 SWKN MAIN0048 LOGICAL DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, DOADEX, MAIN0049 1 DOSPL, DOSTRT, DOUSOU, LUSER, MAIN0050 2 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, MAIN0051 3 SPLINE, WRTBET, WRTCAP MAIN0052 LOGICAL PIVLIN, PIVLAM, STARTV MAIN0053 DIMENSION BOUNDL(2), BOUNDN(2), IHOLER(6) MAIN0054 C***********************************************************************MAIN0055 C THE INSTRUCTIONS SET OFF BY ASTERISKS DESCRIBE ALL POSSIBLE CHANGES MAIN0056 C THAT YOU HAVE TO MAKE IN THIS SUBPROGRAM. (SEE ALSO THE CHANGES MAIN0057 C IN THE BLOCK DATA AND USER SUBPROGRAMS.) THESE CHANGES IN THE MAIN MAIN0058 C PROGRAM ARE ONLY NECESSARY IF YOU CHANGE MAIN0059 C NYMAX, NDMAX, NONLMX, NGAMMX, NLINMX, NLNDMX, NZMAX, OR MTRYMX MAIN0060 C IN THE DATA STATEMENTS BELOW. IF YOU DO, THEN THE MAIN0061 C FOLLOWING DIMENSIONS MUST BE READJUSTED AS DESCRIBED BELOW - MAIN0062 C MAIN0063 DIMENSION Y(1024,3), T(1024,3), YLYFIT(1024,3), MAIN0064 1 SQRTW(1024,3), RESPON(1024) MAIN0065 DIMENSION PLAM(5), PLMTRY(5), EPK(5), DELTAN(5), MAIN0066 1 NOKSET(5), NCOMBI(5), DX(5), IC(5) MAIN0067 DIMENSION BSPL0(4,5), BSPL1(4,5), BSPL2(4,5) MAIN0068 DIMENSION PLMSAV(5,10), DELTAH(5,3) MAIN0069 DIMENSION GTG(5,5) MAIN0070 DIMENSION EK(9), DELTAL(9) MAIN0071 DIMENSION PLIN(9,3), PLNTRY(9,3), PIVLIN(9,3), R(9,3) MAIN0072 DIMENSION ETE(9,9) MAIN0073 DIMENSION EKEK(9,9,3), ASAVE(9,9,3), H(9,9,3) MAIN0074 DIMENSION NY(3), YW(3), ERRFIT(3) MAIN0075 DIMENSION PIVLAM(32), DTOLER(32), BUFF(32) MAIN0076 DIMENSION ATOT(32,32) MAIN0077 DIMENSION ZGRID(52), D(52), U(52), WORK(52), BETA(52) MAIN0078 DIMENSION CCAP(52,52,3), ECAP(52,4,3), YCAP(52,3), MAIN0079 1 SGGCAP(4,4,3), SGYCAP(4,3) MAIN0080 DIMENSION STARTV(50) MAIN0081 C MAIN0082 C THE ABOVE DIMENSION STATEMENTS MUST BE ADJUSTED AS FOLLOWS - MAIN0083 C MAIN0084 C DIMENSION Y(NYMAX,NDMAX), T(NYMAX,NDMAX), MAIN0085 C YLYFIT(NYMAX,NDMAX), SQRTW(NYMAX,NDMAX), MAIN0086 C MAIN0087 C RESPON(NYMAX) MAIN0088 C MAIN0089 C YOU CAN SAVE STORAGE BY DIMENSIONING RESPON AS RESPON(1) MAIN0090 C IF YOU NEVER WORK WITH CONVOLUTIONS (IUSER(9).NE.2) MAIN0091 C MAIN0092 C DIMENSION PLAM(NONLMX), PLMTRY(NONLMX), EPK(NONLMX), MAIN0093 C DELTAN(NONLMX), NOKSET(NONLMX), NCOMBI(NONLMX), MAIN0094 C DX(NONLMX), IC(NONLMX) MAIN0095 C DIMENSION BSPL0(IB,NONLMX), BSPL1(IB,NONLMX), BSPL2(IB,NONLMX) MAIN0096 C DIMENSION PLMSAV(NONLMX,10), DELTAH(NONLMX,NDMAX) MAIN0097 C DIMENSION GTG(NONLMX,NONLMX) MAIN0098 C DIMENSION EK(NLINMX), DELTAL(NLINMX) MAIN0099 C DIMENSION PLIN(NLINMX,NDMAX), PLNTRY(NLINMX,NDMAX), MAIN0100 C PIVLIN(NLINMX,NDMAX), R(NLINMX,NDMAX) MAIN0101 C DIMENSION ETE(NLINMX,NLINMX) MAIN0102 C DIMENSION EKEK(NLINMX,NLINMX,NDMAX), ASAVE(NLINMX,NLINMX,NDMAX) MAIN0103 C H(NLINMX,NLINMX,NDMAX) MAIN0104 C DIMENSION NY(NDMAX), YW(NDMAX), ERRFIT(NDMAX) MAIN0105 C DIMENSION PIVLAM(NLNDMX), DTOLER(NLNDMX), BUFF(NLNDMX) MAIN0106 C DIMENSION ATOT(NLNDMX,NLNDMX)) MAIN0107 C DIMENSION ZGRID(NZMAX), D(NZMAX), U(NZMAX), WORK(NZMAX), MAIN0108 C BETA(NZMAX) MAIN0109 C DIMENSION CCAP(NZMAX,NZMAX,NDMAX), ECAP(NZMAX ,NGAMMX,NDMAX), MAIN0110 C YCAP(NZMAX,NDMAX), SGGCAP(NGAMMX,NGAMMX,NDMAX), MAIN0111 C SGYCAP(NGAMMX,NDMAX) MAIN0112 C DIMENSION STARTV(MTRYMX) MAIN0113 C MAIN0114 C***********************************************************************MAIN0115 COMMON /DBLOCK/ PRECIS, RANGE, DZ, DZINV, VARMIN, VAROLD, MAIN0116 1 VARSAV(10), ZSTART, ZERO, ONE, TWO, THREE, FOUR, MAIN0117 2 SIX, TEN, TWOTHR, FORTHR, SIXTEL MAIN0118 COMMON /SBLOCK/ SRANGE, CONVRG, PLMNMX(2), PNMNMX(2), RUSER(100),MAIN0119 1 DBOX, EXMAX, PNG(10), ZTOTAL MAIN0120 COMMON /IBLOCK/ IRWCAP, IUNIT, IWT, LINEPG, MCONV, METHOD, MAIN0121 1 MIOERR, NABORT, ND, NERFIT, NGAM, NNL, NZ, MAIN0122 2 IFORMT(70), IFORMW(70), IFORMY(70), IPLFIT(2), MAIN0123 3 IPLRES(2), IPRINT(3,2), IPRITR(2), IUSER(50), MAIN0124 4 MTRY(10,2), MXITER(2), NL(10), MAIN0125 5 IALPHA(6), IGAMMA(6), ILAMDA(6), IB, INDEX(10), MAIN0126 6 ISTTXT(28), ITITLE(80), MBOX, MMTRY, NB, NGAMM1, MAIN0127 7 NGAMP1, NLIN, NONL, NONLP1, NPLINE, NYSUM, MAIN0128 8 NIN, NOUT MAIN0129 COMMON /LBLOCK/ DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, MAIN0130 1 DOADEX(10,2), DOSPL(10,2), DOSTRT(10,2), MAIN0131 2 DOUSOU(2), LUSER(30), MAIN0132 3 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, MAIN0133 4 SPLINE, WRTBET, WRTCAP MAIN0134 C***********************************************************************MAIN0135 C YOU CAN SAVE STORAGE BY MAKING THE INTEGERS IN THE FOLLOWING DATA MAIN0136 C STATEMENTS AS SMALL AS THE SIZE OF YOUR PROBLEM WILL ALLOW. (SEE MAIN0137 C SEC 3.5 IN THE USERS MANUAL FOR THE MINIMUM ALLOWABLE VALUES.) MAIN0138 C MAIN0139 DATA NYMAX /1024/, NDMAX /3/, NONLMX /5/, NGAMMX/4/ MAIN0140 C MAIN0141 C NLINMX = NONLMX + NGAMMX MAIN0142 C NLNDMX = NDMAX * NLINMX + NONLMX MAIN0143 C MAIN0144 DATA NLINMX/9/, NLNDMX /32/ MAIN0145 C MAIN0146 DATA NZMAX /52/, MTRYMX/50/ MAIN0147 C MAIN0148 C IF YOU CHANGE THE ABOVE DATA STATEMENTS, THE DIMENSION MAIN0149 C STATEMENTS ABOVE MUST BE READJUSTED, AS DESCRIBED ABOVE. MAIN0150 C MAIN0151 C THIS IS THE END OF ALL POSSIBLE CHANGES THAT YOU MIGHT HAVE TO MAKE MAIN0152 C IN THE MAIN PROGRAM, MAIN0153 C EXECPT THAT IF YOUR SYSTEM DOES NOT AUTOMATICALLY OPEN INPUT AND MAIN0154 C OUTPUT FILES FOR YOU, THEN YOU MIGHT HAVE TO OPEN THEM HERE AND MAIN0155 C GIVE THEM THE NUMBERS NIN (FOR THE INPUT) AND NOUT (FOR THE MAIN0156 C OUTPUT), WHERE NIN AND NOUT HAVE BEEN SET IN THE BLOCK DATA MAIN0157 C SUBPROGRAM. MAIN0158 C IN ADDITION, IF YOU ARE GOING TO INPUT IRWCAP.NE.0, THEN YOU MAY MAIN0159 C HAVE TO OPEN A TEMPORARY SCRATCH FILE NUMBERED IUNIT. DO THIS MAIN0160 C DIRECTLY AFTER STATEMENT 100. THIS IS NOT NECESSARY IF IRWCAP IS MAIN0161 C ZERO OR IF YOUR SYSTEM OPENS FILES AUTOMATICALLY. MAIN0162 C MAIN0163 C***********************************************************************MAIN0164 DATA IHOLER/1HM, 1HA, 1HI, 1HN, 1H , 1H / MAIN0165 C-----------------------------------------------------------------------MAIN0166 C OPEN FILES MAIN0167 C-----------------------------------------------------------------------MAIN0168 C MAIN0169 C MAIN0170 C-----------------------------------------------------------------------MAIN0171 C INITIALIZE CERTAIN VARIABLES MAIN0172 C-----------------------------------------------------------------------MAIN0173 CALL INIT MAIN0174 C-----------------------------------------------------------------------MAIN0175 C READ INPUT DATA MAIN0176 C-----------------------------------------------------------------------MAIN0177 100 CALL INPUT (NYMAX,Y,T,YLYFIT,SQRTW,NY,NONLMX,NLINMX,NGAMMX,NZMAX, MAIN0178 1 NDMAX,NLNDMX,MTRYMX,RESPON,BOUNDL,BOUNDN) MAIN0179 C-----------------------------------------------------------------------MAIN0180 C OPEN FILE FOR INPUT AND OUTPUT OF B-SPLINE COEFFICIENTS AND MAIN0181 C ARRAYS ...CAP IF ONE OF THE LOGICALS REDBET, REDCAP, WRTBET, MAIN0182 C OR WRTCAP IS .TRUE. (I.E. IRWCAP.NE.0) MAIN0183 C-----------------------------------------------------------------------MAIN0184 C IF (IRWCAP .NE. 0) OPEN.......... MAIN0185 C MAIN0186 C-----------------------------------------------------------------------MAIN0187 C CALCULATE SIMULATED DATA MAIN0188 C-----------------------------------------------------------------------MAIN0189 IF (SIMULA) CALL USERSI (NYMAX,Y,T,YLYFIT,NY,RESPON) MAIN0190 IF (SIMULA .AND. PRY) CALL WRITYT ( MAIN0191 1 NYMAX,Y,T,YLYFIT,SQRTW,NY,ND,IWT,SIMULA,NOUT) MAIN0192 C-----------------------------------------------------------------------MAIN0193 C INITIAL COMPUTATIONS. MAIN0194 C-----------------------------------------------------------------------MAIN0195 NGAMM1=NGAM-1 MAIN0196 NGAMP1=NGAM+1 MAIN0197 IF (NGAM .GT. 1) NPLINE=MIN0(8,NGAM) MAIN0198 IWTIST=.TRUE. MAIN0199 SPLINE=.FALSE. MAIN0200 IWTSAV=IWT MAIN0201 I=1 MAIN0202 IF (IWT.EQ.1 .OR. IWT.EQ.4) I=2 MAIN0203 DO 120 J=1,10 MAIN0204 IF (.NOT. DOSPL(J,I)) GO TO 120 MAIN0205 SPLINE=.TRUE. MAIN0206 GO TO 130 MAIN0207 120 CONTINUE MAIN0208 C-----------------------------------------------------------------------MAIN0209 C COMPUTE ACTUAL BOUNDS FOR LINEAR AND NONLINEAR PARAMETERS. MAIN0210 C-----------------------------------------------------------------------MAIN0211 130 AYMAX=0. MAIN0212 DO 134 N=1,ND MAIN0213 MY=NY(N) MAIN0214 DO 132 K=1,MY MAIN0215 132 AYMAX=AMAX1(AYMAX,ABS(Y(K,N))) MAIN0216 134 CONTINUE MAIN0217 PLMNMX(1)=PLMNMX(1)*AYMAX MAIN0218 PLMNMX(2)=PLMNMX(2)*AYMAX MAIN0219 T1 =SRANGE MAIN0220 TN =-SRANGE MAIN0221 DO 140 N=1,ND MAIN0222 MY =NY(N) MAIN0223 DDUM=-SRANGE MAIN0224 DO 138 J=1,5 MAIN0225 DUM=SRANGE MAIN0226 DO 136 K=1,MY MAIN0227 TKN=T(K,N) MAIN0228 IF (J .EQ. 1) TN=AMAX1(TN,TKN) MAIN0229 IF (TKN.GT.DDUM .AND. TKN.LT.DUM) DUM=TKN MAIN0230 136 CONTINUE MAIN0231 DDUM=DUM+1.E-5*ABS(DUM) MAIN0232 IF (J .EQ. 1) TT1=DUM MAIN0233 138 CONTINUE MAIN0234 IF (T1 .LE. TT1) GO TO 140 MAIN0235 T1=TT1 MAIN0236 T5=DUM MAIN0237 140 CONTINUE MAIN0238 DUM =.25*(T5-T1) MAIN0239 DUM =DUM-T1 MAIN0240 PNMNMX(1)=PNMNMX(1)/(TN+DUM) MAIN0241 PNMNMX(2)=PNMNMX(2)/(T1+DUM) MAIN0242 IF (NGAM .EQ. 0) WRITE (NOUT,1000) MAIN0243 1 PLMNMX(1),PLMNMX(2),PNMNMX(1),PNMNMX(2) MAIN0244 IF (NGAM .GT. 0) WRITE (NOUT,2000) MAIN0245 1 PLMNMX(1),PLMNMX(2),PNMNMX(1),PNMNMX(2) MAIN0246 IF (PLMNMX(1).GT.PLMNMX(2) .OR. PNMNMX(1).GT.PNMNMX(2)) MAIN0247 1 CALL ERRMES (1,.TRUE.,IHOLER,NOUT) MAIN0248 C-----------------------------------------------------------------------MAIN0249 C DO ALL EVALUATIONS ON USERTR(PLAM) AXIS MAIN0250 C-----------------------------------------------------------------------MAIN0251 DUM =PNMNMX(1) MAIN0252 PNMNMX(1)=USERTR(DUM,1) MAIN0253 DUM =PNMNMX(2) MAIN0254 PNMNMX(2)=USERTR(DUM,1) MAIN0255 ZTOTAL =PNMNMX(2)-PNMNMX(1) MAIN0256 IF (.NOT. SPLINE) GO TO 200 MAIN0257 C-----------------------------------------------------------------------MAIN0258 C BUILD EQUIDISTANT INTERPOLATION GRID MAIN0259 C-----------------------------------------------------------------------MAIN0260 DZ =ZTOTAL/FLOAT(NZ-3) MAIN0261 DZINV =ONE/DZ MAIN0262 ZSTART =PNMNMX(1)-DZ MAIN0263 ZGRID(1)=ZSTART MAIN0264 DO 150 J=2,NZ MAIN0265 150 ZGRID(J)=ZGRID(J-1)+DZ MAIN0266 C-----------------------------------------------------------------------MAIN0267 C DECOMPOSE SPECIAL TRIDIAGONAL MATRIX FOR SOLVING A SYSTEM OF MAIN0268 C LINEAR EQUATIONS FOR B-SPLINES COEFFICIENTS. MAIN0269 C-----------------------------------------------------------------------MAIN0270 CALL CHOTRD (NZMAX,D,U,NZ,ONE,TWO,FOUR) MAIN0271 IF (REDBET .OR. REDCAP .OR. WRTBET .OR. WRTCAP) REWIND IUNIT MAIN0272 200 IF (IWT.EQ.1 .OR. IWT.EQ.4) GO TO 400 MAIN0273 C-----------------------------------------------------------------------MAIN0274 C INITIALIZE ACCUMULATION ARRAYS ...CAP MAIN0275 C-----------------------------------------------------------------------MAIN0276 IF (SPLINE) CALL INICAP ( MAIN0277 1 NYMAX,Y,T,SQRTW,NY,NZMAX,BETA,ZGRID,D,U,NGAMMX,CCAP,ECAP, MAIN0278 2 YCAP,SGGCAP,SGYCAP,YW,WORK,ZERO,DZ,RESPON) MAIN0279 IF (WRTBET) REDBET=WRTBET MAIN0280 WRTBET=.FALSE. MAIN0281 IF (WRTCAP) REDCAP=WRTCAP MAIN0282 WRTCAP=.FALSE. MAIN0283 VARMIN=ZERO MAIN0284 DO 220 N=1,ND MAIN0285 MY=NY(N) MAIN0286 DO 220 K=1,MY MAIN0287 YKN =Y(K,N) MAIN0288 VARMIN=VARMIN+YKN*YKN MAIN0289 220 CONTINUE MAIN0290 VARMIN=VARMIN*CONVRG*CONVRG MAIN0291 C-----------------------------------------------------------------------MAIN0292 C DO A COMPLTE PRELIMINARY UNWEIGHTED ANALYSIS TO GET A SMOOTH FIT MAIN0293 C TO THE DATA. THIS SMOOTH CURVE IS THEN USED TO CALCULATE THE MAIN0294 C WEIGHTS. MAIN0295 C-----------------------------------------------------------------------MAIN0296 CALL ANALYZ (1, MAIN0297 1 NYMAX,Y,T,YLYFIT,SQRTW,NY,NLINMX,NONLMX,PLIN, MAIN0298 2 PLNTRY,PLAM,PLMTRY,PLMSAV,PIVLIN,PIVLAM,EK,EPK, MAIN0299 3 ETE,R,DELTAL,DELTAN,GTG,H,DELTAH,EKEK,ASAVE, MAIN0300 4 NOKSET,NCOMBI,MTRYMX,STARTV,NLNDMX,ATOT,DTOLER, MAIN0301 5 BUFF,DX,BSPL0,BSPL1,BSPL2,IC,NZMAX,NGAMMX,CCAP, MAIN0302 6 ECAP,YCAP,SGGCAP,SGYCAP,YW,RESPON) MAIN0303 IF (IWT .EQ. 1) GO TO 420 MAIN0304 IWTIST=.FALSE. MAIN0305 C-----------------------------------------------------------------------MAIN0306 C CALCULATE SQRTW (SQUARE ROOT OF LEAST SQUARES WEIGHTS) MAIN0307 C-----------------------------------------------------------------------MAIN0308 CALL SETWT (NYMAX,Y,T,YLYFIT,SQRTW,NY,ND,ERRFIT,NERFIT,IWT, MAIN0309 1 PRWT,SRANGE,NOUT,RESPON) MAIN0310 IF (REDBET) REWIND IUNIT MAIN0311 SPLINE=.FALSE. MAIN0312 DO 300 J=1,10 MAIN0313 IF (.NOT. DOSPL(J,2)) GO TO 300 MAIN0314 SPLINE=.TRUE. MAIN0315 GO TO 400 MAIN0316 300 CONTINUE MAIN0317 400 IF (IWT .EQ. 4) IWTIST=.FALSE. MAIN0318 C-----------------------------------------------------------------------MAIN0319 C INITIALIZE ACCUMULATION ARRAYS ...CAP WITH WEIGHTS MAIN0320 C-----------------------------------------------------------------------MAIN0321 IF (SPLINE) CALL INICAP ( MAIN0322 1 NYMAX,Y,T,SQRTW,NY,NZMAX,BETA,ZGRID,D,U,NGAMMX,CCAP,ECAP, MAIN0323 2 YCAP,SGGCAP,SGYCAP,YW,WORK,ZERO,DZ,RESPON) MAIN0324 420 VARMIN=ZERO MAIN0325 DO 320 N=1,ND MAIN0326 MY=NY(N) MAIN0327 DO 320 K=1,MY MAIN0328 YKN =Y(K,N) MAIN0329 SWKN =SQRTW(K,N) MAIN0330 VARMIN=VARMIN+(YKN*SWKN)**2 MAIN0331 320 CONTINUE MAIN0332 VARMIN=VARMIN*CONVRG*CONVRG MAIN0333 C-----------------------------------------------------------------------MAIN0334 C DO FINAL WEIGHTED ANALYSIS MAIN0335 C-----------------------------------------------------------------------MAIN0336 CALL ANALYZ (2, MAIN0337 1 NYMAX,Y,T,YLYFIT,SQRTW,NY,NLINMX,NONLMX,PLIN,PLNTRY, MAIN0338 2 PLAM,PLMTRY,PLMSAV,PIVLIN,PIVLAM,EK,EPK,ETE,R,DELTAL,MAIN0339 3 DELTAN,GTG,H,DELTAH,EKEK,ASAVE,NOKSET,NCOMBI,MTRYMX, MAIN0340 4 STARTV,NLNDMX,ATOT,DTOLER,BUFF,DX,BSPL0,BSPL1,BSPL2, MAIN0341 5 IC,NZMAX,NGAMMX,CCAP,ECAP,YCAP,SGGCAP,SGYCAP,YW, MAIN0342 6 RESPON) MAIN0343 IF (LAST) STOP MAIN0344 C-----------------------------------------------------------------------MAIN0345 C RESTORE CHANGED CONTROL VARIABLES MAIN0346 C-----------------------------------------------------------------------MAIN0347 PLMNMX(1)=BOUNDL(1) MAIN0348 PLMNMX(2)=BOUNDL(2) MAIN0349 PNMNMX(1)=BOUNDN(1) MAIN0350 PNMNMX(2)=BOUNDN(2) MAIN0351 IWT =IWTSAV MAIN0352 GO TO 100 MAIN0353 1000 FORMAT(///4X,32HBOUNDS FOR LINEAR PARAMETERS:, MAIN0354 1 5X,1PE12.4,20H .LE. ALPHA .LE.,1PE12.4, MAIN0355 2 /4X,32HBOUNDS FOR NONLINEAR PARAMETERS:, MAIN0356 3 5X,1PE12.4,20H .LE. LAMBDA .LE.,1PE12.4) MAIN0357 2000 FORMAT(///4X,32HBOUNDS FOR LINEAR PARAMETERS:, MAIN0358 1 5X,1PE12.4,26H .LE. ALPHA, GAMMA .LE.,1PE12.4, MAIN0359 2 /4X,32HBOUNDS FOR NONLINEAR PARAMETERS:, MAIN0360 3 5X,1PE12.4,26H .LE. LAMBDA .LE.,1PE12.4) MAIN0361 END MAIN0362 C++++++++++++++++ DOUBLE PRECISION VERSION 3DP (JUN 1988) ++++++++++++++BLCKD001 BLOCK DATA BLCKD002 C-----------------------------------------------------------------------BLCKD003 C INITIALIZES COMMON BLOCKS. BLCKD004 C-----------------------------------------------------------------------BLCKD005 C COMMON USED - DBLOCK, SBLOCK, IBLOCK, LBLOCK BLCKD006 C-----------------------------------------------------------------------BLCKD007 DOUBLE PRECISION PRECIS, RANGE, BLCKD008 1 DZ, DZINV, VARMIN, VAROLD, VARSAV, ZSTART, BLCKD009 2 ZERO, ONE, TWO, THREE, FOUR, SIX, TEN, BLCKD010 3 TWOTHR, FORTHR, SIXTEL BLCKD011 LOGICAL DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, DOADEX, BLCKD012 1 DOSPL, DOSTRT, DOUSOU, LUSER, BLCKD013 2 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, BLCKD014 3 SPLINE, WRTBET, WRTCAP BLCKD015 COMMON /DBLOCK/ PRECIS, RANGE, DZ, DZINV, VARMIN, VAROLD, BLCKD016 1 VARSAV(10), ZSTART, ZERO, ONE, TWO, THREE, FOUR, BLCKD017 2 SIX, TEN, TWOTHR, FORTHR, SIXTEL BLCKD018 COMMON /SBLOCK/ SRANGE, CONVRG, PLMNMX(2), PNMNMX(2), RUSER(100),BLCKD019 1 DBOX, EXMAX, PNG(10), ZTOTAL BLCKD020 COMMON /IBLOCK/ IRWCAP, IUNIT, IWT, LINEPG, MCONV, METHOD, BLCKD021 1 MIOERR, NABORT, ND, NERFIT, NGAM, NNL, NZ, BLCKD022 2 IFORMT(70), IFORMW(70), IFORMY(70), IPLFIT(2), BLCKD023 3 IPLRES(2), IPRINT(3,2), IPRITR(2), IUSER(50), BLCKD024 4 MTRY(10,2), MXITER(2), NL(10), BLCKD025 5 IALPHA(6), IGAMMA(6), ILAMDA(6), IB, INDEX(10), BLCKD026 6 ISTTXT(28), ITITLE(80), MBOX, MMTRY, NB, NGAMM1, BLCKD027 7 NGAMP1, NLIN, NONL, NONLP1, NPLINE, NYSUM, BLCKD028 8 NIN, NOUT BLCKD029 COMMON /LBLOCK/ DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, BLCKD030 1 DOADEX(10,2), DOSPL(10,2), DOSTRT(10,2), BLCKD031 2 DOUSOU(2), LUSER(30), BLCKD032 3 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, BLCKD033 4 SPLINE, WRTBET, WRTCAP BLCKD034 C***********************************************************************BLCKD035 C YOU MUST SET THE FOLLOWING 4 VARIABLES TO VALUES APPROPRIATE FOR BLCKD036 C YOUR COMPUTER. (SEE USERS MANUAL) BLCKD037 C BLCKD038 C DATA RANGE/1.E35/, SRANGE/1.E35/, NIN/5/, NOUT/6/ BLCKD039 DATA RANGE/1.D35/, SRANGE/1.E35/, NIN/5/, NOUT/6/ BLCKD040 C***********************************************************************BLCKD041 C BLCKD042 C BLCKD043 C***********************************************************************BLCKD044 C SPLMOD WILL USE THE VALUES OF THE CONTROL VARIABLES THAT ARE IN THE BLCKD045 C DATA STATEMENTS BELOW UNLESS YOU INPUT DIFFERENT VALUES. TO SAVE BLCKD046 C INPUTTING THESE OFTEN, YOU CAN CHANGE THE VALUES BELOW TO THOSE BLCKD047 C THAT YOU WILL USUALLY USE. BLCKD048 C BLCKD049 C THE FOLLOWING DATA STATEMENTS CONTAIN (IN ALPHABETICAL ORDER) THE BLCKD050 C REAL, INTEGER, AND LOGICAL CONTROL VARIABLES, IN THAT ORDER. BLCKD051 C BLCKD052 DATA CONVRG/5.E-5/, PLMNMX/-1.E5, 1.E5/, BLCKD053 1 PNMNMX/.02, 2.08/, RUSER/100*0./ BLCKD054 DATA IFORMT/1H(,1H5,1HE,1H1,1H5,1H.,1H6,1H),62*1H /, BLCKD055 1 IFORMW/1H(,1H5,1HE,1H1,1H5,1H.,1H6,1H),62*1H /, BLCKD056 2 IFORMY/1H(,1H5,1HE,1H1,1H5,1H.,1H6,1H),62*1H /, BLCKD057 3 IPLFIT/2*0/, IPLRES/0, 1/, IPRINT/6*1/, BLCKD058 4 IPRITR/0, 2/, IRWCAP/0/, IUNIT/0/, IUSER/50*0/, BLCKD059 5 IWT/1/, LINEPG/60/, MCONV/3/, METHOD/2/, BLCKD060 6 MIOERR/5/, MTRY/10, 9*20, 20, 9*50/, BLCKD061 7 MXITER/20, 40/, NABORT/3/, ND/1/, NERFIT/10/, BLCKD062 8 NGAM/0/, NL/1, 2, 3, 4, 5, 6, 7, 8, 9, 10/, BLCKD063 9 NNL/1/, NZ/50/ BLCKD064 DATA DOADEX/20*.FALSE./, DOSPL/20*.TRUE./, BLCKD065 1 DOSTRT/20*.FALSE./, DOUSIN/.FALSE./, BLCKD066 2 DOUSOU/2*.FALSE./, LAST/.TRUE./, BLCKD067 3 LUSER/30*.FALSE./, PRWT/.TRUE./, PRY/.TRUE./, BLCKD068 4 SAMET/.FALSE./, SIMULA/.FALSE./ BLCKD069 C***********************************************************************BLCKD070 C DATA ZERO/0.E0/, ONE/1.E0/, TWO/2.E0/, BLCKD071 C 1 THREE/3.E0/, FOUR/4.E0/, SIX/6.E0/, TEN/1.E1/ BLCKD072 DATA ZERO/0.D0/, ONE/1.D0/, TWO/2.D0/, BLCKD073 1 THREE/3.D0/, FOUR/4.D0/, SIX/6.D0/, TEN/1.D1/ BLCKD074 DATA IALPHA/1H , 1HA, 1HL, 1HP, 1HH, 1HA/, BLCKD075 1 IGAMMA/1H , 1HG, 1HA, 1HM, 1HM, 1HA/, BLCKD076 2 ILAMDA/1HL, 1HA, 1HM, 1HB, 1HD, 1HA/, BLCKD077 3 ISTTXT/1H , 1HI, 1HT, 1HR, 1H , 1H , 1H , 1HV, BLCKD078 4 1HA, 1HR, 1HI, 1HA, 1HN, 1HC, 1HE, 1H , BLCKD079 5 1H , 1H , 1HD, 1HA, 1HM, 1HP, 1HI, 1HN, BLCKD080 6 1HG, 1H , 1HQ, 1H /, BLCKD081 7 INDEX/10*0/, IB/4/ BLCKD082 DATA REDBET/.FALSE./, REDCAP/.FALSE./, BLCKD083 1 WRTBET/.FALSE./, WRTCAP/.FALSE./ BLCKD084 END BLCKD085 C++++++++++++++++ DOUBLE PRECISION VERSION 3DP (JUN 1988) ++++++++++++++USREX001 REAL FUNCTION USEREX (NYMAX,T,NY,K,N,RESPON) USREX002 C-----------------------------------------------------------------------USREX003 C USER SUPPLIED SUBROUTINE. (ONLY CALLED IF SIMULA=.TRUE.) USREX004 C TO EVALUATE THE NOISE-FREE VALUE OF THE SIMULATED DATA AT POINT USREX005 C T(K,N). USREX006 C IN CASE OF CONVOLUTION, RESPON CONTAINS THE RESPONSE FUNCTION. USREX007 C BELOW IS ILLUSTRATED THE CASE OF THE DATA BEING COMPOSED USREX008 C OF A SUM OF IUSER(10) (.LE. 10) EXPONENTIALS PLUS A CONSTANT USREX009 C BACKGROUND, WHERE USREX010 C DECAY CONSTANTS ARE STORED IN RUSER( 9+J), J=1,...,IUSER(10) USREX011 C AND AMPLITUDES IN RUSER(19+J), J=1,...,IUSER(10) USREX012 C AND BACKGROUND IN RUSER(30) USREX013 C IT ALSO SHOWS HOW TO PREVENT UNDERFLOWS IN CALCULATING EXP(-X) USREX014 C WHEN X IS LARGE. USREX015 C-----------------------------------------------------------------------USREX016 C SUBPROGRAMS CALLED - NONE USREX017 C COMMON USED - DBLOCK, SBLOCK, IBLOCK, LBLOCK USREX018 C-----------------------------------------------------------------------USREX019 DOUBLE PRECISION PRECIS, RANGE, USREX020 1 DZ, DZINV, VARMIN, VAROLD, VARSAV, ZSTART, USREX021 2 ZERO, ONE, TWO, THREE, FOUR, SIX, TEN, USREX022 3 TWOTHR, FORTHR, SIXTEL USREX023 LOGICAL DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, DOADEX, USREX024 1 DOSPL, DOSTRT, DOUSOU, LUSER, USREX025 2 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USREX026 3 SPLINE, WRTBET, WRTCAP USREX027 DIMENSION T(NYMAX,1), NY(1), RESPON(1), IHOLER(6) USREX028 COMMON /DBLOCK/ PRECIS, RANGE, DZ, DZINV, VARMIN, VAROLD, USREX029 1 VARSAV(10), ZSTART, ZERO, ONE, TWO, THREE, FOUR, USREX030 2 SIX, TEN, TWOTHR, FORTHR, SIXTEL USREX031 COMMON /SBLOCK/ SRANGE, CONVRG, PLMNMX(2), PNMNMX(2), RUSER(100),USREX032 1 DBOX, EXMAX, PNG(10), ZTOTAL USREX033 COMMON /IBLOCK/ IRWCAP, IUNIT, IWT, LINEPG, MCONV, METHOD, USREX034 1 MIOERR, NABORT, ND, NERFIT, NGAM, NNL, NZ, USREX035 2 IFORMT(70), IFORMW(70), IFORMY(70), IPLFIT(2), USREX036 3 IPLRES(2), IPRINT(3,2), IPRITR(2), IUSER(50), USREX037 4 MTRY(10,2), MXITER(2), NL(10), USREX038 5 IALPHA(6), IGAMMA(6), ILAMDA(6), IB, INDEX(10), USREX039 6 ISTTXT(28), ITITLE(80), MBOX, MMTRY, NB, NGAMM1, USREX040 7 NGAMP1, NLIN, NONL, NONLP1, NPLINE, NYSUM, USREX041 8 NIN, NOUT USREX042 COMMON /LBLOCK/ DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, USREX043 1 DOADEX(10,2), DOSPL(10,2), DOSTRT(10,2), USREX044 2 DOUSOU(2), LUSER(30), USREX045 3 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USREX046 4 SPLINE, WRTBET, WRTCAP USREX047 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HE, 1HX/ USREX048 C-----------------------------------------------------------------------USREX049 C THE FOLLOWING STATEMENTS SHOULD BE REPLACED WITH THE ONES USREX050 C APPROPRIATE FOR YOUR SIMULATION. USREX051 C-----------------------------------------------------------------------USREX052 USEREX=RUSER(30) USREX053 NLAM =IUSER(10) USREX054 DO 100 J=1,NLAM USREX055 EZ =RUSER(9+J)*T(K,N) USREX056 EXPEZ=0. USREX057 IF (EZ .LT. EXMAX) EXPEZ=EXP(-EZ) USREX058 USEREX=USEREX+RUSER(19+J)*EXPEZ USREX059 100 CONTINUE USREX060 RETURN USREX061 END USREX062 C++++++++++++++++ DOUBLE PRECISION VERSION 3DP (JUN 1988) ++++++++++++++USRFL001 REAL FUNCTION USERFL (MU,NYMAX,T,NY,K,N,RESPON) USRFL002 C-----------------------------------------------------------------------USRFL003 C USER SUPPLIED SUBROUTINE. USRFL004 C ONLY CALLED WHEN NGAM>0. EVALUATES MU'TH LINEAR FUNCTION AT USRFL005 C T(K,N), I.E. THE LINEAR TERMS IN THE RIGHT SUM OF EQ (3.1-1) USRFL006 C OR EQ (4.1.2-1) IN THE USERS MANUAL. USRFL007 C USRFL008 C IN CASE OF CONVOLUTION, RESPON CONTAINS THE RESPONSE FUNCTION.USRFL009 C USRFL010 C BELOW IS ILLUSTRATED THE CASE WHEN THE FIRST LINEAR FUNCTION IS A USRFL011 C SIMPLE ADDITIVE CONSTANT, OFTEN CALLED "BASELINE" OR "BACKGROUND".USRFL012 C THE OTHER TERMS WITH MU > 1 ARE USED, WHEN ONE IS WORKING WITHUSRFL013 C A CONVOLUTION AND/OR TO CORRECT FOR A SHIFT OF THE RESPONSE USRFL014 C FUNCTION (SEE ALSO REF [6] IN THE USERS MANUAL). USRFL015 C USRFL016 C IN DETAIL THERE ARE THE FOLLOWING POSSIBILITIES: USRFL017 C USRFL018 C NGAM = 1 IF IUSER(1)=1 USERFN(1,...) ACCOUNTS FOR CONSTANT USRFL019 C BACKGROUND USRFL020 C IF IUSER(1)=3 USERFN(1,...) ACCOUNTS FOR CONSTANT USRFL021 C BACKGROUND USRFL022 C NGAM = 2 IF IUSER(1)=2 USERFN(1,...) ACCOUNTS FOR CONSTANT USRFL023 C USERFN(2,...) USRFL024 C USRFL025 C IUSER(1) AND LUSER(1) DEFINE MODEL SWITCHES IN SUBPROGRAM USRFL026 C USERFN, SEE THERE FOR DETAILS. USRFL027 C USRFL028 C THIS VERSION OF USERFL AND THE PROPER SETTINGS OF IUSER(1), USRFL029 C RUSER(1), AND LUSER(1) ARE DESCRIBED IN DETAIL IN SEC 6.1 USRFL030 C OF THE USERS MANUAL. USRFL031 C-----------------------------------------------------------------------USRFL032 C SUBPROGRAMS CALLED - ERRMES USRFL033 C COMMON USED - DBLOCK, SBLOCK, IBLOCK, LBLOCK USRFL034 C-----------------------------------------------------------------------USRFL035 DOUBLE PRECISION PRECIS, RANGE, USRFL036 1 DZ, DZINV, VARMIN, VAROLD, VARSAV, ZSTART, USRFL037 2 ZERO, ONE, TWO, THREE, FOUR, SIX, TEN, USRFL038 3 TWOTHR, FORTHR, SIXTEL USRFL039 LOGICAL DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, DOADEX, USRFL040 1 DOSPL, DOSTRT, DOUSOU, LUSER, USRFL041 2 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRFL042 3 SPLINE, WRTBET, WRTCAP USRFL043 DIMENSION T(NYMAX,1), NY(1), RESPON(1), IHOLER(6) USRFL044 COMMON /DBLOCK/ PRECIS, RANGE, DZ, DZINV, VARMIN, VAROLD, USRFL045 1 VARSAV(10), ZSTART, ZERO, ONE, TWO, THREE, FOUR, USRFL046 2 SIX, TEN, TWOTHR, FORTHR, SIXTEL USRFL047 COMMON /SBLOCK/ SRANGE, CONVRG, PLMNMX(2), PNMNMX(2), RUSER(100),USRFL048 1 DBOX, EXMAX, PNG(10), ZTOTAL USRFL049 COMMON /IBLOCK/ IRWCAP, IUNIT, IWT, LINEPG, MCONV, METHOD, USRFL050 1 MIOERR, NABORT, ND, NERFIT, NGAM, NNL, NZ, USRFL051 2 IFORMT(70), IFORMW(70), IFORMY(70), IPLFIT(2), USRFL052 3 IPLRES(2), IPRINT(3,2), IPRITR(2), IUSER(50), USRFL053 4 MTRY(10,2), MXITER(2), NL(10), USRFL054 5 IALPHA(6), IGAMMA(6), ILAMDA(6), IB, INDEX(10), USRFL055 6 ISTTXT(28), ITITLE(80), MBOX, MMTRY, NB, NGAMM1, USRFL056 7 NGAMP1, NLIN, NONL, NONLP1, NPLINE, NYSUM, USRFL057 8 NIN, NOUT USRFL058 COMMON /LBLOCK/ DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, USRFL059 1 DOADEX(10,2), DOSPL(10,2), DOSTRT(10,2), USRFL060 2 DOUSOU(2), LUSER(30), USRFL061 3 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRFL062 4 SPLINE, WRTBET, WRTCAP USRFL063 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HF, 1HL/ USRFL064 C-----------------------------------------------------------------------USRFL065 IF (N .GT. ND) CALL ERRMES (1,.TRUE.,IHOLER,NOUT) USRFL066 IF (K.LE.0 .OR. K.GT.NY(N)) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) USRFL067 IF (MU.LE.0 .OR. MU.GT.6) CALL ERRMES (3,.TRUE.,IHOLER,NOUT) USRFL068 C-----------------------------------------------------------------------USRFL069 GO TO (100,200,300,400,500,600), MU USRFL070 C-----------------------------------------------------------------------USRFL071 C THE FOLLOWING CONSTANT TERM ACCOUNTS FOR BACKGROUND USRFL072 C-----------------------------------------------------------------------USRFL073 100 USERFL=1. USRFL074 RETURN USRFL075 C-----------------------------------------------------------------------USRFL076 C 1ST CASE (NGAM=2 OR 3) USRFL077 C-----------------------------------------------------------------------USRFL078 200 IF (LUSER(1) .AND. IUSER(1).NE.2) GOTO 210 USRFL079 USERFL=0. USRFL080 IF (NGAM.NE.2 .AND. NGAM.NE.3) CALL ERRMES (4,.TRUE.,IHOLER,NOUT) USRFL081 USERFL=RESPON(K) USRFL082 RETURN USRFL083 C-----------------------------------------------------------------------USRFL084 C 2ND CASE (NGAM=2 AND SHIFT CORRECTION) USRFL085 C-----------------------------------------------------------------------USRFL086 210 USERFL=0. USRFL087 IF (NGAM .NE. 2) CALL ERRMES (5,.TRUE.,IHOLER,NOUT) USRFL088 IF (K .EQ. NY(1)) GOTO 230 USRFL089 IF (K .EQ. 1) GOTO 220 USRFL090 KP1 =K+1 USRFL091 KM1 =K-1 USRFL092 USERFL=(RESPON(KP1)-RESPON(KM1))/(T(KP1,1)-T(KM1,1)) USRFL093 GOTO 240 USRFL094 220 USERFL=(RESPON(2)-RESPON(1))/(T(2,1)-T(1,1)) USRFL095 GOTO 240 USRFL096 230 NY1 =NY(1) USRFL097 NY1M1 =NY1-1 USRFL098 USERFL=(RESPON(NY1)-RESPON(NY1M1))/(T(NY1,1)-T(NY1M1,1)) USRFL099 240 USERFL=RUSER(1)*RESPON(K)+USERFL USRFL100 RETURN USRFL101 C-----------------------------------------------------------------------USRFL102 C NGAM=3 USRFL103 C-----------------------------------------------------------------------USRFL104 300 USERFL=0. USRFL105 IF (NGAM.NE.3 .OR. LUSER(1)) CALL ERRMES (6,.TRUE.,IHOLER,NOUT) USRFL106 IF (K .EQ. NY(1)) GOTO 320 USRFL107 IF (K .EQ. 1) GOTO 310 USRFL108 KP1 =K+1 USRFL109 KM1 =K-1 USRFL110 USERFL=(RESPON(KP1)-RESPON(KM1))/(T(KP1,1)-T(KM1,1)) USRFL111 RETURN USRFL112 310 USERFL=(RESPON(2)-RESPON(1))/(T(2,1)-T(1,1)) USRFL113 RETURN USRFL114 320 NY1 =NY(1) USRFL115 NY1M1 =NY1-1 USRFL116 USERFL=(RESPON(NY1)-RESPON(NY1M1))/(T(NY1,1)-T(NY1M1,1)) USRFL117 RETURN USRFL118 C-----------------------------------------------------------------------USRFL119 C NGAM=4 USRFL120 C-----------------------------------------------------------------------USRFL121 400 USERFL=0. USRFL122 RETURN USRFL123 C-----------------------------------------------------------------------USRFL124 C NGAM=5 USRFL125 C-----------------------------------------------------------------------USRFL126 500 USERFL=0. USRFL127 RETURN USRFL128 C-----------------------------------------------------------------------USRFL129 C NGAM=6 USRFL130 C-----------------------------------------------------------------------USRFL131 600 USERFL=0. USRFL132 RETURN USRFL133 END USRFL134 C++++++++++++++++ DOUBLE PRECISION VERSION 3DP (JUN 1988) ++++++++++++++USRFN001 REAL FUNCTION USERFN (Z,J,NYMAX,T,NY,K,N,IDERIV,RESPON) USRFN002 C-----------------------------------------------------------------------USRFN003 C USER SUPPLIED SUBROUTINE. USRFN004 C EVALUATES THEORETICAL NONLINEAR MODEL FUNCTION OF EQ (3.1-1) OR USRFN005 C EQ (4.1.2-1) IN THE USERS MANUAL FOR A GIVEN PARAMETER USRFN006 C Z=ALOG(PLAM(J)) AND GIVEN T-VALUE T(K,N) (IDERIV=0), OR USRFN007 C ITS 1ST (IDERIV=1), OR 2ND (IDERIV=2) DERIVATIVE W.R.T. Z. USRFN008 C USRFN009 C USES RECURSION FORMULAS IF POSSIBLE. USRFN010 C USRFN011 C IN CASE OF CONVOLUTION, RESPON CONTAINS RESPONSE FUNCTION. USRFN012 C WHEN USING RECURSION FORMULAS YOU WILL NEED AN ADDITIONAL COMMON USRFN013 C BLOCK TO STORE THE FUNCTION VALUE, ITS 1ST, AND 2ND DERIVATIVE, USRFN014 C AND OTHER QUANTITIES FROM PREVIOUS CALCULATIONS. USRFN015 C USRFN016 C BELOW ARE ILLUSTRATED 2 CASES OF RECURSION FORMULAS: USRFN017 C IF IUSER(1) = 1 USE RECURSIONS FOR PURE EXPONENTIALS USRFN018 C IF IUSER(1) = 2 USE RECURSIONS FOR CONVOLUTION USRFN019 C IF IUSER(1) = 3 USE RECURSIONS FOR SPECIAL CONVOLUTION AS USRFN020 C DESCRIBED IN REF. [6] OF USERS MANUAL USRFN021 C USRFN022 C THIS VERSION OF USERFN AND THE PROPER SETTINGS OF IUSER(1), USRFN023 C RUSER(1), AND LUSER(1) ARE DESCRIBED IN DETAIL IN SEC 6.1 USRFN024 C OF THE USERS MANUAL. USRFN025 C-----------------------------------------------------------------------USRFN026 C SUBPROGRAMS CALLED - ERRMES, USERTR USRFN027 C WHICH IN TURN CALL - ERRMES USRFN028 C COMMON USED - DBLOCK, SBLOCK, IBLOCK, LBLOCK, RECURS USRFN029 C-----------------------------------------------------------------------USRFN030 DOUBLE PRECISION PRECIS, RANGE, USRFN031 1 DZ, DZINV, VARMIN, VAROLD, VARSAV, ZSTART, USRFN032 2 ZERO, ONE, TWO, THREE, FOUR, SIX, TEN, USRFN033 3 TWOTHR, FORTHR, SIXTEL USRFN034 LOGICAL DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, DOADEX, USRFN035 1 DOSPL, DOSTRT, DOUSOU, LUSER, USRFN036 2 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRFN037 3 SPLINE, WRTBET, WRTCAP USRFN038 DIMENSION T(NYMAX,1), NY(1), RESPON(1), IHOLER(6) USRFN039 COMMON /DBLOCK/ PRECIS, RANGE, DZ, DZINV, VARMIN, VAROLD, USRFN040 1 VARSAV(10), ZSTART, ZERO, ONE, TWO, THREE, FOUR, USRFN041 2 SIX, TEN, TWOTHR, FORTHR, SIXTEL USRFN042 COMMON /SBLOCK/ SRANGE, CONVRG, PLMNMX(2), PNMNMX(2), RUSER(100),USRFN043 1 DBOX, EXMAX, PNG(10), ZTOTAL USRFN044 COMMON /IBLOCK/ IRWCAP, IUNIT, IWT, LINEPG, MCONV, METHOD, USRFN045 1 MIOERR, NABORT, ND, NERFIT, NGAM, NNL, NZ, USRFN046 2 IFORMT(70), IFORMW(70), IFORMY(70), IPLFIT(2), USRFN047 3 IPLRES(2), IPRINT(3,2), IPRITR(2), IUSER(50), USRFN048 4 MTRY(10,2), MXITER(2), NL(10), USRFN049 5 IALPHA(6), IGAMMA(6), ILAMDA(6), IB, INDEX(10), USRFN050 6 ISTTXT(28), ITITLE(80), MBOX, MMTRY, NB, NGAMM1, USRFN051 7 NGAMP1, NLIN, NONL, NONLP1, NPLINE, NYSUM, USRFN052 8 NIN, NOUT USRFN053 COMMON /LBLOCK/ DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, USRFN054 1 DOADEX(10,2), DOSPL(10,2), DOSTRT(10,2), USRFN055 2 DOUSOU(2), LUSER(30), USRFN056 3 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRFN057 4 SPLINE, WRTBET, WRTCAP USRFN058 C***********************************************************************USRFN059 C IF YOU CHANGE NZMAX OR NONLMX IN THE MAIN PROGRAM, THEN YOU MUST USRFN060 C READJUST THE DIMENSIONS OF THE FOLLOWING ARRAYS AS DESCRIBED BELOW - USRFN061 C USRFN062 COMMON /RECURS/ DERIV0(52), FOLD0(52), PLAMDT(52), USRFN063 1 DERIV1(5), DERIV2(5), USRFN064 2 FOLD1(5), FOLD2(5), USRFN065 3 EZ, DT, DT2, DTOLD0, DTOLD1, DTOLD2, SAVD0, SAVD1USRFN066 C USRFN067 C THE DIMENSIONS OF THE ARRAYS ABOVE MUST BE ADJUSTED AS FOLLOWS - USRFN068 C USRFN069 C COMMON /RECURS/ DERIV0(NZMAX), FOLD0(NZMAX), PLAMDT(NZMAX), USRFN070 C 1 DERIV1(NONLMX), DERIV2(NONLMX), USRFN071 C 2 FOLD1(NONLMX), FOLD2(NONLMX), USRFN072 C 3 EZ, DT, DT2, DTOLD0, DTOLD1, DTOLD2, SAVD0, SAVD1USRFN073 C USRFN074 C***********************************************************************USRFN075 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HF, 1HN/ USRFN076 C-----------------------------------------------------------------------USRFN077 IF (J .GT. NB) CALL ERRMES (1,.TRUE.,IHOLER,NOUT) USRFN078 IF (N .GT. ND) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) USRFN079 IF (K.LE.0 .OR. K.GT.NY(N)) CALL ERRMES (3,.TRUE.,IHOLER,NOUT) USRFN080 IF (IDERIV.LT.0 .OR. IDERIV.GT.2) USRFN081 1 CALL ERRMES (4,.TRUE.,IHOLER,NOUT) USRFN082 C-----------------------------------------------------------------------USRFN083 IF (IUSER(1) .EQ. 3) GO TO 600 USRFN084 IF (IUSER(1) .EQ. 2) GO TO 400 USRFN085 IF (IUSER(1) .EQ. 1) GO TO 100 USRFN086 CALL ERRMES (5,.TRUE.,IHOLER,NOUT) USRFN087 C-----------------------------------------------------------------------USRFN088 C USE RECURSION FORMULAS FOR EXPONENTIALS ONLY USRFN089 C-----------------------------------------------------------------------USRFN090 100 IF (IDERIV-1) 110,140,170 USRFN091 C-----------------------------------------------------------------------USRFN092 C FUNCTION VALUE USRFN093 C-----------------------------------------------------------------------USRFN094 110 IF (K .GT. 1) GO TO 120 USRFN095 PLAM =USERTR(Z,2) USRFN096 DTOLD0 =T(2,N)-T(1,N) USRFN097 EZ =PLAM*DTOLD0 USRFN098 FOLD0(J)=0. USRFN099 IF (EZ .LT. EXMAX) FOLD0(J)=EXP(-EZ) USRFN100 PLAMDT(J)=EZ USRFN101 EZ =PLAM*T(1,N) USRFN102 USERFN =0. USRFN103 IF (EZ .LT. EXMAX) USERFN=EXP(-EZ) USRFN104 DERIV0(J)=USERFN USRFN105 RETURN USRFN106 120 DT=T(K,N)-T(K-1,N) USRFN107 IF (DT .EQ. DTOLD0) GO TO 130 USRFN108 PLAM =USERTR(Z,2) USRFN109 EZ =PLAM*DT USRFN110 FOLD0(J)=0. USRFN111 IF (EZ .LT. EXMAX) FOLD0(J)=EXP(-EZ) USRFN112 PLAMDT(J)=EZ USRFN113 IF (J .EQ. NONL) DTOLD0=DT USRFN114 130 USERFN =DERIV0(J)*FOLD0(J) USRFN115 SAVD0 =DERIV0(J) USRFN116 DERIV0(J)=USERFN USRFN117 RETURN USRFN118 C-----------------------------------------------------------------------USRFN119 C 1ST DERIVATIVE USRFN120 C-----------------------------------------------------------------------USRFN121 140 IF (K .GT. 1) GO TO 150 USRFN122 DTOLD1 =DTOLD0 USRFN123 USERFN =-EZ*DERIV0(J) USRFN124 DERIV1(J)=USERFN USRFN125 FOLD1(J) =-PLAMDT(J)*FOLD0(J) USRFN126 RETURN USRFN127 150 IF (DT .EQ. DTOLD1) GO TO 160 USRFN128 FOLD1(J)=-PLAMDT(J)*FOLD0(J) USRFN129 IF (J .EQ. NONL) DTOLD1=DT USRFN130 160 USERFN =DERIV1(J)*FOLD0(J)+SAVD0*FOLD1(J) USRFN131 SAVD1 =DERIV1(J) USRFN132 DERIV1(J)=USERFN USRFN133 RETURN USRFN134 C-----------------------------------------------------------------------USRFN135 C 2ND DERIVATIVE USRFN136 C-----------------------------------------------------------------------USRFN137 170 IF (K .GT. 1) GO TO 180 USRFN138 DTOLD2 =DTOLD1 USRFN139 USERFN =DERIV1(J)*(1.0-EZ) USRFN140 DERIV2(J)=USERFN USRFN141 FOLD2(J) =FOLD1(J)*(1.0-PLAMDT(J)) USRFN142 RETURN USRFN143 180 IF (DT .EQ. DTOLD2) GO TO 190 USRFN144 FOLD2(J)=FOLD1(J)*(1.0-PLAMDT(J)) USRFN145 IF (J .EQ. NONL) DTOLD2=DT USRFN146 190 USERFN =DERIV2(J)*FOLD0(J)+2.0*SAVD1*FOLD1(J)+SAVD0*FOLD2(J) USRFN147 DERIV2(J)=USERFN USRFN148 RETURN USRFN149 C-----------------------------------------------------------------------USRFN150 C USE RECURSION FORMULAS FOR CONVOLUTIONS USRFN151 C-----------------------------------------------------------------------USRFN152 400 IF (IDERIV-1) 410,450,500 USRFN153 C-----------------------------------------------------------------------USRFN154 C FUNCTION VALUE USRFN155 C-----------------------------------------------------------------------USRFN156 410 IF (K .GT. 2) GO TO 430 USRFN157 IF (K .GT. 1) GO TO 420 USRFN158 PLAM =USERTR(Z,2) USRFN159 DTOLD0 =T(2,N)-T(1,N) USRFN160 EZ =PLAM*DTOLD0 USRFN161 FOLD0(J)=0. USRFN162 IF (EZ .LT. EXMAX) FOLD0(J)=EXP(-EZ) USRFN163 PLAMDT(J)=EZ USRFN164 USERFN =RESPON(1)*DTOLD0/2.0 USRFN165 SAVD0 =USERFN USRFN166 DERIV0(J)=USERFN USRFN167 RETURN USRFN168 420 USERFN =DERIV0(J)*FOLD0(J)+RESPON(2)*DTOLD0/2.0 USRFN169 DERIV0(J)=USERFN USRFN170 RETURN USRFN171 430 DT =T(K,N)-T(K-1,N) USRFN172 DT2=DT/2.0 USRFN173 IF (DT .EQ. DTOLD0) GO TO 440 USRFN174 PLAM =USERTR(Z,2) USRFN175 EZ =PLAM*DT USRFN176 FOLD0(J)=0. USRFN177 IF (EZ .LT. EXMAX) FOLD0(J)=EXP(-EZ) USRFN178 PLAMDT(J)=EZ USRFN179 IF (J .EQ. NONL) DTOLD0=DT USRFN180 440 USERFN =(DERIV0(J)+RESPON(K-1)*DT2)*FOLD0(J)+RESPON(K)*DT2 USRFN181 SAVD0 =DERIV0(J) USRFN182 DERIV0(J)=USERFN USRFN183 RETURN USRFN184 C-----------------------------------------------------------------------USRFN185 C 1ST DERIVATIVE USRFN186 C-----------------------------------------------------------------------USRFN187 450 IF (K .GT. 2) GO TO 470 USRFN188 IF (K .GT. 1) GO TO 460 USRFN189 DTOLD1 =DTOLD0 USRFN190 USERFN =0.0 USRFN191 DERIV1(J)=USERFN USRFN192 FOLD1(J) =-PLAMDT(J)*FOLD0(J) USRFN193 RETURN USRFN194 460 USERFN =SAVD0*FOLD1(J) USRFN195 DERIV1(J)=USERFN USRFN196 RETURN USRFN197 470 IF (DT .EQ. DTOLD1) GO TO 480 USRFN198 FOLD1(J)=-PLAMDT(J)*FOLD0(J) USRFN199 IF (J .EQ. NONL) DTOLD1=DT USRFN200 480 USERFN =DERIV1(J)*FOLD0(J)+(SAVD0+RESPON(K-1)*DT2)*FOLD1(J) USRFN201 SAVD1 =DERIV1(J) USRFN202 DERIV1(J)=USERFN USRFN203 RETURN USRFN204 C-----------------------------------------------------------------------USRFN205 C 2ND DERIVATIVE USRFN206 C-----------------------------------------------------------------------USRFN207 500 IF (K .GT. 2) GO TO 520 USRFN208 IF (K .GT. 1) GO TO 510 USRFN209 DTOLD2 =DTOLD1 USRFN210 USERFN =0.0 USRFN211 DERIV2(J)=USERFN USRFN212 FOLD2(J) =FOLD1(J)*(1.0-PLAMDT(J)) USRFN213 RETURN USRFN214 510 USERFN =SAVD0*FOLD2(J) USRFN215 DERIV2(J)=USERFN USRFN216 RETURN USRFN217 520 IF (DT .EQ. DTOLD2) GO TO 530 USRFN218 FOLD2(J)=FOLD1(J)*(1.0-PLAMDT(J)) USRFN219 IF (J .EQ. NONL) DTOLD2=DT USRFN220 530 USERFN =DERIV2(J)*FOLD0(J)+2.0*SAVD1*FOLD1(J) USRFN221 USERFN =USERFN+(SAVD0+RESPON(K-1)*DT2)*FOLD2(J) USRFN222 DERIV2(J)=USERFN USRFN223 RETURN USRFN224 C-----------------------------------------------------------------------USRFN225 C USE RECURSION FORMULAS FOR SPECIAL CONVOLUTION (REFERENCE LIFETIME USRFN226 C KNOWN AND STORED IN RUSER(1), SEE ALSO REF [6] OF USERS MANUAL) USRFN227 C-----------------------------------------------------------------------USRFN228 600 IF (IDERIV-1) 620,700,800 USRFN229 C-----------------------------------------------------------------------USRFN230 C FUNCTION VALUE USRFN231 C-----------------------------------------------------------------------USRFN232 620 IF (K .GT. 2) GO TO 660 USRFN233 IF (K .GT. 1) GO TO 640 USRFN234 PLAM =USERTR(Z,2) USRFN235 PLAMDT(J)=RUSER(1)-PLAM USRFN236 DT =T(2,N)-T(1,N) USRFN237 DT2 =DT/2. USRFN238 EZ =PLAM*DT USRFN239 FOLD0(J) =0. USRFN240 IF (EZ .LT. EXMAX) FOLD0(J)=EXP(-EZ) USRFN241 USERFN =PLAMDT(J)*RESPON(1)*DT2 USRFN242 DERIV0(J)=USERFN USRFN243 USERFN =USERFN+RESPON(1) USRFN244 RETURN USRFN245 640 USERFN =DERIV0(J)*FOLD0(J)+PLAMDT(J)*RESPON(2)*DT2 USRFN246 DERIV0(J)=USERFN USRFN247 USERFN =USERFN+RESPON(2) USRFN248 RETURN USRFN249 660 USERFN =(DERIV0(J)+PLAMDT(J)*RESPON(K-1)*DT2)*FOLD0(J)+ USRFN250 1 PLAMDT(J)*RESPON(K) *DT2 USRFN251 SAVD0 =DERIV0(J) USRFN252 DERIV0(J)=USERFN USRFN253 USERFN =USERFN+RESPON(K) USRFN254 RETURN USRFN255 C-----------------------------------------------------------------------USRFN256 C 1ST DERIVATIVE USRFN257 C-----------------------------------------------------------------------USRFN258 700 IF (K .GT. 2) GO TO 740 USRFN259 IF (K .GT. 1) GO TO 720 USRFN260 USERFN =-RESPON(1)*DT2 USRFN261 DERIV1(J)=USERFN USRFN262 FOLD1(J) =-DT*FOLD0(J) USRFN263 GO TO 760 USRFN264 720 USERFN =((PLAMDT(J)*FOLD1(J)-FOLD0(J))*RESPON(1)- USRFN265 1 RESPON(2))*DT2 USRFN266 DERIV1(J)=USERFN USRFN267 GO TO 760 USRFN268 740 USERFN =DERIV1(J)*FOLD0(J)+ USRFN269 1 (SAVD0+PLAMDT(J)*RESPON(K-1)*DT2)*FOLD1(J) USRFN270 USERFN =USERFN-(RESPON(K-1)*FOLD0(J)+RESPON(K))*DT2 USRFN271 DERIV1(J)=USERFN USRFN272 760 USERFN =(RUSER(1)-PLAMDT(J))*USERFN USRFN273 RETURN USRFN274 C-----------------------------------------------------------------------USRFN275 C 2ND DERIVATIVE USRFN276 C-----------------------------------------------------------------------USRFN277 800 CALL ERRMES (6,.TRUE.,IHOLER,NOUT) USRFN278 RETURN USRFN279 END USRFN280 C++++++++++++++++ DOUBLE PRECISION VERSION 3DP (JUN 1988) ++++++++++++++USRIN001 SUBROUTINE USERIN (NYMAX,Y,T,SQRTW,NY,RESPON) USRIN002 C-----------------------------------------------------------------------USRIN003 C THIS IS A USER-SUPPLIED ROUTINE (ONLY CALLED WHEN DOUSIN=.TRUE.) USRIN004 C THAT IS CALLED RIGHT AFTER THE INITIALIZATION AND INPUT OF THE USRIN005 C COMMON VARIABLES, T, Y, AND THE LEAST-SQUARES WEIGHTS. USRIN006 C THEREFORE, IT CAN BE USED TO MODIFY ANY OF THESE. USRIN007 C SEE THE USERS MANUAL (SEC 4.2) FOR THE POSITION IN THE INPUT USRIN008 C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. USRIN009 C USRIN010 C BELOW IS ILLUSTRATED THE CASE OF A CONVOLUTION, WHERE USRIN011 C IUSER(1) > 1 HAS TO BE SET IN THE INPUT DECK, AND ND=1. USRIN012 C USERIN READS THE RESPONSE/REFERENCE FUNCTION (FOR DETAILS USRIN013 C SEE SEC 6.1 AND REF. [6] OF THE USERS MANUAL), STORES IT IN USRIN014 C ARRAY RESPON, WRITES IT OUT (IF PRY=.TRUE.), AND MODIFIES USRIN015 C THE BOUNDS FOR THE LINEAR PARAMETERS. USRIN016 C USERIN HAS ALSO THE OPTION TO SUBTRACT BACKGROUND FROM THE USRIN017 C RESPONSE/REFERENCE FUNCTION. USRIN018 C USRIN019 C RESPON - REAL ARRAY OF DIMENSION AT LEAST NY. USRIN020 C LUSER(10) - LOGICAL, IF .TRUE. SUBTRACT BACKGROUND (STORED USRIN021 C IN RUSER(10)) FROM ARRAY RESPON. USRIN022 C RUSER(10) - REAL, CONTAINS BACKGROUND OF RESPON (ONLY USRIN023 C NEEDED IF LUSER(10)=.TRUE.). USRIN024 C-----------------------------------------------------------------------USRIN025 C SUBPROGRAMS CALLED - NONE USRIN026 C COMMON USED - DBLOCK, SBLOCK, IBLOCK, LBLOCK USRIN027 C-----------------------------------------------------------------------USRIN028 DOUBLE PRECISION PRECIS, RANGE, USRIN029 1 DZ, DZINV, VARMIN, VAROLD, VARSAV, ZSTART, USRIN030 2 ZERO, ONE, TWO, THREE, FOUR, SIX, TEN, USRIN031 3 TWOTHR, FORTHR, SIXTEL USRIN032 LOGICAL DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, DOADEX, USRIN033 1 DOSPL, DOSTRT, DOUSOU, LUSER, USRIN034 2 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRIN035 3 SPLINE, WRTBET, WRTCAP USRIN036 DIMENSION Y(NYMAX,1), T(NYMAX,1), SQRTW(NYMAX,1), NY(1), USRIN037 1 RESPON(1), IHOLER(6) USRIN038 COMMON /DBLOCK/ PRECIS, RANGE, DZ, DZINV, VARMIN, VAROLD, USRIN039 1 VARSAV(10), ZSTART, ZERO, ONE, TWO, THREE, FOUR, USRIN040 2 SIX, TEN, TWOTHR, FORTHR, SIXTEL USRIN041 COMMON /SBLOCK/ SRANGE, CONVRG, PLMNMX(2), PNMNMX(2), RUSER(100),USRIN042 1 DBOX, EXMAX, PNG(10), ZTOTAL USRIN043 COMMON /IBLOCK/ IRWCAP, IUNIT, IWT, LINEPG, MCONV, METHOD, USRIN044 1 MIOERR, NABORT, ND, NERFIT, NGAM, NNL, NZ, USRIN045 2 IFORMT(70), IFORMW(70), IFORMY(70), IPLFIT(2), USRIN046 3 IPLRES(2), IPRINT(3,2), IPRITR(2), IUSER(50), USRIN047 4 MTRY(10,2), MXITER(2), NL(10), USRIN048 5 IALPHA(6), IGAMMA(6), ILAMDA(6), IB, INDEX(10), USRIN049 6 ISTTXT(28), ITITLE(80), MBOX, MMTRY, NB, NGAMM1, USRIN050 7 NGAMP1, NLIN, NONL, NONLP1, NPLINE, NYSUM, USRIN051 8 NIN, NOUT USRIN052 COMMON /LBLOCK/ DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, USRIN053 1 DOADEX(10,2), DOSPL(10,2), DOSTRT(10,2), USRIN054 2 DOUSOU(2), LUSER(30), USRIN055 3 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRIN056 4 SPLINE, WRTBET, WRTCAP USRIN057 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HI, 1HN/ USRIN058 C-----------------------------------------------------------------------USRIN059 C READ IN RESPONSE/REFERENCE FUNCTION AND WRITE IT (IF PRY=.TRUE.) USRIN060 C (IUSER(1) > 1 (I.E. CONVOLUTION) HAS TO BE SET IN THE INPUT DECK). USRIN061 C-----------------------------------------------------------------------USRIN062 MY=NY(1) USRIN063 READ (NIN,IFORMY) (RESPON(K),K=1,MY) USRIN064 IF (.NOT.PRY) GO TO 100 USRIN065 WRITE (NOUT,1000) USRIN066 WRITE (NOUT,1100) (RESPON(K),K=1,MY) USRIN067 C-----------------------------------------------------------------------USRIN068 C IF LUSER(10)=.TRUE. SUBTRACT BACKGROUND STORED IN RUSER(10). USRIN069 C-----------------------------------------------------------------------USRIN070 100 IF (.NOT. LUSER(10)) GO TO 300 USRIN071 DO 200 K=1,MY USRIN072 200 RESPON(K)=RESPON(K)-RUSER(10) USRIN073 C-----------------------------------------------------------------------USRIN074 C RECOMPUTE PLMNMX(*), THE FACTORS TO SET BOUNDS FOR THE LINEAR USRIN075 C PARAMETERS (SEE USERS MANUAL SEC 3.4.3 FOR MORE INFORMATION). USRIN076 C-----------------------------------------------------------------------USRIN077 300 REFMAX=RESPON(1)*(T(2,1)-T(1,1)) USRIN078 MYM1 =MY-1 USRIN079 DO 400 K=2,MYM1 USRIN080 400 REFMAX=AMAX1(REFMAX,RESPON(K)*.5*(T(K+1,1)-T(K-1,1))) USRIN081 REFMAX =AMAX1(REFMAX,RESPON(MY)*(T(MY,1)-T(MYM1,1))) USRIN082 PLMNMX(1)=PLMNMX(1)/REFMAX USRIN083 PLMNMX(2)=PLMNMX(2)/REFMAX USRIN084 RETURN USRIN085 1000 FORMAT(31H1 REFERENCE/RESPONSE - FUNCTION/) USRIN086 1100 FORMAT(1X,1P10E13.4) USRIN087 END USRIN088 C++++++++++++++++ DOUBLE PRECISION VERSION 3DP (JUN 1988) ++++++++++++++USROU001 SUBROUTINE USEROU (IROUTE,NYMAX,Y,T,YLYFIT,SQRTW,NY,NLINMX,NONLMX,USROU002 1 PLIN,PLAM,PLMHLP,STDDEV,VAR,ISTAGE,RESPON) USROU003 C-----------------------------------------------------------------------USROU004 C USER SUPPLIED SUBROUTINE. (ONLY CALLED WHEN DOUSOU=.TRUE.) USROU005 C TO PRODUCE YOUR OWN EXTRA OUTPUT. USROU006 C USROU007 C BELOW IS THE ILLUSTRATED THE CASE IROUTE = 5, USROU008 C TO OUTPUT THE RECIPROCAL PLAM AND ITS STANDARD DEVIATION. USROU009 C ONLY DONE IF ISTAGE=2. USROU010 C-----------------------------------------------------------------------USROU011 C USROU012 C DESCRIPTION OF PARAMETERS: USROU013 C USROU014 C IROUTE - TELLS YOU FROM WHICH SUBROUTINE CALL WAS, SO YOU USROU015 C HAVE AN IDEA WHAT IS COMPUTED ALREADY, AND WHAT CANUSROU016 C BE USED FOR DOING YOUR OWN OUTPUT. USROU017 C IROUTE = 1 FROM ANALYZ, TO OUTPUT YOUR OWN USROU018 C STARTING VALUES FOR A COMPLETE ANALYSISUSROU019 C IROUTE = 2 FROM LSTSQR, AFTER A NOT NECESSARILY USROU020 C SUCCESFULL ITERATION USROU021 C IROUTE = 3 FROM OUTALL, TO OUTPUT PARAMETERS USROU022 C AND STD.DEV. USROU023 C IROUTE = 4 FROM OUTALL, TO OUTPUT RESIDUALS USROU024 C AND FITTED DATA. USROU025 C E.G. CALL A PLOT ROUTINE. USROU026 C IROUTE = 5 FROM OUTALL, TO OUTPUT BEST PARAMETERS USROU027 C AND STD.DEV. USROU028 C IROUTE = 6 FROM OUTALL, TO OUTPUT RESIDUALS AND USROU029 C FITTED DATA FOR BEST SOLUTION USROU030 C E.G. CALL A PLOT ROUTINE. USROU031 C NYMAX - FIRST DIMENSION OF ARRAYS Y, T, YLYFIT, AND SQRTW, USROU032 C I.E. MAXIMUM NUMBER OF DATA POINTS USROU033 C (NYMAX = MAX(NY(I)),I=1,...,ND) USROU034 C Y - REAL ARRAY OF DIMENSION (NYMAX,ND), CONTAINS DATA USROU035 C POINTS FOR ND DIFFERENT DATA SETS USROU036 C T - REAL ARRAY OF DIMENSION (NYMAX,ND), CONTAINS USROU037 C ABSICCA-VALUES FOR ND DIFFERENT DATA SETS USROU038 C YLYFIT - REAL ARRAY OF DIMENSION (NYMAX,ND), MAY CONTAIN USROU039 C RESIDUALS FOR ND DIFFERENT DATA SETS, IT DEPENDS USROU040 C ON IROUTE AND IPRINT (IN COMMON) USROU041 C SQRTW - REAL ARRAY OF DIMENSION (NYMAX,ND), CONTAINS USROU042 C WEIGHTS USROU043 C NY - ARRAY OF DIMENSION ND, CONTAINS NUMBER OF DATA USROU044 C POINTS FOR THE ND DIFFERENT DATA SETS USROU045 C NLINMX - FIRST DIMENSION OF ARRAY PLIN, I.E. MAXIMUM NUMBER USROU046 C OF LINEAR PARAMETERS USROU047 C NONLMX - DIMENSION OF ARRAY PLAM, I.E. MAXIMUM NUMBER OF USROU048 C NONLINEAR PARAMETERS USROU049 C PLIN - REAL ARRAY OF DIMENSION (NLINMX,ND), CONTAINS USROU050 C LINEAR PARAMETERS FOR ND DIFFERENT DATA SETS USROU051 C PLAM - REAL ARRAY OF DIMENSION NONLMX, CONTAINS NONLINEAR USROU052 C PARAMETERS USROU053 C PLMHLP - REAL ARRAY OF DIMENSION NONLMX, USED AS WORKING USROU054 C SPACE USROU055 C STDDEV - REAL ARRAY OF DIMENSION NLNDNL, CONTAINS DIAGONAL USROU056 C ELEMENTS OF INVERTED FULL LEAST SQUARES MATRIX, USROU057 C FOR IROUTE = 3 OR 5 ONLY. USROU058 C VAR - REAL, CONTAINS SUM OF SQUARED WEIGHTED RESIDUALS, USROU059 C FOR IROUTE > 1 ONLY. USROU060 C ISTAGE - IF ISTAGE = 1, THEN YOU ARE IN THE PRELIMINARY USROU061 C ANALYSIS TO DETERMINE WEIGHTS USROU062 C IF ISTAGE = 2, THEN YOU ARE IN THE FINAL ANALYSIS USROU063 C RESPON - REAL ARRAY, CONTAINS RESPONSE FUNCTION IN CASE OF USROU064 C CONVOLUTION (IUSER(1)=2). USROU065 C USROU066 C-----------------------------------------------------------------------USROU067 C SUBPROGRAMS CALLED - ERRMES, USERTR USROU068 C WHICH IN TURN CALL - ERRMES USROU069 C COMMON USED - DBLOCK, SBLOCK, IBLOCK, LBLOCK USROU070 C-----------------------------------------------------------------------USROU071 DOUBLE PRECISION PRECIS, RANGE, USROU072 1 DZ, DZINV, VARMIN, VAROLD, VARSAV, ZSTART, USROU073 2 ZERO, ONE, TWO, THREE, FOUR, SIX, TEN, USROU074 3 TWOTHR, FORTHR, SIXTEL USROU075 DOUBLE PRECISION VAR USROU076 LOGICAL DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, DOADEX, USROU077 1 DOSPL, DOSTRT, DOUSOU, LUSER, USROU078 2 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USROU079 3 SPLINE, WRTBET, WRTCAP USROU080 DIMENSION Y(NYMAX,1), T(NYMAX,1), YLYFIT(NYMAX,1), USROU081 1 SQRTW(NYMAX,1), NY(1), PLIN(NLINMX,1), USROU082 2 PLAM(NONLMX), PLMHLP(NONLMX), STDDEV(1), USROU083 3 RESPON(1) USROU084 DIMENSION IHOLER(6) USROU085 COMMON /DBLOCK/ PRECIS, RANGE, DZ, DZINV, VARMIN, VAROLD, USROU086 1 VARSAV(10), ZSTART, ZERO, ONE, TWO, THREE, FOUR, USROU087 2 SIX, TEN, TWOTHR, FORTHR, SIXTEL USROU088 COMMON /SBLOCK/ SRANGE, CONVRG, PLMNMX(2), PNMNMX(2), RUSER(100),USROU089 1 DBOX, EXMAX, PNG(10), ZTOTAL USROU090 COMMON /IBLOCK/ IRWCAP, IUNIT, IWT, LINEPG, MCONV, METHOD, USROU091 1 MIOERR, NABORT, ND, NERFIT, NGAM, NNL, NZ, USROU092 2 IFORMT(70), IFORMW(70), IFORMY(70), IPLFIT(2), USROU093 3 IPLRES(2), IPRINT(3,2), IPRITR(2), IUSER(50), USROU094 4 MTRY(10,2), MXITER(2), NL(10), USROU095 5 IALPHA(6), IGAMMA(6), ILAMDA(6), IB, INDEX(10), USROU096 6 ISTTXT(28), ITITLE(80), MBOX, MMTRY, NB, NGAMM1, USROU097 7 NGAMP1, NLIN, NONL, NONLP1, NPLINE, NYSUM, USROU098 8 NIN, NOUT USROU099 COMMON /LBLOCK/ DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, USROU100 1 DOADEX(10,2), DOSPL(10,2), DOSTRT(10,2), USROU101 2 DOUSOU(2), LUSER(30), USROU102 3 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USROU103 4 SPLINE, WRTBET, WRTCAP USROU104 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HO, 1HU/ USROU105 C-----------------------------------------------------------------------USROU106 IF (ISTAGE.LT.1.OR.ISTAGE.GT.2) CALL ERRMES (1,.TRUE.,IHOLER,NOUT)USROU107 IF (IROUTE.LT.1.OR.IROUTE.GT.6) CALL ERRMES (2,.TRUE.,IHOLER,NOUT)USROU108 C-----------------------------------------------------------------------USROU109 IF (ISTAGE .EQ. 1) RETURN USROU110 GO TO (100,200,300,400,500,600), IROUTE USROU111 C-----------------------------------------------------------------------USROU112 C IROUTE = 1 USROU113 C-----------------------------------------------------------------------USROU114 100 CONTINUE USROU115 RETURN USROU116 C-----------------------------------------------------------------------USROU117 C IROUTE = 2 USROU118 C-----------------------------------------------------------------------USROU119 200 CONTINUE USROU120 RETURN USROU121 C-----------------------------------------------------------------------USROU122 C IROUTE = 3 USROU123 C-----------------------------------------------------------------------USROU124 300 CONTINUE USROU125 RETURN USROU126 C-----------------------------------------------------------------------USROU127 C IROUTE = 4 USROU128 C-----------------------------------------------------------------------USROU129 400 CONTINUE USROU130 RETURN USROU131 C-----------------------------------------------------------------------USROU132 C IROUTE = 5 USROU133 C-----------------------------------------------------------------------USROU134 500 CONTINUE USROU135 IF (ISTAGE .EQ. 1) RETURN USROU136 WRITE (NOUT,5000) USROU137 NLNDNL=NLIN*ND+NONL USROU138 STDVAR=VAR/FLOAT(NYSUM-NLNDNL) USROU139 STDVAR=SQRT(STDVAR) USROU140 WRITE (NOUT,5300) USROU141 DO 520 J=1,NONL USROU142 JD =J+NLIN*ND USROU143 PLMINV=1./USERTR(PLAM(J),2) USROU144 ERRPLM=STDVAR*SQRT(STDDEV(JD))*ABS(PLMINV) USROU145 PRCPLM=100.*ERRPLM/ABS(PLMINV) USROU146 WRITE (NOUT,5200) PLIN(J+NGAM,1),PLMINV,ERRPLM,PRCPLM USROU147 520 CONTINUE USROU148 IF (ND.LE.1) RETURN USROU149 DO 530 K = 2,ND USROU150 WRITE (NOUT,5300) USROU151 DO 531 J = 1,NONL USROU152 WRITE (NOUT,5200) PLIN(J+NGAM,K) USROU153 531 CONTINUE USROU154 530 CONTINUE USROU155 RETURN USROU156 5000 FORMAT(///15X,5HALPHA,10X,36H1/(LAMBDA) +- STD. ERROR PERCENT)USROU157 5200 FORMAT(8X,1PE12.4,E20.4,4H +-,E12.4,0PF10.3) USROU158 5300 FORMAT(1H ) USROU159 C-----------------------------------------------------------------------USROU160 C IROUTE = 6 USROU161 C-----------------------------------------------------------------------USROU162 600 CONTINUE USROU163 RETURN USROU164 END USROU165 C++++++++++++++++ DOUBLE PRECISION VERSION 3DP (JUN 1988) ++++++++++++++USRSI001 SUBROUTINE USERSI (NYMAX,Y,T,EXACT,NY,RESPON) USRSI002 C-----------------------------------------------------------------------USRSI003 C USER SUPPLIED SUBROUTINE. (ONLY CALLED IF SIMULA=.TRUE.) USRSI004 C FOR CALCULATING EXACT (THE SIMULATED DATA BEFORE NOISE IS ADDED) USRSI005 C AND Y(K,N) (THE SIMULATED NOISY DATA) FOR K=1,..,NY(N), N=1,..,ND.USRSI006 C IUSER(3) = STARTING INTEGER FOR RANDOM NUMBER GENERATOR RANDOM. USRSI007 C IUSER(3) AND RUSER(3) MUST BE SUPPLIED BY THE USER. USRSI008 C IUSER(3) MUST BE BETWEEN 1 AND 2147483646, IF IT IS NOT, THEN USRSI009 C IT IS SET TO THE DEFAULT VALUE OF 30171. USRSI010 C USAGE OF RUSER(3) IS DESCRIBED BELOW. USRSI011 C SEE THE USERS MANUAL (SEC 4.2) FOR THE POSITION IN THE INPUT USRSI012 C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. USRSI013 C-----------------------------------------------------------------------USRSI014 C SUBPROGRAMS CALLED - RGAUSS, USEREX USRSI015 C WHICH IN TURN CALL - RANDOM USRSI016 C COMMON USED - DBLOCK, SBLOCK, IBLOCK, LBLOCK USRSI017 C-----------------------------------------------------------------------USRSI018 DOUBLE PRECISION PRECIS, RANGE, USRSI019 1 DZ, DZINV, VARMIN, VAROLD, VARSAV, ZSTART, USRSI020 2 ZERO, ONE, TWO, THREE, FOUR, SIX, TEN, USRSI021 3 TWOTHR, FORTHR, SIXTEL USRSI022 DOUBLE PRECISION DUB USRSI023 LOGICAL DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, DOADEX, USRSI024 1 DOSPL, DOSTRT, DOUSOU, LUSER, USRSI025 2 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRSI026 3 SPLINE, WRTBET, WRTCAP USRSI027 DIMENSION Y(NYMAX,1), T(NYMAX,1), EXACT(NYMAX,1), USRSI028 1 NY(1), RN(2), RESPON(1), IHOLER(6) USRSI029 COMMON /DBLOCK/ PRECIS, RANGE, DZ, DZINV, VARMIN, VAROLD, USRSI030 1 VARSAV(10), ZSTART, ZERO, ONE, TWO, THREE, FOUR, USRSI031 2 SIX, TEN, TWOTHR, FORTHR, SIXTEL USRSI032 COMMON /SBLOCK/ SRANGE, CONVRG, PLMNMX(2), PNMNMX(2), RUSER(100),USRSI033 1 DBOX, EXMAX, PNG(10), ZTOTAL USRSI034 COMMON /IBLOCK/ IRWCAP, IUNIT, IWT, LINEPG, MCONV, METHOD, USRSI035 1 MIOERR, NABORT, ND, NERFIT, NGAM, NNL, NZ, USRSI036 2 IFORMT(70), IFORMW(70), IFORMY(70), IPLFIT(2), USRSI037 3 IPLRES(2), IPRINT(3,2), IPRITR(2), IUSER(50), USRSI038 4 MTRY(10,2), MXITER(2), NL(10), USRSI039 5 IALPHA(6), IGAMMA(6), ILAMDA(6), IB, INDEX(10), USRSI040 6 ISTTXT(28), ITITLE(80), MBOX, MMTRY, NB, NGAMM1, USRSI041 7 NGAMP1, NLIN, NONL, NONLP1, NPLINE, NYSUM, USRSI042 8 NIN, NOUT USRSI043 COMMON /LBLOCK/ DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, USRSI044 1 DOADEX(10,2), DOSPL(10,2), DOSTRT(10,2), USRSI045 2 DOUSOU(2), LUSER(30), USRSI046 3 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRSI047 4 SPLINE, WRTBET, WRTCAP USRSI048 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HS, 1HI/ USRSI049 C-----------------------------------------------------------------------USRSI050 TWOPI=6.2831853072D0 USRSI051 DUB =DBLE(FLOAT(IUSER(3))) USRSI052 IF (DUB.LT.1.D0 .OR. DUB.GT.2147483646.D0) DUB=30171.D0 USRSI053 DO 200 N=1,ND USRSI054 MY=NY(N) USRSI055 NN=N USRSI056 DO 100 K=1,MY USRSI057 KK=K USRSI058 C-----------------------------------------------------------------------USRSI059 C USEREX DELIVERES EXACT FUNCTION VALUE USRSI060 C-----------------------------------------------------------------------USRSI061 EXACT(K,N)=USEREX(NYMAX,T,NY,KK,NN,RESPON) USRSI062 C-----------------------------------------------------------------------USRSI063 C RGAUSS DELIVERS TWO NORMAL DEVIATES. THEREFORE IT IS ONLY USRSI064 C CALLED FOR ODD K. USRSI065 C-----------------------------------------------------------------------USRSI066 J=2-MOD(K,2) USRSI067 IF (J .EQ. 1) CALL RGAUSS (RN(1),RN(2),TWOPI,DUB) USRSI068 C-----------------------------------------------------------------------USRSI069 C THE FOLLOWING (COMMENT) STATEMENT IS APPROPRIATE WHEN THE EXPECTED USRSI070 C ERROR IS PROPORTIONAL TO SQRT(EXACT), AS IN POISSON PROCESSES. USRSI071 C FOR POISSON STATISTICS, RUSER(3) IS SIMPLY TO CORRECT FOR ANY USRSI072 C PREVIOUS SCALING IN EXACT, USRSI073 C I.E., RUSER(3)=SQRT(EXACT/(TOTAL COUNTS IN CHANNEL K)). USRSI074 C IF EXACT ALREADY IS IN COUNTS, THEN RUSER(3)=1. USRSI075 C-----------------------------------------------------------------------USRSI076 C Y(K,N)=EXACT(K,N)+RUSER(3)*RN(J)*SQRT(EXACT(K,N)) USRSI077 C-----------------------------------------------------------------------USRSI078 C THE FOLLOWING STATEMENT WOULD SIMULATE ZERO-MEAN GAUSSIAN NOISE USRSI079 C WITH STANDARD DEVIATION RUSER(3). USRSI080 C-----------------------------------------------------------------------USRSI081 Y(K,N)=EXACT(K,N)+RUSER(3)*RN(J) USRSI082 100 CONTINUE USRSI083 200 CONTINUE USRSI084 RETURN USRSI085 END USRSI086 C++++++++++++++++ DOUBLE PRECISION VERSION 3DP (JUN 1988) ++++++++++++++USRST001 SUBROUTINE USERST (ISTAGE,NYMAX,T,NY,NONLMX,PLAM) USRST002 C-----------------------------------------------------------------------USRST003 C USER SUPPLIED SUBROUTINE. USRST004 C (ONLY CALLED WHEN DOSTRT(*,ISTAGE)=.TRUE.) USRST005 C IF YOU WANT TO START FIRST TRY OF MTRY(*,ISTAGE) ANALYSES ASSUMINGUSRST006 C NONL=NL(*) COMPONENTS WITH YOUR OWN STARTING VALUES FOR NONLINEAR USRST007 C PARAMETERS PLAM. USRST008 C BELOW IS ILLUSTRATED THE VERY SIMPLE CASE, THAT YOU HAVE STORED USRST009 C IN RUSER(50+I), I=1,...,NONL, THE STARTING VALUES FOR PLAM. USRST010 C-----------------------------------------------------------------------USRST011 C USRST012 C DESCRIPTION OF PARAMETERS: USRST013 C USRST014 C ISTAGE - IF ISTAGE = 1 YOU ARE IN PRELIMINARY ANALYSIS TO USRST015 C DETERMINE WEIGHTS USRST016 C IF ISTAGE = 2 YOU ARE IN FINAL ANALYSIS USRST017 C NONLMX - DIMENSION OF ARRAY PLAM, I.E. MAXIMUM NUMBER USRST018 C OF NONLINEAR PARAMETERS USRST019 C PLAM - REAL ARRAY OF DIMENSION NONLMX, CONTAINS ON USRST020 C OUTPUT STARTING VALUES FOR NONLINEAR PARAMETERS USRST021 C USRST022 C-----------------------------------------------------------------------USRST023 C SUBPROGRAMS CALLED - ERRMES, USERTR USRST024 C WHICH IN TURN CALL - ERRMES USRST025 C COMMON USED - DBLOCK, SBLOCK, IBLOCK, LBLOCK USRST026 C-----------------------------------------------------------------------USRST027 DOUBLE PRECISION PRECIS, RANGE, USRST028 1 DZ, DZINV, VARMIN, VAROLD, VARSAV, ZSTART, USRST029 2 ZERO, ONE, TWO, THREE, FOUR, SIX, TEN, USRST030 3 TWOTHR, FORTHR, SIXTEL USRST031 LOGICAL DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, DOADEX, USRST032 1 DOSPL, DOSTRT, DOUSOU, LUSER, USRST033 2 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRST034 3 SPLINE, WRTBET, WRTCAP USRST035 DIMENSION T(NYMAX,1), NY(1), PLAM(NONLMX), IHOLER(6) USRST036 COMMON /DBLOCK/ PRECIS, RANGE, DZ, DZINV, VARMIN, VAROLD, USRST037 1 VARSAV(10), ZSTART, ZERO, ONE, TWO, THREE, FOUR, USRST038 2 SIX, TEN, TWOTHR, FORTHR, SIXTEL USRST039 COMMON /SBLOCK/ SRANGE, CONVRG, PLMNMX(2), PNMNMX(2), RUSER(100),USRST040 1 DBOX, EXMAX, PNG(10), ZTOTAL USRST041 COMMON /IBLOCK/ IRWCAP, IUNIT, IWT, LINEPG, MCONV, METHOD, USRST042 1 MIOERR, NABORT, ND, NERFIT, NGAM, NNL, NZ, USRST043 2 IFORMT(70), IFORMW(70), IFORMY(70), IPLFIT(2), USRST044 3 IPLRES(2), IPRINT(3,2), IPRITR(2), IUSER(50), USRST045 4 MTRY(10,2), MXITER(2), NL(10), USRST046 5 IALPHA(6), IGAMMA(6), ILAMDA(6), IB, INDEX(10), USRST047 6 ISTTXT(28), ITITLE(80), MBOX, MMTRY, NB, NGAMM1, USRST048 7 NGAMP1, NLIN, NONL, NONLP1, NPLINE, NYSUM, USRST049 8 NIN, NOUT USRST050 COMMON /LBLOCK/ DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, USRST051 1 DOADEX(10,2), DOSPL(10,2), DOSTRT(10,2), USRST052 2 DOUSOU(2), LUSER(30), USRST053 3 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRST054 4 SPLINE, WRTBET, WRTCAP USRST055 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HS, 1HT/ USRST056 C-----------------------------------------------------------------------USRST057 IF (ISTAGE.LT.1.OR.ISTAGE.GT.2) CALL ERRMES (1,.TRUE.,IHOLER,NOUT)USRST058 DO 100 I=1,NONL USRST059 DUM =RUSER(50+I) USRST060 DUM =USERTR(DUM,1) USRST061 PLAM(I)=AMAX1(PNMNMX(1),AMIN1(DUM,PNMNMX(2))) USRST062 IF (PLAM(I) .NE. DUM) CALL ERRMES (2,.FALSE.,IHOLER,NOUT) USRST063 100 CONTINUE USRST064 RETURN USRST065 END USRST066 C++++++++++++++++ DOUBLE PRECISION VERSION 3DP (JUN 1988) ++++++++++++++USRTR001 FUNCTION USERTR (X,IFUNCT) USRTR002 C-----------------------------------------------------------------------USRTR003 C USER SUPPLIED SUBROUTINE. USRTR004 C FOR COMPUTING THE INTERPOLATION-GRID TRANSFORMATION (CALL IT H) USRTR005 C H(G) WHEN IFUNCT=1, THE INVERSE TRANSFORMATION WHEN IFUNCT=2, USRTR006 C AND THE DERIVATIVE OF THE TRANSFORMATION WHEN IFUNCT=3. USRTR007 C BELOW IS ILLUSTRATED THE CASE H(G)=ALOG(G). FOR ANOTHER H, YOU USRTR008 C CAN REPLACE THE STATEMENTS NUMBERED 110, 120, AND 130. USRTR009 C THESE ARE THE ONLY STATEMENTS THAT CAN BE REPLACED. ALSO NOTE USRTR010 C THAT ONLY AN H THAT IS MONOTONIC IN THE RANGE OF INTERPOLATION USRTR011 C MAKES SENSE. USRTR012 C-----------------------------------------------------------------------USRTR013 C SUBPROGRAMS CALLED - ERRMES USRTR014 C COMMON USED - DBLOCK, SBLOCK, IBLOCK, LBLOCK USRTR015 C-----------------------------------------------------------------------USRTR016 DOUBLE PRECISION PRECIS, RANGE, USRTR017 1 DZ, DZINV, VARMIN, VAROLD, VARSAV, ZSTART, USRTR018 2 ZERO, ONE, TWO, THREE, FOUR, SIX, TEN, USRTR019 3 TWOTHR, FORTHR, SIXTEL USRTR020 LOGICAL DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, DOADEX, USRTR021 1 DOSPL, DOSTRT, DOUSOU, LUSER, USRTR022 2 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRTR023 3 SPLINE, WRTBET, WRTCAP USRTR024 DIMENSION IHOLER(6) USRTR025 COMMON /DBLOCK/ PRECIS, RANGE, DZ, DZINV, VARMIN, VAROLD, USRTR026 1 VARSAV(10), ZSTART, ZERO, ONE, TWO, THREE, FOUR, USRTR027 2 SIX, TEN, TWOTHR, FORTHR, SIXTEL USRTR028 COMMON /SBLOCK/ SRANGE, CONVRG, PLMNMX(2), PNMNMX(2), RUSER(100),USRTR029 1 DBOX, EXMAX, PNG(10), ZTOTAL USRTR030 COMMON /IBLOCK/ IRWCAP, IUNIT, IWT, LINEPG, MCONV, METHOD, USRTR031 1 MIOERR, NABORT, ND, NERFIT, NGAM, NNL, NZ, USRTR032 2 IFORMT(70), IFORMW(70), IFORMY(70), IPLFIT(2), USRTR033 3 IPLRES(2), IPRINT(3,2), IPRITR(2), IUSER(50), USRTR034 4 MTRY(10,2), MXITER(2), NL(10), USRTR035 5 IALPHA(6), IGAMMA(6), ILAMDA(6), IB, INDEX(10), USRTR036 6 ISTTXT(28), ITITLE(80), MBOX, MMTRY, NB, NGAMM1, USRTR037 7 NGAMP1, NLIN, NONL, NONLP1, NPLINE, NYSUM, USRTR038 8 NIN, NOUT USRTR039 COMMON /LBLOCK/ DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, USRTR040 1 DOADEX(10,2), DOSPL(10,2), DOSTRT(10,2), USRTR041 2 DOUSOU(2), LUSER(30), USRTR042 3 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRTR043 4 SPLINE, WRTBET, WRTCAP USRTR044 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HT, 1HR/ USRTR045 C-----------------------------------------------------------------------USRTR046 IF (IFUNCT.LT.1.OR.IFUNCT.GT.3) CALL ERRMES (1,.TRUE.,IHOLER,NOUT)USRTR047 GO TO (110,120,130),IFUNCT USRTR048 C-----------------------------------------------------------------------USRTR049 C COMPUTE TRANSFORMATION. USRTR050 C-----------------------------------------------------------------------USRTR051 110 IF (X .LE. 0.) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) USRTR052 USERTR=ALOG(X) USRTR053 RETURN USRTR054 C-----------------------------------------------------------------------USRTR055 C COMPUTE INVERSE TRANSFORMATION. USRTR056 C-----------------------------------------------------------------------USRTR057 120 USERTR=0. USRTR058 IF (-X .GT. EXMAX) RETURN USRTR059 IF ( X .GT. EXMAX) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) USRTR060 USERTR=EXP(X) USRTR061 RETURN USRTR062 C-----------------------------------------------------------------------USRTR063 C COMPUTE DERIVATIVE OF TRANSFORMATION. USRTR064 C-----------------------------------------------------------------------USRTR065 130 IF (X .EQ. 0.) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) USRTR066 USERTR=1./X USRTR067 RETURN USRTR068 END USRTR069 C++++++++++++++++ DOUBLE PRECISION VERSION 3DP (JUN 1988) ++++++++++++++USRWT001 SUBROUTINE USERWT (NYMAX,Y,T,YLYFIT,SQRTW,NY,ERRFIT,RESPON) USRWT002 C-----------------------------------------------------------------------USRWT003 C USER SUPPLIED SUBROUTINE. (ONLY NEEDED WHEN IWT=5) USRWT004 C FOR CALCULATING SQRTW (SQUARE ROOTS OF THE LEAST SQUARES WEIGHTS) USRWT005 C FROM Y, YLYFIT, AND ERRFIT, AS EXPLAINED IN DETAIL IN SEC 4.1.3.4 USRWT006 C OF THE USERS MANUAL. USRWT007 C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT DATA DECK OF USRWT008 C ANY INPUT FOR THIS USER SUBPROGRAM. USRWT009 C BELOW IS ILLUSTRATED A CASE FROM PHOTON CORRELATION SPECTROSCOPY, USRWT010 C WHERE THE VARIANCE OF Y IS PROPORTIONAL TO (Y**2+1)/(4*Y**2). USRWT011 C-----------------------------------------------------------------------USRWT012 C SUBPROGRAMS CALLED - NONE USRWT013 C COMMON USED - DBLOCK, SBLOCK, IBLOCK, LBLOCK USRWT014 C-----------------------------------------------------------------------USRWT015 DOUBLE PRECISION PRECIS, RANGE, USRWT016 1 DZ, DZINV, VARMIN, VAROLD, VARSAV, ZSTART, USRWT017 2 ZERO, ONE, TWO, THREE, FOUR, SIX, TEN, USRWT018 3 TWOTHR, FORTHR, SIXTEL USRWT019 LOGICAL DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, DOADEX, USRWT020 1 DOSPL, DOSTRT, DOUSOU, LUSER, USRWT021 2 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRWT022 3 SPLINE, WRTBET, WRTCAP USRWT023 DIMENSION Y(NYMAX,1), T(NYMAX,1), YLYFIT(NYMAX,1), USRWT024 1 SQRTW(NYMAX,1), NY(1), ERRFIT(1), RESPON(1) USRWT025 DIMENSION IHOLER(6) USRWT026 COMMON /DBLOCK/ PRECIS, RANGE, DZ, DZINV, VARMIN, VAROLD, USRWT027 1 VARSAV(10), ZSTART, ZERO, ONE, TWO, THREE, FOUR, USRWT028 2 SIX, TEN, TWOTHR, FORTHR, SIXTEL USRWT029 COMMON /SBLOCK/ SRANGE, CONVRG, PLMNMX(2), PNMNMX(2), RUSER(100),USRWT030 1 DBOX, EXMAX, PNG(10), ZTOTAL USRWT031 COMMON /IBLOCK/ IRWCAP, IUNIT, IWT, LINEPG, MCONV, METHOD, USRWT032 1 MIOERR, NABORT, ND, NERFIT, NGAM, NNL, NZ, USRWT033 2 IFORMT(70), IFORMW(70), IFORMY(70), IPLFIT(2), USRWT034 3 IPLRES(2), IPRINT(3,2), IPRITR(2), IUSER(50), USRWT035 4 MTRY(10,2), MXITER(2), NL(10), USRWT036 5 IALPHA(6), IGAMMA(6), ILAMDA(6), IB, INDEX(10), USRWT037 6 ISTTXT(28), ITITLE(80), MBOX, MMTRY, NB, NGAMM1, USRWT038 7 NGAMP1, NLIN, NONL, NONLP1, NPLINE, NYSUM, USRWT039 8 NIN, NOUT USRWT040 COMMON /LBLOCK/ DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, USRWT041 1 DOADEX(10,2), DOSPL(10,2), DOSTRT(10,2), USRWT042 2 DOUSOU(2), LUSER(30), USRWT043 3 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, USRWT044 4 SPLINE, WRTBET, WRTCAP USRWT045 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HW, 1HT/ USRWT046 C-----------------------------------------------------------------------USRWT047 C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE USRWT048 C FOR YOUR APPLICATION. USRWT049 C-----------------------------------------------------------------------USRWT050 DO 100 N=1,ND USRWT051 MY=NY(N) USRWT052 DO 100 K=1,MY USRWT053 DUM =AMAX1(ABS(Y(K,N)-YLYFIT(K,N)),ERRFIT(N)) USRWT054 SQRTW(K,N)=2.*DUM/SQRT(DUM*DUM+1.) USRWT055 100 CONTINUE USRWT056 RETURN USRWT057 END USRWT058 C++++++++++++++++ DOUBLE PRECISION VERSION 3DP (JUN 1988) ++++++++++++++ANLYZ001 SUBROUTINE ANALYZ (ISTAGE,NYMAX,Y,T,YLYFIT,SQRTW,NY,NLINMX,NONLMX,ANLYZ002 1 PLIN,PLNTRY,PLAM,PLMTRY,PLMSAV,PIVLIN,PIVLAM, ANLYZ003 2 EK,EPK,ETE,R,DELTAL,DELTAN,GTG,H,DELTAH,EKEK, ANLYZ004 3 ASAVE,NOKSET,NCOMBI,MTRYMX,STARTV,NLNDMX,ATOT, ANLYZ005 4 DTOLER,BUFF,DX,BSPL0,BSPL1,BSPL2,IC,NZMAX, ANLYZ006 5 NGAMMX,CCAP,ECAP,YCAP,SGGCAP,SGYCAP,YW,RESPON) ANLYZ007 C-----------------------------------------------------------------------ANLYZ008 C FIND STARTING VALUES FOR PLAM, PERFORMS LEAST SQUARES ANALYSIS, AND ANLYZ009 C OUTPUT FINAL RESULTS AS SPECIFIED THROUGH IPRINT, IPLFIT, AND IPLRES ANLYZ010 C-----------------------------------------------------------------------ANLYZ011 C SUBPROGRAMS CALLED - ERRMES, FISHNI, LSTSQR, NOKSUB, NPASCL, ANLYZ012 C OUTALL, USEROU, USERST, USERTR, VARIAN ANLYZ013 C WHICH IN TURN CALL - BETAIN, BSPLIN, ERRMES, FULLSM, GAMLN, ANLYZ014 C GETPRU, HESEXA, HESSEN, HESSPL, INTERP, ANLYZ015 C INTERV, LINEQS, LINPAR, NPASCL, OUTCOR, ANLYZ016 C OUTITR, OUTPAR, PGAUSS, PLPRIN, PLRES, ANLYZ017 C PIVOT, PIVOT1, USERFL, USERFN, USEROU, ANLYZ018 C USERTR, VARIAN ANLYZ019 C COMMON USED - DBLOCK, SBLOCK, IBLOCK, LBLOCK ANLYZ020 C-----------------------------------------------------------------------ANLYZ021 DOUBLE PRECISION PRECIS, RANGE, ANLYZ022 1 DZ, DZINV, VARMIN, VAROLD, VARSAV, ZSTART, ANLYZ023 2 ZERO, ONE, TWO, THREE, FOUR, SIX, TEN, ANLYZ024 3 TWOTHR, FORTHR, SIXTEL ANLYZ025 DOUBLE PRECISION EK, EPK, ETE, R, DELTAL, DELTAN, GTG, H, ANLYZ026 1 DELTAH, EKEK, ASAVE, ATOT, DTOLER, DX, BSPL0, ANLYZ027 2 BSPL1, BSPL2, CCAP, ECAP, YCAP, SGGCAP, ANLYZ028 3 SGYCAP, YW, VAR, VARBES ANLYZ029 LOGICAL DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, DOADEX, ANLYZ030 1 DOSPL, DOSTRT, DOUSOU, LUSER, ANLYZ031 2 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, ANLYZ032 3 SPLINE, WRTBET, WRTCAP ANLYZ033 LOGICAL PIVLIN, PIVLAM, STARTV, OWNST, MORE, BETTER, ANLYZ034 1 PRBEST, PRHEAD, ERRONE, ADEXCT ANLYZ035 DIMENSION Y(NYMAX,1), T(NYMAX,1), YLYFIT(NYMAX,1), ANLYZ036 1 SQRTW(NYMAX,1), NY(1), RESPON(1), PLIN(NLINMX,1),ANLYZ037 2 PLNTRY(NLINMX,1), PLAM(NONLMX), PLMTRY(NONLMX), ANLYZ038 3 PLMSAV(NONLMX,1), PIVLIN(NLINMX,1), ANLYZ039 4 PIVLAM(NLNDMX), EK(NLINMX), EPK(NONLMX), ANLYZ040 5 ETE(NLINMX,NLINMX), R(NLINMX,1), DELTAL(NLINMX), ANLYZ041 6 DELTAN(NONLMX), GTG(NONLMX,NONLMX) ANLYZ042 DIMENSION H(NLINMX,NLINMX,1), DELTAH(NONLMX,1), ANLYZ043 1 EKEK(NLINMX,NLINMX,1), ASAVE(NLINMX,NLINMX,1), ANLYZ044 2 NOKSET(NONLMX), NCOMBI(NONLMX), STARTV(MTRYMX), ANLYZ045 3 ATOT(NLNDMX,NLNDMX), DTOLER(NLNDMX), DX(NONLMX), ANLYZ046 4 BUFF(NLNDMX), BSPL0(4,NONLMX), BSPL1(4,NONLMX), ANLYZ047 5 BSPL2(4,NONLMX), IC(NONLMX), CCAP(NZMAX,NZMAX,1),ANLYZ048 6 ECAP(NZMAX,NGAMMX,1), YCAP(NZMAX,1), ANLYZ049 7 SGGCAP(NGAMMX,NGAMMX,1), SGYCAP(NGAMMX,1), YW(1) ANLYZ050 DIMENSION IHOLER(6) ANLYZ051 COMMON /DBLOCK/ PRECIS, RANGE, DZ, DZINV, VARMIN, VAROLD, ANLYZ052 1 VARSAV(10), ZSTART, ZERO, ONE, TWO, THREE, FOUR, ANLYZ053 2 SIX, TEN, TWOTHR, FORTHR, SIXTEL ANLYZ054 COMMON /SBLOCK/ SRANGE, CONVRG, PLMNMX(2), PNMNMX(2), RUSER(100),ANLYZ055 1 DBOX, EXMAX, PNG(10), ZTOTAL ANLYZ056 COMMON /IBLOCK/ IRWCAP, IUNIT, IWT, LINEPG, MCONV, METHOD, ANLYZ057 1 MIOERR, NABORT, ND, NERFIT, NGAM, NNL, NZ, ANLYZ058 2 IFORMT(70), IFORMW(70), IFORMY(70), IPLFIT(2), ANLYZ059 3 IPLRES(2), IPRINT(3,2), IPRITR(2), IUSER(50), ANLYZ060 4 MTRY(10,2), MXITER(2), NL(10), ANLYZ061 5 IALPHA(6), IGAMMA(6), ILAMDA(6), IB, INDEX(10), ANLYZ062 6 ISTTXT(28), ITITLE(80), MBOX, MMTRY, NB, NGAMM1, ANLYZ063 7 NGAMP1, NLIN, NONL, NONLP1, NPLINE, NYSUM, ANLYZ064 8 NIN, NOUT ANLYZ065 COMMON /LBLOCK/ DOUSIN, LAST, PRWT, PRY, SAMET, SIMULA, ANLYZ066 1 DOADEX(10,2), DOSPL(10,2), DOSTRT(10,2), ANLYZ067 2 DOUSOU(2), LUSER(30), ANLYZ068 3 ALPLAM, ALPLIN, IWTIST, NSTORD, REDBET, REDCAP, ANLYZ069 4 SPLINE, WRTBET, WRTCAP ANLYZ070 DATA IHOLER/1HA, 1HN, 1HA, 1HL, 1HY, 1HZ/ ANLYZ071 C-----------------------------------------------------------------------ANLYZ072 C INITIAL COMPUTATIONS ANLYZ073 C-----------------------------------------------------------------------ANLYZ074 PRBEST=.FALSE. ANLYZ075 PRHEAD=.FALSE. ANLYZ076 DO 100 J=1,3 ANLYZ077 IF (IPRINT(J,ISTAGE).EQ.1 .OR. IPRINT(J,ISTAGE).EQ.3) ANLYZ078 1 PRBEST=.TRUE. ANLYZ079 IF (IPRINT(J,ISTAGE) .GE. 2) PRHEAD=.TRUE. ANLYZ080 100 CONTINUE ANLYZ081 IF (IPRITR(ISTAGE) .GE. 1) PRHEAD=.TRUE. ANLYZ082 IAPFIT=IABS(IPLFIT(ISTAGE)) ANLYZ083 IF (IAPFIT .EQ.1 .OR. IAPFIT .EQ.3) PRBEST=.TRUE. ANLYZ084 IF (IPLRES(ISTAGE).EQ.1 .OR. IPLRES(ISTAGE).EQ.3) PRBEST=.TRUE. ANLYZ085 IF (IAPFIT .GE. 2) PRHEAD=.TRUE. ANLYZ086 IF (IPLRES(ISTAGE) .GE. 2) PRHEAD=.TRUE. ANLYZ087 C-----------------------------------------------------------------------ANLYZ088 C START OF NNL DIFFERENT ANALYSIS, ASSUMING NL(L) NONLINEAR ANLYZ089 C COMPONENTS, L=1,...,NNL ANLYZ090 C-----------------------------------------------------------------------ANLYZ091 DO 700 L=1,NNL ANLYZ092 LOOP =L ANLYZ093 NONL =NL(L) ANLYZ094 NONLP1=NONL+1 ANLYZ095 NLIN =NONL+NGAM ANLYZ096 NLND =NLIN*ND ANLYZ097 NLNDNL=NLND+NONL ANLYZ098 IF (NGAM .LE. 1) NPLINE=MIN0(4,NONL) ANLYZ099 VARSAV(L)=RANGE ANLYZ100 SPLINE =DOSPL(L,ISTAGE) ANLYZ101 ADEXCT =DOADEX(L,ISTAGE) .AND. SPLINE ANLYZ102 ERRONE =.NOT.PRHEAD ANLYZ103 C-----------------------------------------------------------------------ANLYZ104 C COMPUTE MAXIMUM NUMBER OF TRIALS MMTRY, SO THAT ANLYZ105 C MBOX ANLYZ106 C MMTRY := ( ) .LE. MTRY(L,ISTAGE) ANLYZ107 C NONL ANLYZ108 C-----------------------------------------------------------------------ANLYZ109 MMTRY =MTRY(L,ISTAGE) ANLYZ110 MBOX =NONL ANLYZ111 NOVERK=1 ANLYZ112 200 NTEMP =NOVERK ANLYZ113 NOVERK=NPASCL(NTEMP,MBOX,NONL,1) ANLYZ114 IF (NOVERK .EQ. MMTRY) GO TO 220 ANLYZ115 IF (NOVERK .LT. MMTRY) GO TO 200 ANLYZ116 MBOX =MBOX-1 ANLYZ117 MMTRY=NTEMP ANLYZ118 220 DBOX=ZTOTAL/FLOAT(MBOX) ANLYZ119 DO 240 J=1,MMTRY ANLYZ120 240 STARTV(J)=.FALSE. ANLYZ121 OWNST=DOSTRT(L,ISTAGE) ANLYZ122 MORE =.FALSE. ANLYZ123 M =0 ANLYZ124 C-----------------------------------------------------------------------ANLYZ125 C START OF MMTRY ANALYSIS, ASSUMING NONL:=NL(L) NONLINEAR COMPONENTS ANLYZ126 C-----------------------------------------------------------------------ANLYZ127 IF (ISTAGE.EQ.1 .AND. PRHEAD) WRITE (NOUT,1000) ITITLE,NONL ANLYZ128 IF (ISTAGE.EQ.2 .AND. PRHEAD) WRITE (NOUT,1500) ITITLE,NONL ANLYZ129 300 M=M+1 ANLYZ130 IF (M .GT. MMTRY) GO TO 500 ANLYZ131 C-----------------------------------------------------------------------ANLYZ132 C LOOK FOR STARTING VALUES FOR PARAMETER PLAM ANLYZ133 C-----------------------------------------------------------------------ANLYZ134 IF (.NOT. OWNST) GO TO 330 ANLYZ135 OWNST=.FALSE. ANLYZ136 M =M-1 ANLYZ137 CALL USERST (ISTAGE,NYMAX,T,NY,NONLMX,PLAM) ANLYZ138 IF (DOUSOU(ISTAGE)) CALL USEROU (1,NYMAX,Y,T,YLYFIT,SQRTW, ANLYZ139 1 NY,NLINMX,NONLMX,PLIN,PLAM,PLMTRY,BUFF,VAR,ISTAGE,RESPON)ANLYZ140 C-----------------------------------------------------------------------ANLYZ141 C CHECK IF USER'S PLAM(J), J=1,...,NONL, ARE FEASIBLE STARTING VALUES ANLYZ142 C-----------------------------------------------------------------------ANLYZ143 DO 310 J=1,NONL ANLYZ144 IF (PLAM(J).LT.PNMNMX(1).OR.PLAM(J).GT.PNMNMX(2)) GO TO 300 ANLYZ145 310 CONTINUE ANLYZ146 DO 320 J=1,NONL ANLYZ147 320 NCOMBI(J)=0 ANLYZ148 GO TO 350 ANLYZ149 330 CALL NOKSUB (MBOX,NONL,NOKSET,LTEMPL,MTEMPM,MORE) ANLYZ150 IF (STARTV(M)) GO TO 300 ANLYZ151 DO 340 J=1,NONL ANLYZ152 I =NOKSET(J) ANLYZ153 NCOMBI(J)=I ANLYZ154 PLAM(J) =(FLOAT(I)-.5)*DBOX+PNMNMX(1) ANLYZ155 340 CONTINUE ANLYZ156 350 IERROR=0 ANLYZ157 ALPLAM=.TRUE. ANLYZ158 C-----------------------------------------------------------------------ANLYZ159 C SET FEASIBLE STARTING VALUES FOR PLIN ANLYZ160 C-----------------------------------------------------------------------ANLYZ161 DO 360 J=1,NLIN ANLYZ162 360 PLIN(J,1)=AMIN1(PLMNMX(2),AMAX1(0.E0,PLMNMX(1))) ANLYZ163 IF (ND .EQ. 1) GO TO 400 ANLYZ164 DO 370 N=2,ND ANLYZ165 DO 370 J=1,NLIN ANLYZ166 370 PLIN(J,N)=PLIN(J,1)