Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/nonlpr/nonlps.for
There are 2 other files named nonlps.for in the archive. Click here to see a list.
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