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
	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(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
C ARGUMENTS COME IN IN ARGUMENTS IN LINE
C RESULTS GO INTO PERCENT (XAC) AND WHEREVER ELSE DESIRED...
	LOGICAL*1 FNAMS(6,10)
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
     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
	N=1
C **** BE SURE THE 2ND BOUND ON N IS THE SAME AS THE DIMENSION OF
C ****  FNAMS   **************************
	DO 100 N=1,10
	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),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.
	XAC=VAC*10000.+TAC*100.+UAC
	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
	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
