C++++++++++++++++ CD PACKAGE - VERSION 2DP (MAR 1984) ++++++++++++++++++ 0001 C CONTIN. MAIN SUBPROGRAM. 0002 C FOR THE REGULARIZED SOLUTION OF LINEAR ALGEBRAIC AND 0003 C LINEAR FREDHOLM INTEGRAL EQUATIONS OF THE FIRST KIND, WITH 0004 C OPTIONS FOR PEAK CONSTRAINTS AND LINEAR EQUALITY AND INEQUALITY 0005 C CONSTRAINTS. 0006 C REFERENCES - S.W. PROVENCHER (1982) COMPUT. PHYS. COMMUN., VOL. 27, 0007 C PAGES 213-227, 229-242. 0008 C (1984) CONTIN USERS MANUAL (EMBL 0009 C TECHNICAL REPORT DA07). 0010 C C AS A NEW USER YOU SHOULD NOTIFY S.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) 0017 C AND HAVE YOUR NAME PUT ON A MAILING LIST FOR UPDATES, ETC. 0018 C 0025 C----------------------------------------------------------------------- 0026 C CALLS SUBPROGRAMS - BLOCK DATA, INIT, INPUT, SETGRD, USERSI, 0027 C WRITYT, USERSX, USERNQ, SETNNG, ANALYZ, SETWT 0028 C WHICH IN TURN CALL - STORIN, READYT, ERRMES, WRITIN, USERIN, 0029 C USERGR, CQTRAP, USERTR, USEREX, RGAUSS, RANDOM, SETSCA, SEQACC, 0030 C H12, GETROW, USERK, USERLF, SETREG, USERRG, USEREQ, ELIMEQ, 0031 C LH1405, SVDRS2, QRBD, G1,G2, DIFF, DIAREG, DIAGA, SETGA1, 0032 C SETVAL, LDPETC, LDP, NNLS, CVNEQ, FISHNI, GAMLN, 0033 C BETAIN, PLPRIN, USEROU, MOMENT, MOMOUT, RUNRES, GETPRU, GETYLY, 0034 C PGAUSS, PLRES, SETSGN, ANPEAK, UPDSGN, UPDDON, FFLAT, UPDLLS, 0035 C USERWT 0036 C----------------------------------------------------------------------- 0037 DOUBLE PRECISION PRECIS, RANGE 0038 DOUBLE PRECISION A, AA, AEQ, AINEQ, PIVOT, REG, RHSNEQ, 0039 1 S, SOLBES, SOLUTN, SSCALE, VALPCV, VALPHA, VK1Y1, WORK 0040 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0041 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0042 LOGICAL LBIND 0043 C 0044 C*********************************************************************** 0045 C THE INSTRUCTIONS SET OFF BY ASTERISKS DESCRIBE ALL POSSIBLE CHANGES 0046 C THAT YOU MAY HAVE TO MAKE IN THIS MAIN SUBPROGRAM. (SEE ALSO THE 0047 C CHANGES IN THE BLOCK DATA AND USER SUBPROGRAMS.) THESE CHANGES 0048 C IN THE MAIN SUBPROGRAM ARE ONLY NECESSARY IF YOU CHANGE MY, MA, 0049 C MG, MREG, MINEQ, MEQ, MDONE, OR MWORK IN THE DATA STATEMENT 0050 C BELOW. IF YOU DO, THEN THE FOLLOWING DIMENSIONS MUST BE 0051 C READJUSTED AS DESCRIBED BELOW - 0052 C 0053 DIMENSION T(252), SQRTW(252), Y(252), EXACT(252), YLYFIT(252) 0054 DIMENSION G(18), CQUAD(18), VK1Y1(18), S(18,3), VALPHA(18), 0055 1 VALPCV(18), SOLUTN(18), IISIGN(18), SOLBES(18), 0056 2 AA(18,18), SSCALE(18) 0057 DIMENSION AINEQ(3,18), RHSNEQ(3), LBIND(3) 0058 DIMENSION A(18,18), IWORK(18) 0059 DIMENSION REG(17,18) 0060 DIMENSION AEQ(1,18), PIVOT(1) 0061 DIMENSION WORK(288) 0062 DIMENSION LSDONE(1,3,2), VDONE(1) 0063 C 0064 C THE ABOVE DIMENSION STATEMENTS MUST BE ADJUSTED AS FOLLOWS - 0065 C 0066 C DIMENSION T(MY), SQRTW(MY), Y(MY), EXACT(MY), YLYFIT(MY) 0067 C DIMENSION G(MG), CQUAD(MG), VK1Y1(MG), S(MG,3), VALPHA(MG), 0068 C 1 VALPCV(MG), SOLUTN(MG), IISIGN(MG), SOLBES(MG), 0069 C 2 AA(MG,MG), SSCALE(MG) 0070 C DIMENSION AINEQ(MINEQ,MG), RHSNEQ(MINEQ), LBIND(MINEQ) 0071 C DIMENSION A(MA,MG), IWORK(MA) 0072 C DIMENSION REG(MREG,MG) 0073 C DIMENSION AEQ(MEQ,MG), PIVOT(MEQ) 0074 C DIMENSION WORK(MWORK) 0075 C DIMENSION LSDONE(MDONE,3,2), VDONE(MDONE) 0076 C*********************************************************************** 0077 C 0078 COMMON /DBLOCK/ PRECIS, RANGE 0079 COMMON /SBLOCK/ DFMIN, SRMIN, 0080 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 0081 2 SRANGE 0082 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0083 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0084 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0085 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0086 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0087 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0088 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0089 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0090 2 LUSER(30) 0091 COMMON /SCD/ AINTRP(251,16), PSINV(6,15), WLEND, WLSTRT 0092 COMMON /ICD/ MIDEG1, MIPTS, MSTORE, NSPECT, NWAVEL 0093 C 0094 C*********************************************************************** 0095 C YOU CAN SAVE STORAGE BY MAKING THE INTEGERS IN THE FOLLOWING DATA 0096 C STATEMENT AS SMALL AS THE SIZE OF YOUR PROBLEM WILL ALLOW. (SEE 0097 C USERS MANUAL FOR THE MINIMUM ALLOWABLE VALUES.) 0098 C 0099 DATA MY/252/, MA/18/, MG/18/, MREG/17/, MINEQ/3/, MEQ/1/, 0100 1 MDONE/1/, MWORK/288/ 0101 C 0102 C IF YOU CHANGE THE ABOVE DATA STATEMENT, THEN THE DIMENSION 0103 C STATEMENTS ABOVE MUST BE READJUSTED, AS DESCRIBED ABOVE. 0104 C 0105 C THIS IS THE END OF ALL POSSIBLE CHANGES THAT YOU MIGHT HAVE TO MAKE 0106 C IN THE MAIN PROGRAM, 0107 C EXCEPT THAT IF YOUR SYSTEM DOES NOT AUTOMATICALLY 0108 C OPEN INPUT AND OUTPUT FILES FOR YOU, THEN YOU MIGHT HAVE TO OPEN 0109 C THEM HERE AND GIVE THEM THE NUMBERS NIN (FOR THE INPUT) AND 0110 C NOUT (FOR THE OUTPUT), WHERE NIN AND NOUT HAVE BEEN SET IN THE 0111 C BLOCK DATA SUBPROGRAM. 0112 C IN ADDITION, IF YOU ARE GOING TO INPUT IUNIT.GE.0, THEN YOU MAY 0113 C HAVE TO OPEN A TEMPORARY SCRATCH FILE NUMBERED IUNIT. DO THIS 0114 C DIRECTLY AFTER STATEMENT 100. THIS IS NOT NECESSARY IF IUNIT 0115 C IS NEGATIVE OR IF YOUR SYSTEM OPENS FILES AUTOMATICALLY. 0116 C*********************************************************************** 0117 C 0118 C----------------------------------------------------------------------- 0119 C INITIALIZE VARIABLES 0120 C----------------------------------------------------------------------- 0121 CALL INIT 0122 C----------------------------------------------------------------------- 0123 C READ INPUT DATA 0124 C----------------------------------------------------------------------- 0125 100 CALL INPUT (EXACT,G,MA,MEQ,MG,MINEQ,MREG,MWORK,MY,SQRTW,T,Y) 0126 C----------------------------------------------------------------------- 0127 C SET UP QUADRATURE GRID 0128 C----------------------------------------------------------------------- 0129 CALL SETGRD (CQUAD,G,GMNMX,IGRID,IQUAD,MG,NG,NOUT) 0130 C----------------------------------------------------------------------- 0131 C CALCULATE SIMULATED DATA 0132 C----------------------------------------------------------------------- 0133 IF (SIMULA) CALL USERSI (EXACT,G,MG,MY,SQRTW,T,Y) 0134 C----------------------------------------------------------------------- 0135 C WRITE OUT SIMULATED DATA 0136 C----------------------------------------------------------------------- 0137 IF (SIMULA) CALL WRITYT (EXACT,G,IPRINT,IUSROU,IWT,MG,NOUT,NY, 0138 1 PRY,SIMULA,SQRTW,T,Y) 0139 C----------------------------------------------------------------------- 0140 C PUT SECOND CURVE TO BE PLOTTED WITH SOLUTION IN EXACT. 0141 C----------------------------------------------------------------------- 0142 IF (.NOT.ONLY1) CALL USERSX (EXACT,G,MG) 0143 NINEQ=0 0144 NGL=NG+NLINF 0145 NGLP1=NGL+1 0146 C----------------------------------------------------------------------- 0147 C SET SPECIAL USER-SUPPLIED INEQUALITY CONSTRAINTS 0148 C----------------------------------------------------------------------- 0149 IF (DOUSNQ) CALL USERNQ (AINEQ,MG,MINEQ) 0150 C----------------------------------------------------------------------- 0151 C SET NG NONNEGATIVITY CONSTRAINTS AT ALL NG GRID POINTS 0152 C----------------------------------------------------------------------- 0153 IF (NONNEG) CALL SETNNG (AINEQ,MINEQ,NG,NGLP1,NINEQ,NOUT) 0154 IF (IWT.EQ.1 .OR. IWT.EQ.4) GO TO 200 0155 C----------------------------------------------------------------------- 0156 C DO A COMPLETE PRELIMINARY UNWEIGHTED ANALYSIS TO GET A SMOOTH FIT 0157 C TO THE DATA. THIS SMOOTH CURVE IS THEN USED TO CALCULATE THE 0158 C WEIGHTS. 0159 C----------------------------------------------------------------------- 0160 CALL ANALYZ (1, 0161 1 A,AA,AEQ,AINEQ,CQUAD,EXACT,G,IISIGN,IWORK,LBIND,LSDONE,MA, 0162 2 MDONE,MEQ,MG,MINEQ,MREG,MWORK,MY,PIVOT,REG,RHSNEQ,S,SOLBES, 0163 3 SOLUTN,SQRTW,SSCALE,T,VALPCV,VALPHA,VDONE,VK1Y1,WORK, 0164 4 Y,YLYFIT) 0165 C----------------------------------------------------------------------- 0166 C CALCULATE SQRTW (SQUARE ROOT OF LEAST SQUARES WEIGHTS). 0167 C----------------------------------------------------------------------- 0168 CALL SETWT ( 0169 1 CQUAD,G,IUNIT,IWT,MWORK,MY,NERFIT,NG,NGL,NLINF,NOUT,NY,PRWT, 0170 2 SOLBES,SQRTW,SRANGE,SSCALE,T,WORK,Y,YLYFIT) 0171 C----------------------------------------------------------------------- 0172 C DO FINAL WEIGHTED ANALYSIS. 0173 C----------------------------------------------------------------------- 0174 200 CALL ANALYZ (2, 0175 1 A,AA,AEQ,AINEQ,CQUAD,EXACT,G,IISIGN,IWORK,LBIND,LSDONE,MA, 0176 2 MDONE,MEQ,MG,MINEQ,MREG,MWORK,MY,PIVOT,REG,RHSNEQ,S,SOLBES, 0177 3 SOLUTN,SQRTW,SSCALE,T,VALPCV,VALPHA,VDONE,VK1Y1,WORK, 0178 4 Y,YLYFIT) 0179 IF (.NOT.LAST) GO TO 100 0180 STOP 0181 END 0182 C++++++++++++++++ CD PACKAGE - VERSION 2DP (MAR 1984) ++++++++++++++++++ 0183 BLOCK DATA 0184 DOUBLE PRECISION PRECIS, RANGE 0185 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0186 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0187 COMMON /DBLOCK/ PRECIS, RANGE 0188 COMMON /SBLOCK/ DFMIN, SRMIN, 0189 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 0190 2 SRANGE 0191 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0192 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0193 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0194 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0195 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0196 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0197 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0198 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0199 2 LUSER(30) 0200 COMMON /SCD/ AINTRP(251,16), PSINV(6,15), WLEND, WLSTRT 0201 COMMON /ICD/ MIDEG1, MIPTS, MSTORE, NSPECT, NWAVEL 0202 C*********************************************************************** 0203 C YOU MUST SET THE FOLLOWING 4 VARIABLES TO VALUES APPROPRIATE FOR 0204 C YOUR COMPUTER. (SEE USERS MANUAL.) 0205 C DATA RANGE/1.E35/, SRANGE/1.E35/, NIN/5/, NOUT/6/!SP 0206 DATA RANGE/1.D35/, SRANGE/1.E35/, NIN/5/, NOUT/6/ 0207 C*********************************************************************** 0208 C 0209 C 0210 C*********************************************************************** 0211 C CONTIN WILL USE THE VALUES OF THE CONTROL VARIABLES THAT ARE IN THE 0212 C DATA STATEMENTS BELOW UNLESS YOU INPUT DIFFERENT VALUES. TO SAVE 0213 C INPUTTING THESE OFTEN, YOU CAN CHANGE THE VALUES BELOW TO THOSE 0214 C THAT YOU WILL USUALLY USE. 0215 C THE FOLLOWING DATA STATEMENTS CONTAIN (IN ALPHABETICAL ORDER) THE 0216 C REAL, INTEGER, AND LOGICAL CONTROL VARIABLES, IN THAT ORDER. 0217 DATA ALPST/2*0./, DFMIN/3./, GMNMX/1., 16./, PLEVEL/4*.5/, 0218 1 RSVMNX/2*1., 2*0./, RUSER/13*0., 1., .001, 500., 84*0./, 0219 2 SRMIN/.01/ 0220 DATA ICRIT/2*1/, 0221 1 IFORMT/1H(, 1H5, 1HE, 1H1, 1H5, 1H., 1H6, 1H), 62*1H /, 0222 2 IFORMW/1H(, 1H5, 1HE, 1H1, 1H5, 1H., 1H6, 1H), 62*1H /, 0223 3 IFORMY/1H(, 1H5, 1HE, 1H1, 1H5, 1H., 1H6, 1H), 62*1H /, 0224 4 IGRID/1/, IPLFIT/2*2/, IPLRES/2*2/, IPRINT/2, 3/, IQUAD/1/, 0225 5 IUNIT/-1/, IUSER/15*0, 4, 7, 33*0/, IUSROU/2*3/, IWT/1/, 0226 6 LINEPG/60/,LSIGN/16*0/, MIOERR/5/, MOMNMX/2*0/, MPKMOM/0/, 0227 7 MQPITR/35/, NENDZ/2*0/, NEQ/0/, NERFIT/0/, NFLAT/8*0/, NG/16/, 0228 8 NINTT/1/, NLINF/0/, NNSGN/2*0/, NORDER/-1/, NQPROG/2*6/, 0229 9 NSGN/4*0/ 0230 DATA DOCHOS/.TRUE./, DOMOM/.FALSE./, DOUSIN/.TRUE./, 0231 1 DOUSNQ/.TRUE./, LAST/.TRUE./, 0232 2 LUSER/30*.FALSE./, NEWPG1/.FALSE./, NONNEG/.FALSE./, 0233 3 ONLY1/.TRUE./, PRWT/.TRUE./, PRY/.TRUE./, 0234 4 SIMULA/.FALSE./ 0235 C*********************************************************************** 0236 C 0237 C*********************************************************************** 0238 C IF YOU MAKE ANY CHANGES TO CONTIN, YOU SHOULD PUT A NAME (UP TO 6 0239 C CHARACTERS) IN IAPACK TO UNIQUELY IDENTIFY YOUR VERSION OF 0240 C CONTIN. THIS WILL BE PRINTED IN THE HEADING OF VARIOUS PARTS OF 0241 C THE OUTPUT. YOU MUST SPECIFY THE NAME IN THE FOLLOWING STATEMENT 0242 DATA IAPACK/1H ,1HC, 1HD, 1H-, 1H1, 1H / 0243 C*********************************************************************** 0244 C 0245 C 0246 C----------------------------------------------------------------------- 0247 C INITIALIZATION OF SPECIAL COMMON BLOCKS FOR CD PACKAGE - 0248 C SEE THE USERS MANUAL (SEC 6.5) FOR A DESCRIPTION OF THE CD 0249 C PACKAGE. 0250 C----------------------------------------------------------------------- 0251 DATA WLEND/190./, WLSTRT/240./ 0252 DATA MIDEG1/6/, MIPTS/15/, MSTORE/251/, NSPECT/16/, NWAVEL/51/ 0253 END 0254 C++++++++++++++++ CD PACKAGE - VERSION 2DP (MAR 1984) ++++++++++++++++++ 0255 C SUBROUTINE USERIN. THIS IS A USER-SUPPLIED ROUTINE (ONLY CALLED 0256 C WHEN DOUSIN=.TRUE.) THAT IS CALLED RIGHT AFTER THE INITIALIZATION 0257 C AND INPUT OF THE COMMON VARIABLES, T, Y, AND THE LEAST-SQUARES 0258 C WEIGHTS. THEREFORE, IT CAN BE USED TO MODIFY ANY OF THESE. 0259 C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 0260 C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 0261 C SEE THE USERS MANUAL FOR A DESCRIPTION OF THE CD PACKAGE. 0262 SUBROUTINE USERIN (T,Y,SQRTW,MY) 0263 DOUBLE PRECISION PRECIS, RANGE 0264 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0265 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0266 DIMENSION T(MY), Y(MY), SQRTW(MY) 0267 COMMON /DBLOCK/ PRECIS, RANGE 0268 COMMON /SBLOCK/ DFMIN, SRMIN, 0269 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 0270 2 SRANGE 0271 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0272 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0273 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0274 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0275 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0276 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0277 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0278 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0279 2 LUSER(30) 0280 C----------------------------------------------------------------------- 0281 C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 0282 C FOR YOUR APPLICATION. 0283 C----------------------------------------------------------------------- 0284 NY=NY+1 0285 IF (NY.LE.MY .AND. AMIN1(ABS(RUSER(14)),RUSER(15),RUSER(16)) 0286 1.GT.0.) GO TO 110 0287 5110 FORMAT (45H0NY OR RUSER(J), J=14, 15, OR 16 OUT OF RANGE, 0288 1 10H IN USERIN,I10,1P3E13.3) 0289 WRITE (NOUT,5110) NY,(RUSER(J),J=14,16) 0290 STOP 0291 110 Y(NY)=RUSER(14) 0292 SQRTW(NY)=(RUSER(16)/RUSER(15))**2 0293 T(NY)=0. 0294 LUSER(2)=.FALSE. 0295 LUSER(3)=.FALSE. 0296 800 RETURN 0297 END 0298 C++++++++++++++++ CD PACKAGE - VERSION 2DP (MAR 1984) ++++++++++++++++++ 0299 C FUNCTION USERK. THIS IS A USER-SUPPLIED ROUTINE (ALWAYS NEEDED) 0300 C TO COMPUTE THE FREDHOLM KERNEL, USERK, WHICH DEPENDS ON T(JT) 0301 C (THE INDEPENDENT VARIABLE AT DATA POINT JT) AND G(JG) (THE 0302 C VALUE OF THE JG TH GRID POINT. 0303 C AUTHORS: J. GLOECKNER AND S.W. PROVENCHER 0304 C SEE THE USERS MANUAL FOR A DESCRIPTION OF THE CD PACKAGE. 0305 C----------------------------------------------------------------------- 0306 C CALLS SUBPROGRAMS - ERRMES, SVDRS2 0307 C WHICH IN TURN CALL - H12, QRBD, G1, G2, DIFF 0308 C----------------------------------------------------------------------- 0309 FUNCTION USERK (JT,T,JG,G) 0310 DOUBLE PRECISION PRECIS, RANGE 0311 DOUBLE PRECISION XD, ED, SD 0312 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0313 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0314 DIMENSION T(JT), G(JG) 0315 DIMENSION IHOLER(6) 0316 DIMENSION XD(15,6), ED(15,15), SD(6,3), AREF(51,16) 0317 DIMENSION SPEC01(51),SPEC02(51),SPEC03(51),SPEC04(51) 0318 DIMENSION SPEC05(51),SPEC06(51),SPEC07(51),SPEC08(51) 0319 DIMENSION SPEC09(51),SPEC10(51),SPEC11(51),SPEC12(51) 0320 DIMENSION SPEC13(51),SPEC14(51),SPEC15(51),SPEC16(51) 0321 EQUIVALENCE (SPEC01(1),AREF(1,1)),(SPEC02(1),AREF(1,2)) 0322 EQUIVALENCE (SPEC03(1),AREF(1,3)),(SPEC04(1),AREF(1,4)) 0323 EQUIVALENCE (SPEC05(1),AREF(1,5)),(SPEC06(1),AREF(1,6)) 0324 EQUIVALENCE (SPEC07(1),AREF(1,7)),(SPEC08(1),AREF(1,8)) 0325 EQUIVALENCE (SPEC09(1),AREF(1,9)),(SPEC10(1),AREF(1,10)) 0326 EQUIVALENCE (SPEC11(1),AREF(1,11)),(SPEC12(1),AREF(1,12)) 0327 EQUIVALENCE (SPEC13(1),AREF(1,13)),(SPEC14(1),AREF(1,14)) 0328 EQUIVALENCE (SPEC15(1),AREF(1,15)),(SPEC16(1),AREF(1,16)) 0329 COMMON /DBLOCK/ PRECIS, RANGE 0330 COMMON /SBLOCK/ DFMIN, SRMIN, 0331 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 0332 2 SRANGE 0333 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0334 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0335 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0336 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0337 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0338 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0339 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0340 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0341 2 LUSER(30) 0342 COMMON /SCD/ AINTRP(251,16), PSINV(6,15), WLEND, WLSTRT 0343 COMMON /ICD/ MIDEG1, MIPTS, MSTORE, NSPECT, NWAVEL 0344 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HK, 1H / 0345 C 1 MYOGLOBIN 0346 DATA SPEC01 / 0347 1 -2290., -2860., -3840., -4780., -5890., -7140., -8800., 0348 2 -10400.,-12300.,-14100.,-15700.,-17900.,-19800.,-21400., 0349 3 -22500.,-23800.,-24300.,-24900.,-24900.,-24700.,-24400., 0350 4 -24000.,-23500.,-22800.,-22300.,-21700.,-21700.,-21900., 0351 5 -22100.,-22600.,-22800.,-22900.,-21800.,-20000.,-17600., 0352 6 -14000.,-10000., -3860., 1170., 8660., 16200., 25100., 0353 7 32900., 40100., 47000., 51100., 53100., 54900., 55100., 0354 8 54400., 52900./ 0355 C 2 LYSOZYME 0356 DATA SPEC02 / 0357 1 -1400., -1760., -2210., -2710., -3290., -4060., -4910., 0358 2 -5820., -6540., -7170., -7660., -8110., -8430., -8750., 0359 3 -9020., -9290., -9510., -9650., -9830., -9960.,-10050., 0360 4 -10100.,-10190.,-10230.,-10410.,-10680.,-10950.,-11360., 0361 5 -11900.,-12400.,-12850.,-13210.,-13430.,-13250.,-12620., 0362 6 -11630.,-10190., -8300., -5230., -1000., 2040., 4910., 0363 7 7360., 9820., 12680., 15140., 16360., 17390., 17590., 0364 8 16770., 15960./ 0365 C 3 RIBONUCLEASE A 0366 DATA SPEC03 / 0367 1 180., 130., 50., -100., -380., -670., -1050., 0368 2 -1540., -2030., -2590., -3290., -4000., -4750., -5491., 0369 3 -6240., -6970., -7640., -8310., -8850., -9340., -9720., 0370 4 -10000.,-10260.,-10410.,-10540.,-10650.,-10830.,-10950., 0371 5 -11030.,-11110.,-11080.,-10900.,-10520., -9880., -8980., 0372 6 -7900., -6520., -5380., -3310., -1390., 470., 2160., 0373 7 3500., 4720., 5660., 6300., 6470., 6410., 6180., 0374 8 5770., 5730./ 0375 C 4 PAPAIN 0376 DATA SPEC04 / 0377 1 -1050., -1670., -1970., -2430., -2870., -3460., -4160., 0378 2 -5040., -6050., -7140., -8150., -9290.,-10100.,-10780., 0379 3 -11290.,-11640.,-11770.,-11830.,-11700.,-11550.,-11370., 0380 4 -11200.,-11070.,-11040.,-11040.,-11130.,-11310.,-11590., 0381 5 -11900.,-12340.,-12580.,-12620.,-12320.,-11480.,-10470., 0382 6 -9270., -8784., -8385., -8185., -8490., -8880., -8390., 0383 7 -6990., -4940., -2300., 900., 4290., 7590., 10480., 0384 8 12930., 14820./ 0385 C 5 LACTATE DEHYDROGENASE 0386 DATA SPEC05 / 0387 1 -1790., -2220., -2730., -3370., -4010., -4710., -5540., 0388 2 -6370., -7270., -8210., -9170.,-10110.,-10960.,-11660., 0389 3 -12270.,-12780.,-13080.,-13210.,-13240.,-13130.,-12940., 0390 4 -12730.,-12490.,-12220.,-12060.,-11930.,-11870.,-11900., 0391 5 -11980.,-12140.,-12330.,-12030.,-11340.,-10350., -8990., 0392 6 -7110., -4870., -2350., -160., 3780., 7070., 10960., 0393 7 14740., 18150., 20890., 22530., 23750., 24120., 24000., 0394 8 23390., 22170./ 0395 C 6 ALPHA-CHYMOTRYPSIN 0396 DATA SPEC06 / 0397 1 -670., -806., -920., -1090., -1400., -1840., -2400., 0398 2 -3170., -3990., -4570., -4780., -4600., -4390., -4260., 0399 3 -4140., -4110., -4110., -4160., -4260., -4360., -4490., 0400 4 -4660., -4850., -5060., -5300., -5540., -5830., -6200., 0401 5 -6520., -6890., -7270., -7670., -8080., -8490., -8960., 0402 6 -9320., -9670., -9820., -9730., -9260., -8830., -7250., 0403 7 -6204., -4980., -3670., -2530., -1220., 870., 2270., 0404 8 3410., 4110./ 0405 C 7 CONCANAVALIN A 0406 DATA SPEC07 / 0407 1 -790., -910., -1080., -1300., -1560., -1870., -2210., 0408 2 -2660., -3170., -3880., -4590., -5240., -5780., -6170., 0409 3 -6480., -6680., -6710., -6600., -6370., -6170., -5890., 0410 4 -5550., -5180., -4760., -4300., -3850., -3340., -2720., 0411 5 -2040., -1250., -450., 260., 1100., 2000., 2970., 0412 6 3870., 4902., 6000., 7220., 8390., 9680., 10970., 0413 7 11680., 11930., 11870., 11480., 10450., 8640., 6450., 0414 8 4640., 2770./ 0415 C 8 CYTOCHROME C 0416 DATA SPEC08 / 0417 1 -710., -1030., -1330., -1730., -2250., -2920., -3680., 0418 2 -4620., -5780., -6850., -8080., -9060.,-10000.,-10640., 0419 3 -11480.,-12010.,-12090.,-12240.,-12240.,-12160.,-11940., 0420 4 -11640.,-11050.,-10650.,-10270.,-10030., -9820., -9720., 0421 5 -9730., -9830., -9880., -9620., -9070., -8230., -7000., 0422 6 -5390., -3020., 0., 2880., 5350., 7540., 9600., 0423 7 11240., 12610., 13030., 12890., 12340., 11790., 10970., 0424 8 10280., 9460./ 0425 C 9 ELASTASE 0426 DATA SPEC09 / 0427 1 -420., -500., -630., -750., -800., -840., -840., 0428 2 -840., -750., -840., -920., -1000., -1170., -1340., 0429 3 -1550., -1760., -2010., -2260., -2510., -2800., -3140., 0430 4 -3470., -3810., -4090., -4360., -4830., -5310., -5750., 0431 5 -6200., -6700., -7210., -7650., -8080., -8480., -8710., 0432 6 -8890., -9530.,-10490.,-11250.,-12010.,-12390.,-12580., 0433 7 -12110.,-11440.,-10490., -9340., -8010., -6670., -4770., 0434 8 -3430., -2100./ 0435 C 10 NUCLEASE 0436 DATA SPEC10 / 0437 1 -1160., -1410., -1780., -2140., -2610., -3090., -3700., 0438 2 -4370., -5140., -6000., -6770., -7540., -8290., -8910., 0439 3 -9450., -9960.,-10300.,-10600.,-10700.,-10800.,-10800., 0440 4 -10800.,-10700.,-10500.,-10300.,-10200.,-10200.,-10300., 0441 5 -10500.,-10700.,-11110.,-11300.,-11100.,-10700.,-10100., 0442 6 -9240., -8200., -6880., -5100., -3380., -1070., 1650., 0443 7 4120., 6090., 7980., 9540., 10600., 11300., 11900., 0444 8 12100., 12400./ 0445 C 11 INSULIN 0446 DATA SPEC11 / 0447 1 -930., -1280., -1580., -1940., -2460., -3200., -3920., 0448 2 -4520., -5290., -6130., -6970., -7940., -8680., -9530., 0449 3 -10300.,-10900.,-11100.,-11300.,-11300.,-11300.,-11200., 0450 4 -11100.,-11100.,-11200.,-11400.,-11500.,-11800.,-12200., 0451 5 -12800.,-13100.,-13300.,-13900.,-13600.,-12400.,-10800., 0452 6 -8260., -5570., -1470., 2200., 6610., 11800., 17600., 0453 7 22000., 25000., 25700., 25800., 23500., 20600., 16900., 0454 8 15300., 12200./ 0455 C 12 PARVALBUMIN 0456 DATA SPEC12 / 0457 1 -1769., -2239., -2739., -3364., -4066., -4892., -5748., 0458 2 -6757., -7758., -8918.,-10009.,-11042.,-12003.,-12815., 0459 3 -13483.,-13997.,-14368.,-14535.,-14511.,-14377.,-14199., 0460 4 -13982.,-13714.,-13435.,-13151.,-12892.,-12853.,-12978., 0461 5 -13256.,-13712.,-14192.,-14376.,-14102.,-13318.,-11779., 0462 6 -9595., -6655., -2858., 457., 5432., 10141., 14251., 0463 7 17149., 18808., 19606., 19616., 19057., 17516., 14801., 0464 8 10226., 6005./ 0465 C 13 CARBOXYPEPTIDASE A 0466 DATA SPEC13 / 0467 1 -1445., -1826., -2282., -2764., -3296., -3905., -4487., 0468 2 -5045., -5502., -5856., -6136., -6389., -6643., -6790., 0469 3 -6972., -7124., -7302., -7429., -7631., -7859., -8113., 0470 4 -8417., -8721., -8950., -9229., -9431., -9634., -9761., 0471 5 -9888., -9938., -9508., -8772., -7708., -5983., -3651., 0472 6 -534., 2368., 5255., 7853., 10568., 12646., 13917., 0473 7 14783., 14725., 13917., 12820., 11549., 9875., 8604., 0474 8 7218., 5775./ 0475 C 14 TRYPSIN INHIBITOR 0476 DATA SPEC14 / 0477 1 -876., -986., -1150., -1388., -1625., -1954., -2300., 0478 2 -2702., -3140., -3670., -4200., -4765., -5350., -5952., 0479 3 -6518., -7066., -7559., -8033., -8398., -8782., -9074., 0480 4 -9312., -9530., -9713., -9859.,-10060.,-10297.,-10735., 0481 5 -11174.,-11612.,-12540.,-13724.,-14680.,-16052.,-17050., 0482 6 -17840.,-18464.,-18713.,-18380.,-17549.,-16468.,-15304., 0483 7 -12434., -9731., -7153., -4492., -2661., -416., 2412., 0484 8 4076., 4991./ 0485 C 15 ADENYLATE KINASE 0486 DATA SPEC15 / 0487 1 -1445., -1772., -2208., -2780., -3434., -4224., -5069., 0488 2 -5968., -6840., -7794., -8775., -9729.,-10492.,-11255., 0489 3 -11990.,-12481.,-12971.,-13326.,-13462.,-13435.,-13407., 0490 4 -13353.,-13298.,-13271.,-13244.,-13244.,-13271.,-13407., 0491 5 -13570.,-13843.,-14007.,-13898.,-13244.,-11936.,-10192., 0492 6 -8548., -6104., -3488., 3104., 6518., 10429., 13905., 0493 7 17629., 20858., 23403., 25637., 26940., 27934., 28306., 0494 8 27996., 27500./ 0495 C 16 RIBONUCLEASE S 0496 DATA SPEC16 / 0497 1 173., 96., 0., -211., -481., -788., -1153., 0498 2 -1615., -2173., -2854., -3364., -4056., -4787., -5555., 0499 3 -6267., -7036., -7670., -8324., -8900., -9304., -9727., 0500 4 -10073.,-10304.,-10457.,-10573.,-10688.,-10841.,-10976., 0501 5 -11034.,-11092.,-10995.,-10842.,-10419., -9688., -8727., 0502 6 -7498., -5998., -4268., -2364., -884., 2014., 3809., 0503 7 5035., 5779., 6348., 6567., 6567., 6260., 5867., 0504 8 5298., 4597./ 0505 IF (JT.GT.NY .OR. JG.GT.NG .OR. MIN0(JT,JG).LE.0) CALL ERRMES (1, 0506 1.TRUE.,IHOLER,NOUT) 0507 C----------------------------------------------------------------------- 0508 C THE FOLLOWING STATEMENTS SHOULD BE REPLACED BY THOSE 0509 C APPROPRIATE FOR YOUR KERNEL. 0510 C FOR EXAMPLE, FOR THE LAPLACE INTEGRAL EQUATION, YOU WOULD 0511 C SIMPLY WRITE - 0512 C USERK=EXP(-T(JT)*G(JG)) 0513 C----------------------------------------------------------------------- 0514 IF (JT .NE. NY) GO TO 110 0515 C----------------------------------------------------------------------- 0516 C THE LAST DATA POINT IS SIMPLY THE SUM OF THE GAMMAS SET EQUAL TO 0517 C RUSER(14). 0518 C----------------------------------------------------------------------- 0519 USERK=1. 0520 LUSER(3)=.TRUE. 0521 RETURN 0522 110 IF (LUSER(3)) GO TO 300 0523 IDEG1=IUSER(16) 0524 IF (IDEG1.LE.1 .OR. IDEG1.GT.MIDEG1) CALL ERRMES (2,.TRUE., 0525 1IHOLER, NOUT) 0526 IPTS=IUSER(17) 0527 IF (IPTS.LT.IDEG1 .OR. IPTS.GT.MIPTS) CALL ERRMES (3,.TRUE., 0528 1 IHOLER,NOUT) 0529 IF (LUSER(2)) GO TO 200 0530 LUSER(2)=.TRUE. 0531 IF (NY-1 .GT. MSTORE) CALL ERRMES (4,.TRUE.,IHOLER,NOUT) 0532 IF (NG .GT. NSPECT) CALL ERRMES (5,.TRUE.,IHOLER,NOUT) 0533 C----------------------------------------------------------------------- 0534 C COMPUTE PSEUDOINVERSE FOR LEAST SQUARES SMOOTHING AND INTERPOLATION. 0535 C----------------------------------------------------------------------- 0536 DUM=-FLOAT(IPTS/2) 0537 DO 120 IROW=1,IPTS 0538 XD(IROW,1)=1. 0539 DO 130 ICOL=2,IDEG1 0540 XD(IROW,ICOL)=XD(IROW,ICOL-1)*DUM 0541 130 CONTINUE 0542 DUM=DUM+1. 0543 DO 140 ICOL=1,IPTS 0544 ED(IROW,ICOL)=0. 0545 140 CONTINUE 0546 ED(IROW,IROW)=1. 0547 120 CONTINUE 0548 CALL SVDRS2 (XD,MIPTS,IPTS,IDEG1,ED,MIPTS,IPTS,SD,IERR,RANGE) 0549 IF (IERR .NE. 1) CALL ERRMES (6,.TRUE.,IHOLER,NOUT) 0550 DO 150 ICOL=1,IDEG1 0551 DUM=SD(ICOL,1) 0552 IF (DUM .LE. 0.) CALL ERRMES (7,.TRUE.,IHOLER,NOUT) 0553 DUM=1./DUM 0554 DO 160 IROW=1,IDEG1 0555 XD(IROW,ICOL)=XD(IROW,ICOL)*DUM 0556 160 CONTINUE 0557 150 CONTINUE 0558 DO 170 IROW=1,IDEG1 0559 DO 180 ICOL=1,IPTS 0560 PSINV(IROW,ICOL)=0. 0561 DO 190 J=1,IDEG1 0562 PSINV(IROW,ICOL)=PSINV(IROW,ICOL)+XD(IROW,J)*ED(J,ICOL) 0563 190 CONTINUE 0564 180 CONTINUE 0565 170 CONTINUE 0566 C----------------------------------------------------------------------- 0567 C DO INTERPOLATION. 0568 C----------------------------------------------------------------------- 0569 200 TINT=(T(JT)-WLSTRT)*FLOAT(NWAVEL-1)/(WLEND-WLSTRT)+1. 0570 IF (TINT.LT..9 .OR. TINT.GT.FLOAT(NWAVEL)+.1) CALL ERRMES (8, 0571 1 .TRUE.,IHOLER,NOUT) 0572 ISTART=MIN0(NWAVEL+1-IPTS,MAX0(1,INT(TINT+.5)-IPTS/2)) 0573 TREL=TINT-FLOAT(ISTART+IPTS/2) 0574 SUM=0. 0575 POWER=1. 0576 DO 220 J=1,IDEG1 0577 DUM=0. 0578 I=ISTART 0579 DO 230 II=1,IPTS 0580 DUM=DUM+PSINV(J,II)*AREF(I,JG) 0581 I=I+1 0582 230 CONTINUE 0583 SUM=SUM+DUM*POWER 0584 POWER=POWER*TREL 0585 220 CONTINUE 0586 AINTRP(JT,JG)=SUM 0587 300 USERK=AINTRP(JT,JG) 0588 RETURN 0589 END 0590 C++++++++++++++++ CD PACKAGE - VERSION 2DP (MAR 1984) ++++++++++++++++++ 0591 C SUBROUTINE USERNQ. THIS IS A USER-SUPPLIED ROUTINE (ONLY 0592 C CALLED WHEN DOUSNQ IS INPUT AS .TRUE.) TO SET NINEQ 0593 C (THE NO. OF INEQUALITY CONSTRAINTS) AND TO PUT THE 0594 C INEQUALITY-CONSTRAINT MATRIX IN THE FIRST NINEQ 0595 C ROWS OF AINEQ AND TO PUT THE RIGHT HAND SIDES OF THE 0596 C INEQUALITIES IN COLUMN NGLP1 OF AINEQ. 0597 C I.E., (SUM FROM J=1 TO NGL OF (AINEQ(I,J)*SOLUTION(J))) .GE. 0598 C AINEQ(I,NGLP1), I=1,NINEQ. 0599 C NOTE - IF THE SOLUTION IS TO BE NONNEGATIVE AT ALL NG GRID 0600 C POINTS, DO NOT USE USERNQ TO SET THESE CONSTRAINTS - 0601 C SIMPLY INPUT NONNEG=.TRUE. INSTEAD. 0602 C SEE THE USERS MANUAL (SEC 4.2) FOR THE POSITION IN THE INPUT 0603 C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 0604 C SEE THE USERS MANUAL (SEC 6.5) FOR A DESCRIPTION OF THE CD 0605 C PACKAGE. 0606 C----------------------------------------------------------------------- 0607 C CALLS SUBPROGRAMS - ERRMES 0608 C----------------------------------------------------------------------- 0609 SUBROUTINE USERNQ (AINEQ,MG,MINEQ) 0610 DOUBLE PRECISION PRECIS, RANGE 0611 DOUBLE PRECISION AINEQ, ZERO 0612 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0613 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0614 DIMENSION AINEQ(MINEQ,MG) 0615 DIMENSION FCAP(3,16) 0616 DIMENSION IHOLER(6) 0617 COMMON /DBLOCK/ PRECIS, RANGE 0618 COMMON /SBLOCK/ DFMIN, SRMIN, 0619 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 0620 2 SRANGE 0621 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0622 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0623 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0624 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0625 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0626 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0627 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0628 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0629 2 LUSER(30) 0630 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HN, 1HQ/ 0631 DATA NCLASS/3/, FCAP/ 0632 1 .79,.00,.21 , .41,.16,.43 , .23,.40,.37 , 0633 2 .28,.14,.58 , .45,.24,.31 , .09,.34,.57 , 0634 3 .02,.51,.47 , .39,.00,.61 , .07,.52,.41 , 0635 4 .24,.15,.61 , .51,.24,.25 , .62,.05,.33 , 0636 5 .37,.15,.48 , .28,.33,.39 , 0637 6 .54,.12,.34 , .26,.44,.30 / 0638 C ZERO=0.E0!SP 0639 ZERO=0.D0 0640 C----------------------------------------------------------------------- 0641 C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 0642 C FOR YOUR APPLICATION. 0643 C----------------------------------------------------------------------- 0644 NINEQ=NCLASS 0645 IF (NINEQ .GT. MINEQ) CALL ERRMES (1,.TRUE.,IHOLER,NOUT) 0646 DO 110 J=1,NINEQ 0647 DO 120 K=1,NG 0648 AINEQ(J,K)=FCAP(J,K) 0649 120 CONTINUE 0650 AINEQ(J,NGLP1)=ZERO 0651 110 CONTINUE 0652 RETURN 0653 END 0654 C++++++++++++++++ CD PACKAGE - VERSION 2DP (MAR 1984) ++++++++++++++++++ 0655 C SUBROUTINE USERRG. THIS IS A USER-SUPPLIED ROUTINE (NEEDED 0656 C WHEN NORDER IS INPUT NEGATIVE) TO SET NREG AND TO PUT A 0657 C SPECIAL USER-DEFINED REGULARIZOR IN THE FIRST NREG COLUMNS 0658 C AND ROWS OF REG AND THE RIGHT-HAND-SIDE (R.H.S.) OF THE 0659 C REGULARIZOR IN COLUMN NGLP1 OF REG. 0660 C NOTE - IF WEIGHTS ARE TO BE CALCULATED (I.E., IF IWT=2,3, OR 5), 0661 C THEN USERRG IS CALLED TWICE. THEREFORE IT IS BEST TO HAVE ANY 0662 C INPUT DATA NEEDED BY USERRG READ IN ONLY ONCE AND STORED, E.G., 0663 C IN RUSER, AS ILLUSRATED BELOW. OTHERWISE THIS DATA WOULD HAVE 0664 C TO BE READ TWICE. 0665 C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 0666 C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 0667 C SEE THE USERS MANUAL FOR A DESCRIPTION OF THE CD PACKAGE. 0668 SUBROUTINE USERRG (REG,MREG,MG,NREG) 0669 DOUBLE PRECISION PRECIS, RANGE 0670 DOUBLE PRECISION REG 0671 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0672 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0673 DIMENSION REG(MREG,MG) 0674 COMMON /DBLOCK/ PRECIS, RANGE 0675 COMMON /SBLOCK/ DFMIN, SRMIN, 0676 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 0677 2 SRANGE 0678 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0679 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0680 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0681 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0682 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0683 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0684 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0685 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0686 2 LUSER(30) 0687 C----------------------------------------------------------------------- 0688 C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 0689 C FOR YOUR APPLICATION. 0690 C----------------------------------------------------------------------- 0691 NREG=NG 0692 DO 110 J=1,NREG 0693 DO 120 K=1,NGL 0694 REG(J,K)=0. 0695 120 CONTINUE 0696 REG(J,J)=1. 0697 REG(J,NGLP1)=RUSER(14)/FLOAT(NG-NEQ) 0698 110 CONTINUE 0699 RETURN 0700 END 0701 C++++++++++++++++ CD PACKAGE - VERSION 2DP (MAR 1984) ++++++++++++++++++ 0702 C SUBROUTINE USERWT. THIS IS A USER-SUPPLIED ROUTINE (ONLY 0703 C NEEDED WHEN IWT IS INPUT AS 5) FOR CALCULATING SQRTW (SQUARE 0704 C ROOTS OF THE LEAST SQUARES WEIGHTS) FROM Y, YLYFIT, AND 0705 C ERRFIT, AS EXPLAINED IN DETAIL IN THE USERS MANUAL. 0706 C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 0707 C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 0708 C SEE THE USERS MANUAL FOR A DESCRIPTION OF THE CD 0709 C PACKAGE. 0710 SUBROUTINE USERWT (Y,YLYFIT,MY,ERRFIT,SQRTW) 0711 DOUBLE PRECISION PRECIS, RANGE 0712 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0713 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0714 DIMENSION Y(MY), YLYFIT(MY), SQRTW(MY) 0715 COMMON /DBLOCK/ PRECIS, RANGE 0716 COMMON /SBLOCK/ DFMIN, SRMIN, 0717 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 0718 2 SRANGE 0719 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0720 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0721 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0722 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0723 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0724 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0725 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0726 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0727 2 LUSER(30) 0728 C----------------------------------------------------------------------- 0729 C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 0730 C FOR YOUR APPLICATION. 0731 C----------------------------------------------------------------------- 0732 NNY=NY-1 0733 IF (IUSER(14).GE.1 .AND. IUSER(14).LE.NNY) GO TO 120 0734 5120 FORMAT (/34H IUSER(14) OUT OF RANGE IN USERWT.) 0735 WRITE (NOUT,5120) 0736 STOP 0737 120 S1=0. 0738 S2=0. 0739 DO 130 J=1,NNY 0740 IF (J .LE. IUSER(14)) S1=S1+YLYFIT(J)**2 0741 IF (J .GT. IUSER(14)) S2=S2+YLYFIT(J)**2 0742 130 CONTINUE 0743 E1=SQRT(S1/FLOAT(IUSER(14))) 0744 E2=SRANGE 0745 IF (IUSER(14) .LT. NNY) E2=SQRT(S2/FLOAT(NNY-IUSER(14))) 0746 5130 FORMAT (//27H RMS RESIDUAL FOR PTS. 1 TO,I3,2H =,1PE11.2/ 0747 1 34H RMS RESIDUAL FOR REMAINING PTS. =,E9.2) 0748 WRITE (NOUT,5130) IUSER(14),E1,E2 0749 IF (FLOAT(IUSER(15))*(E1-E2).GE.0. .AND. AMIN1(E1,E2).GT.0.) 0750 1 GO TO 150 0751 E1=SQRT((S1+S2)/FLOAT(NNY)) 0752 E2=E1 0753 150 IF (AMAX1(E1,E2) .GT. 0.) GO TO 170 0754 5150 FORMAT (54H0ALL RESIDUALS SOMEHOW ZERO. CANNOT ESTIMATE WEIGHTS.) 0755 WRITE (NOUT,5150) 0756 STOP 0757 170 E1=1./E1 0758 E2=1./E2 0759 DO 180 J=1,NNY 0760 IF (J .LE. IUSER(14)) SQRTW(J)=E1 0761 IF (J .GT. IUSER(14)) SQRTW(J)=E2 0762 180 CONTINUE 0763 SQRTW(NY)=1./RUSER(15) 0764 RETURN 0765 END 0766 C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0248 C SUBROUTINE USEREQ. THIS IS A USER-SUPPLIED ROUTINE (NEEDED 0249 C WHEN NEQ IS INPUT POSITIVE) TO PUT THE EQUALITY-CONSTRAINT 0250 C MATRIX IN THE FIRST NEQ ROWS AND NGL COLUMNS OF AEQ, AND 0251 C THE RIGHT-HAND-SIDE OF THE EQUALITIES IN COLUMN NGLP1=NGL+1. 0252 C CQUAD CONTAINS THE COEFFICIENTS OF THE QUADRATURE FORMULA THAT WERE 0253 C SET IN SETGRD. 0254 C NOTE - IF WEIGHTS ARE TO BE CALCULATED (I.E., IF IWT=2,3, OR 5), THEN 0255 C USEREQ WILL BE CALLED TWICE. THEREFORE IT IS BEST TO HAVE ANY 0256 C DATA NEEDED BY USEREQ READ IN ONLY ONCE AND STORED, E.G., IN 0257 C RUSER, AS ILLUSTRATED BELOW. 0258 C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 0259 C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 0260 C BELOW IS ILLUSTRATED THE CASE WHERE THE END POINTS AND THE 0261 C INTEGRAL OVER THE SOLUTION CAN BE CONSTRAINED - 0262 C IF NEQ=1, THEN THE SOLUTION IS CONSTRAINED TO BE RUSER(1) AT 0263 C THE GRID POINT NG (THE LAST POINT). 0264 C IF NEQ=2, THEN, IN ADDITION, THE SOLUTION IS CONSTRAINED TO 0265 C RUSER(2) AT GRID POINT 1 0266 C IF NEQ=3, THEN, IN ADDITION, THE INTEGRAL OVER THE SOLUTION 0267 C (USING THE QUADRATURE APPROXIMATION) IS CONSTRAINED TO BE 0268 C RUSER(6). NEQ = 1, 2 OR 3 AND RUSER(J), J=1, (AND 2 AND 6 0269 C IF NEQ=2 OR 3) WOULD HAVE TO HAVE BEEN PREVIOUSLY INPUT. 0270 SUBROUTINE USEREQ (AEQ,CQUAD,MEQ,MG) 0271 DOUBLE PRECISION PRECIS, RANGE 0272 DOUBLE PRECISION AEQ 0273 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0274 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0275 DIMENSION AEQ(MEQ,MG), CQUAD(MG) 0276 COMMON /DBLOCK/ PRECIS, RANGE 0277 COMMON /SBLOCK/ DFMIN, SRMIN, 0278 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 0279 2 SRANGE 0280 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0281 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0282 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0283 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0284 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0285 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0286 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0287 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0288 2 LUSER(30) 0289 C ZERO=0.E0!SP 0290 ZERO=0.D0 0291 C ONE=1.E0!SP 0292 ONE=1.D0 0293 C----------------------------------------------------------------------- 0294 C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 0295 C FOR YOUR APPLICATION. 0296 C----------------------------------------------------------------------- 0297 IF (NEQ.GE.1 .AND. NEQ.LE.3) GO TO 105 0298 5105 FORMAT (/6H NEQ =,I3,28H IS NOT 1, 2 OR 3 IN USEREQ.) 0299 WRITE (NOUT,5105) NEQ 0300 STOP 0301 105 L=MIN0(NEQ,2) 0302 DO 110 J=1,L 0303 DO 120 K=1,NGL 0304 AEQ(J,K)=ZERO 0305 120 CONTINUE 0306 AEQ(J,NGLP1)=RUSER(J) 0307 110 CONTINUE 0308 AEQ(1,NG)=ONE 0309 IF (NEQ .GT. 1) AEQ(2,1)=ONE 0310 IF (NEQ .NE. 3) GO TO 800 0311 DO 130 K=1,NGL 0312 AEQ(3,K)=ZERO 0313 IF (K .LE. NG) AEQ(3,K)=CQUAD(K) 0314 130 CONTINUE 0315 AEQ(3,NGLP1)=RUSER(6) 0316 800 RETURN 0317 END 0318 C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0319 C FUNCTION USEREX. THIS IS A USER-SUPPLIED FUNCTION (ONLY USED 0320 C IF SIMULA=.TRUE.) TO EVALUATE THE NOISE-FREE VALUE OF THE 0321 C SIMULATED DATA AT T(IROW) (THE INDEPENDENT VARIABLE AT THE 0322 C DATA POINT IROW). 0323 C BELOW IS ILLUSTRATED THE CASE OF THE SOLUTION BEING COMPOSED OF 0324 C A SUM OF IUSER(11) DIRAC DELTA FUNCTIONS LOCATED AT 0325 C RUSER(25), RUSER(27), ... , RUSER(23+2*IUSER(11)) 0326 C WITH RESPECTIVE AMPLITUDES 0327 C RUSER(26), RUSER(28), ... , RUSER(24+2*IUSER(11)). 0328 C IN ADDITION, THE NLINF FUNCTIONS IN USERLF ARE ADDED WITH 0329 C RESPECTIVE AMPLITUDES 0330 C RUSER(41), RUSER(42), ... , RUSER(40+NLINF). 0331 C IUSER(11) AND NLINF CANNOT BOTH BE LESS THAN ONE. 0332 C IUSER(11) CANNOT EXCEED 7. 0333 C NLINF CANNOT EXCEED 10. 0334 C WHEN IWT=5, THE NORMALIZATION CONSTANT, RUSER(14), IS COMPUTED 0335 C FOR THE SPECIAL CASE OF PHOTON CORRELATION SPECTROSCOPY SO THAT 0336 C GAMMA = RUSER(26)+RUSER(28)+...+RUSER(24+2*IUSER(11)), WHERE 0337 C GAMMA IS GIVEN IN EQ.(6) IN MAKROMOL. CHEM., VOL. 180, P. 201 0338 C (1979). 0339 C----------------------------------------------------------------------- 0340 C CALLS SUBPROGRAMS - USERLF, USERK, ERRMES 0341 C WHICH IN TURN CALL - USERTR 0342 C----------------------------------------------------------------------- 0343 FUNCTION USEREX (IROW,T,MY,G,MG) 0344 DOUBLE PRECISION PRECIS, RANGE 0345 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0346 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0347 DIMENSION T(MY), IHOLER(6), ADUM(1), G(MG) 0348 COMMON /DBLOCK/ PRECIS, RANGE 0349 COMMON /SBLOCK/ DFMIN, SRMIN, 0350 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 0351 2 SRANGE 0352 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0353 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0354 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0355 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0356 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0357 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0358 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0359 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0360 2 LUSER(30) 0361 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HE, 1HX/ 0362 C----------------------------------------------------------------------- 0363 C THE FOLLOWING STATEMENTS SHOULD BE REPLACED WITH THE ONES 0364 C APPROPRIATE FOR YOUR SIMULATION. 0365 C----------------------------------------------------------------------- 0366 IF (NLINF.GT.10 .OR. IUSER(11).GT.8 .OR. 0367 1 MAX0(NLINF,IUSER(11)).LE.0) CALL ERRMES (1,.TRUE.,IHOLER,NOUT) 0368 USEREX=0. 0369 IF (NLINF .LE. 0) GO TO 150 0370 DO 110 J=1,NLINF 0371 JJ=J 0372 IF (ABS(RUSER(J+40)) .GT. 0.) USEREX=USEREX+RUSER(J+40)* 0373 1 USERLF(IROW,JJ,T,MY) 0374 110 CONTINUE 0375 150 IF (IUSER(11) .LE. 0) GO TO 800 0376 JJ=IUSER(11) 0377 IF (IROW .GT. 1) GO TO 158 0378 RUSER(14)=1. 0379 IF (IWT .NE. 5) GO TO 158 0380 RUSER(14)=0. 0381 AMPSUM=0. 0382 DO 155 J=1,JJ 0383 AMPSUM=AMPSUM+RUSER(2*J+24) 0384 ADUM(1)=0. 0385 G(NG+1)=RUSER(2*J+23) 0386 RUSER(14)=RUSER(14)+RUSER(2*J+24)*USERK(1,ADUM,NG+1,G) 0387 155 CONTINUE 0388 IF (RUSER(14) .LE. 0.) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) 0389 RUSER(14)=AMPSUM/RUSER(14) 0390 158 DO 160 J=1,JJ 0391 G(NG+1)=RUSER(2*J+23) 0392 USEREX=USEREX+RUSER(2*J+24)*USERK(IROW,T,NG+1,G)*RUSER(14) 0393 160 CONTINUE 0394 800 RETURN 0395 END 0396 C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0397 C SUBROUTINE USERGR. THIS IS A USER-SUPPLIED ROUTINE (ONLY CALLED 0398 C WHEN IGRID IS INPUT AS 3) FOR COMPUTING A NONSTANDARD GRID G 0399 C AND THE QUADRATURE WEIGHTS CQUAD. 0400 C IN THE EXAMPLE BELOW, G IS SIMPLY READ IN. CQUAD IS 0401 C THEN COMPUTED FOR THE TRAPEZOIDAL RULE. 0402 C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 0403 C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 0404 C----------------------------------------------------------------------- 0405 C CALLS SUBPROGRAMS - CQTRAP 0406 C----------------------------------------------------------------------- 0407 SUBROUTINE USERGR (G,CQUAD,MG) 0408 DOUBLE PRECISION PRECIS, RANGE 0409 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0410 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0411 DIMENSION G(MG), CQUAD(MG) 0412 COMMON /DBLOCK/ PRECIS, RANGE 0413 COMMON /SBLOCK/ DFMIN, SRMIN, 0414 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 0415 2 SRANGE 0416 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0417 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0418 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0419 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0420 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0421 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0422 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0423 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0424 2 LUSER(30) 0425 C----------------------------------------------------------------------- 0426 C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 0427 C FOR YOUR APPLICATION. 0428 C----------------------------------------------------------------------- 0429 READ (NIN,5100) (G(J),J=1,NG) 0430 5100 FORMAT (5E15.6) 0431 C----------------------------------------------------------------------- 0432 C CQTRAP USES G TO PUT TRAPEZOIDAL RULE WEIGHTS IN CQUAD. 0433 C----------------------------------------------------------------------- 0434 CALL CQTRAP (G,CQUAD,NG) 0435 RETURN 0436 END 0437 C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0846 C FUNCTION USERLF. THIS IS A USER-SUPPLIED ROUTINE (NEEDED IF 0847 C NLINF IS INPUT POSITIVE) TO CALCULATE THE NLINF ADDITIONAL 0848 C KNOWN FUNCTIONS THAT ARE PRESENT IN UNKNOWN AMOUNTS 0849 C IN THE DATA. 0850 C JY = THE SUBSCRIPT OF THE DATA POINT. 0851 C JLINF = INDEX TELLING WHICH OF THE NLINF FUNCTIONS IS TO BE 0852 C EVALUATED. (1 .LE. JLINF .LE. NLINF) 0853 C NYDIM = THE DIMENSION OF THE ARRAY T. IT CAN BE MY OR NY. 0854 C BELOW IS ILLUSTRATED THE CASE WHERE 2 SETS OF DATA CAN BE 0855 C COMBINED AND EACH CAN HAVE A DIFFERENT CONSTANT ADDITIVE 0856 C BACKGROUND (BASELINE). THE FIRST IUSER(2) LOCATIONS IN 0857 C Y, T, ETC. ARE OCCUPIED BY THE FIRST DATA SET. THE SECOND SET 0858 C OCCUPIES LOCATIONS IUSER(2)+1 THRU NY. THEREFORE JLINF=1 OR 2, 0859 C AND USERLF=1. OR 0. DEPENDING ON WHETHER THE (JY)TH DATA POINT 0860 C BELONGS TO DATA SET JLINF OR NOT. 0861 C IF, HOWEVER, IUSER(2)=0 AND NLINF=1, THEN IT IS ASSUMED THAT 0862 C THERE IS ONLY ONE SET OF DATA AND ONLY ONE BASELINE. (THIS 0863 C SAVES YOU THE TROUBLE OF SETTING IUSER(2)=NY EVERY TIME NY 0864 C CHANGES.) 0865 C IF NLINF=0, THEN THERE WILL BE NO BASELINE. 0866 C----------------------------------------------------------------------- 0867 C CALLS SUBPROGRAMS - ERRMES 0868 C----------------------------------------------------------------------- 0869 FUNCTION USERLF (JY,JLINF,T,NYDIM) 0870 DOUBLE PRECISION PRECIS, RANGE 0871 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0872 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0873 DIMENSION T(NYDIM), IHOLER(6) 0874 COMMON /DBLOCK/ PRECIS, RANGE 0875 COMMON /SBLOCK/ DFMIN, SRMIN, 0876 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 0877 2 SRANGE 0878 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0879 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0880 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0881 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0882 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0883 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0884 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0885 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0886 2 LUSER(30) 0887 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HL, 1HF/ 0888 IF (JY.GT.NY .OR. JY.LE.0) CALL ERRMES (1,.TRUE.,IHOLER,NOUT) 0889 C----------------------------------------------------------------------- 0890 C THE FOLLOWING STATEMENTS MUST BE REPLACED BY THE APPROPRIATE 0891 C ONES FOR YOUR ADDITIVE FUNCTIONS. 0892 C----------------------------------------------------------------------- 0893 IF (JLINF.LT.1 .OR. JLINF.GT.2) CALL ERRMES (2,.TRUE.,IHOLER, 0894 1 NOUT) 0895 USERLF=0. 0896 IF ((JY.LE.IUSER(2) .AND. JLINF.EQ.1) .OR. (JY.GT.IUSER(2) .AND. 0897 1 JLINF.EQ.2) .OR. IUSER(2).LE.0) USERLF=1. 0898 RETURN 0899 END 0900 C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0962 C SUBROUTINE USEROU. THIS IS A USER-SUPPLIED ROUTINE (ONLY CALLED 0963 C WHEN DOUSOU=.TRUE.) THAT YOU CAN USE TO PRODUCE YOUR OWN EXTRA 0964 C OUTPUT. 0965 C G CONTAINS THE NG GRID POINTS. 0966 C SOL CONTAINS NG+NLINF VALUES - FIRST THE NG SOLUTION VALUES AND THEN 0967 C THE NLINF LINEAR COEFFICIENTS. 0968 C EXACT CONTAINS THE NG VALUES OF THE SECOND CURVE TO BE PLOTTED THAT 0969 C HAS BEEN COMPUTED IN USERSX (ONLY IF ONLY1=.FALSE.). 0970 C AA IS ( SQUARE ROOT OF THE COVARIANCE MATRIX OF THE SOLUTION)/STDDEV. 0971 C YOU CAN USE AA AND THE LAST 4 ARGUMENTS TO COMPUTE ERROR 0972 C ESTIMATES, AS SHOWN BELOW. 0973 C (CONTIN SETS DOERR=.FALSE. IF IT WAS NOT ABLE TO COMPUTE AA OR 0974 C STDDEV. THE ERROR ESTIMATES THEN CANNOT BE COMPUTED AND ARE SET 0975 C TO ZERO.) 0976 C BELOW IS ILLUSTRATED THE CASE FROM THE ANALYSIS OF CIRCULAR DICHROIC 0977 C SPECTRA, WHERE THE NSPECT-BY-1 VECTOR SOL (THE FRACTIONS OF EACH 0978 C REFERENCE PROTEIN SPECTRUM IN THE SPECTRUM BEING ANALYZED) IS 0979 C RIGHT-MULTIPLIED BY THE NCLASS-BY-NSPECT MATRIX FCAP (THE X-RAY 0980 C VALUES OF THE FRACTIONS OF THE REFERENCE PROTEINS IN EACH CLASS) 0981 C TO YIELD THE NCLASS-BY-1 VECTOR F (THE FRACTION OF THE PROTEIN 0982 C BEING ANALYZED IN EACH CLASS). F AND THE ERROR ESTIMATES ARE 0983 C THEN NORMALIZED BY DIVIDING EACH ELEMENT BY THE SUM OF THE F'S 0984 C (WHICH IS ALSO PRINTED AS THE SCALE FACTOR, SUMF). 0985 C NSPECT = THE NUMBER OF REFERENCE PROTEIN CD SPECTRA. 0986 C NCLASS = THE NUMBER OF CONFORMATIONAL CLASSES. 0987 SUBROUTINE USEROU (CQUAD,G,SOL,EXACT,AA,MG,STDDEV,DOERR,NGLEY) 0988 DOUBLE PRECISION AA 0989 DOUBLE PRECISION PRECIS, RANGE 0990 LOGICAL DOCHOS, DOERR, DOMOM, DOUSIN, DOUSNQ, LAST, 0991 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0992 DIMENSION CQUAD(MG), G(MG), SOL(MG), EXACT(MG), AA(MG,MG) 0993 C----------------------------------------------------------------------- 0994 C IF YOU CHANGE NSPECT OR NCLASS IN THE DATA STATEMENT BELOW, THEN YOU 0995 C MUST READJUST THE DIMENSIONS IN THE FOLLOWING STATEMENT TO 0996 C DIMENSION FCAP(NCLASS,NSPECT), F(NCLASS), ERROR(NCLASS) 0997 C AND, OF COURSE, MODIFY FCAP (THE X-RAY FRACTIONS) AND NSPECT 0998 C AND NCLASS IN THE DATA STATEMENT BELOW. 0999 C----------------------------------------------------------------------- 1000 DIMENSION FCAP(3,16), F(3), ERROR(3) 1001 COMMON /DBLOCK/ PRECIS, RANGE 1002 COMMON /SBLOCK/ DFMIN, SRMIN, 1003 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 1004 2 SRANGE 1005 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 1006 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 1007 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 1008 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 1009 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 1010 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 1011 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 1012 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 1013 2 LUSER(30) 1014 DATA NSPECT/16/, NCLASS/3/, FCAP/ 1015 1 .79,.00,.21 , .41,.16,.43 , .23,.40,.37 , 1016 2 .28,.14,.58 , .45,.24,.31 , .09,.34,.57 , 1017 3 .02,.51,.47 , .39,.00,.61 , .07,.52,.41 , 1018 4 .24,.15,.61 , .51,.24,.25 , .62,.05,.33 , 1019 5 .37,.15,.48 , .28,.33,.39 , 1020 6 .54,.12,.34 , .26,.44,.30 / 1021 SUMF=0. 1022 DO 110 J=1,NCLASS 1023 F(J)=0. 1024 DO 120 K=1,NSPECT 1025 F(J)=F(J)+SOL(K)*FCAP(J,K) 1026 120 CONTINUE 1027 SUMF=SUMF+F(J) 1028 110 CONTINUE 1029 DO 140 J=1,NCLASS 1030 F(J)=F(J)/SUMF 1031 140 CONTINUE 1032 FACTOR=STDDEV/SUMF 1033 C----------------------------------------------------------------------- 1034 C RUSER(14) IS THE SCALE FACTOR BY WHICH THE DATA SPECTRUM IS 1035 C DIVIDED (SEE USERS MANUAL). NORMALLY, FACTOR=STDDEV ABOVE. 1036 IF (ABS(RUSER(14)) .GT. 0.) SUMF=SUMF/RUSER(14) 1037 C----------------------------------------------------------------------- 1038 C ABOVE, THE NCLASS-BY-1 VECTOR F WAS PRODUCED BY MULTIPLYING THE 1039 C NSPECT-BY-1 SOLUTION SOL BY THE NCLASS-BY-NSPECT MATRIX FCAP. 1040 C THIS IS A VERY COMMON TYPE OF TRANSFORMATION, THAT OCCURS IN 1041 C MANY APPPLICATIONS, E.G., WHEN THE SOLUTION IS REPRESENTED BY 1042 C A SUM OF BASIS (E.G., SPLINE) FUNCTIONS. IN THIS LATTER CASE, 1043 C SOL WOULD BE THE EXPANSION COEFFICIENTS AND F WOULD BE THE 1044 C SOLUTION VALUES AT NCLASS GRID POINTS (SEE USERS MANUAL.) 1045 C THE COMPUTATION OF THE NCLASS-BY-1 VECTOR ERROR, THE 1046 C STANDARD ERROR ESTIMATES FOR F, BELOW IS COMPLETELY GENERAL 1047 C FOR ANY FCAP. SO YOU DO NOT HAVE TO CHANGE ANYTHING BELOW IF 1048 C FCAP, NCLASS, OR NSPECT ARE CHANGED. (USUALLY NSPECT=NGL, BUT 1049 C THIS IS NOT NECESSARY.) 1050 C----------------------------------------------------------------------- 1051 DO 210 J=1,NCLASS 1052 ERROR(J)=0. 1053 IF (.NOT.DOERR) GO TO 210 1054 DO 220 K=1,NGLEY 1055 SUM=0. 1056 DO 230 L=1,NSPECT 1057 SUM=SUM+FCAP(J,L)*AA(L,K) 1058 230 CONTINUE 1059 ERROR(J)=ERROR(J)+SUM**2 1060 220 CONTINUE 1061 ERROR(J)=FACTOR*SQRT(ERROR(J)) 1062 210 CONTINUE 1063 C----------------------------------------------------------------------- 1064 C IF YOU CHANGE THE CLASSES, THEN YOU MUST MODIFY THE FOLLOWING 1065 C FORMAT AND WRITE STATEMENTS. 1066 C----------------------------------------------------------------------- 1067 5140 FORMAT (/23X,5HHELIX,3X,10HBETA-SHEET,4X, 1068 1 9HREMAINDER,16X,12HSCALE FACTOR/ 1069 2 9H FRACTION,F19.2,2F13.2,F28.3/ 1070 3 15H STANDARD ERROR,1P3E13.1) 1071 WRITE (NOUT,5140) F,SUMF,ERROR 1072 RETURN 1073 END 1074 C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 1141 C SUBROUTINE USERSI. THIS IS A USER-SUPPLIED ROUTINE (ONLY 1142 C CALLED WHEN SIMULA IS .TRUE.) FOR CALCULATING EXACT(J) 1143 C (THE SIMULATED DATA BEFORE NOISE IS ADDED) AND 1144 C Y(J) (THE SIMULATED NOISY DATA) FOR J=1,NY. 1145 C EXACT(J) MUST BE COMPUTED IN USEREX. 1146 C IUSER(3) = STARTING INTEGER FOR RANDOM NUMBER GENERATOR RANDOM. 1147 C IUSER(3) AND RUSER(3) MUST BE SUPPLIED BY THE USER. 1148 C IUSER(3) MUST BE BETWEEN 1 AND 2147483646. IF IT IS NOT, THEN 1149 C IT IS SET TO THE DEFAULT VALUE OF 30171. 1150 C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 1151 C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 1152 C BELOW IS ILLUSTRATED THE CASE WHERE ZERO-MEAN 1153 C NORMALLY DISTRIBUTED PSEUDORANDOM NOISE IS ADDED TO EXACT(J). 1154 C SD(J), THE STANDARD DEVIATION OF THE NOISE AT POINT J, IS 1155 C DETERMINED BY IWT AND RUSER(3) AS FOLLOWS - 1156 C IWT = 1 CAUSES SD(J)=RUSER(3) FOR ALL J. 1157 C IWT = 2 CAUSES SD(J)=RUSER(3)*SQRT(EXACT(J)), AS APPROPRIATE 1158 C FOR POISSON STATISTICS. IN THE POISSON CASE, 1159 C RUSER(3) IS JUST A SCALE FACTOR FOR THE CASE THAT 1160 C EXACT(J) IS NOT IN NUMBER OF EVENTS, I.E., 1161 C RUSER(3)=SQRT(EXACT(J)/(NO. OF EVENTS IN CHANNEL J)). 1162 C THUS, EXACT(J)/RUSER(3)**2 IS THE POISSON VARIABLE. 1163 C IF EXACT(J) IS ALREADY IN NO. OF EVENTS, THEN YOU 1164 C SHOULD SET RUSER(3)=1. 1165 C IWT = 3 CAUSES SD(J)=RUSER(3)*EXACT(J). 1166 C IWT = 4 CAUSES SD(J)=RUSER(3)/SQRTW(J). 1167 C IWT = 5 IS FOR THE SPECIAL CASE OF A POISSON SECOND ORDER 1168 C CORRELATION FUNCTION IN PHOTON CORRELATION SPECTROSCOPY. 1169 C YOU SHOULD SET RUSER(3)=1/SQRT(B), AND 1170 C USEREX MUST COMPUTE THE SIMULATED (1ST ORDER CORRELATION 1171 C FUNCTION)*GAMMA, WHERE GAMMA AND B ARE GIVEN IN EQ.(6) 1172 C OF MAKROMOL. CHEM., VOL. 180, P. 201 (1979). 1173 C (SEE ALSO THE USERS MANUAL.) 1174 C----------------------------------------------------------------------- 1175 C CALLS SUBPROGRAMS - USEREX, RGAUSS, ERRMES 1176 C WHICH IN TURN CALL - RANDOM, USERK, USERLF, USERTR 1177 C----------------------------------------------------------------------- 1178 SUBROUTINE USERSI (EXACT,G,MG,MY,SQRTW,T,Y) 1179 DOUBLE PRECISION PRECIS, RANGE 1180 DOUBLE PRECISION DUB 1181 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 1182 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 1183 DIMENSION T(MY), EXACT(MY), Y(MY), SQRTW(MY), G(MG) 1184 DIMENSION RN(2), IHOLER(6) 1185 COMMON /DBLOCK/ PRECIS, RANGE 1186 COMMON /SBLOCK/ DFMIN, SRMIN, 1187 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 1188 2 SRANGE 1189 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 1190 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 1191 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 1192 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 1193 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 1194 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 1195 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 1196 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 1197 2 LUSER(30) 1198 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HS, 1HI/ 1199 TWOPI=6.2831853072D0 1200 DUB=DBLE(FLOAT(IUSER(3))) 1201 IF (DUB.LT.1.D0 .OR. DUB.GT.2147483646.D0) DUB=30171.D0 1202 DO 150 J=1,NY 1203 JJ=J 1204 C----------------------------------------------------------------------- 1205 C RGAUSS DELIVERS TWO NORMAL DEVIATES WITH ZERO MEANS AND STANDARD 1206 C DEVIATIONS 1. THEREFORE IT IS ONLY CALLED FOR ODD J. 1207 C----------------------------------------------------------------------- 1208 K=2-MOD(J,2) 1209 IF (K .EQ. 1) CALL RGAUSS (RN(1),RN(2),TWOPI,DUB) 1210 C----------------------------------------------------------------------- 1211 C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 1212 C FOR YOUR APPLICATION. 1213 C----------------------------------------------------------------------- 1214 EXACT(J)=USEREX(JJ,T,MY,G,MG) 1215 C----------------------------------------------------------------------- 1216 C THE NEXT STATEMENT TEMPORARILY SETS THE ERROR TO A NORMAL DEVIATE 1217 C WITH MEAN = ZERO AND STANDARD DEVIATION = RUSER(3). 1218 C----------------------------------------------------------------------- 1219 ERROR=RUSER(3)*RN(K) 1220 IF (IWT .NE. 2) GO TO 160 1221 IF (EXACT(J) .LT. 0.) CALL ERRMES (1,.TRUE.,IHOLER,NOUT) 1222 ERROR=ERROR*SQRT(EXACT(J)) 1223 GO TO 190 1224 160 IF (IWT .EQ. 3) ERROR=ERROR*EXACT(J) 1225 IF (IWT .NE. 4) GO TO 170 1226 IF (SQRTW(J) .LE. 0.) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) 1227 ERROR=ERROR/SQRTW(J) 1228 170 IF (IWT .NE. 5) GO TO 190 1229 C----------------------------------------------------------------------- 1230 C SPECIAL CASE OF A POISSON SECOND ORDER CORRELATION FUNCTION IN 1231 C PHOTON CORRELATION SPECTROSCOPY. 1232 C G2 = NOISE-FREE NORMALIZED 2ND ORDER CORRELATION FUNCTION. 1233 C G2N1 = NOISY 2ND ORDER CORRELATION FUNCTION LESS 1. 1234 C Y(J) = NOISY FIRST ORDER CORRELATION FUNCTION. IT IS EVALUATED 1235 C BELOW WITH THE FUNCTION SIGN() SO THAT NEGATIVE DATA POINTS 1236 C REMAIN NEGATIVE. THIS PREVENTS A BIAS TOWARD POSITIVE 1237 C NOISE COMPONENTS. 1238 C----------------------------------------------------------------------- 1239 G2=1.+EXACT(J)**2 1240 G2N1=G2-1.+ERROR*SQRT(G2) 1241 Y(J)=SIGN(SQRT(ABS(G2N1)),G2N1) 1242 GO TO 150 1243 190 Y(J)=EXACT(J)+ERROR 1244 150 CONTINUE 1245 RETURN 1246 END 1247 C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 1248 C SUBROUTINE USERSX. THIS IS A USER-SUPPLIED ROUTINE (ONLY 1249 C CALLED WHEN YOU HAVE INPUT ONLY1=.FALSE.) 1250 C TO COMPUTE EXACT(J),J=1,NG, WHICH WILL BE PLOTTED WITH 1251 C EACH SOLUTION. USUALLY EXACT IS THE EXACT THEORETICAL 1252 C SOLUTION USED TO SIMULATE YOUR DATA IN USERSI EVALUATED AT 1253 C THE GRID POINTS G(J),J=1,NG. 1254 C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 1255 C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 1256 C BELOW IS ILLUSTRATED THE CASE WHERE EXACT IS 1257 C G(J)**RUSER(8)*EXP(-G(J))/(FACTORIAL OF RUSER(8)), WHERE 1258 C 1. .LE. RUSER(8) .LE. 20. AND THE G(J) ARE NONNEGATIVE. 1259 C----------------------------------------------------------------------- 1260 C CALLS SUBPROGRAMS - GAMLN 1261 C----------------------------------------------------------------------- 1262 SUBROUTINE USERSX (EXACT,G,MG) 1263 DOUBLE PRECISION PRECIS, RANGE 1264 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 1265 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 1266 DIMENSION EXACT(MG), G(MG) 1267 COMMON /DBLOCK/ PRECIS, RANGE 1268 COMMON /SBLOCK/ DFMIN, SRMIN, 1269 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 1270 2 SRANGE 1271 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 1272 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 1273 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 1274 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 1275 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 1276 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 1277 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 1278 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 1279 2 LUSER(30) 1280 C----------------------------------------------------------------------- 1281 C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 1282 C FOR YOUR EVALUATION OF EXACT. 1283 C----------------------------------------------------------------------- 1284 IF (RUSER(8).GE.1. .AND. RUSER(8).LE.20.) GO TO 120 1285 5120 FORMAT (/11H RUSER(8) =,E12.4,27H IS OUT OF RANGE IN USERSX.) 1286 WRITE (NOUT,5120) RUSER(8) 1287 STOP 1288 120 EXMIN=-ALOG(SRANGE) 1289 FACTL=GAMLN(RUSER(8)+1.) 1290 DO 150 J=1,NG 1291 EXACT(J)=0. 1292 IF (G(J)) 160,150,180 1293 160 WRITE (NOUT,5160) 1294 5160 FORMAT (/22H NEGATIVE G IN USEREX.) 1295 STOP 1296 180 EX=RUSER(8)*ALOG(G(J))-G(J)-FACTL 1297 IF (EX .GE. EXMIN) EXACT(J)=EXP(EX) 1298 150 CONTINUE 1299 RETURN 1300 END 1301 C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 1302 C FUNCTION USERTR. THIS IS A USER-SUPPLIED ROUTINE 1303 C FOR COMPUTING THE QUADRATURE-GRID TRANSFORMATION 1304 C (CALL IT H) H(G) WHEN IFUNCT=1, THE INVERSE TRANSFORMATION 1305 C WHEN IFUNCT=2, AND THE DERIVATIVE OF THE TRANSFORMATION 1306 C WHEN IFUNCT=3. WHEN IGRID=2, G (THE QUADRATURE GRID) WILL 1307 C BE IN EQUAL INTERVALS OF H(G) RATHER THAN IN EQUAL 1308 C INTERVALS OF G (AS IT IS WHEN IGRID=1 AND H IS THE 1309 C IDENTITY TRANSFORMATION). 1310 C BELOW IS ILLUSTRATED THE CASE H(G)=ALOG(G). FOR ANOTHER H, 1311 C YOU CAN REPLACE THE STATEMENTS NUMBERED 210, 220, AND 230. 1312 C THESE ARE THE ONLY STATEMENTS THAT CAN BE REPLACED. ALSO 1313 C NOTE THAT ONLY AN H THAT IS MONOTONIC IN THE RANGE OF 1314 C INTEGRATION MAKES SENSE. 1315 C----------------------------------------------------------------------- 1316 C CALLS SUBPROGRAMS - ERRMES 1317 C----------------------------------------------------------------------- 1318 FUNCTION USERTR (X,IFUNCT) 1319 DOUBLE PRECISION PRECIS, RANGE 1320 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 1321 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 1322 COMMON /DBLOCK/ PRECIS, RANGE 1323 COMMON /SBLOCK/ DFMIN, SRMIN, 1324 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 1325 2 SRANGE 1326 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 1327 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 1328 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 1329 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 1330 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 1331 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 1332 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 1333 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 1334 2 LUSER(30) 1335 DIMENSION IHOLER(6) 1336 DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HT, 1HR/ 1337 IF (IFUNCT.LT.1 .OR. IFUNCT.GT.3) CALL ERRMES (1,.TRUE., 1338 1 IHOLER,NOUT) 1339 IF (IGRID .NE. 1) GO TO 200 1340 C----------------------------------------------------------------------- 1341 C COMPUTE TRANSFORMATION, INVERSE, AND DERIVATIVE FOR IDENTITY 1342 C TRANSFORMATION. 1343 C----------------------------------------------------------------------- 1344 USERTR=1. 1345 IF (IFUNCT .NE. 3) USERTR=X 1346 RETURN 1347 200 IF (IGRID .NE. 2) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) 1348 GO TO (210,220,230),IFUNCT 1349 C----------------------------------------------------------------------- 1350 C COMPUTE TRANSFORMATION. 1351 C----------------------------------------------------------------------- 1352 210 USERTR=ALOG(X) 1353 RETURN 1354 C----------------------------------------------------------------------- 1355 C COMPUTE INVERSE TRANSFORMATION. 1356 C----------------------------------------------------------------------- 1357 220 USERTR=EXP(X) 1358 RETURN 1359 C----------------------------------------------------------------------- 1360 C COMPUTE DERIVATIVE OF TRANSFORMATION. 1361 C----------------------------------------------------------------------- 1362 230 USERTR=1./X 1363 RETURN 1364 END 1365 C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 1403 C SUBROUTINE ANALYZ. DOES COMPLETE CONSTRAINED REGULARIZED 1404 C ANALYSIS. 1405 C----------------------------------------------------------------------- 1406 C CALLS SUBPROGRAMS - SETSCA, SEQACC, SETREG, USEREQ, ELIMEQ, SVDRS2, 1407 C ERRMES, DIAREG, DIAGA, SETGA1, SETVAL, LDPETC, RUNRES, 1408 C ANPEAK, SETSGN, SETNNG 1409 C WHICH IN TURN CALL - H12, GETROW, USERK, USERLF, USERRG, LH1405, 1410 C QRBD, G1, G2, DIFF, LDP, CVNEQ, FISHNI, BETAIN, GAMLN, 1411 C NNLS, PLRES, UPDSGN, UPDDON, FFLAT, UPDLLS, 1412 C GETPRU, GETYLY, PGAUSS, MOMENT, MOMOUT, PLPRIN, USEROU, USERTR 1413 C----------------------------------------------------------------------- 1414 SUBROUTINE ANALYZ (ISTAGE, 1415 1 A,AA,AEQ,AINEQ,CQUAD,EXACT,G,IISIGN,IWORK,LBIND,LSDONE,MA, 1416 2 MDONE,MEQ,MG,MINEQ,MREG,MWORK,MY,PIVOT,REG,RHSNEQ,S,SOLBES, 1417 3 SOLUTN,SQRTW,SSCALE,T,VALPCV,VALPHA,VDONE,VK1Y1,WORK, 1418 4 Y,YLYFIT) 1419 DOUBLE PRECISION PRECIS, RANGE 1420 DOUBLE PRECISION A, AA, ABS, AEQ, AINEQ, ALPBES, ALPHA, 1421 1 ALPOLD, DUB, ONE, PIVOT, RALPFL, REG, RHSNEQ, 1422 2 S, SOLBES, SOLUTN, SSCALE, VALPCV, VALPHA, VK1Y1, WORK, ZERO 1423 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 1424 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 1425 LOGICAL LDUM, LBIND, FLAT, PPLTPR, HEADNG, NEWPAG, NOPAGE 1426 DIMENSION A(MA,MG), T(MY), Y(MY), SQRTW(MY), G(MG), CQUAD(MG), 1427 1 REG(MREG,MG), AEQ(MEQ,MG), PIVOT(MEQ), VK1Y1(MG), S(MG,3), 1428 2 AINEQ(MINEQ,MG), VALPHA(MG), VALPCV(MG), RHSNEQ(MINEQ), 1429 3 WORK(MWORK), IWORK(MA), EXACT(MG), SOLUTN(MG), LBIND(MINEQ), 1430 4 IISIGN(MG), SOLBES(MG), LSDONE(MDONE,3,2), VDONE(MDONE), 1431 5 YLYFIT(MY), AA(MG,MG), SSCALE(MG) 1432 DIMENSION IHOLER(6), PREJ(2), LLSIGN(5), RS2MNX(2) 1433 COMMON /DBLOCK/ PRECIS, RANGE 1434 COMMON /SBLOCK/ DFMIN, SRMIN, 1435 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(100), 1436 2 SRANGE 1437 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 1438 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 1439 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 1440 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 1441 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 1442 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 1443 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 1444 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 1445 2 LUSER(30) 1446 DATA IHOLER/1HA, 1HN, 1HA, 1HL, 1HY, 1HZ/, RALPHA/1./ 1447 ABS(ALPHA)=DABS(ALPHA) 1448 SINGLE(ONE)=SNGL(ONE) 1449 C SINGLE(VAR)=VAR!SP 1450 C ZERO=0.E0!SP 1451 ZERO=0.D0 1452 C ONE=1.E0!SP 1453 ONE=1.D0 1454 NGLE=NGL-NEQ 1455 C----------------------------------------------------------------------- 1456 C PUT SCALE FACTORS FOR SOLUTION IN SSCALE AND SCALE INEQUALITY 1457 C CONSTRAINTS. 1458 C----------------------------------------------------------------------- 1459 MG1=MG+1 1460 MG2=MG1+MG 1461 MG3=MG2+MG 1462 CALL SETSCA (WORK,WORK(MG1),WORK(MG2),WORK(MG3), 1463 1 AINEQ,CQUAD,G,ISTAGE,MG,MINEQ,MREG,MY,NGLE,REG,S,SQRTW,SSCALE,T, 1464 2 Y) 1465 C----------------------------------------------------------------------- 1466 C DO SEQUENTIAL HOUSEHOLDER TRANSFORMATIONS TO COMPRESS 1467 C COEFFICIENT MATRIX INTO AN UPPER TRIANGLE. 1468 C----------------------------------------------------------------------- 1469 CALL SEQACC ( 1470 1 A,CQUAD,G,IUNIT,IWT,MA,MG,NG,NGL,NGLP1,NLINF,NY, 1471 2 RANGE,SQRTW,SSCALE,T,Y) 1472 NGLY=MIN0(NGL,NY) 1473 C----------------------------------------------------------------------- 1474 C SET UP REGULARIZOR. 1475 C----------------------------------------------------------------------- 1476 CALL SETREG (MG,MREG,NENDZ,NG,NGL,NGLE,NGLP1,NORDER, 1477 1 NOUT,NREG,REG,SSCALE) 1478 IF (NEQ .LE. 0) GO TO 200 1479 C----------------------------------------------------------------------- 1480 C SET EQUALITY CONSTRAINTS IN USER-SUPPLIED PROGRAM USEREQ AND THEN 1481 C SCALE THEM AND NORMALIZE THEIR ROWS (IN THE L1 METRIC). 1482 C----------------------------------------------------------------------- 1483 CALL USEREQ (AEQ,CQUAD,MEQ,MG) 1484 DO 150 ICOL=1,NGL 1485 DO 160 IROW=1,NEQ 1486 AEQ(IROW,ICOL)=AEQ(IROW,ICOL)*SSCALE(ICOL) 1487 160 CONTINUE 1488 150 CONTINUE 1489 DO 170 IROW=1,NEQ 1490 DUB=ZERO 1491 DO 180 ICOL=1,NGL 1492 DUB=DUB+ABS(AEQ(IROW,ICOL)) 1493 180 CONTINUE 1494 IF (DUB .LE. ZERO) CALL ERRMES (0,.TRUE.,IHOLER,NOUT) 1495 DUB=ONE/DUB 1496 DO 190 ICOL=1,NGLP1 1497 AEQ(IROW,ICOL)=AEQ(IROW,ICOL)*DUB 1498 190 CONTINUE 1499 170 CONTINUE 1500 C----------------------------------------------------------------------- 1501 C ELIMINATE EQUALITY CONSTRAINTS USING A BASIS OF NULL SPACE OF AEQ. 1502 C----------------------------------------------------------------------- 1503 200 CALL ELIMEQ (AEQ,MEQ,MG,PIVOT,NEQ,NGL,A,MA,REG,MREG,NREG,NGLP1, 1504 1 NGLY,VK1Y1,RANGE) 1505 C----------------------------------------------------------------------- 1506 C DO SINGULAR VALUE DECOMPOSITION AND DIAGONALIZATION OF REGULARIZOR. 1507 C----------------------------------------------------------------------- 1508 CALL SVDRS2 (REG(1,NEQ+1),MREG,NREG,NGLE,REG(1,NGLP1),MREG, 1509 1 1,S,IERROR,RANGE) 1510 IF (IERROR .NE. 1) CALL ERRMES (1,.TRUE.,IHOLER,NOUT) 1511 CALL DIAREG (A,AEQ,MA,MEQ,MG,MREG,NEQ,NGL,NGLE,NGLY,NOUT, 1512 1 NUNREG,PIVOT,PRECIS,RANGE,REG,S) 1513 C----------------------------------------------------------------------- 1514 C DO SINGULAR VALUE DECOMPOSITION AND DIAGONALIZATION OF MODIFIED 1515 C COMPRESSED COEFFICIENT MATRIX A. 1516 C----------------------------------------------------------------------- 1517 CALL SVDRS2 (A,MA,NGLY,NGLE,A(1,NGLP1),MA,1,S,IERROR,RANGE) 1518 IF (IERROR .NE. 1) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) 1519 CALL DIAGA (A,MA,MG,MREG,NEQ,NGL,NGLE,NGLP1,NGLY,REG,S) 1520 DO 320 J=1,NGLE 1521 YLYFIT(J)=S(J,1) 1522 320 CONTINUE 1523 5320 FORMAT (/16H SINGULAR VALUES/(1X,1P10E13.3)) 1524 WRITE (NOUT,5320) (YLYFIT(J),J=1,NGLE) 1525 C----------------------------------------------------------------------- 1526 C LEFT MULTIPLY BY INEQUALITY MATRIX. 1527 C----------------------------------------------------------------------- 1528 IF (NINEQ .GT. 0) CALL SETGA1 (NINEQ, 1529 1 A,AINEQ,MA,MG,MINEQ,MREG,NGL,NGLE,REG) 1530 IF (MAX0(NQPROG(1),NQPROG(2)).LE.0 .AND. ALPST(ISTAGE).LE.0.) CALL 1531 1 ERRMES (3,.TRUE.,IHOLER,NOUT) 1532 RSVM2J=SRANGE 1533 IF (DFMIN.LT.0. .OR. DFMIN+FLOAT(NUNREG).GE.FLOAT(NGLE)) GO TO 390 1534 C----------------------------------------------------------------------- 1535 C GUARANTEE AT LEAST NUNREG+DFMIN DEGREES OF FREEDOM (IGNORING 1536 C INEQUALITY CONSTRAINTS) BY USING RSVM2J INSTEAD OF RSVMNX(2,J), 1537 C IF NECESSARY. 1538 C----------------------------------------------------------------------- 1539 ALP2=S(1,1)**2 1540 DO 350 J=1,500 1541 IF (ALP2 .LE. 0.) GO TO 390 1542 DF=0. 1543 DO 360 K=1,NGLE 1544 DUM=S(K,1)**2 1545 DF=DF+DUM/(DUM+ALP2) 1546 360 CONTINUE 1547 IF (DF .GT. DFMIN+FLOAT(NUNREG)) GO TO 370 1548 ALP2=ALP2*.1 1549 350 CONTINUE 1550 370 RSVM2J=SQRT(10.*ALP2)/S(1,1) 1551 390 BTEST=SRANGE 1552 ALPBES=ZERO 1553 VARZ=SRANGE 1554 RS2MNX(1)=-1. 1555 RS2MNX(2)=-1. 1556 C----------------------------------------------------------------------- 1557 C START OF NQPROG(1) REGULARIZED SOLUTIONS 1558 C----------------------------------------------------------------------- 1559 LDUM=.TRUE. 1560 PPLTPR=MAX0(IPRINT(ISTAGE),IUSROU(ISTAGE)) .GE. 3 1561 NOPAGE=.NOT.PPLTPR 1562 LHEDNG=0 1563 IF (NOPAGE) LHEDNG=4 1564 K=NQPROG(1) 1565 IF (K .LE. 1) GO TO 410 1566 RTOT=AMIN1(RSVM2J,RSVMNX(2,1))/(RSVMNX(1,1)*PRECIS) 1567 IF (RTOT .LE. 1.) RTOT=RSVMNX(2,1)/(RSVMNX(1,1)*PRECIS) 1568 IF (RTOT .LE. 1.) CALL ERRMES (4,.TRUE.,IHOLER,NOUT) 1569 RALPHA=RTOT**(1./FLOAT(K-1)) 1570 410 ALPHA=RSVMNX(1,1)*PRECIS*S(1,1) 1571 DO 420 J=1,K 1572 CALL SETVAL (ALPHA,LDUM,NINEQ, 1573 1 A,AINEQ,MA,MG,MINEQ,MREG,NGL,NGLE,REG,RHSNEQ,S,VALPCV,VALPHA, 1574 2 VK1Y1) 1575 NEWPAG=MAX0(IPRINT(ISTAGE),IUSROU(ISTAGE), 1576 1 IPLRES(ISTAGE)+1).GE.4 .OR. LDUM 1577 HEADNG=PPLTPR .OR. LDUM .OR. 1578 1 MAX0(IPLRES(ISTAGE),IPLFIT(ISTAGE)).GE.3 1579 CALL LDPETC (3,.TRUE.,NINEQ,.TRUE.,ICRIT(ISTAGE),DOMOM,PPLTPR, 1580 1 .TRUE.,ALPHA,HEADNG,NEWPAG,ALPBES,VAR, 1581 2 A,AA,AINEQ,BTEST,CQUAD,DEGFRE,DEGFRZ,EXACT,G,IERROR, 1582 3 ISTAGE,IWORK,LBIND,MA,MG,MINEQ,MREG,MWORK,MY, 1583 4 NGLE,NGLY,PREJ,REG,RHSNEQ,RS2MNX,S,SOLBES, 1584 5 SOLUTN,SQRTW,SSCALE,T,VALPCV,VALPHA,VARREG,VARZ,WORK,Y,YLYFIT) 1585 IF (IERROR.EQ.1) CALL RUNRES (3,SOLUTN,.FALSE., 1586 A SINGLE(ALPHA/S(1,1)),NOPAGE, 1587 1 CQUAD,G,IPLFIT,IPLRES,ISTAGE,ITITLE,IUNIT,IWT,LINEPG-LHEDNG, 1588 2 MWORK,NG,NGL,NLINF,NOUT,NY,SQRTW,SRANGE,SSCALE,T,WORK,Y,YLYFIT) 1589 ALPHA=ALPHA*RALPHA 1590 LDUM=.FALSE. 1591 420 CONTINUE 1592 C----------------------------------------------------------------------- 1593 C START OF NQPROG(2) REGULARIZED SOLUTIONS. 1594 C----------------------------------------------------------------------- 1595 450 IF (ALPST(ISTAGE) .LE. 0.) GO TO 455 1596 K=1 1597 ALPHA=ALPST(ISTAGE) 1598 BTEST=SRANGE 1599 GO TO 475 1600 455 IF (NQPROG(2) .LE. 0) GO TO 490 1601 K=NQPROG(2) 1602 IF (AMIN1(RSVMNX(1,2),RSVMNX(2,2)) .LE. 0.) GO TO 465 1603 ALPHA=RSVMNX(1,2)*PRECIS*S(1,1) 1604 IF (K .LE. 1) GO TO 475 1605 RTOT=AMIN1(RSVM2J,RSVMNX(2,2))/(RSVMNX(1,2)*PRECIS) 1606 IF (RTOT .LE. 1.) RTOT=RSVMNX(2,2)/(RSVMNX(1,2)*PRECIS) 1607 IF (RTOT .LE. 1.) CALL ERRMES (4,.TRUE.,IHOLER,NOUT) 1608 RALPHA=RTOT**(1./FLOAT(K-1)) 1609 GO TO 475 1610 465 NABUT=2 1611 IF (RS2MNX(1) .GT. 0.) NABUT=NABUT-1 1612 IF (RS2MNX(2) .GT. 0.) NABUT=NABUT-1 1613 L=K-NABUT+1 1614 IF (L .LE. 0) GO TO 470 1615 RTOT=ABS(AMIN1(RSVM2J,RS2MNX(2))/(RS2MNX(1)*PRECIS)) 1616 IF (RTOT .LE. 1.) RTOT=ABS(RS2MNX(2)/(RS2MNX(1)*PRECIS)) 1617 RALPHA=RTOT**(1./FLOAT(L)) 1618 IF (RS2MNX(1) .GT. 0.) RS2MNX(1)=RS2MNX(1)*RALPHA 1619 470 ALPHA=ABS(RS2MNX(1)*PRECIS)*S(1,1) 1620 475 DO 480 J=1,K 1621 CALL SETVAL (ALPHA,LDUM,NINEQ, 1622 1 A,AINEQ,MA,MG,MINEQ,MREG,NGL,NGLE,REG,RHSNEQ,S,VALPCV,VALPHA, 1623 2 VK1Y1) 1624 NEWPAG=MAX0(IPRINT(ISTAGE),IUSROU(ISTAGE), 1625 1 IPLRES(ISTAGE)+1).GE.4 .OR. LDUM 1626 HEADNG=PPLTPR .OR. LDUM .OR. 1627 1 MAX0(IPLRES(ISTAGE),IPLFIT(ISTAGE)).GE.3 1628 CALL LDPETC (3,.TRUE.,NINEQ,.TRUE.,ICRIT(ISTAGE),DOMOM,PPLTPR, 1629 1 .TRUE.,ALPHA,HEADNG,NEWPAG,ALPBES,VAR, 1630 2 A,AA,AINEQ,BTEST,CQUAD,DEGFRE,DEGFRZ,EXACT,G,IERROR, 1631 3 ISTAGE,IWORK,LBIND,MA,MG,MINEQ,MREG,MWORK,MY, 1632 4 NGLE,NGLY,PREJ,REG,RHSNEQ,RS2MNX,S,SOLBES, 1633 5 SOLUTN,SQRTW,SSCALE,T,VALPCV,VALPHA,VARREG,VARZ,WORK,Y,YLYFIT) 1634 IF (IERROR.EQ.1) CALL RUNRES (3,SOLUTN,.FALSE., 1635 A SINGLE(ALPHA/S(1,1)),NOPAGE, 1636 1 CQUAD,G,IPLFIT,IPLRES,ISTAGE,ITITLE,IUNIT,IWT,LINEPG-LHEDNG, 1637 2 MWORK,NG,NGL,NLINF,NOUT,NY,SQRTW,SRANGE,SSCALE,T,WORK,Y,YLYFIT) 1638 ALPHA=ALPHA*RALPHA 1639 LDUM=.FALSE. 1640 480 CONTINUE 1641 490 IF (BTEST .GE. SRANGE) CALL ERRMES (5,.TRUE.,IHOLER,NOUT) 1642 IF (NNSGN(ISTAGE) .LE. 0) GO TO 700 1643 C----------------------------------------------------------------------- 1644 C START PEAK-CONSTRAINED SOLUTION BY SETTING UP NEW INEQUALITY 1645 C MATRIX IN AINEQ AND INITIALIZING SO THAT UNSCALED SOLUTION IS 1646 C MONOTONICALLY DECREASING. 1647 C----------------------------------------------------------------------- 1648 NGM1=NG-1 1649 NNINEQ=NINEQ-1 1650 IF (NONNEG) GO TO 510 1651 NNINEQ=NINEQ+NGM1 1652 IF (NNINEQ .LE. MINEQ) GO TO 510 1653 CALL ERRMES (6,.FALSE.,IHOLER,NOUT) 1654 GO TO 790 1655 510 IROW=NNINEQ-NGM1 1656 DO 520 J=1,NGM1 1657 IROW=IROW+1 1658 DO 525 ICOL=1,NGLP1 1659 AINEQ(IROW,ICOL)=ZERO 1660 525 CONTINUE 1661 AINEQ(IROW,J)=ONE*SSCALE(J) 1662 AINEQ(IROW,J+1)=-ONE*SSCALE(J+1) 1663 IISIGN(J)=1 1664 520 CONTINUE 1665 CALL SETGA1 (NNINEQ, 1666 1 A,AINEQ,MA,MG,MINEQ,MREG,NGL,NGLE,REG) 1667 NNSGNI=MIN0(NNSGN(ISTAGE),4) 1668 ALPOLD=ZERO 1669 IF (ISTAGE .EQ. 1) BTEST=SRANGE 1670 C----------------------------------------------------------------------- 1671 C START LOOP TO DO NNSGN(ISTAGE) DIFFERENT PEAK-CONSTRAINED 1672 C ANALYSES. 1673 C----------------------------------------------------------------------- 1674 DO 600 INSGN=1,NNSGNI 1675 NSGNI=NSGN(INSGN) 1676 LNINEQ=NNINEQ 1677 IF (NONNEG) LNINEQ=LNINEQ+(NSGNI-(1+LSIGN(1,INSGN))/2)/2+1 1678 IF (LNINEQ .LE. MINEQ) GO TO 610