C	WESTERN MICHIGAN UNIVERSITY
C	NONLPS. FOR (FILE NAME ON LIBRARY DECTAPE)
C	THIS PROGRAM IS VERSION 2 OF SUMT PROGRAMMED AT THE RESEARCH
C	 ANALYSIS CORPORATION, MCLEAN, VIRGINIA.  THE MAIN PROGRAM
C	 AND SUBR. OUTPUT WAS SUBSTANTIALLY MODIFIED FOR WMU BY
C	 R.R. BARR.
C	FORWMU SUBR. USED:  TTYPTY, PRINTS
C	INTERNAL SUBR. USED:  BODY, ESTIM, EVALU, FEAS, FINAL,
C	 GRAD, INVERS, OUTPUT, OPT, STORE, PUNCH, REJECT,
C	 RHOCOM, PEVALU, SECORD, CONVRG
C	EXTERNAL SUBR. USED:  EQUATE, IDENTF, (THESE ARE GENERATED
C	 BY NONLPR.FOR)
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS  PUT IN BY WG
C
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M,MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON/OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON/VALUE/ F,G,P0,RSIGMA, RJ(200),RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO,EPSI,THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
	COMMON/IOBLK/IDLG,IRSP,INP,IOUT,IDEVI,IDEVO,ICODE,IB,NAME(2)
	COMMON/IOBLKA/NAM(2),IPJ,IPG,NCOPYS,ITYCH
	DIMENSION DAY(2)
	DATA EPSI,RHOIN,THETA0,RATIO,MMAX,IMAX
     1	/1.E-7,1.,1.E-4,16.,100,10/
	DATA NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10/
     1	3,2,1,1,1,1,1,1,1,1/
	IRERUN=0
	INP=4
	IOUT=6
	IDLG=-1
	IRSP=-4
	ITYCH=7
	CALL TTYPTY(ICODE)
C---------------SUBR. EQUATE IS IN NONLPP.FOR WHEN USER CHOOSES
C--------------- POLYNOM. COEFF. METHOD.  SEE NONLPR.FOR
C--------------- ST. 302+5.  NOTE AT THIS POINT NONLPX.FOR
C--------------- CONTAINS IDENT, PARAMS.  ALSO SEE NONLPP.FOR ST. 10+6.
C---------------WHEN USER CHOOSES FUNCTION EQUATION METHOD OF INPUT,
C--------------- SUBR. EQUATE IS IN NONLPX.FOR.  SEE NONLPR.FOR ST. 110,
C---------------COMMENTS IN FRONT OF ST. 120-1 ADN ST. 186+1.
C---------------M, N ARE OUTPUT FROM EQUATE THRU COMMON /SHARE/.
C--------------- MZ OUTPUT FROM EQUATE THRU COMMON /EQAL/.
	CALL EQUATE(1,IDUM,DUM,KDUM)
412	IERR=0
	WRITE(IDLG,305)
305	FORMAT(' ENTER INITIAL VALUES OF THE X VECTOR(6 PER LINE,'
     1	' SEPARATED BY COMMAS)',/)
      READ (IRSP,102) (X(I),I=1,N)
102	FORMAT(6F)
      NTCTR=0
      NP1=N+1
      NM1=N-1
	IF(ICODE.NE.0)IMAX=MMAX+1
	IF(IRERUN.EQ.1)GO TO 315
	WRITE(IDLG,400)
400	FORMAT(' DO YOU WISH SIMPLIFIED DIALOGUE AND PRINTOUT?',
     1	'(YES OR NO) ',$)
	READ(IRSP,301)ANS
	IF(ANS.NE.'YES')GO TO 315
	NT3=3
	GO TO 313
315	WRITE(IDLG,300)EPSI,RHOIN,THETA0,RATIO,MMAX,IMAX
300	FORMAT(' THE PARAMETER VALUES ARE AS FOLLOWS:',/,
     1	' EPSILON=',F13.7,2X,'R=',F,2X,'THETA ZERO=',F,/,
     1	' RATIO=',F,2X,'M MAX=',3X,I4,2X,'IM AX=',I4,/,
     1	' DO YOU WISH TO CHANGE THEM?',
     1	'(YES OR NO) ',$)
	READ(IRSP,301,END=999)ANS
301	FORMAT(A3)
	IF(ANS.NE.'YES')GO TO 306
	WRITE(IDLG,302)
302	FORMAT(' TYPE EPSILON,R,THETA ZERO,RATIO,M MAX,I MAX',/)
	READ(IRSP,4)EPSI,RHOIN,THETA0,RATIO,MMAX,IMAX
4	FORMAT(4F,2I)
304	FORMAT(5F)
306	WRITE(IDLG,310)NT1,NT4,NT5,NT7,NT9
310	FORMAT(' THE STANDARD OPTIONS ARE: ',5I3,//,
     1	' DO YOU WISH TO CHANGE THEM?(YES OR NO) ',$)
	READ(IRSP,301)ANS
	IF(ANS.NE.'YES')GO TO 312
	WRITE(IDLG,307)
307	FORMAT(' ENTER OPTIONS',/)
      READ (IRSP,103)NT1,NT4,NT5,NT7,NT9
103	FORMAT(11I)
312	WRITE(IDLG,320)
320	FORMAT(' PRINTOUT OPTION:',/'  TYPE 1 FOR STANDARD, 2 FOR ',
     1	'EXTENDED OR 3 FOR ABBREVIATED: ',$)
	READ(IRSP,103)NT3
	IF(NT3.GE.1.AND.NT3.LE.3)GO TO 313
	WRITE(IDLG,331)
	GO TO 312
313	WRITE(IDLG,322)
322	FORMAT(' TYPE:',/,' 1 IF TRIVIAL CONSTRAINTS: X(I)>=0 ,',
     1	' I=1,...,N',/,'   SHOULD BE INCLUDED IN PROBLEM',/,
     2	' 2 IF NOT',/)
	READ(IRSP,103)NT2
	IF(NT2.GE.1.AND.NT2.LE.2)GO TO 323
	WRITE(IDLG,331)
	GO TO 313
323	WRITE(IDLG,324)
324	FORMAT(' TYPE:',/,
     1	'	1 IF AT LEAST ONE NON-LINEAR CONSTRAINT',/,
     2	'	2 IF ALL LINEAR CONSTRAINTS',/,
     3	'	3 IF LINEAR FUNCTION AND CONSTRAINTS',/)
	READ(IRSP,103)NT10
	IF(NT10.GT.0.AND.NT10.LE.3)GO TO 330
	WRITE(IDLG,331)
331	FORMAT(' ?RESPONSE OUT OF LIMITS',/)
	GO TO 323
C---------------IOUT, IDLG ARE INPUT THRU COMMON /IOBLK/ OF
C--------------- SUBR. IDENTF.
330	CALL IDENTF
C	USER ID FROM STARTER SEQMENT
	CALL TIME(XT)
	CALL DATE(DAY)
	WRITE(IOUT,104)DAY,XT
104   FORMAT (1X,2A5,1X,A5)
	WRITE(IOUT,105) N,M,MZ,EPSI,RHOIN,THETA0,RATIO,MMAX,IMAX
105   FORMAT (1H0,5X,2HN=I3,6X,2HM=I3,6X,3HMZ=I3,//,
     1	' EPSILON=',F13.7,2X,'R=',F,2X,'THETA ZERO=',F,/,
     1	' RATIO=',F,2X,'M MAX=',3X,I4,2X,'IMAX=',I4,/)
	WRITE(IOUT,106) NT1,NT4,NT5,NT7,NT9,NT3,NT2,NT10

106	FORMAT(//,' STANDARD OPTIONS: '5I3,/,' SPECIAL OPTIONS:  ',3I3)
      NPHASE = 4
	IERR=0
C --- JUST TO GET AN INITIAL PRINTOUT
      CALL EVALU
      P0=0.0
      G=0.0
      H=0.0
      RSIGMA=0.0
      CALL OUTPUT(1)
      CALL STORE
      CALL FEAS
C     NPHASE=5 IS USED TO INDICATE NO FEASIBLE POINT EXISTS
	IF(NPHASE.EQ.5)GO TO 500
	NPHASE=2
      NTCTR=0
	CALL BODY
500	IF(NT3.NE.3)GO TO 510
	NT3=1
	CALL OUTPUT(1)
510	IF(IERR.LT.0)WRITE(IOUT,512)
512	FORMAT(/,' NO SOLUTION ENCOUNTERED',/)
	IF(IERR.EQ.0)WRITE(IOUT,518)
518	FORMAT(/,' SOLUTION ENCOUNTERED',/)
316	WRITE(IDLG,514)
514	FORMAT(//,' TYPE:',/,' 1 TO RESTART CURRENT PROBLEM',/,
     1	' 2 TO EXIT',/)
	IRERUN=1
	READ(IRSP,103)IBR
	IF(IBR.EQ.1)GO TO 412
	IF(IBR.EQ.2)GO TO 520
	WRITE(IDLG,331)
	GO TO 316
520	CALL RELEAS(IOUT)
	IF(IDEVO.EQ.'LPT')CALL PRINTS(NAM,2,1,1)
	CALL EXIT
999      END
      SUBROUTINE BODY
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M,MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON/OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON/VALUE/ F,G,P0,RSIGMA, RJ(200),RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO,EPSI,THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
      COMMON/CONPAR/ NF1, NF2, NF3
	COMMON/IOBLK/IDLG,IRSP,INP,IOUT,IDEVI,IDEVO,ICODE,IB,NAME(2)
      NF2=2
      NF3=2
      MN=0
      NUMINI=0
C     OPTION OF GETTING INITIAL RHO
      CALL RHOCOM
2     CALL PEVALU
C     (1)  MEANS ACCUMULATE MATRIX OF SECOND PARTIALS
3     CALL GRAD(1)
      CALL SECORD(1)
666   DO 10 I=1,N
10    DELX(I)=DELX0(I)
      CALL INVERS(1)
      CALL STORE
15    CALL OPT
	IF(IERR.LT.0)RETURN
17      GO TO (18,171,18),NT3
171    CALL OUTPUT(1)
18    IF(MMAX - NTCTR) 105,20,20
C     IN FEASIBILITY PHASE 4 MEANS FEAS ACHIEVED
20    GO TO (21,21,21,105), NSATIS
21    CALL CONVRG(N1)
	IF(IERR.LT.0)RETURN
      GO TO (25,3),N1
C     MINIMUM ACHIEVED IF N1=1
26	CONTINUE
25    GO TO (261,27,27),NT3
261    CALL OUTPUT(1)
C --- NUMBER OF MINIMA ACHIEVED INCREASED BY 1
27    NUMINI=NUMINI+1
      MN=0
      GO TO (300,28,28),NPHASE
28    CALL ESTIM
C     FINAL MIGHT HAVE BEEN CALLED BY ESTIM-CONVERGED IF N2=1
      GO TO (40,402,403),NT4
C     NT4=1 FINAL CONVERGENCE ON 0 ORDER ESTIMATES, NT4=2 CONVERGE ON FI
C     ORDER ESTIMATES, NT4=3 CONVERGE ON SECOND ORDER ESTIMATES.
40    CALL FINAL(NF1)
      GO TO (41,45),NF1
402   GO TO (41,45),NF2
403   GO TO (41,45),NF3
41    RETURN
45    RHO=RHO/RATIO
      CALL PEVALU
C     A VECTOR IS LEFT IN DELX(I) BY ESTIM
      IF(NUMINI-2) 3,46,46
46    GO TO (3,47,47),NT7
47    CALL GRAD(2)
      CALL OPT
	IF(IERR.LT.0)RETURN
48      GO TO (49,4811,49),NT3
4811	WRITE(IOUT,481)
481   FORMAT (6X,30HMOVED ON EXTRAPOLATION VECTOR )
      CALL OUTPUT(1)
49    GO TO 21
300   IF(G) 28,28,105
105	IERR=-1
	NT3=1
	CALL OUTPUT(1)
	RETURN	
C --- DUAL VALUE GREATER THAN 0 MEANS NO FEASIBLE POINT EXISTS
      END
      SUBROUTINE CONVRG(N1)
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M,MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON/OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON/VALUE/ F,G,P0,RSIGMA, RJ(200),RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO,EPSI,THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
      COMMON/TSW/NSWW
C     LEAVES A 1 IF MIN HAS BEEN FOUND, 2 IF NOT
C     DOT PRODUCT OF GRAD AND INVERS PRODUCT LEFT (DURING OPT) AT DOTT
	COMMON/IOBLK/IDLG,IRSP,INP,IOUT,IDEVI,IDEVO,ICODE,IB,NAME(2)
      N1=2
      IF(MN-1) 95,95,2
95    Q1=P0
2     GO TO (5,10,15),NT9
5     IF(ABS(DOTT)-EPSI) 25,25,50
10    IF(ABS(DOTT)-(P1-P0)/5.) 25,25,50
15    IF(ADELX.GT.EPSI) GO TO 50
25    N1=1
      GO TO 555
50    GO TO (55,60),NSWW
C     NSWW IS SET TO 2 IN SUBROOTIN TIME WHEN THE TIME LIMIT HAS BEEN
C     EXCEEDED.
55    IF(P0.GT.Q1) GO TO 56
555   Q1=P0
      RETURN
56    N1=1
      Q1=P0
	WRITE(IOUT,90)
      RETURN
60    CALL PUNCH
C     PUNCH IS ALWAYS CALLED IF PROGRAM QUITS BECAUSE TIME IS UP, EXCEPT
C     IN CASE OF AN INFINITE LOOP.
	WRITE(IOUT,99)
99    FORMAT (48H0**** TIME IS UP, CALLING EXIT FROM CONVRG. ****)
	NT3=1
	CALL OUTPUT(1)
	IERR=-1
	RETURN
90    FORMAT (/,' APPARENTLY ROUNDOFF ERRORS PREVENT A MORE ACCURATE',/,
     1	' DETERMINATION OF THE MINIMUM OF THIS SUBPROBLEM.')
      GO TO 555
      END
      SUBROUTINE ESTIM
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M,MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON/OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON/VALUE/ F,G,P0,RSIGMA, RJ(200),RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO,EPSI,THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
      COMMON/CONPAR/ NF1, NF2, NF3
	COMMON/IOBLK/IDLG,IRSP,INP,IOUT,IDEVI,IDEVO,ICODE,IB,NAME(2)
      CALL STORE
      Z9=SQRT(RATIO)
      Z1=(1./Z9+1./RATIO)
      Z2=Z1+1./Z9**3
      Z3=1./Z9**3
      Z4=RATIO+Z9
      Z5=Z9**3
      Z6=1./((RATIO-1.)*(Z9-1.))
      Z7=1./Z9
      Z8=1./(Z9-1.)
      RQ=1.0/SQRT(RHO)
      IF(NUMINI-2) 105,20,30
30      P0=(PR2-Z4*PR1+Z5*P1)*Z6
      G=(RATIO*G1-GR1)/(RATIO-1.)
      DO 35 I=1,N
35    X(I)=(XR2(I)-Z4*XR1(I)+Z5*X1(I))*Z6
      NP=NPHASE
      NPHASE=4
      CALL EVALU
      NPHASE=NP
	IF(NT3.EQ.3)GO TO 2000
	WRITE(IOUT,102)
102   FORMAT (/26H0    2ND  ORDER ESTIMATES )
      CALL OUTPUT(2)
C     CHECK TO SEE IF ESTIMATES HAVE CONVERGED
2000      GO TO (39,37,39),NPHASE
37    DO 38 J=1,M
      IF(RJ(J)) 375,38,38
375   IF(THETA0+RJ(J)) 39,38,38
38    CONTINUE
      GO TO (39,39,385),NT4
385   CALL FINAL(NF3)
39    CONTINUE
20      P0=(Z9*P1-PR1)*Z8
      G=(RATIO*G1-GR1)/(RATIO-1.)
      DO 25 I=1,N
25    X(I)=(Z9*X1(I)-XR1(I))*Z8
      NP=NPHASE
      NPHASE=4
      CALL EVALU
      NPHASE=NP
	IF(NT3.EQ.3)GO TO 2001
	WRITE(IOUT,101)
101   FORMAT (/26H0    1ST  ORDER ESTIMATES )
      CALL OUTPUT(2)
C     CHECK TO SEE IF ESTIMATES HAVE CONVERGED
2001      GO TO (49,47,49),NPHASE
47    DO 48 J=1,M
      IF(RJ(J)) 475,48,48
475   IF(RJ(J)+THETA0) 49,48,48
48    CONTINUE
      GO TO (49,485,49),NT4
485   CALL FINAL(NF2)
49    CONTINUE
105      IF(M) 405,405,404
404   DO 40 J=1,M
40    RJ(J)=RHO/RJ1(J)**2
405   IF(MZ) 43,43,41
41    DO 42 J=1,MZ
      MNJ=M+J
42    RJ(MNJ)=2.*RJ1(MNJ)*RQ
43    GO TO (44,46),NT2
44    DO 45 I=1,N
45    X(I)=RHO/X1(I)**2
46	CONTINUE
	IF(NT3.EQ.3)GO TO 2002
	WRITE(IOUT,103)
103   FORMAT (/25H0  LAGRANGE MULTIPLIERS  )
	CALL OUTPUT(2)
2002      CALL REJECT
      IF(NUMINI-2) 63,71,50
50    GO TO (63,715,60),NT7
C     SECOND ORDER MOVE FOR NEXT MINIMUM
60    DO 62 I=1,N
62    DELX(I)=Z1*X1(I)-Z2*XR1(I)+Z3*XR2(I)
63    PR2=PR1
      GR2=GR1
      PR1=P1
      GR1=G1
      DO 65 I=1,N
      XR2(I)=XR1(I)
65    XR1(I)=X1(I)
66    RETURN
71    GO TO (63,715,715),NT7
715   DO 72 I=1,N
72    DELX(I)=(X1(I)-XR1(I))*Z7
      GO TO 63
      END
      SUBROUTINE EVALU
C---------------NT2 INPUT THRU COMMON /OPTNS/ OF MAIN PROG. OF
C--------------- NONLPS.FOR.  M, N ARE INPUT THRU COMMON /SHARE/
C--------------- FROM SUBR. EQUATE IN MAIN PROG.  MZ IS INPUT THRU
C--------------- COMMON /EQAL/ FROM SUBR. EQUATE.
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M,MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON/OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON/VALUE/ F,G,P0,RSIGMA, RJ(200),RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO,EPSI,THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
C     COMBINES EVALU-PEVALU- AND SPECFE
      H=0.0
      RSIGMA=0.0
      F=0.0
      NSATIS=2
      GO TO (100,200,300,400),NPHASE
C     =1  FEASIBILITY
C     =2  NORMAL
C     =3  GUESS
C     =4  ALL FUNCTIONS ARE TO BE EVALUATED
C     FEASIBILITY
100   GO TO (105,111),NT2
C     NON-NEGATIVES INCLUDED
105   DO 106 I=1,N
      IF(X(I)) 555,555,106
106   RSIGMA=RSIGMA+RHO/X(I)
111   IF(M.EQ.0) GO TO 116
      DO 115 J=1,M
	CALL EQUATE(2,J,RJ(J),KDUM)
      IF(RJ1(J).LE.0.0) GO TO 112
      IF(RJ(J).GT.0.0) GO TO 113
C     VIOLATION OF A PREVIOUSLY SATISFIED CONSTRAINT
      GO TO 555
112   IF(RJ(J).GT.0.0) GO TO 114
C     ALL VIOLATED CONSTRAINTS ADDED INTO OBJECTIVE FUNCTION
      F=F-RJ(J)
      GO TO 115
113   RSIGMA=RSIGMA+RHO/RJ(J)
      GO TO 115
C     INDICATES SATISFACTION OF CONSTRAINT (1 OR MORE)
114   NSATIS=1
      RSIGMA=RSIGMA+RHO/RJ(J)
115   CONTINUE
116   CONTINUE
C     EQUALITIES NOT COMPUTED IN FEAS. PHASE
      P0=F+RSIGMA
      G=F-RSIGMA
      RETURN
C     REGULAR PHASE
200   GO TO (205,211),NT2
C     NON NEGATIVITIES INCLUDED
205   DO 206 I=1,N
      IF(X(I)) 555,555,206
206   RSIGMA=RSIGMA+RHO/X(I)
211   IF(M.EQ.0) GO TO 216
      DO 215 J=1,M
	CALL EQUATE(2,J,RJ(J),KDUM)
      IF(RJ(J).LE.0.0) GO TO 555
      RSIGMA=RSIGMA+RHO/RJ(J)
215   CONTINUE
C     EVALUATE AND ADD IN EQUALITY CONSTRAINTS
216	CALL EQUATE(2,0,F,KDUM)
      IF(MZ) 250,250,220
220   DO 240 I=1,MZ
      J=I+M
	CALL EQUATE(2,J,RJ(J),KDUM)
C     ADD INTO THIRD TERM OF P FUNCTION
      H=H+(RJ(J))**2
240   CONTINUE
      RS=SQRT(RHO)
      H=H/RS
250   P0=RSIGMA+H
      P0=F+P0
      G=-RSIGMA+2.0*H
      G=G+F
C     DUAL VALUE
      RETURN
C     GUESS PHASE NOT CODED
300   RETURN
C --- STRAIGHT FUNCTION EVALUATION  ( MAIN + FEAS ONLY)
400   CONTINUE
411   IF(M.EQ.0) GO TO 416
      DO 415 I=1,M
	CALL EQUATE(2,I,RJ(I),KDUM)
415   CONTINUE
416	CALL EQUATE(2,0,F,KDUM)
C     EQUALITY RESTRAINTS
      IF(MZ) 450,450,420
420   DO 440 I=1,MZ
      KZ=M+I
440	CALL EQUATE(2,KZ,RJ(KZ),KDUM)
450   RETURN
C     CONSTRAINTS VIOLATED NOT SO BEFORE
555   NSATIS=3
      P0=10.E35
      RETURN
      END
      SUBROUTINE FEAS
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M,MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON/OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON/VALUE/ F,G,P0,RSIGMA, RJ(200),RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO,EPSI,THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
	COMMON/IOBLK/IDLG,IRSP,INP,IOUT,IDEVI,IDEVO,ICODE,IB,NAME(2)
      NPHASE=1
      GO TO (10,15),NT2
10    NFIX=1
      DO 12 I=1,N
      IF(X(I)) 11,11,12
11    NFIX=2
      X(I)=1.E-05
12    CONTINUE
      GO TO (15,13),NFIX
13    NPHASE=4
      CALL EVALU
C     JUST GET ALL CONSTRAINTS EVALUATED
      NPHASE=1
	IF(NT3.EQ.3)GO TO 15
	WRITE(IOUT,998)
998   FORMAT (1H0,'  MADE VIOLATED NON-NEGATIVITIES SLIGHTLY POSITIVE')
      CALL OUTPUT(2)
15    IF(M) 18,18,16
16    DO 17 I=1,M
      IF(RJ(I)) 25,25,17
17    CONTINUE
175      G=0.0
	CALL EQUATE(2,0,F,KDUM)
	IF(NT3.EQ.3)GO TO 18
	WRITE(IOUT,997)
997   FORMAT (51H0*****THE FEASIBLE STARTING POINT TO BE USED IS ...)
      CALL OUTPUT(2)
18    RETURN
25    CALL BODY
      DO 26 I=1,M
      IF(RJ(I)) 28,28,26
26    CONTINUE
      GO TO 175
28	WRITE(IOUT,999)
999   FORMAT ('   THIS PROBLEM POSSESSES NO FEASIBLE STARTING POINT')  
C     TO INDICATE TO MAIN TO START ON NEXT PROBLEM.
      NPHASE=5
	NT3=3
	IERR=-1
      GO TO 18
      END
      SUBROUTINE FINAL(N2)
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M,MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON/OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON/VALUE/ F,G,P0,RSIGMA, RJ(200),RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO,EPSI,THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
C     FINAL CONVERGENCE CRITERION
C     1 LEFT IF CONVERGENCE MET   2 IF NOT
      GO TO (5,15,20),NT5
5     EPSIL=ABS(F/G-1.)
      IF(EPSIL-THETA0) 30,30,50
15    IF(RSIGMA-THETA0) 30,30,50
20    IF(NUMINI-1) 30,22,22
22    PEST=PR1-(PR1-P0)/(1.-1./SQRT(RATIO))
      EPSIL=ABS(PEST/G-1.)
      IF(EPSIL-THETA0) 30,50,50
30    N2=1
      GO TO (100,200),NT6
200   CALL PUNCH
      GO TO 100
50    N2=2
100   RETURN
      END
      SUBROUTINE GRAD(IS)
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M,MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON/OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON/VALUE/ F,G,P0,RSIGMA, RJ(200),RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO,EPSI,THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
C     IF (IS=1) ACCUM. MATRIX OF 2ND PARTIALS.  IF (IS=2) DONT.
      GO TO (100,1055),IS
100   DO 105 I=1,N
      DO 105 J=1,I
105   A(I,J)=0.
1055  DO 106 I=1,N
106   DELX0(I)=0.
C    THIS SECTION WORKS CORRECTLY IN FEASIBILITY PHASE AS WELL AS NORMAL
      GO TO (110,120),NT2
110   DO 115 I=1,N
112   DELX0(I)=-RHO/X(I)**2
      GO TO (113,115),IS
113   A(I,I)=-2.*DELX0(I)/X(I)
115   CONTINUE
120   CONTINUE
      IF(M.LE.0) GO TO 141
121   DO 130 K=1,M
	CALL EQUATE(3,K,DUM,KDUM)
      IF(RJ(K)) 1130,1130,122
122   TT=RHO/RJ(K)**2
      DO 125 I=1,N
      IF(DEL(I)) 123,125,123
C     IF DEL(I)=0 SKIP ALL THE FOLLOWING COMPUTATION INVOLVING * BY DEL(
123   T=TT*DEL(I)
      DELX0(I)=DELX0(I)-T
      GO TO (124,125), IS
124   T=2.*T/RJ(K)
      DO 1245 JJ=1,I
      IF(DEL(JJ)) 1244,1245,1244
1244  A(I,JJ)=A(I,JJ)+T*DEL(JJ)
1245  CONTINUE
125   CONTINUE
130   CONTINUE
C     EQUALITY CHANGES FOR GRAD
141   IF(MZ.LE.0) GO TO 131
      GO TO (131,140,131),NPHASE
140   RQ=2./SQRT(RHO)
      DO 150 J=1,MZ
      K=M+J
	CALL EQUATE(3,K,DUM,KDUM)
      TT=RQ*RJ(K)
      DO 145 I=1,N
      IF(DEL(I).EQ.0.0) GO TO 145
      DELX0(I)=DELX0(I)+DEL(I)*TT
      GO TO (144,145),IS
144   T=RQ*DEL(I)
      DO 1445 JJ=1,I
      IF(DEL(JJ)) 1444,1445,1444
1444  A(I,JJ)=A(I,JJ)+T*DEL(JJ)
1445  CONTINUE
145   CONTINUE
150   CONTINUE
131   GO TO (132,135),IS
132   DO 133 I=1,N
133   DIAG(I)=A(I,I)
135   GO TO (136,138,136),NPHASE
C     LEAVES NEGATIVE GRADIENT IN DELP
136   DO 137 I=1,N
137   DELX0(I)=-DELX0(I)
200   ADELX=0.
      DO 580 I=1,N
580   ADELX=ADELX+DELX0(I)**2
      ADELX=SQRT(ADELX)
      RETURN
138	CALL EQUATE(3,0,DUM,KDUM)
      DO 139 I=1,N
139   DELX0(I)=-DELX0(I)-DEL(I)
C     LEAVES THE NEG. GRAD OF P IN DELX0
      GO TO 200
C     ALL VIOLATED CONSTRAINT GRADS ADDED TO OBJ. FUNCTION
1130  DO 1135 I=1,N
      IF(DEL(I)) 1131,1135,1131
1131  DELX0(I)=DELX0(I)-DEL(I)
1135  CONTINUE
      GO TO 130
      END
      SUBROUTINE INVERS(NSME)
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M,MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON/OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON/VALUE/ F,G,P0,RSIGMA, RJ(200),RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO,EPSI,THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
      DIMENSION B(50)
      EQUIVALENCE (B(1),DELX)
C     GETTING THE PRODUCT INVERSE USING THE CROUT METHOD AND TAKING ADVA
C     TAGE OF THE SYMMETRY OF THE A MATRIX
C     IF A DIVISOR OF ZERO IS GENERATED, THEN DELX IS SET EQUAL TO DELX0
C     IF NSME=1 HAVE NEW A MATRIX  IF NSME=2  COMPUTE INVERSE PRODUCT US
C     A AS LEFT FROM LAST TIME.
	COMMON/IOBLK/IDLG,IRSP,INP,IOUT,IDEVI,IDEVO,ICODE,IB,NAME(2)
      GO TO (11,2),NSME
2     GO TO (20,2070),KSW
11    KSW=1
      IF(A(1,1)) 66,66,150
150   A(1,1)=1./A(1,1)
200   DO 1 I=2,N
1     A(1,I)=A(1,I)*A(1,1)
      DO 50 J=2,N
      JM1=J-1
4     T=0.
      DO 5 I=1,JM1
      IF(A(I,J)) 44,5,44
44    T=T+A(J,I)*A(I,J)
5     CONTINUE
      A(J,J)=A(J,J)-T
      IF(A(J,J)) 77,77,8
8     A(J,J)=1./A(J,J)
      IF(J.EQ.N) GO TO 20
      JP1=J+1
      DO 10 L=JP1,N
      T=0.
      DO 15 I=1,JM1
      IF(A(I,J)) 14,15,14
14    T=T+A(L,I)*A(I,J)
15    CONTINUE
      A(L,J)=A(L,J)-T
      A(J,L)=A(L,J)*A(J,J)
10    CONTINUE
50    CONTINUE
20    B(1)=B(1)*A(1,1)
      DO 55 J=2,N
      T=0.
      JM1=J-1
      DO 25 I=1,JM1
      IF(A(J,I)) 24,25,24
24    T=T+A(J,I)*B(I)
25    CONTINUE
      B(J)=(B(J)-T)*A(J,J)
55    CONTINUE
      DO 100 I=1,NM1
      NMK=N-I
      DO 75 J=1,I
      L=NP1-J
      IF(A(NMK,L)) 74,75,74
74    B(NMK)=B(NMK)-A(NMK,L)*B(L)
75    CONTINUE
100   CONTINUE
102   FORMAT (7E17.8)
2070	GO TO (2072,2071,2072),NT3
2071	WRITE(IOUT,103)
	WRITE(IOUT,102)(DELX0(I),I=1,N)
      GO TO (2073,2072),KSW
2073	WRITE(IOUT,104)
	WRITE(IOUT,102)(DELX(I),I=1,N)
103   FORMAT (1H0,6X,12HDEL P VECTOR)
104   FORMAT (1H0,6X,24HSECOND ORDER MOVE VECTOR)
2072  RETURN
66    J=1
77	WRITE(IOUT,999) J,A(J,J)
999   FORMAT (/,' POSSIBLE NON-CONVEX PROBLEM,',
     1	' THE',I3,'TH DISC.=',E12.4)
      KSW=2
      DO 78 I=1,N
78    DELX(I)=DELX0(I)
      GO TO 2070
      END
      SUBROUTINE OUTPUT(K)
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M,MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON/OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON/VALUE/ F,G,P0,RSIGMA, RJ(200),RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO,EPSI,THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
	COMMON/IOBLK/IDLG,IRSP,INP,IOUT,IDEVI,IDEVO,ICODE,IB,NAME(2)
	ANS=0
	ANSA=0
	ANSB=0
      NZ=M+MZ
      GO TO (1,2,1),K
1	IF(NTCTR.EQ.0)GO TO 5
	WRITE(IOUT,1000)
1000	FORMAT(//,2(1X,8('*'),/))
	WRITE(IOUT,999)NTCTR,F,DOTT,RHO
999   FORMAT (6X,6HPOINT=I4,5X,2HF=E12.4,/,
     1	6X,5HDOTT=E12.4,6X,4HRHO=E12.4)
	IF(K.NE.3)GO TO 5
	IF(IDEVO.NE.'TTY')GO TO 5
	WRITE(IDLG,980)
980	FORMAT(/,' DO YOU WANT A SUMMARY?(YES OR NO) ',$)
	READ(IRSP,982)ANS
982	FORMAT(A3)
	IF(ANS.NE.'YES')GO TO 7
5	IF(NT3.NE.3)WRITE(IOUT,998)P0,G,RSIGMA,H,ADELX,NPHASE
998   FORMAT (/,6X,2HP=E12.4,5X,2HG=E12.4,5X,7HRSIGMA=,
     1	E12.4,/,6X,2HH=E12.4,11H MAGNITUDE=E12.4,6X,6HPHASE=I2)
2	WRITE(IOUT,997)(X(I),I=1,N)
997   FORMAT (/,6X,25HTHE CURRENT VALUE OF X IS/(4E17.4))
	IF(M+MZ.EQ.0)GO TO 12
	WRITE(IOUT,994)
994   FORMAT (/,6X,21HTHE CONSTRAINT VALUES)
      GO TO (3,4),NT2
3	WRITE(IOUT,993)
993   FORMAT (28X,34HNOT INCLUDING THE NON-NEGATIVITIES)
4	WRITE(IOUT,996)(RJ(I),I=1,NZ)
996   FORMAT (4E17.4)
12	IF(NT3.EQ.0)GO TO 11
	IF(K.NE.3)RETURN
7	IF(IDEVO.NE.'TTY')RETURN
	WRITE(IDLG,984)
984	FORMAT(/,' DO YOU WISH TO CHANGE THE PARAMETERS?(YES OR NO) ',$)
	READ(IRSP,982)ANSA
	IF(ANSA.NE.'YES')GO TO 8
	WRITE(IDLG,986)EPSI,RHOIN,THETA0,RATIO,MMAX,IMAX
986	FORMAT(' EPSILON=',F,' R=',F,' THETA ZERO=',F,/,
     1	' RATIO=',F,'    M MAX=',I4,'    I MAX=',I4,//,
     1	' ENTER EPSILON,R,THETA ZERO,RATIO,M MAX,I MAX',/)
	READ(IRSP,988)EPSI,RHOIN,THETA0,RATIO,MMAX,IMAX
988	FORMAT(4F,2I)
	RETURN
8	WRITE(IDLG,9)
9	FORMAT(' DO YOU WISH TO TERMINATE SOLUTION?(YES OR NO) ',$)
	READ(IRSP,982)ANSB
	IF(ANSB.NE.'YES')RETURN
	IF(ANS.EQ.'YES')GO TO 11
	NT3=0
	GO TO 5
11	NT3=1
	IERR=-1
	RETURN
	END
      SUBROUTINE OPT
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M,MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON/OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON/VALUE/ F,G,P0,RSIGMA, RJ(200),RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO,EPSI,THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
	COMMON/IOBLK/IDLG,IRSP,INP,IOUT,IDEVI,IDEVO,ICODE,IB,NAME(2)
      KSW=1
      N405=1
      DOTT=0.
      DO 2027 J=1,N
 2027 DOTT=DOTT+DELX(J)*DELX0(J)
      IF (DOTT)2065,2065,2070
 2065 DO 2067 I=1,N
 2067 DELX(I)=DELX0(I)
 2070 CONTINUE
      N404=0
      MN=MN+1
C MN IS NOW NUMB. OF POINTS AFTER MIN ACHIEVED
      NTCTR=NTCTR+1
	IF(NTCTR.EQ.(NTCTR/IMAX)*IMAX)CALL OUTPUT(3)
      DO 10 I=1,N
   10 X2(I)=X(I)
      PX1=P0
   11 N401=0
   15 N401=N401+1
      DO 20 I=1,N
   20 X(I)=X2(I)+DELX(I)
      CALL EVALU
C 1 MEANS SATIS. OF CONSTRAINT NT. PREV.  2 MEANS NO CHANGE.  3 MEANS VI
C     IF POINT IS NOT FEASIBLE GIVE IT AN ARBITRARILY HIGH VALUE.
      GO TO (472,25,26),NSATIS
   26 PX2=10.E35
      P0=10.E35
      GO TO 30
   25 CONTINUE
      PX2=P0
      IF (PX1-PX2)30,30,40
 30   IF (N401-2)35,33,33
 33   DO 34 I=1,N
   34 X1(I)=X(I)
      P1=PX2
      GO TO 529
C ONLY ONE POINT SO FAR COMPUTED
   35 DO 37 I=1,N
   37 X3(I)=X2(I)
      PREV3=PX1
      GO TO 430
   40 DO 45 I=1,N
      X3(I)=X2(I)
      X2(I)=X(I)
   45 DELX(I)=1.61803399E0*DELX(I)
      PREV3=PX1
      PX1=PX2
      GO TO 15
C FIBONACCI   METHOD
C B VECTOR  GOES TO X1(I)
  475 P0=1.E36
      N404=N404+1
  430 DO 431 I=1,N
  431 X1(I)=X(I)
      P1=P0
      DO 500 I=1,N
      X(I)=.38196601E0*(X1(I)-X3(I))+X3(I)
  500 X2(I)=X(I)
      CALL EVALU
      GO TO (472,505,4750),NSATIS
 4750 IF (N404.LT.30) GO TO 475
C--IT IS POSSIBLE NO FEASIBLE POINT EXIST, IF NOT TRY MOVING ON DELX0.
C--  IF IT IS NOT POSSIBLE TO MOVE ON DELX0 THEN WE MUST BE AT A
C--  S3LUTION OF NLP PROBLEM.
      IF (N404.GT.100) GO TO 4752
 4755 DO 4751 I=1,N
      IF (ABS(ABS(X3(I)/X1(I))-1.)  .GT.1.E-7) GO TO 475
 4751 CONTINUE
 4752 GO TO (4753,4754),N405
 4753 N405=2
C--TRY TO MOVE ON GRADIENT.
      NTCTR=NTCTR-1
      MN=MN-1
      GO TO 2065
 4754	WRITE(IOUT,9999)
9999  FORMAT (106H  OPT CAN'T FIND A FEASIBLE POINT, THAT GIVES A LOWER
     1VALUE OF THE P-FUNCTION, FINAL CONVERGENCE ASSUMED.  )
	NT3=1
      CALL OUTPUT(2)
C	DONT SAY SOLUTION ENCOUTERED TWICE
	IERR=1
	RETURN
C--IN 7040 IBSYS THIS CAUSES CONTROL TO REVERT TO MAIN.
  505 CONTINUE
      N404=0
      PX1=P0
      DO 510 I=1,N
  510 X(I)=0.38196601E0*(X1(I)-X2(I))+X2(I)
      CALL EVALU
      GO TO (472,515,4755),NSATIS
  515 PX2=P0
      N401=1
  516 N401=N401+1
 5162 IF (N401-25)5183,518,518
  518 KSW=2
      IF (N401-40)5181,535,535
 5181 DO 5182 I=1,N
      IF (ABS(X2(I)/X(I)-1.0)  .GE.1.E-7) GO TO 5183
 5182 CONTINUE
      GO TO 535
 5183 IF (ABS(PX1/PX2-1.)  .LE.1.E-7) GO TO 535
      IF (PX1-PX2)520,535,525
C FROM LEFTORIGHT  X3(I)(PREV3)X2(I)(PX1)X(I)PX2 X1(I)P1
  520 DO 5205 I=1,N
 5205 X1(I)=X(I)
C THROW AWAY RIGHT PART
      P1=PX2
  522 DO 523 I=1,N
C POINTXP1 BECOMES XP2
  523 X(I)=.38196601E0*(X1(I)-X3(I))+X3(I)
C TEMPORARILY IN  X STORAGE
      CALL EVALU
      GO TO (472,524,475),NSATIS
  524 CONTINUE
      PX2=PX1
C SWITCH VECTORS TO PROPER POSITION
      PX1=P0
      DO 5245 I=1,N
      XX=X2(I)
      X2(I)=X(I)
 5245 X(I)=XX
      GO TO 516
C LEFT SIDE  TOSSED  AWAY
  525 DO 526 I=1,N
      X3(I)=X2(I)
  526 X2(I)=X(I)
      PREV3=PX1
      PX1=PX2
  529 DO 530 I=1,N
  530 X(I)=0.38196601E0*(X1(I)-X2(I))+X2(I)
      CALL EVALU
      GO TO (472,531,475),NSATIS
  531 CONTINUE
      PX2=P0
      GO TO 516
C FIBONNACCI  POINTS HAVE EQUAL VALUE NOW COMPUTE MIDPOINT
  535 DO 536 I=1,N
      DELX0(I)=X(I)
  536 X(I)=(DELX0(I)+X2(I))*0.5
      PP1=PX2
      CALL EVALU
      GO TO (5361,537),KSW
 5361 IF (ABS(P0/PX1-1.)  .GT.1.E-7) GO TO 538
  537 RETURN
  538 DO 539 I=1,N
  539 X(I)=DELX0(I)
      GO TO 520
C ARE  WE NOW IN FEASIBILITY PHASE
  472 DO 1450 I=1,M
      IF (RJ(I))1460,1460,1450
 1450 CONTINUE
      NSATIS=4
      RETURN
C---  PROBLEM HAS BECOME FEASIBLE
C ---  P - FUNCTION CHANGES IF A CONSTRAINT BECOMES FEASIBLE
 1460 MN=0
      DO 1461 I=1,M
 1461 RJ1(I)=RJ(I)
      RETURN
      END
      SUBROUTINE STORE
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M, MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON /OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON /VALUE/ F,G,P0,RSIGMA, RJ(200), RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO, EPSI, THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT,PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
  605 DO 610 I=1,N
  610 X1(I)=X(I)
      MMZ=M+MZ
      DO 615 J=1,MMZ
  615 RJ1(J)=RJ(J)
      P1=P0
      F1=F
      G1=G
      RSIG1=RSIGMA
      H1=H
      RETURN
      END
      SUBROUTINE PUNCH
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M, MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON /OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON /VALUE/ F,G,P0,RSIGMA, RJ(200), RHO
      COMMON/CRST/ DELX(50), DELX0(50),RHOIN,RATIO,EPSI, THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX

C  A ROUTINE TO PUNCH OUT  PARAMETER CARD, STARTING P0INT, AND OPTION
C   CARD FOR RESTARTS.
C  THIS ROUTINE IS CALLED IF NT6=2.
	COMMON/IOBLK/IDLG,IRSP,INP,IOUT,IDEVI,IDEVO,ICODE,IB,NAME(2)
	WRITE(IOUT,4)
    4 FORMAT('1  VALUES FOR RESTART ')
	WRITE(IOUT,1)EPSI,RHO,THETA0,RATIO,MMAX,M,N
    1 FORMAT (E12.4,E12.6,0P2E12.4,I12  ,2I4)
	WRITE(IOUT,2)(X(I),I=1,N)
    2 FORMAT (6E12.5)
      NT1=3
C SET  RHO OPTION SO THIS VALUE OF RHO WILL BE USE FOR THE RESTART.
	WRITE(IOUT,3)NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10
    3 FORMAT (10I7)
      RETURN
      END
      SUBROUTINE REJECT
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M, MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON /OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON /VALUE/ F,G,P0,RSIGMA, RJ(200), RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO, EPSI, THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
  505 DO 510 I=1,N
  510 X(I)=X1(I)
      MMZ=M+MZ
      DO 515 J=1,MMZ
  515 RJ(J)=RJ1(J)
      P0=P1
      RSIGMA=RSIG1
      G=G1
      F=F1
      H=H1
      RETURN
      END
      SUBROUTINE RHOCOM
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M, MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON /OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON /VALUE/ F,G,P0,RSIGMA, RJ(200), RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO, EPSI, THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT,PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
C     SUBROUTINE TO COMPUTE INITIAL RHO VALUE
C     CONTROLLED BY COL. 7 ON OPTION CARD
      GO TO (9015,9004,8000,9050),NT1
 8000 RHO=RHOIN
 9000 IF (RHO)9001,9001,9002
 9001 RHO=1.
 9002 RETURN
 9004 NPAR1=1
 9005 RHO=1.
C     2 MEANS  RHD WHICH  MINIMIZES GRADIENT  MAG.
      CALL GRAD(2)
      DO 9007 I=1,N
 9007 PGRAD(I)=DELX0(I)
      RHO=2.
      CALL GRAD(2)
      DO 9009 I=1,N
      DELX0(I)=DELX0(I)-PGRAD(I)
 9009 PGRAD(I)=PGRAD(I)-DELX0(I)
      GO TO (9019,9029),NPAR1
 9019 DOT1=0.
      DOT2=0.
      DO 9010 I=1,N
      DOT1=DOT1+DELX0(I)*PGRAD(I)
 9010 DOT2=DOT2+DELX0(I)**2
      RHO=ABS(DOT1/DOT2)
      GO TO 9000
C  3 MEANS COMPUTE RHO SO AS TO MINIMINIZE DEL P(/DDP/1.)DEL P
 9015 NPAR2=1
 9028 NPAR1=2
C     USE DF  AND DR  SUBROOTIN 
      GO TO 9005
 9029 RHO=1.
C     ASSUME   SIGMA TERM IS CONSID. GRTER THAN  F TERM
      CALL SECORD(2)
      DO 9032 I=1,N
 9032 DELX(I)=PGRAD(I)
      CALL INVERS(1)
      DO 9034 I=1,N
      X1(I)=DELX(I)
 9034 DELX(I)=DELX0(I)
      CALL SECORD(2)
      CALL INVERS(1)
      DO 9035 I=1,N
 9035 XR2(I)=DELX(I)
      GO TO (9040,9060),NPAR2
 9040 DOT1=0.
      DOT2=0.
      DO 9043 I=1,N
      DOT1=DOT1+PGRAD(I)*X1(I)
 9043 DOT2=DOT2+DELX0(I)*XR2(I)
      RHO=SQRT(ABS(DOT1/DOT2))
      GO TO 9000
 9050 NPAR2=2
C     RHO MINIMIZES  2ND ORDER MOVE
      GO TO 9028
C     USES INTERNAL  SUB. TO COM  /DDP/-1 DF AND /DDP/- DR
 9060 DOT1=0.
      DOT2=0
      DO 9063 I=1,N
      DOT1=X1(I)**2+DOT1
 9063 DOT2=X1(I)*XR2(I)+DOT2
      RHO=ABS(DOT1/DOT2)
      GO TO 9000
      END
      SUBROUTINE PEVALU
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M, MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON/OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON /VALUE/ F,G,P0,RSIGMA, RJ(200), RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO,EPSI, THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT, PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
      H=0.0
      RSIGMA=0.0
C NONNEGS IF INCLUDED ARE ADDED TO P--ARE P0S. IN ALL PHASES
      GO TO (50,55),NT2
   50 DO 52 I=1,N
   52  RSIGMA=RSIGMA+RHO/X(I)
   55 GO TO (60,100,300),NPHASE
C OBJECTIVE FUNCTION   -SIGMA VIOL. CONSTS.
  60  F=0.0
 100  IF (M)160,160,105
  105  DO 150 J=1,M
      IF (RJ(J))120,120,110
  110  RSIGMA=RSIGMA+RHO/RJ(J)
      GO TO 150
  120  F=F-RJ(J)
  150  CONTINUE
C EQUALITIES NOT ADDED IN FEAS. PHASE
 160  IF (MZ)180,180,163
  163  GO TO (180,165,300),NPHASE
  165  DO 170 I=1,MZ
      K=M+I
 170  H=H+RJ(K)**2
      H=H/SQRT(RHO)
  180  HS=H+RSIGMA
      P0=F+HS
      HMS=-RSIGMA+2.*H
      G=F+HMS
      RETURN
  300  RETURN
      END
      SUBROUTINE SECORD(IS)
      COMMON/SHARE/ X(50), DEL(50), A(50,50),N,M, MN,NP1,NM1
      COMMON/EQAL/ MZ, H, H1, H3, NA, NC,IERR
      COMMON /OPTNS/ NT1,NT2,NT3,NT4,NT5,NT6,NT7,NT8,NT9,NT10,NT11,NT12
      COMMON /VALUE/ F,G,P0,RSIGMA, RJ(200), RHO
      COMMON/CRST/ DELX(50), DELX0(50), RHOIN,RATIO, EPSI, THETA0,
     1NTCTR, NUMINI, X1(50), X2(50), X3(50), XR2(50), XR1(50), PR1,
     2PR2, P1, F1, RJ1(200), DOTT,PGRAD(50), DIAG(50), RJ3(200),
     3PREV3, F3, ADELX, RSIG1, G1, G3, RSIG3,NVI,NPHASE,NSATIS,MMAX,IMAX
C- (1) MEANS DONT  COMPUTE GRAD. OUTER PRODUCT(IN SECORD)
  99  DO 102 I=1,N
      DO 102 J=I,N
  102 A(I,J)=0.
  110 GO TO (135,111),IS
CGRAD .TERM NOT PREV. COMPUTED
  111 GO TO (112,115),NT2
  112 DO 113 I=1,N
  113 A(I,I)=2.*RHO/X(I)**3
  115 CONTINUE
  116 IF (M.LE.0) GO TO 131
      DO 130 K=1,M
      IF (RJ(K))130,130,117
  117	CALL EQUATE(3,K,DUM,KDUM)
      TT=2.*RHO/RJ(K)**3
      DO 125 I=1,N

      IF (DEL(I))120,125,120
  120 T=TT*DEL(I)
      DO 123 J=1,I
      IF (DEL(J))122,123,122
  122 A(I,J)=A(I,J)+T*DEL(J)
  123 CONTINUE
  125 CONTINUE
  130 CONTINUE
C EQUALITY  CONSTRAINTS
  131 IF (MZ)1131,1131,133
  133 GO TO (1131,134,135),NPHASE
  134 RQ=2./SQRT(RHO)
      DO 1350 JJ=1,MZ
      K=M+JJ
	CALL EQUATE(3,K,DUM,KDUM)
      DO 1325 I=1,N
      IF (DEL(I))1320,1325,1320
 1320  T=RQ*DEL(I)
      DO 1323 J=1,I
      IF (DEL(J))1322,1323,1322
 1322 A(I,J)=A(I,J)+T*DEL(J)
 1323 CONTINUE
 1325 CONTINUE
 1350 CONTINUE
 1131 DO 1132 I=1,N
      DIAG(I)=A(I,I)
 1132 A(I,I)=0.
C  READY NOW FOR MATRIX OF 2ND PARTIALS OF RESTRAINTS
  135 GO TO (136,1000,1010),NT10
  136 IF (M.LE.0) GO TO 2151
      DO 150 K=1,M
  137 LORN=2
C CONSTRAINT ASSUMED NONLINEAR
	CALL EQUATE(4,K,DUM,LORN)
      IF (LORN.LT.2) GO TO 150
      IF (RJ(K))1150,1150,138
  138 T=-RHO/RJ(K)**2
      DO 145 I=2,N
      IM1=I-1
      DO 145 J=1,IM1
      IF (A(J,I))143,145,143
  143 A(I,J)=A(I,J)+T*A(J,I)
      A(J,I)=0.
  145 CONTINUE
      DO 151 I=1,N
      IF (A(I,I))147,151,147
  147  DIAG(I)=DIAG(I)+T*A(I,I)
      A(I,I)=0.
  151 CONTINUE
  150 CONTINUE
 2151 CONTINUE
      GO TO (1010,165,1010),NPHASE
C--   EQUALITY  SECOND PARTIALS HERE
C  GET MATRIX OF 2ND PARTIAL OF OBJECTIVE FUNCTION
165	CALL EQUATE(4,0,DUM,0)
      DO 170 I=2,N
      IM1=I-1
      DO 170 J=1,IM1
      IF (A(J,I))167,170,167
  167 A(I,J)=A(I,J)+A(J,I)
  170 A(J,I)=A(I,J)
      DO 174 I=1,N
      IF (A(I,I))172,173,172
  172 A(I,I)=DIAG(I)+A(I,I)
      GO TO 174
  173 A(I,I)=DIAG(I)
  174 CONTINUE
  175 RETURN
 1000 GO TO (1010,165,165),NPHASE
 1010 DO 1020 I=2,N
      IM1=I-1
      DO 1020 J=1,IM1
 1020 A(J,I)=A(I,J)
      DO 1030 I=1,N
 1030 A(I,I)=DIAG(I)
      GO TO 175
 1150 DO 1145 I=2,N
      IM1=I-1
      DO 1145 J=1,IM1
      IF (A(J,I))1143,1145,1143
 1143 A(I,J)=A(I,J)-A(J,I)
      A(J,I)=0.
 1145 CONTINUE
      DO 1151 I=1,N
      DIAG(I)=DIAG(I)-A(I,I)
 1151 A(I,I)=0.0
      GO TO 150
      END

