C COPYRIGHT (C) 1983 GLENN EVERHART C PERMISSION IS GIVEN TO ANYONE TO USE, DISTRIBUTE, OR COPY THIS C PROGRAM FREELY BUT NOT TO SELL IT COMMERICALLY. C USER FUNCTION ROUTINE C GENERATES PARSING AND EXECUTION OF ROUTINE CALLS OF FORM C *U FNAME (ARGUMENTS) C WHERE LINE (80 BYTES) CONTAINS COMMAND LINE AND ALL C ARGUMENTS MAY BE PARSED. C CALLED FROM CMND C C VAX VERSION: MOST MATRIX ROUTINES AVAILABLE C BUT ASSUMES SUBSTANTIAL SPACE AVAILABLE. C c available parsing aid: c call varscn(line,ibgn,lend,lstchr,id1,id2,ivalid) c where line(ibgn... lend) is scanned. If variable found c ivalid=1 else ivalid=0. id1,id2 are dims in xvbls for c variable found if any. lstchr is last char found+1... C OTHER USEFUL ROUTINES IN THE SHEET: C GN(LAST,LEND,NUMBER,LINE) C LOOKS FROM LINE(LAST) THRU LINE(LEND) FOR A NUMBER AND C RETURNS ANY NUMBER IN "NUMBER" ARG. ASSUMES "LINE" IS A C BYTE ARRAY. (NO INDICATION OF WHERE THE NUMBER WAS FOUND C HOWEVER). THROWS OUT LEADING SPACES, TERMINATES ON A NON C NUMERIC. C INDEX(LINE,CHAR) C EXPECTS LINE TO BE NULL TERMINATED AND RETURNS EITHER C THE SUBSCRIPT (COUNTING FROM 1) OF CHAR IN LINE OR THE C MAX SUBSCRIPT IN LINE (I.E., WHERE IT HIT THE NULL TERMINATOR). C NOTE THIS DIFFERS FROM THE "STANDARD" VERSION OF INDEX WHICH C RETURNS 0 FOR "NOT FOUND" -- THIS VERSION RETURNS MAX LENGTH C FOR "NOT FOUND". STOPS AT 512 BYTES HOWEVER... C PARSING IS UP TO USER. NOTE VARSCN MAY BE CALLED TO PARSE C VARIABLE NAMES. SUPPLIED VERSION CALLS IDATE WHICH RETURNS C SYSTEM DATE IN RSX OR VMS AS INTEGER DAY, MONTH, AND YEAR. C THIS RETURNS HERE IN AC T, U, AND V SUBROUTINE USRFCT(LINE,RETCD) INCLUDE 'VKLUGPRM.FTN' BYTE LINE(80) INTEGER RETCD LOGICAL*1 AVBLS(100,27),WRK(128),VBLS(8,RRW,RCL) INTEGER*2 TYPE(RRW,RCL),VLEN(9) REAL*8 XAC,XVBLS(RRW,RCL) REAL*8 TAC,UAC,VAC,TMP,WAC,YAC INTEGER*4 JVBLS(2,RRW,RCL) EQUIVALENCE(XAC,AVBLS(1,27)) EQUIVALENCE(TAC,AVBLS(1,20)) EQUIVALENCE(UAC,AVBLS(1,21)) EQUIVALENCE(VAC,AVBLS(1,22)) EQUIVALENCE(WAC,AVBLS(1,23)) EQUIVALENCE(YAC,AVBLS(1,25)) EQUIVALENCE(VBLS(1,1,1),JVBLS(1,1,1)) EQUIVALENCE(VBLS(1,1,1),XVBLS(1,1)) COMMON/V/TYPE,AVBLS,VBLS,VLEN INTEGER*2 XTNCNT,XTCFG,IPSET LOGICAL*1 XTNCMD(80) INTEGER*2 FORMFG,RCFGX,PZAP,RCONE INTEGER*2 IDOL1,IDOL2,IDOL3,IDOL4,IDOL5,IDOL6 INTEGER*2 RRWACT,RCLACT COMMON/RCLACT/RRWACT,RCLACT COMMON/DOLLR/IDOL1,IDOL2,IDOL3,IDOL4,IDOL5,IDOL6 COMMON/FFGG/FORMFG,RCFGX,PZAP,RCONE COMMON/XCMD/XTNCNT,XTNCMD,XTCFG,IPSET C LOOP CONTROL FOR VARY FUNCTION. SET ZERO IN SPREDSHT AND C MUST BE SET POSITIVE HERE IF WE NEED ITERATIONS. C (IMPLEMENT FOR VAX ONLY) INTEGER KALKIT COMMON/VARYIT/KALKIT C ARGUMENTS COME IN IN ARGUMENTS IN LINE C RESULTS GO INTO PERCENT (XAC) AND WHEREVER ELSE DESIRED... INTEGER*2 PROW,PCOL,DROW,DCOL,DRWV,DCLV COMMON/DCTL/PROW,PCOL,DROW,DCOL,DRWV,DCLV DIMENSION NRDSP(DRW,DCL),NCDSP(DRW,DCL) COMMON/D2R/NRDSP,NCDSP LOGICAL*1 FNAMS(6,17) C FNAMS IS NAME OF FUNCTION CALLED. DATA FNAMS /'I','D','A','T','E',0, 1 'M','T','X','E','Q',0, 2 'M','O','V','E','V',0, 3 'M','D','E','T',0,0, 4 'M','P','R','O','D',0, 5 'M','A','D','D','V',0, 6 'M','S','U','B','V',0, 7 'M','M','P','Y','T',0, 8 'M','M','P','Y','C',0, 9 'V','A','R','Y',0,0, 1 'X','Q','T','C','M',0, 2 'S','T','R','V','L',0, 3 'H','E','R','E',0,0, 4 'Y','R','M','O','D',0, 5 'J','D','A','T','E',0, 6 'J','T','O','C','H',0, 7 'D','A','T','E',0,0 9 / C NULL TERMINATE ANY NAMES (ALLOWS 5 CHARACTERS) C START LOOKING PAST THE *U C GET FUNCTION NAME AND GO TO PROCESS EACH FUNCTION SEPARATELY C GET NONBLANK CHAR FOR FUNCTION NAME START K=3 30 IF(LINE(K).NE.' ')GOTO 40 K=K+1 IF(K.LT.60)GOTO 30 40 CONTINUE C UNCOMMENT THE DO 100 STMT IF DIM 2 OF FNAMS > 1 400 N=1 C **** BE SURE THE 2ND BOUND ON N IS THE SAME AS THE DIMENSION OF C **** FNAMS ************************** DO 100 N=1,17 KF=N DO 110 NN=1,6 IF(LINE(K+NN-1).NE.FNAMS(NN,N).AND.FNAMS(NN,N).GT.0) 1 GOTO 100 110 CONTINUE GOTO 200 100 CONTINUE C UNRECOGNIZED FUNCTION... IGNORE 300 RETCD=3 RETURN 200 CONTINUE C NOW HAVE FOUND FUNCTION IDENTIFIED BY KF. CALL IT AND ALLOW TO WORK GOTO (1100,1200,1300,1400,1500,1600,1700,1800, 1 1900,2000,2100,2200,2300,2400,2500,2600,2700),KF GOTO 300 1100 CONTINUE C IDATE FUNCTION C RETURNS MONTH, DAY, YEAR IN AC'S T,U,V CALL IDATE(IMO,IDA,IYR) TAC=IMO UAC=IDA VAC=IYR C RETURN A FLOATING VALUE OF DATE FORM AS YYMMDD SO IT CAN BE C USED FOR SORTING AND SIMILAR APPLICATIONS. COULD BE USED ALSO C FOR INTERVALS IF A JULIAN DATE WERE RETURNED, BUT THIS WILL DO C FOR COMPARISONS AND ORDERING. C XAC=VAC*10000.+TAC*100.+UAC XAC=JULMDY(IYR,IMO,IDA) C RETURN JULIAN DATE FOR DATE ARITHMETIC IF NEEDED. C BASED ON 1/1/80. SINCE THIS IS DONE ON 4/24/1984, THAT SHOULD C BE FAR BACK ENOUGH. RETURN 1200 CONTINUE C MATRIX EQUATION. NOTE WE MUST NOW START SCAN FOR ARGUMENTS... C K+5 IS START OF ARG LIST. START AT K+6 TO ALLOW ( TO BE THERE... C FORMAT DESIRED: C *U MTXEQ(A1:A2,X1:X2,B1:B2) GENERATING SOLUTION MATRIX X1:X2 C FROM MATRICES A,B AND SOLVING EQUATION AX=B WHERE A IS AN N BY C N SQUARE MATRIX, AND X AND B ARE N BY M MATRICES. RETCD=1 C COLLECT ARGUMENTS. NOTE THAT VARSCN AND GN TRASH POINTERS PASSED C TO THEM IN IBGN, LEND, SO MAKE UP EVERY TIME. USE VARSCN TO C COLLECT POINTERS TO THE SHEET ARRAY FIRST OFF COMMAND LINE, C THEN PROCESS IN OUR MAGICAL MYSTICAL ROUTINE... IBGN=K+6 LEND=IBGN+20 C GET LOCATIONS OF MATRICES A, X, AND B (FOR AX=B EQN) C A MUST BER N BY N, SQUARE. X,B ARE N BY M. CALL PMTX2(RETCD,3,LINE,IBGN,ID1A,ID2A,ID1B,ID2B, 1 IDXA,IDXB,IDYA,IDYB,IDBA,IDBB,IDCA,IDCB) N=IABS(ID1B-ID1A)+1 C CHECK THAT MATRIX A IS SQUARE IF(N.NE.(IABS(ID2B-ID2A)+1))GOTO 300 C CHECK THAT MATRIX X AND B HAVE THE SAME DIMENSIONS IF((IDYA-IDXA).NE.(IDCA-IDBA))GOTO 300 IF((IDYB-IDXB).NE.(IDCB-IDBB))GOTO 300 M=IABS(IDYA-IDXA)+1 C CHECK THAT THE X AND B MATRIX DIMENSIONS ARE N BY M C WHERE THE N IS THE SAME AS FOR THE A MATRIX NN=IABS(IDYB-IDXB)+1 IF(NN.NE.N)GOTO 300 C NOW HAVE DIMENSIONS FOR ALL THIS STUFF... C SINCE MTXEQU TRASHES ITS' B MATRIX, COPY IT INTO X MATRIX C AND THEN CALL... DO 1210 NN=IDBA,IDCA DO 1210 MM=IDBB,IDCB XVBLS(NN-IDBA+IDXA,MM-IDBB+IDXB)=XVBLS(NN,MM) 1210 CONTINUE C NOW ALL THE ARGUMENTS ARE SET UP... GO DO THE WORK. C CALL UTILITY ROUTINE, THEN DONE... CALL MTXEQU(XVBLS(ID1A,ID2A),XVBLS(IDXA,IDXB),N,M,XAC) RETURN 1300 CONTINUE C MOVEV MTX1 MTX2 MOVE MTX1 VALUES TO MTX2 RETCD=1 IBGN=K+6 CALL PMTX2(RETCD,2,LINE,IBGN,IR1T,IC1T,IR1B,IC1B,IR2T,IC2T, 1 IR2B,IC2B) C CHECK FOR SAME SIZE MATRICES IF((IC1T-IC1B).NE.(IC2T-IC2B))GOTO 300 IF((IR1T-IR1B).NE.(IR2T-IR2B))GOTO 300 C DO THE COPY HERE (EASIER THAN CALLING SOMETHING...) DO 1301 NN=IR1T,IR1B DO 1301 MM=IC1T,IC1B XVBLS(NN-IR1T+IR2T,MM-IC1T+IC2T)=XVBLS(NN,MM) 1301 CONTINUE RETURN 1400 CONTINUE C MDET - DETERMINANT OF SQUARE MATRIX C 1 ARGUMENT, VIZ., MATRIX COORDS RETCD=1 C ACCOUNT FOR "MDET" BEING 4 CHARS NOT 5 IBGN=K+5 CALL PMTX2(RETCD,1,LINE,IBGN,IR1T,IC1T,IR1B,IC1B) C CALL A DETERMINANT ROUTINE TO DO THE WORK C NOTE IT CHECKS FOR SQUARE MATRIX INTERNALLY AND RETURNS 0 IF NOT C SQUARE... CALL MDET(XVBLS,IR1T,IC1T,IR1B,IC1B,XAC) RETURN 1500 CONTINUE C MPROD A,B,C C=A*B MATRIX WISE IBGN=K+6 RETCD=1 IMXX=3 CALL PMTX2(RETCD,IMXX,LINE,IBGN,ID1A,ID2A,ID1B,ID2B, 1 IDXA,IDXB,IDYA,IDYB,IDBA,IDBB,IDCA,IDCB) C A=N BY M C B=M BY L C C=N BY L M=1+ID1B-ID1A N=1+ID2B-ID2A IF(N.NE.(1+IDYB-IDXB))GOTO 300 L=1+IDYA-IDXA IF(M.NE.(1+IDCB-IDBB))GOTO 300 IF(L.NE.(1+IDCA-IDBA))GOTO 300 C DIMENSIONS LOOK OK NOW SO DO THE WORK C USE SLIGHTLY MODIFIED GMPRD CALL GMPRD(XVBLS(ID1A,ID2A),XVBLS(IDXA,IDXB), 1 XVBLS(IDBA,IDBB),N,M,L) RETURN 1600 CONTINUE C MADDV A,B,C C=A+B IMXX=3 IBGN=K+6 RETCD=1 CALL PMTX2(RETCD,IMXX,LINE,IBGN,ID1A,ID2A,ID1B,ID2B, 1 IDXA,IDXB,IDYA,IDYB,IDBA,IDBB,IDCA,IDCB) N=1+ID1B-ID1A M=1+ID2B-ID2A IF(N.NE.(1+IDYA-IDXA))GOTO 300 IF(N.NE.(1+IDCA-IDBA))GOTO 300 IF(M.NE.(1+IDYB-IDXB))GOTO 300 IF(M.NE.(1+IDCB-IDBB))GOTO 300 C USE MODIFIED GMADD CALL GMADD(XVBLS(ID1A,ID2A),XVBLS(IDXA,IDXB), 1 XVBLS(IDBA,IDBB),M,N) RETURN 1700 CONTINUE C MSUBV A,B,C C=A-B IMXX=3 IBGN=K+6 RETCD=1 CALL PMTX2(RETCD,IMXX,LINE,IBGN,ID1A,ID2A,ID1B,ID2B, 1 IDXA,IDXB,IDYA,IDYB,IDBA,IDBB,IDCA,IDCB) N=1+ID1B-ID1A M=1+ID2B-ID2A IF(N.NE.(1+IDYA-IDXA))GOTO 300 IF(N.NE.(1+IDCA-IDBA))GOTO 300 IF(M.NE.(1+IDYB-IDXB))GOTO 300 IF(M.NE.(1+IDCB-IDBB))GOTO 300 CALL GMSUB(XVBLS(ID1A,ID2A),XVBLS(IDXA,IDXB), 1 XVBLS(IDBA,IDBB),M,N) RETURN 1800 CONTINUE C MMPYT A,B,C C=AT*B C GET 3 MATRICES IMXX=3 IBGN=K+6 RETCD=1 CALL PMTX2(RETCD,IMXX,LINE,IBGN,ID1A,ID2A,ID1B,ID2B, 1 IDXA,IDXB,IDYA,IDYB,IDBA,IDBB,IDCA,IDCB) C TRANSPOSE DIMENSIONS OF A... M=1+ID1B-ID1A N=1+ID2B-ID2A IF(N.NE.(1+IDYB-IDXB))GOTO 300 L=1+IDYA-IDXA IF(M.NE.(1+IDCB-IDBB))GOTO 300 IF(L.NE.(1+IDCA-IDBA))GOTO 300 CALL GTPRD(XVBLS(ID1A,ID2A),XVBLS(IDXA,IDXB), 1 XVBLS(IDBA,IDBB),N,M,L) RETURN 1900 CONTINUE C MMPYC A,B,K B=A*K (K=CONSTANT) C FOR MPY BY CONSTANT WE GET MATRICES IN ORDER A,C, THEN AC WITH CONST C IN IT LAST... IBGN=K+6 RETCD=1 IMXX=2 CALL PMTX2(RETCD,IMXX,LINE,IBGN,ID1A,ID2A,ID1B,ID2B, 1 IDXA,IDXB,IDYA,IDYB,IDBA,IDBB,IDCA,IDCB) IF(LINE(IBGN-1).NE.',')GOTO 300 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,IDCA,IDCB,IVALID) IF(IVALID.EQ.0)GOTO 300 C NOW HAVE EVERYTHING OF ARGS... CHECK DIMENSIONS OF MATRICES.... N=1+ID1B-ID1A M=1+ID2B-ID2A IF(N.NE.(1+IDYA-IDXA))GOTO 300 IF(N.NE.(1+IDCA-IDBA))GOTO 300 DO 1901 NN=ID1A,ID1B DO 1901 MM=ID2A,ID2B XVBLS(NN-ID1A+IDXA,MM-ID2A+IDXB)=XVBLS(NN,MM) 1 *XVBLS(IDCA,IDCB) 1901 CONTINUE RETURN C *U VARY X,A,W,I,P;Q;R;S;T C REPEATEDLY COMPUTE SHEET FOR I ITERATIONS (DEFAULTS TO 1 C IF NONE GIVEN) AND VARY AC P,Q,R,S, T (POSITIONAL...WHATEVER C IS NAMED) UNTIL CONDITION THAT AC X (WHATEVER IS NAMED THERE) C IS MADE EQUAL TO AC A AS CLOSELY AS POSSIBLE. DOES MULTI-DIMENSIONAL C STEPPING SEARCH SAVING AC'S AND MODIFYING. ACTUALLY WILL HANDLE ANY C CELL. UP TO 8 DIMENSIONS PERMITTED (ARBITRARY LIMIT). C NOTE THAT RECALCULATE SPECIAL VARY FLAG WILL BE SET HERE IF C VARYING MORE THAN ONCE... C WILL VARY ONE OF THE AC'S IN THE LIST P,Q,R,S,T... BY INITIAL C FRACTION W (AN ARBITRARY "STEP SIZE" FRACTION) AND COMPUTE THE C GRADIENT OF (X-A) WRT THAT AC, THEN WILL REPLACE ALL AC'S AND C VARY THAT AC BY W * THE GRADIENT, MEANING THAT AS THE GRADIENT C DECREASES, THE VARIANCE DOES ALSO. LAST GRADIENTS ARE SAVED AND C USED AS INITIAL VARIANCES, SO THAT THE W FRACTION IS AN INITIAL C GUESS. HOWEVER IT ALSO IS A LIMIT SO NO STEP VARIES AN AC BY C MORE FRACTIONALLY THAN W. C ONCE THIS IS DONE ANOTHER ONE OF THE P,Q,R,S,T,... LIST IS C CHOSEN CIRCULARLY AND THE PROCESS REPEATS. THIS MAY CONTINUE C INDEFINITELY TO LOOK FOR CONVERGENCE. C NOTE THAT X AND A MAY BE ANY CELL AND NEED NOT BE ACCUMULATORS. C HOWEVER ALL OTHER CELLS TO VARY MUST BE AC'S AND MUST BE THE C INDEPENDENT VARIABLES. CALCULATIONS ELSEWHERE ON THE SHEET C (PERHAPS LATER IN THE SAME CELL...)MUST ESTABLISH DEPENDENT C VARIABLES OR BOUNDARY OR NORMALIZATION CONDITIONS. 2000 CONTINUE RETCD=1 C SPLIT OFF THESE FUNCTIONS INTO A COMMON SUBROUTINE CALL VVARY(LINE,RETCD,K) RETURN 2100 CONTINUE C EXECUTE COMMAND. FILL IN COMMAND FROM GIVEN FUNCTION AND C CALL XQTCMD TO DO IT. SETS UP NECESSARY VARIABLES FIRST. C ASSUME THE COMMAND LINE MUST BE ALONE ON LINE AFTER THIS CALL... KK=1 KKK=K+6 DO 2101 NN=KKK,80 XTNCMD(KK)=LINE(NN) IF(XTNCMD(KK).LE.0)GOTO 2102 KK=KK+1 2101 CONTINUE 2102 CONTINUE XTNCMD(KK+1)=0 XTNCMD(KK+2)=0 XTNCNT=KK XTCFG=1 IPSET=1 CALL XQTCMD(ICODE) RETURN 2200 CONTINUE C RETURN PACKED FORMULA STRING TO EXTRACT UP TO 8 CHARS OF C FORMULA. C START AT K+6 XAC=0. IBGN=K+6 IEND=IBGN+20 CALL VARSCN(LINE,IBGN,IEND,LSTC,I1,I2,IVLD) IF(IVLD.LE.0)RETURN C GET START, LENGTH NOW IN FORMULA... IBGN=LSTC+1 IEND=IBGN+20 CALL GN(IBGN,IEND,ISTART,LINE) IBGN=INDEX(LINE,';') C LOOK FOR ';' CHAR AS START OF 2ND NUMBER IF(IBGN.GT.50.OR.ISTART.LE.0.OR.ISTART.GT.80)RETURN C BUMP IBGN PAST THE ; CHAR IBGN=IBGN+1 IEND=80 CALL GN(IBGN,IEND,ILN,LINE) ILN=MIN0(ILN,8) IF(ILN.LE.0)RETURN C READ IN FORMULA INTO WRK ARRAY IRX=(I2-1)*RRW+I1 CALL WRKFIL(IRX,WRK,0) KZ=0 DO 991 NN=1,ILN K=WRK(ISTART+NN-1) K=K.AND.127 IF(K.EQ.0)KZ=1 IF(KZ.EQ.1)K=0 C STOP THE ENCODE ON SEEING ANY NULLS TMP=K XAC=XAC*128.D0+TMP 991 CONTINUE C XAC RETURNS WITH ENCODED VALUE. RETURN 2300 CONTINUE C RETURN PRESENT LOCATION IN THE MATRIX. TAC=PROW UAC=PCOL XAC=(PCOL-1)*RRW+PROW C RETURN ENCODED VM FLAG, RM FLAG, AND RF FLAG C IN V ACCUMULATOR. 4 BIT IS VM FLAG IF 1 C 2 BIT IS RM FLAG IF 1. 1 BIT IS FORCE FLAG IF ONE (SELDOM IF EVER) VAC=4*FORMFG+2*RCFGX+RCONE C VAC=(DROW-1)*DRW+DCOL C RESULT IN % IS PHYS SHEET HASHCODE C RESULT IN V ACCUMULATOR IS DISPLAY SHEET LOC HASHCODE C T AND U ACCUMULATORS GET PHYS COL, ROW OFFSET. WAC=RRWACT YAC=RCLACT RETURN 2400 CONTINUE C YRMOD RETCD=1 IBGN=K+6 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ID1A,ID2A,IVALID) IF(IVALID.EQ.0)GOTO 9300 IF(LINE(LSTCHR).NE.',')GOTO 9300 IBGN=LSTCHR+1 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ID1B,ID2B,IVALID) IF(IVALID.EQ.0)GOTO 9300 IF(LINE(LSTCHR).NE.',')GOTO 9300 IBGN=LSTCHR+1 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ID1C,ID2C,IVALID) IF(IVALID.EQ.0)GOTO 9300 C C V1, V2, V3 ARE YR, MONTH, DAY FOR RETURN OF JULIAN DATE C IYR=XVBLS(ID1A,ID2A) IMO=XVBLS(ID1B,ID2B) IDA=XVBLS(ID1C,ID2C) C RETURN JULIAN DATE FROM Y, M, D GIVEN XAC=JULMDY(IYR,IMO,IDA) RETURN 2500 CONTINUE C JDATE RETCD=1 IBGN=K+6 LEND=IBGN+20 C GET V1 WHICH HAS VARIABLE WITH THE STRING IN IT CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ID1A,ID2A,IVALID) IF(IVALID.EQ.0)GOTO 9300 C RETURN JULIAN DATE NOW AFTER FETCHING FORMULA. IRX=(ID2A-1)*RRW+ID1A CALL WRKFIL(IRX,WRK,0) XAC=JULIAN(WRK) RETURN 2600 CONTINUE C JTOCH RETCD=1 IBGN=K+6 LEND=IBGN+20 C V1 = JULIAN DATE C V2 IS WHERE TO STORE ASCII DATE STRING AS FORMULA. CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ID1A,ID2A,IVALID) IF(IVALID.EQ.0)GOTO 9300 IF(LINE(LSTCHR).NE.',')GOTO 9300 IBGN=LSTCHR+1 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ID1B,ID2B,IVALID) IF(IVALID.EQ.0)GOTO 9300 IJUL=XVBLS(ID1A,ID2A) IRX=(ID2B-1)*RRW+ID1B CALL WRKFIL(IRX,WRK,0) DO 2502 N=1,110 2502 WRK(N)=0 CALL JULASC(IJUL,WRK,IYR,IMO,IDA) CALL WRKFIL(IRX,WRK,1) C WRITE THE FORMULA BACK OUT C RETURN T = MONTH C U =DAY C V = YEAR TAC=IMO UAC=IDA VAC=IYR RETURN 2700 CONTINUE C DATE RETCD=1 IBGN=K+5 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ID1A,ID2A,IVALID) IF(IVALID.EQ.0)GOTO 9300 IF(LINE(LSTCHR).NE.',')GOTO 9300 IBGN=LSTCHR+1 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ID1B,ID2B,IVALID) IF(IVALID.EQ.0)GOTO 9300 IF(LINE(LSTCHR).NE.',')GOTO 9300 IBGN=LSTCHR+1 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ID1C,ID2C,IVALID) IF(IVALID.EQ.0)GOTO 9300 IF(LINE(LSTCHR).NE.',')GOTO 9300 IBGN=LSTCHR+1 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ID1D,ID2D,IVALID) IF(IVALID.EQ.0)GOTO 9300 IYR=XVBLS(ID1A,ID2A) IMO=XVBLS(ID1B,ID2B) IDA=XVBLS(ID1C,ID2C) IRX=(ID2D-1)*RRW+ID1D CALL WRKFIL(IRX,WRK,0) DO 2702 N=1,110 2702 WRK(N)=0 IJUL=JULMDY(IYR,IMO,IDA) CALL JULASC(IJUL,WRK,IYR,IMO,IDA) CALL WRKFIL(IRX,WRK,1) 9300 RETURN END C SPLIT OFF MATRIX PARSING LOGIC HERE TO PICK OFF UP TO 3 MATRICES C COORDINATES C THIS ALLOWS US TO CALL ONE ROUTINE TO LOCATE UP TO 3 MATRIX C SPECIFICATIONS SEPARATED BY COMMAS. SUBROUTINE PMTX2(IRTCD,IMXX,LINE,IBGN,ID1A,ID2A,ID1B,ID2B, 1 IDXA,IDXB,IDYA,IDYB,IDBA,IDBB,IDCA,IDCB) LOGICAL*1 LINE(80) CALL GMTX(LINE,IBGN,LSTCHR,ID1A,ID2A,ID1B, 1 ID2B,RETCD) C GET LOC OF MATRIX A (MUST BE SQUARE) IBGN=LSTCHR+1 IF(RETCD.NE.0.OR.IMXX.LE.1)GOTO 1000 IF(LINE(LSTCHR).NE.',')GOTO 300 CALL GMTX(LINE,IBGN,LSTCHR,IDXA,IDXB,IDYA, 1 IDYB,RETCD) C GET LOC OF MATRIX X (RESULT). IBGN=LSTCHR+1 IF(RETCD.NE.0.OR.IMXX.LE.2)GOTO 1000 IF(LINE(LSTCHR).NE.',')GOTO 300 CALL GMTX(LINE,IBGN,LSTCHR,IDBA,IDBB,IDCA, 1 IDCB,RETCD) IBGN=LSTCHR+1 C GET LOC OF MATRIX B (AX=B), THE OTHER HALF OF OUR GIVENS C IF WE FALL TO HERE, ALL LOOKS OK, SO LEAVE RETCD ALONE. C HOWEVER IF ANY ERRS HAVE OCCURRED, RETCD IS ALREADY SET TO 3 C FOR ERROR... 1000 RETURN 300 CONTINUE RETCD=3 RETURN END C GET SPECS FOR A MATRIX (2 VARS SEPARATED BY COLONS) SUBROUTINE GMTX(LINE,IBGN,LSTCHR,ID1A,ID2A,ID1B, 1 ID2B,RETCD) LOGICAL*1 LINE(80) C REQUIRE END OF MATRIX NAME WITHIN 20 CHARS OF START. C SHOULD BE OK IN ALL REASONABLE CASES. LEND=IBGN+20 C GET LOC OF MATRIX A (MUST BE SQUARE) CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ID1A,ID2A,IVALID) IF(IVALID.EQ.0)GOTO 300 IF(LINE(LSTCHR).NE.':')GOTO 300 IBGN=LSTCHR+1 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ID1B,ID2B,IVALID) IF(IVALID.EQ.0)GOTO 300 1000 RETURN 300 RETCD=3 RETURN END C C VARY CONTROL ROUTINE C NOTE: THIS ROUTINE RELIES UPON HAVING ITS DATA AREAS REMAIN INTACT C ACROSS CALLS. IT MUST NOT BE IN AN OVERLAY SEGMENT OR THAT WILL FAIL C AND IT WILL NOT WORK. SPECIFICALLY IT EXPECTS THE AC ARRAY TO BE C SET CORRECTLY. SUBROUTINE VVARY(LINE,RETCD,K) INCLUDE 'VKLUGPRM.FTN' BYTE LINE(80) INTEGER RETCD LOGICAL*1 AVBLS(100,27),WRK(128),VBLS(8,RRW,RCL) INTEGER*2 TYPE(RRW,RCL),VLEN(9) REAL*8 XAC,XVBLS(RRW,RCL) EQUIVALENCE(XAC,AVBLS(1,27)) INTEGER*4 JVBLS(2,RRW,RCL) EQUIVALENCE(VBLS(1,1,1),JVBLS(1,1,1)) EQUIVALENCE(VBLS(1,1,1),XVBLS(1,1)) COMMON/V/TYPE,AVBLS,VBLS,VLEN C LOOP CONTROL FOR VARY FUNCTION. SET ZERO IN SPREDSHT AND C MUST BE SET POSITIVE HERE IF WE NEED ITERATIONS. C (IMPLEMENT FOR VAX ONLY) INTEGER KALKIT COMMON/VARYIT/KALKIT EXTERNAL SIGN REAL*8 SIGN LOGICAL*1 LAC(8) REAL*8 XVAC,VW EQUIVALENCE(LAC(1),XVAC) REAL *8 AC(26) REAL*8 DERIV(8) REAL*8 DEL(8) REAL*8 OLDVV INTEGER ACV(8) INTEGER CAC INTEGER CCNT(8) C UNCOMMENT THIS COMMON DECLARATION AND MOVE DATA STMTS INTO BLOCK C IN ORDER TO OVERLAY THIS... C COMMON/VRYDAT/AC,DERIV,DEL,CAC,CCNT,OLDVV C C ACV POINTS TO AC'S VARYING C CAC = CURRENT INDEX INTO ACV TO FIND AC BEING VARIED C AC IS LAST SET OF ACCUMULATORS SEEN C IF ACV ENTRY IS 0, MEANS NO AC TO VARY THERE. INTEGER LW,LX,LI ! LOGICAL W,X,I AC'S INTEGER LA ! LOGICAL A AC C DATA DERIV/8*1./,DEL/8*0./ DATA CAC/1/,CCNT/8*0/ DATA ACV/8*0/ DATA OLDVV/1./ C C PARSE ARGUMENTS FIRST C FIRST 2 ARGS ARE X AND A AC'S (OR GENERAL CELLS) C DEFAULT NO REDOING THIS... KALKIT=0 IBGN=K+5 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,LX,ID2A,IVALID) IF (IVALID.EQ.0)GOTO 9900 IF(LINE(LSTCHR).NE.',')GOTO 9900 IBGN=LSTCHR+1 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,LA,ID2B,IVALID) IF (IVALID.EQ.0)GOTO 9900 IF(LINE(LSTCHR).NE.',')GOTO 9900 IBGN=LSTCHR+1 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,LW,ID3B,IVALID) IF (IVALID.EQ.0)GOTO 9900 IF(LINE(LSTCHR).NE.',')GOTO 9900 IF(ID3B.NE.1)GOTO 9900 IBGN=LSTCHR+1 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,LI,ID3B,IVALID) IF (IVALID.EQ.0)GOTO 9900 IF(LINE(LSTCHR).NE.',')GOTO 9900 IF(ID3B.NE.1)GOTO 9900 C IBGN=LSTCHR+1 C LEND=IBGN+20 C LOOP OVER VALUES TO VARY NOW DO 99 N=1,8 99 ACV(N)=0 DO 100 N=1,8 C ALLOW UP TO 8 DIMENSIONS VARIATION IBGN=LSTCHR+1 LEND=IBGN+20 CALL VARSCN(LINE,IBGN,LEND,LSTCHR,ACV(N),ID3B,IVALID) IF (IVALID.EQ.0)GOTO 9900 IF(LINE(LSTCHR).NE.';')GOTO 110 IF(ID3B.NE.1)GOTO 9900 IBGN=LSTCHR+1 LEND=IBGN+20 100 CONTINUE 110 CONTINUE C NOW HAVE ALL AC POINTERS SET UP. C IF I IS NOW 0 OR NEGATIVE (ITER COUNT), REINITIALIZE. ASSIGN 111 TO LGET LLL=LI GOTO 500 111 CONTINUE IF(XVAC.GT.0.)GOTO 112 C INITIALIZE COUNTS LLL=LW C GET VALUE OF W FRACTION ASSIGN 114 TO LGET GOTO 500 114 CONTINUE VW=XVAC OLDVV=1. DO 113 N=1,8 CCNT(N)=0 DERIV(N)=1. DEL(N)=VW 113 CONTINUE CAC=1 C COPY CURRENT AC'S INTO SAVED ONES NOW. DO 117 N=1,26 LLL=N ASSIGN 118 TO LGET GOTO 500 118 AC(N)=XVAC 117 CONTINUE C AFTER THE INIT, JUST RETURN SINCE WE DON'T WANT TO TRY ANY ITERATIONS C WHEN ITER COUNT EXPIRES. KALKIT=0 RETURN C HERE WHEN ITER COUNT IS POSITIVE. 112 CONTINUE XVAC=XVAC-1. C UPDATE ITERATION COUNT NOW... KALKIT=XVAC ASSIGN 120 TO LPUT GOTO 600 120 CONTINUE C C NOW PROCEED WITH VARIATIONS... IF(CAC.LT.1.OR.CAC.GT.8)CAC=1 IF(CCNT(CAC).GE.1)GOTO 200 C CCNT WAS 0 SO WE DIDN'T GET OUR PARTIAL YET. VARY THE C AC WE'RE LOOKING AT (CAC [= CURRENT AC]) AND USE NEW VALUE OF C (X-A) FOR A NUMERICAL DERIVATIVE RESULT AFTER A RECALC OF SCREEN... CCNT(CAC)=1 C JUST STARTED THIS AC SO VARY BY THE APPROPRIATE DELTA AND C EXIT, ALLOWING PARTIAL TO BE COMPUTED NEXT TIME. LLL=LW ASSIGN 400 TO LGET GOTO 500 400 CONTINUE C GET W ACC. VALUE VW=XVAC IF(VW.EQ.0.)VW=.5 C GET CURRENT AC, FIND HOW TO UPDATE IT. LLL=ACV(CAC) IF(LLL.LE.0)GOTO 9900 ASSIGN 121 TO LGET GOTO 500 121 CONTINUE C NOW XVAC HAS CURRENT AC FOR THE ONE WE'RE VARYING C ADD DEL TO IT AND GET NEW ONE... C SAVE OLD X AC VALUE FOR NEXT ITERATION. C NOTE LLL IS STILL SET AT CURRENTLY VARYING AC C SAVE CURRENT (UNVARIED) VALUE TOO FOR NEXT TIME AROUND. OLDVV=XVAC IF(OLDVV.EQ.0.)OLDVV=1. IF(DEL(CAC).EQ.0.)DEL(CAC)=VW XVAC=XVAC*(1.+DEL(CAC)) C NOW ALL SET... JUST SAVE CURRENT AC'S AND CURRENT X,A C SO WE CAN GET DIFFERENCE NEXT TIME AROUND. C AC(ACV(CAC))=XVAC C STORE XVAC INTO REAL ACCUMULATORS TOO, SO IT'LL WORK C WHEN ALL AC'S ARE RELOADED BELOW. ASSIGN 412 TO LPUT GOTO 600 412 CONTINUE C AT 1000, RELOAD AC ARRAY FROM REAL AC'S... BUT GET OUR MODIFIED C ONE WE JUST STORED TOO. GOTO 1000 200 CONTINUE C COUNT HERE IS 1 SO WE ALREADY HAVE INFO NOW FOR OUR PARITAL C DERIVATIVE. COMPUTE IT AND VARY THE SELECTED AC USING IT C THEN STORE IT AND RESET CCNT(CAC) TO 0 CCNT(CAC)=0 C MUST GET NEW X AND A VALUES NOW. XVAC=XVBLS(LX,ID2A) IF(ID2A.NE.1)GOTO 201 LLL=LX ASSIGN 201 TO LGET C EXTRACT CURRENT X FROM AVBLS GOTO 500 201 CONTINUE XCURR=XVAC XVAC=XVBLS(LA,ID2B) IF(ID2B.NE.1)GOTO 202 LLL=LA ASSIGN 202 TO LGET GOTO 500 202 CONTINUE ACURR=XVAC C NOW WE HAVE ENOUGH TO COMPUTE PARTIAL DERIVATIVE WE NEED. IF(ACV(CAC).LE.0)GOTO 9900 IF(OLDVV.EQ.0.)OLDVV=AC(ACV(CAC)) IF(OLDVV.EQ.0.)OLDVV=1. DERIV(CAC)=((XCURR-ACURR)-(OLDX-OLDA))/(DEL(CAC)*OLDVV) C NEGATIVE FEEDBACK: IF GOING POSITIVE, MAKE IT NEGATIVE... C THIS IS NOT AN ANALYTICAL PROCEDURE ... JUST STEPS IN RIGHT DIRECTION C BY APPROPRIATE AMOUNT AND CONTINUES... C CLAMP VARIATION TO INITIAL PERCENTAGE IN W ACCUMULATOR LLL=LW C OBTAIN VALUE OF W VARIATION NOW...IN CASE USER SETS IT UP TO VARY ASSIGN 203 TO LGET GOTO 500 203 CONTINUE VW=XVAC C C TO ATTEMPT TO GET TO THE ZERO OF (X-A), WE REALLY NEED TO C DIVIDE BY THE DERIVATIVE. HOWEVER, IN CASES WHERE THE FUNCTION C IS NEAR ITS LOCAL MINIMUM AND SLOWLY VARYING, WE REALLY DON'T WANT C TO STEP FAR AWAY (IT MAY NEVER REACH THE ZERO). THEREFORE, TEST C TO SEE IF THE DERIVATIVE IS LARGE AND ALLOW DIVISION WHERE IT IS C OVER A SOMEWHAT ARBITRARY THRESHOLD (USED 1.0 BELOW), BUT C MULTIPLY BY DERIVATIVE OTHERWISE, SO THAT AS THE FUNCTION APPROACHES C ZERO SLOPE, THE STEPS GET FINER TO GET INTO THE LOCAL MINIMUM (IF ANY). C C FORCE NONZERO VARIATION JUST SO WE DON'T GET STUCK. IF(DERIV(CAC).EQ.0.)DERIV(CAC)=.01 IF(DABS(DERIV(CAC)).GT.1.)GOTO 405 DEL(CAC)=-(OLDX-OLDA)*VW*DERIV(CAC) GOTO 406 405 CONTINUE DEL(CAC)=-(OLDX-OLDA)*VW/DERIV(CAC) 406 CONTINUE C VERY IMPORTANT TO CLAMP SIZE OF STEPS HERE SO WE DON'T WILDLY DIVERGE C IN EARLY GOING. SMALL STEPS TAKE LONGER BUT GET TO MINIMA; LARGER ONES C WHERE WE DON'T KNOW FUNCTION SHAPE CAN BE DISASTERS. IF(DABS(DEL(CAC)).GT.VW)DEL(CAC)=VW*SIGN(DEL(CAC)) C NOW RESTORE AC'S TO OLD ONES AND VARY CURRENT ONE BY C THE NEW DELTA. IF(ACV(CAC).LE.0)GOTO 9900 C NEXT LINE MAKES ADJUSTMENT NEEDED TO OUR VARYING AC. AC(ACV(CAC))=OLDVV*(1.+DEL(CAC)) C NOW COPY SAVED OLD AC'S ONTO NEW ONES SO WE START WITH AC'S ALL AS THEY C WERE IN FIRST STEP SO WE VARY FROM INITIAL X, NOT FROM FIRST VARIED X C LOCATION... DO 204 N=1,26 XVAC=AC(N) LLL=N ASSIGN 205 TO LPUT GOTO 600 205 CONTINUE 204 CONTINUE C MOVE ON TO THE NEXT CAC VALUE CAC=CAC+1 IF(ACV(CAC).LE.0.OR.CAC.GT.8)CAC=1 1000 CONTINUE C SAVE OLD AC'S NOW FOR NEXT TIME DO 1100 N=1,26 LLL=N ASSIGN 1101 TO LGET GOTO 500 1101 AC(N)=XVAC 1100 CONTINUE C REMEMBER OLD X AND A VALUES SINCE WE LOOK FOR X=A AS C A SEARCH CONDITION. WE MUST ASSUME THAT SOME SORT OF C VARIATION OF ACCUMULATORS GIVEN WILL ALLOW US TO SATISFY C THE EQUATION (X-A)=0. OLDX=AC(LX) IF(ID2A.NE.1)OLDX=XVBLS(LX,ID2A) OLDA=AC(LA) IF(ID2B.NE.1)OLDA=XVBLS(LA,ID2B) RETURN 9900 CONTINUE RETCD=3 RETURN C PROC TO LOAD XVAC WITH VBLS(LLL) 500 CONTINUE DO 501 KKKKN=1,8 501 LAC(KKKKN)=AVBLS(KKKKN,LLL) GOTO LGET C PROC TO STORE XVAC INTO VBLS(LLL) 600 CONTINUE DO 601 KKKKN=1,8 601 AVBLS(KKKKN,LLL)=LAC(KKKKN) GOTO LPUT END REAL *8 FUNCTION SIGN(VAR) REAL*8 VAR C ALWAYS RETURN 1. OR -1. FOR THIS PROGRAM ... NEVER 0. SIGN=1. IF(VAR.LT.0.)SIGN=-1. RETURN END