C C THE MAIN PROGRAM IS A DUMMY ROUTINE THAT SETS UP THE WORK C VECTORS DP AND INT AFTER READING THE DIMENSION N OF THE C PROBLEM. DP AND INT ARE DIMENSIONED FOR PROBLEMS WITH C N .LE. 50 WITH NO SPECIAL SUBROUTINE FNEVL1, N .LE. 40 WITH C AT MOST 50 CONSTRAINTS, OR N .LE. 30 WITH AT MOST 170 CONSTRAINTS. C FOR LARGER N, DP SHOULD BE OF LENGTH (7*N**2 + 29*N) / 2 + 1020 C IF THERE ARE NO CONSTRAINTS, AND OF LENGTH (9*N**2 + 37*N) / 2 C + 1010 + (N + 1) * NUMBER OF CONSTRAINTS OTHERWISE. INT SHOULD C BE OF LENGTH N**2 + 11*N + 1010 + NUMBER OF CONSTRAINTS. C C MAIN THEN PASSES CONTROL TO THE MAIN SUBROUTINE PLALGO. C CHARACTER CHA(72) INTEGER N,NP1,NP12,LENINT,LENDP,LIWA,M,MIN0,INT(1700) INTEGER I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14 INTEGER I15,I16,I17,I18,I19,I20,I21,I22,I23 INTEGER J,J2,J3,J4,J5,J6,J7,J8,J9,J10,J11,J12 DOUBLE PRECISION DP(2710) LOGICAL QADRIV C C BEGIN PROGRAM C OPEN (5,FILE='FI5') OPEN (6,FILE='FI6') LENINT = 1700 LENDP = 2710 READ (5,100) (CHA(J),J = 1,72) WRITE (6,200) (CHA(J),J = 1,72) READ (5,300) N IF (N.LE.0) WRITE (6,400) IF (N.LE.0) GO TO 20 NP1 = N + 1 N2 = N*N LIWA = 2*NP1 NP12 = NP1*NP1 C C CALCULATE MAXIMUM NUMBER M OF CONSTRAINTS C NN = 1 M = 1 IF ((9*N2 + 37*N + 2022) .GT. 2*LENDP) GO TO 10 NN = N M =MIN0((LENINT-N2-11*N-1008),(2*LENDP-9*N2-37*N-2020)/(2*NP1)) 10 CONTINUE NN2 = NN*NN I2 = N2 + 1 I3 = I2 + NP12 I4 = I3 + NP12 I5 = I4 + (N*NP1)/2 I6 = I5 + NP1 I7 = I6 + NP1 I8 = I7 + N I9 = I8 + N I10 = I9 + NP1 I11 = I10 + NP1 I12 = I11 + NP1 I13 = I12 + N I14 = I13 + NP1 I15 = I14 + NP1 I16 = I15 + 1000 I17 = I16 + NN I18 = I17 + NN I19 = I18 + NN I20 = I19 + NN2 I21 = I20 + NN I22 = I21 + 1 I23 = I22 + M J2 = 1 + NP1 J3 = J2 + NP1 J4 = J3 + NP1 J5 = J4 + NP1 J6 = J5 + NP1 J7 = J6 + N*NP1 J8 = J7 + N J9 = J8 + N J10 = J9 + NP1 J11 = J10 +LIWA J12 = J11 + 1000 CALL PLALGO(N,NP1,DP(1),DP(I2),DP(I3),DP(I4),DP(I5),DP(I6),DP(I7), * DP(I8),DP(I9),DP(I10),DP(I11),DP(I12),DP(I13),DP(I14), * DP(I15),DP(I16),DP(I17),DP(I18),DP(I19),DP(I20), * DP(I21),DP(I22),DP(I23), * INT(1),INT(J2),INT(J3),INT(J4),INT(J5),INT(J6), * INT(J7),INT(J8),INT(J9),INT(J10),INT(J11),INT(J12), * LIWA,QADRIV) C C IF TERMINATION ROUTINE REQUIRED REMOVE THE C'S BELOW C C CALL TERMIN(N,NP1,DP(1),DP(I2),DP(I3),DP(I4),DP(I5),DP(I6),DP(I7), C * DP(I8),DP(I9),DP(I10),DP(I11),DP(I12),DP(I13),DP(I14), C * DP(I15),DP(I16),DP(I17),DP(I18),DP(I19),DP(I20), C * DP(I21),DP(I22),DP(I23), C * INT(1),INT(J2),INT(J3),INT(J4),INT(J5),INT(J6), C * INT(J7),INT(J8),INT(J9),INT(J10),INT(J11),INT(J12), C * LIWA,QADRIV) C 20 STOP 100 FORMAT (72A1) 200 FORMAT (6X,72A1) 300 FORMAT (I5) 400 FORMAT (50H THE PROBLEM DIMENSION IS NONPOSITIVE ) END SUBROUTINE PLALGO(N,NP1,ADERIV,BAS,R,FR,WEIGHT,VTWT,XORIG,XNEW, * NEWCOL,UPDCOL,RCOL,DELTA,SIGMA,WA, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA, * LPERM,FLPERM,IPOINT,KPOINT,MOSTV, * VERT,PI,SIGN,DIM,IWA,IDAT,IFACT,LIWA,QADRIV) INTEGER N,NP1,LIWA INTEGER LPERM(NP1),FLPERM(NP1),IPOINT(NP1),KPOINT(NP1),MOSTV(NP1) INTEGER VERT(N,NP1),PI(N),SIGN(N),IWA(LIWA),DIM(NP1) LOGICAL QADRIV DOUBLE PRECISION ADERIV(N,N),BAS(NP1,NP1),R(NP1,NP1),FR(*) DOUBLE PRECISION WEIGHT(NP1),VTWT(NP1),XORIG(N),XNEW(N) DOUBLE PRECISION NEWCOL(NP1),UPDCOL(NP1),RCOL(NP1),DELTA(N) DOUBLE PRECISION SIGMA(NP1),WA(NP1) INTEGER IDAT(1000),IFACT(*) DOUBLE PRECISION XDAT(1000),FX(*),FXORIG(*),FF(*),FM(*) DOUBLE PRECISION FDIAG(*),FDELTA,FB(*),FA(*) C C PLALGO COMPUTES AN APPROXIMATE ZERO OF A FUNCTION FROM R**N C INTO ITSELF (OR OF AN UPPER SEMI-CONTINUOUS MAPPING FROM R**N C INTO NONEMPTY CONVEX COMPACT SUBSETS OF R**N). IT USES A C PIECEWISE LINEAR RESTART ALGORITHM WITH INTERPOSED QUASI-NEWTON C STEPS IF THE FUNCTION IS SMOOTH. C C THE USER MUST PROVIDE A C C SUBROUTINE FNEVAL(N,X,FUN,JOB,INFO,XDAT,IDAT) C C THAT PERFORMS THE FOLLOWING TASKS. C C IF JOB.EQ.1, ANY INITIAL DATA FOR THE COMPUTATION OF THE FUNCTION C SHOULD BE READ INTO THE DOUBLE PRECISION ARRAY XDAT C OF LENGTH 1000 AND THE INTEGER ARRAY IDAT OF THE SAME C LENGTH. C C IF JOB.EQ.2, THE SUBROUTINE SHOULD PERFORM ANY INTERMEDIATE C (3) (OR FINAL) PRINTING DESIRED BEYOND THE C APPROXIMATE ZERO AND ITS FUNCTION VALUE (E.G., C FEASIBILITIES AND OBJECTIVE VALUE IN OPTIMIZATION C PROBLEMS, PRICE VECTOR AND DEMANDS AND PROFITS C IN ECONOMIC EQUILIBRIUM PROBLEMS). C EACH SUCH CALL FOLLOWS A CALL WITH THE SAME X C AND JOB.EQ.0, SO THAT FURTHER COMPUTATION MAY C BE UNNECESSARY. C C IF JOB.EQ.0, THE SUBROUTINE SHOULD RETURN IN THE VECTOR FUN THE C VALUE OF THE FUNCTION F (OR A SELECTION FROM THE C MAPPING F) AT THE POINT X. HERE X AND FUN ARE C DOUBLE PRECISION LINEAR ARRAYS OF LENGTH N. C INFO IS AN INTEGER TO BE EXPLAINED BELOW. C C EACH APPROXIMATE ZERO OBTAINED IS GENERATED BY PUTTING WEIGHTS C ON NEARBY POINTS AT WHICH F HAS BEEN EVALUATED. FOR CERTAIN C APPLICATIONS, ONLY SOME OF THESE POINTS SHOULD BE INCLUDED. C FOR EXAMPLE, IN OPTIMIZATION PROBLEMS WITH CONVEX CONSTRAINT C SETS ONLY FEASIBLE POINTS, AND IN EQUILIBRIUM PROBLEMS ONLY C PRICE VECTORS FOR WHICH NO ACTIVITY IS PROFITABLE, SHOULD BE C INCLUDED. TO USE THIS OPTION, EACH TIME FNEVAL IS CALLED WITH C JOB.EQ.0 INFO SHOULD BE SET TO 0 IF X IS TO BE INCLUDED AND C 1 OTHERWISE C C PLALGO OPERATES BY MOVING THROUGH PIECES OF A SUBDIVISION OF C R**N X [0,1] UNTIL AN APPROXIMATE ZERO IS FOUND. THEN, C AFTER POSSIBLE ACCELERATION, THE PROCESS IS REPEATED C WITH A FINER SUBDIVISION GIVING A BETTER APPROXIMATION, C AND SO ON. A MAJOR CYCLE IS ONE SUCH TRAVERSAL. C C ON ENTRY, NONE OF THE PARAMETERS EXCEPT N AND NP1 (=N+1) NEED C BE DEFINED. ALL PARAMETERS OF THE ALGORITHM ARE SET C IN SUBROUTINE SETUP, WHERE THEY ARE DESCRIBED. C C ON RETURN, XNEW CONTAINS THE APPROXIMATE ZERO FOUND AND (THE C FIRST N ENTRIES OF) NEWCOL CONTAINS ITS FUNCTION C VALUE. IF QADRIV IS .TRUE., ADERIV CONTAINS AN C APPROXIMATE DERIVATIVE MATRIX OF F AT XNEW. C C THE SUBROUTINES CALLED BY PLALGO ARE: C C SETUP -- READS INITIAL PARAMETERS (DEFAULT VALUES ARE C PROVIDED FOR MOST); C C A COLLECTION OF SUBROUTINES CALLED SUBDIV THAT DEAL WITH C THE PIECES OF THE SUBDIVISION; C C A COLLECTION OF SUBROUTINES CALLED BASIS THAT DEAL WITH C THE LINEAR SYSTEMS THAT DETERMINE THE PROGRESSION C OF PIECES; C C A COLLECTION OF SUBROUTINES CALLED LINALG THAT PERFORM C LINEAR ALGEBRA TASKS; C C ENDCYC -- CALCULATES THE APPROXIMATE ZERO AT THE END OF C EACH MAJOR CYCLE; AND C C ACCEL -- TAKES QUASI-NEWTON (BROYDEN) STEPS AND ACCELERATES C THE REFINEMENT OF THE GRID SIZE IF JUDGED WORTHWHILE. C FUNCTIONS CALLED -- DBLE C C THE PROGRAM IS WRITTEN IN 'VANILLA' FORTRAN AND FEW CHANGES C SHOULD BE NECESSARY FOR ANY COMPILER. THE CODE HAS BEEN C VERIFIED TO CONFORM WITH PFORT, A PORTABLE SUBSET OF ANSI C FORTRAN, WITH THE EXCEPTION OF THE FOLLOWING FEATURES. THE C MAIN PROGRAM HAS OPEN STATEMENTS WHICH MAY BE UNNECESSARY OR C UNDESIRABLE IN YOUR ENVIRONMENT. IT ALSO HAS CHARACTER VARIABLES, C WHICH SHOULD BE CHANGED TO INTEGER FOR FORTRAN 66 COMPILERS; C CHARACTER VARIABLES ARE ALSO USED IN SUBROUTINE PRTPCE IN SUBDIV. C IN ADDITION, THE INTRINSIC FUNCTION DBLE IS USED INSTEAD OF C DFLOAT TO CONVERT INTEGER VARIABLES TO DOUBLE PRECISION. THIS C FUNCTION IS USED IN SUBROUTINES PLALGO AND ACCEL, IN SUBROUTINES C NEWPCE AND LSTPCE IN SUBDIV, AND IN SUBROUTINES FSTBAS, NEWBAS, C AND MINRAT IN BASIS. C C **** VERSION OF AUGUST 15, 1983 C INTEGER NEFFP1,PRTLEV,MAXDEP,MAXIT,TYPJAC,JACFAC,TYPFUN INTEGER NTDEP,NTIT,VDLT,NDEP INTEGER NTOTIT,NTOTFN,NCYCIT,NCYCFN,NBOUND,NTRIV INTEGER IIN,KIN,KINSG,IOUT,KOUT,KOUTSG,ITOP,FLAG,INFO,NWTDIM,I,J LOGICAL OUTTOP DOUBLE PRECISION SCALE,STPDWN,TOLFUN,NMFUN,BIG,DBLE C C BEGIN PROGRAM. INITIALIZE PARAMETERS. C FLAG = 0 CALL SETUP(N,NP1,NEFFP1,PRTLEV,ADERIV,DELTA,XORIG,NEWCOL,SCALE, * STPDWN,MAXDEP,MAXIT,TOLFUN,TYPJAC,JACFAC,TYPFUN,NTDEP, * NTIT,VDLT,BIG,NTOTIT,NTOTFN,QADRIV,FLAG,ITOP,MOSTV, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT,IFACT) NBOUND = 0 NTRIV = 0 IF (FLAG.EQ.2) GO TO 80 C DO 10 I = 1,NP1 DIM(I) = 0 10 CONTINUE C DO 40 NDEP = 1,MAXDEP C C START MAJOR CYCLE. INITIALIZE BASIS AND PIECE. C CALL FSTBAS(N,NP1,NEFFP1,PRTLEV,BAS,R,FR,WEIGHT,ADERIV,DELTA, * XORIG,SCALE,LPERM,FLPERM,IPOINT,KPOINT,TYPJAC, * JACFAC,STPDWN,SIGMA,WA,FLAG,VDLT) IF (FLAG.EQ.2) GO TO 80 CALL FSTPCE(N,NP1,PRTLEV,VDLT,VERT,ITOP,PI,SIGN, * DELTA,XORIG,IOUT,KOUT,KOUTSG,VTWT,XNEW, * IIN,KIN,KINSG) NCYCIT = 0 NCYCFN = 0 C C START ITERATION. C C 20 CALL LABEL(N,NP1,NEFFP1,IIN,KIN,XNEW,MOSTV,ADERIV,DELTA, * SCALE,NEWCOL,BIG,TOLFUN,TYPFUN,FLAG,ITOP,VDLT, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT, * IFACT) DIM(ITOP) = DIM(ITOP) + 1 NCYCIT = NCYCIT + 1 NTOTIT = NTOTIT + 1 IF (IIN.GT.0) NCYCFN = NCYCFN + 1 IF (IIN.GT.0) NTOTFN = NTOTFN + 1 IF (FLAG.EQ.1) GO TO 50 IF (FLAG.EQ.2) GO TO 80 CALL NEWBAS(N,NP1,PRTLEV,BAS,R,WEIGHT,ADERIV,DELTA,SCALE, * NEWCOL,UPDCOL,RCOL,LPERM,IPOINT,KPOINT,MOSTV, * IIN,KIN,KINSG,IOUT,KOUT,KOUTSG,OUTTOP, * NBOUND,NTRIV,TYPFUN,SIGMA,WA,IWA,LIWA) IF (OUTTOP) GO TO 30 IF (NTOTIT.GT.MAXIT.OR.NTOTFN.GT.MAXIT) GO TO 80 CALL NEWPCE(N,NP1,PRTLEV,VDLT,VERT,ITOP,PI,SIGN,IPOINT, * DELTA,XORIG,IOUT,KOUT,KOUTSG,VTWT,XNEW, * IIN,KIN,KINSG,MOSTV,FLAG) IF (FLAG.EQ.2) GO TO 80 GO TO 20 C C MAJOR CYCLE COMPLETED. C 30 CALL ENDCYC(N,NP1,NEFFP1,PRTLEV,XNEW,NEWCOL,WEIGHT,VTWT,NDEP, * IPOINT,MOSTV,DELTA,SCALE,STPDWN,VDLT,VERT,ITOP, * PI,SIGN,XORIG,IOUT,KOUT,KOUTSG,IIN,KIN,KINSG, * FLAG,NMFUN,TOLFUN,TYPFUN,TYPJAC, * NCYCIT,NCYCFN,NTOTIT,NTOTFN,MAXDEP, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT,IFACT) IF (FLAG.EQ.1) GO TO 50 IF (PRTLEV.GE.3) CALL PRTPCE(N,NP1,PRTLEV,VDLT,VERT,ITOP, * PI,SIGN,DELTA,XORIG,IOUT,KOUT, * KOUTSG,VTWT,XNEW,IIN,KIN,KINSG) IF (PRTLEV.GE.3) CALL PRTBAS(N,NP1,BAS,R,WEIGHT,ADERIV, * DELTA,SCALE,WA,LPERM,IPOINT, * KPOINT,IIN,KIN,KINSG,IOUT,KOUT, * KOUTSG) CALL ACCEL(N,NP1,NEFFP1,PRTLEV,BAS,R,ADERIV,DELTA,SCALE, * STPDWN,LPERM,IPOINT,KPOINT,NDEP,NTDEP,NTIT,MAXIT, * NTOTFN,TOLFUN,NMFUN,FLAG,TYPFUN,TYPJAC,JACFAC, * QADRIV,VERT,PI,SIGN,ITOP,NCYCIT,XNEW,XORIG,NEWCOL, * UPDCOL,RCOL,WEIGHT,VTWT,SIGMA,WA,WEIGHT,VTWT, * VDLT,MOSTV, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT,IFACT) IF (FLAG.EQ.1) GO TO 50 IF (FLAG.EQ.2) GO TO 80 40 CONTINUE C C SUCCESSFUL TERMINATION. C 50 WRITE (6,1100) NTOTIT,NTOTFN WRITE (6,1300) NBOUND IF (NTRIV.EQ.0) GO TO 60 WRITE (6,1400) NTRIV 60 CALL FNEVL1(N,XNEW,NEWCOL,3,INFO,NP1,NEFFP1,ITOP,VDLT,MOSTV, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT,IFACT) WRITE (6,1500) WRITE (6,1600) (J,XNEW(J),NEWCOL(J),J = 1,N) WRITE (6,1700) WRITE (6,1800) NMFUN IF (PRTLEV.GT.1) CALL PRTPCE(N,NP1,PRTLEV,VDLT,VERT,ITOP, * PI,SIGN,DELTA,XORIG,IOUT,KOUT, * KOUTSG,VTWT,XNEW,IIN,KIN,KINSG) IF (PRTLEV.GT.1)CALL PRTBAS(N,NP1,BAS,R,WEIGHT,ADERIV,DELTA,SCALE, * WA,LPERM,IPOINT,KPOINT,IIN,KIN,KINSG, * IOUT,KOUT,KOUTSG) C IF (PRTLEV .EQ. 0) GO TO 100 WRITE (6,1900) NWTDIM = 0 J = 0 DO 70 I = 1,NP1 IPOINT(I) = J NWTDIM = NWTDIM + J*DIM(I) J = J + 1 70 CONTINUE WRITE (6,2000) (IPOINT(I),I = 1,NP1) WRITE (6,2000) (DIM(I),I = 1,NP1) SCALE = DBLE(NWTDIM)/ DBLE(N*NTOTIT) WRITE (6,2100) SCALE C GO TO 100 C C UNSUCCESSFUL TERMINATION. C 80 WRITE (6,1200) NTOTIT,NTOTFN WRITE (6,1300) NBOUND IF (NTRIV.EQ.0) GO TO 90 WRITE (6,1400) NTRIV 90 CALL FNEVL1(N,XNEW,NEWCOL,2,INFO,NP1,NEFFP1,ITOP,VDLT,MOSTV, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT,IFACT) WRITE (6,1500) WRITE (6,1600) (J,XNEW(J),NEWCOL(J),J = 1,N) CALL PRTPCE(N,NP1,PRTLEV,VDLT,VERT,ITOP,PI,SIGN,DELTA,XORIG, * IOUT,KOUT,KOUTSG,VTWT,XNEW,IIN,KIN,KINSG) CALL PRTBAS(N,NP1,BAS,R,WEIGHT,ADERIV,DELTA,SCALE,WA, * LPERM,IPOINT,KPOINT,IIN,KIN,KINSG,IOUT,KOUT,KOUTSG) 100 RETURN C C FORMATS C 1100 FORMAT(33H0**** SUCCESSFUL TERMINATION WITH ,I5,11H ITERATIONS * ,4H AND,I5,25H FUNCTION EVALUATIONS ) 1200 FORMAT(35H0**** UNSUCCESSFUL TERMINATION WITH ,I5, * 15H ITERATIONS AND ,I5,25H FUNCTION EVALUATIONS ) 1300 FORMAT(6X,33HNUMBER OF BETWEEN BOUND PIVOTS IS ,I5) 1400 FORMAT(6X,33HNUMBER OF TRIVIAL PIVOTS IS ,I5) 1500 FORMAT(53H0 THE FINAL SOLUTION AND ITS FUNCTION VALUE ARE /) 1600 FORMAT(6X,I2,4X,2D26.10) 1700 FORMAT(45H0 THE FUNCTION VALUE SUP-NORM IS ) 1800 FORMAT(6X,D20.10) 1900 FORMAT(6X,35HDIMENSION OF SIMPLICES ENCOUNTERED ) 2000 FORMAT(20I4) 2100 FORMAT (6X,30HTHE AVERAGE DIMENSION RATIO IS ,F10.3) C C LAST CARD OF PLALGO FOLLOWS C END SUBROUTINE SETUP(N,NP1,NEFFP1,PRTLEV,ADERIV,DELTA,XORIG,NEWCOL, * SCALE,STPDWN,MAXDEP,MAXIT,TOLFUN,TYPJAC,JACFAC, * TYPFUN,NTDEP,NTIT,VDLT,BIG,NTOTIT,NTOTFN,QADRIV, * FLAG,ITOP,MOSTV, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT, * IFACT) INTEGER N,NP1,NEFFP1,PRTLEV,MAXDEP,MAXIT,TYPJAC,JACFAC,TYPFUN INTEGER NTDEP,NTIT,VDLT,NTOTIT,NTOTFN,FLAG,ITOP,MOSTV(NP1) LOGICAL QADRIV DOUBLE PRECISION ADERIV(N,N),DELTA(N),XORIG(N),NEWCOL(NP1) DOUBLE PRECISION SCALE,STPDWN DOUBLE PRECISION TOLFUN,BIG INTEGER IDAT(1000),IFACT(*) DOUBLE PRECISION XDAT(1000),FX(*),FXORIG(*),FF(*),FM(*) DOUBLE PRECISION FDIAG(*),FDELTA,FB(*),FA(*) C C SETUP READS THE PARAMETERS DETERMINING THE ALGORITHM AND C INITIALIZES THE COUNTERS. C C PRTLEV DETERMINES THE AMOUNT OF OUTPUT PRODUCED. FOR MINIMUM C OUTPUT SET PRTLEV TO 0. FOR MAXIMUM OUTPUT SET TO 5. C C TYPJAC DETERMINES THE TYPE OF INITIAL LINEAR MAP USED. C IF TYPJAC <=0 THIS IS A FIRST ORDER TAYLOR EXPANSION C USING A FINITE DIFFERENCE APPROXIMATION TO THE JACOBIAN C IF TYPJAC=1 THIS MAP HAS MATRIX A SCALED C IDENTITY, APPROPRIATE FOR OPTIMIZATION PROBLEMS C IN MINIMIZATION FORM (F = THE GRADIENT OF THE OBJECTIVE). C IF TYPJAC=2 THE MATRIX IS A NEGATIVE SCALED IDENTITY, C APPROPRIATE FOR FIXED POINT PROBLEMS C (F(X) = G(X) - X, G(R**N) COMPACT, E.G.). C IF TYPJAC.GE.3 THE MATRIX IS READ AS ADERIV. C C ADERIV IS AN N BY N ARRAY, THE MATRIX OF THE INITIAL LINEAR C MAP. IF TYPJAC.GE.3 THIS IS READ IN, OTHERWISE IT IS SET. C C TYPFUN IS 0 FOR A SET-VALUED MAPPING, 1 FOR A NONSMOOTH FUNCTION, C AND 2 FOR A SMOOTH FUNCTION. C C XORIG IS A LINEAR ARRAY OF LENGTH N, THE INITIAL GUESSED C SOLUTION TO THE PROBLEM AND THE ZERO OF THE INITIAL C LINEAR MAP. C C DELTA IS A LINEAR ARRAY OF LENGTH N GIVING SCALINGS FOR C EACH COORDINATE. DELTA(J) SHOULD BE OF THE ORDER OF C ONE FIFTH THE UNCERTAINTY IN XORIG(J). C IF ONLY ONE DELTA IS SET ALL DELTA(J) ARE SET TO THIS C VALUE. C IF TYPJAC=0, THE DELTA'S ARE SET AUTOMATICALLY AND ANY C VALUE MAY BE READ. C C THE REMAINING PARAMETERS HAVE DEFAULT VALUES IF READ AS 0. C C STPDWN IS THE RATIO BETWEEN SUCCESSIVE GRID SIZES BEFORE C ACCELERATION. DEFAULT VALUE 0.5. C C MAXDEP IS THE MAXIMUM NUMBER OF MAJOR CYCLES TO BE PERFORMED. C DEFAULT VALUE 10. C C MAXIT IS THE MAXIMUM NUMBER OF ITERATIONS AND OF FUNCTION C EVALUATIONS TO BE PERFORMED. DEFAULT VALUE 20*N**2. C C TOLFUN IS THE TOLERANCE FOR THE FUNCTION USED AS A STOPPING C CRITERION. IF SOME X IS FOUND WITH NORM(F(X)).LT.TOLFUN C THEN THE ALGORITHM TERMINATES. DEFAULT VALUE 1.0D-10. C C BIG IS A LIMIT ON THE VECTORS CONSIDERED. IF SOME X WITH C NORM(X).GT.BIG IS FOUND, THE ALGORITHM TERMINATES. C DEFAULT VALUE 1.0D+10. C C NTDEP IS THE NUMBER OF MAJOR CYCLES BEFORE ACCELERATION IS C ATTEMPTED. DEFAULT VALUE 0. C C NTIT IS THE NUMBER OF BROYDEN STEPS TAKEN BEFORE RESTARTING THE C PL ALGORITHM. IF ZERO, SET TO MAXIT. DEFAULT VALUE 0. C C VDLT IS 0 FOR THE NORMAL SUBDIVISION, 1 FOR VAN DER LAAN AND C TALMAN'S. FOR DETAILS SEE SUBDIV. DEFAULT VALUE 0. C C FACTOR IS THE NUMBER USED TO CALCULATE APPROPRIATE COORDINATE C SCALING FACTORS WHEN TYPJAC IS READ AS ZERO. DEFAULT VALUE C ONE HALF. C C NTOTIT IS THE NUMBER OF ITERATIONS TAKEN. C C NTOTFN IS THE NUMBER OF FUNCTION EVALUATIONS PERFORMED. C C SCALE IS THE APPROXIMATE SIZE OF ENTRIES IN THE BASIS MATRIX. C C SUBROUTINES AND FUNCTIONS USED -- FNEVL1, VSNORM C C **** VERSION OF AUGUST 15, 1983 C INTEGER INFO,I,J,NMIN8 DOUBLE PRECISION TEMP,TEMP2,DEL,ZERO,EIGHTH,FACTOR,HALF,ONE,FNORM DOUBLE PRECISION VSNORM,V2NORM,DMAX1 DATA ZERO,EIGHTH,HALF,ONE /0.0D0,0.125D0,0.5D0,1.0D0/ C C BEGIN PROGRAM C FLAG=0 NEFFP1=NP1 READ (5,200) PRTLEV READ (5,200) TYPJAC SCALE=ONE QADRIV=.FALSE. IF (TYPJAC.LE.0) QADRIV=.TRUE. IF (TYPJAC.LT.3) GO TO 30 C C READ STARTING MATRIX ADERIV C TYPJAC=4 JACFAC=0 SCALE=ZERO DO 20 J=1,N READ (5,300) (ADERIV(I,J),I=1,N) TEMP=VSNORM(N,ADERIV(1,J)) IF (TEMP.EQ.ZERO) FLAG=2 IF (TEMP.GT.SCALE) SCALE=TEMP 20 CONTINUE IF (FLAG.EQ.2) GO TO 1200 30 CONTINUE READ (5,200) TYPFUN C C READ STARTING POINT AND COORDINATE SCALINGS. C READ (5,300) (XORIG(J),J=1,N) NMIN8=N IF (N.GT.8) NMIN8=8 READ (5,300) (DELTA(J),J=1,NMIN8) IF (DELTA(1).LE.ZERO) DELTA(1)=ONE IF (N.EQ.1) GO TO 50 IF (DELTA(2).GT.ZERO) GO TO 45 TEMP=DELTA(1) DO 40 J=2,N DELTA(J)=TEMP 40 CONTINUE GO TO 50 45 IF (N.LE.8) GO TO 50 READ (5,300) (DELTA(J),J=9,N) 50 CONTINUE IF (TYPJAC.EQ.4) SCALE=SCALE*VSNORM(N,DELTA) IF (TYPJAC.EQ.4.OR.TYPJAC.EQ.0) GO TO 80 TEMP=ONE IF (TYPJAC.EQ.2) TEMP=-ONE DO 70 J=1,N DO 60 I=1,N ADERIV(I,J)=ZERO 60 CONTINUE ADERIV(J,J)=TEMP/DELTA(J) 70 CONTINUE C C THE REMAINING PARAMETERS HAVE DEFAULT VALUES IF READ AS ZERO C 80 CONTINUE READ (5,300) STPDWN IF (STPDWN.EQ.ZERO) STPDWN=HALF READ (5,200) MAXDEP IF (MAXDEP.EQ.0) MAXDEP=10 READ (5,200) MAXIT IF (MAXIT.EQ.0) MAXIT=20*N*N READ (5,300) TOLFUN IF (TOLFUN.EQ.ZERO) TOLFUN=1.0D-10 READ (5,300) BIG IF (BIG.EQ.ZERO) BIG=1.0D+10 READ (5,200) NTDEP READ (5,200) NTIT IF (NTIT.EQ.0) NTIT=MAXIT READ (5,200) VDLT READ (5,300) FACTOR IF (FACTOR.LE.ZERO) FACTOR=HALF C C SET COUNTERS C NTOTIT=0 NTOTFN=0 C C READ ANY DATA REQUIRED FOR FUNCTION EVALUATION C CALL FNEVL1(N,XORIG,DELTA,1,INFO,NP1,NEFFP1,ITOP,VDLT,MOSTV, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT,IFACT) IF (TYPJAC.GT.0) GO TO 140 C C COMPUTE APPROXIMATE JACOBIAN BY FINITE DIFFERENCES C TEMP=EIGHTH DO 90 I=1,50 IF (ONE+TEMP*TEMP.EQ.ONE) GO TO 100 TEMP=TEMP*EIGHTH 90 CONTINUE 100 CONTINUE CALL FNEVL1(N,XORIG,NEWCOL,0,INFO,NP1,NEFFP1,ITOP,VDLT,MOSTV, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT,IFACT) SCALE=FACTOR*V2NORM(N,NEWCOL) DO 130 J=1,N DEL=TEMP*XORIG(J) DEL=DMAX1(DEL,TEMP) XORIG(J)=XORIG(J)+DEL CALL FNEVL1(N,XORIG,ADERIV(1,J),0,INFO,NP1,NEFFP1,ITOP,VDLT, * MOSTV, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT,IFACT) DO 110 I=1,N ADERIV(I,J)=(ADERIV(I,J)-NEWCOL(I))/DEL 110 CONTINUE XORIG(J)=XORIG(J)-DEL IF (TYPJAC.LT.0) GO TO 120 TEMP2=V2NORM(N,ADERIV(1,J)) IF (TEMP2.EQ.ZERO) FLAG=2 IF (TEMP2.GT.ZERO) DELTA(J)=SCALE/TEMP2 120 CONTINUE 130 CONTINUE IF (FLAG.EQ.2) GO TO 1200 TYPJAC=4 JACFAC=0 NTOTFN=N+1 WRITE (6,1100) FACTOR C C PRINT PARAMETERS C 140 CONTINUE WRITE (6,400) WRITE (6,500) N,PRTLEV,TYPJAC,TYPFUN,MAXDEP,MAXIT,NTDEP,NTIT, * VDLT IF (TYPJAC.LE.2) GO TO 160 WRITE (6,600) DO 150 I=1,N WRITE (6,700) (ADERIV(I,J),J=1,N) 150 CONTINUE 160 CONTINUE WRITE (6,800) WRITE (6,700) (XORIG(J),J=1,N) WRITE (6,900) WRITE (6,700) (DELTA(J),J=1,N) WRITE (6,1000) STPDWN,TOLFUN,BIG C C FORMATS C 200 FORMAT (I5) 300 FORMAT (8D10.3) 400 FORMAT (20H0 PARAMETERS ARE /7X,1HN,2X,6HPRTLEV,2X,6HTYPJAC, * 2X,6HTYPFUN,3X,6HMAXDEP,2X,5HMAXIT,3X,5HNTDEP,4X,4HNTIT,4X, * 4HVDLT,3X) 500 FORMAT (9I8) 600 FORMAT (35H0 THE INITIAL MATRIX IS /) 700 FORMAT (6X,7D10.3) 800 FORMAT (35H0 THE STARTING POINT IS /) 900 FORMAT (45H0 THE COORDINATE SCALING FACTORS ARE /) 1000 FORMAT (14H0 STPDWN= ,D10.3,10H TOLFUN= , * D10.3,6H BIG= ,D10.3) 1100 FORMAT (6X,43HFACTOR FOR INITIAL COORDINATE SCALINGS = ,D10.3) 1200 RETURN C C LAST CARD OF SETUP FOLLOWS C END SUBROUTINE ACCEL(N,NP1,NEFFP1,PRTLEV,BAS,R,ADERIV,DELTA,SCALE, * STPDWN,LPERM,IPOINT,KPOINT,NDEP,NTDEP,NTIT,MAXIT, * NTOTFN,TOLFUN,NMFUN,FLAG,TYPFUN,TYPJAC,JACFAC, * QADRIV,VERT,PI,SIGN,ITOP,NCYCIT,XOLD,XORIG,FOLD, * XNEW,FNEW,WEIGHT,VTWT,S,PS,U,SCALPS,VDLT,MOSTV, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT, * IFACT) INTEGER N,NP1,NEFFP1,PRTLEV,LPERM(NP1),IPOINT(NP1),KPOINT(NP1) INTEGER NDEP,NTDEP,NTIT,MAXIT,NTOTFN,FLAG,TYPFUN,TYPJAC,JACFAC INTEGER VERT(N,NP1),PI(N),SIGN(N),ITOP,NCYCIT,VDLT,MOSTV(NP1) LOGICAL QADRIV DOUBLE PRECISION BAS(NP1,NP1),R(NP1,NP1),ADERIV(N,N),DELTA(N) DOUBLE PRECISION SCALE,TOLFUN,NMFUN,XOLD(N),XORIG(N),FOLD(NP1) DOUBLE PRECISION XNEW(NP1),FNEW(NP1),WEIGHT(NP1),VTWT(NP1),STPDWN DOUBLE PRECISION S(NP1),PS(NP1),U(NP1),SCALPS(NP1) INTEGER IDAT(1000),IFACT(*) DOUBLE PRECISION XDAT(1000),FX(*),FXORIG(*),FF(*),FM(*) DOUBLE PRECISION FDIAG(*),FDELTA,FB(*),FA(*) C C ACCEL DETERMINES WHETHER ACCELERATION CAN BE PERFORMED. C C ACCELERATION CONSISTS OF A SEQUENCE OF BROYDEN SECANT C STEPS, STARTING WITH THE CURRENT APPROXIMATE ZERO AND C AN APPROXIMATE DERIVATIVE MATRIX OBTAINED FROM THE FINAL C BASIS. IF THE FUNCTION AND THE DEPTH ARE SUITABLE AN C INITIAL STEP CAN BE TAKEN WITHOUT OBTAINING THIS APPROXIMATE C DERIVATIVE. THIS FIRST STEP IS USED TO DETERMINE WHETHER C SUCH STEPS SEEM WORTHWHILE. BROYDEN STEPS ARE REJECTED IF C EITHER THE FIRST STEP IS LONG COMPARED TO THE SIMPLEX SIZE C OR THE NEW FUNCTION VALUE DOES NOT DECREASE SUFFICIENTLY C COMPARED TO THE CURRENT VALUE. IF THIS TEST IS PASSED, A C MAXIMUM OF NTIT SUCH STEPS ARE TAKEN, PROVIDED THE STEPS C DECREASE IN SIZE, THE FUNCTION VALUES DECREASE IN SIZE, AND C THE APPROXIMATE DERIVATIVE MATRIX DOES NOT CHANGE THE SIGN C OF ITS DETERMINANT. C C FUNCTIONS AND SUBROUTINES USED -- DABS, DBLE ,V2NORM,VSNORM, C DOTPRD,SOLVE,LSTPCE,FNEVL1, C QRHESS,QRGEN,QRUPDT. C C **** VERSION OF AUGUST 15, 1983 C INTEGER ITEST,IABS,INFO,I,IJ,IT,J,JP1,K,L,LJ,LL,M,NEFF,N2 INTEGER NACTP1,NTITM1 DOUBLE PRECISION SIG,TEMP,NMFNEW,NMSTEP,ZERO,ONE,P05,P5 DOUBLE PRECISION STS,DABS, DBLE ,V2NORM,VSNORM,DOTPRD DATA ZERO,ONE,P05,P5 /0.0D0,1.0D0,.05D0,0.5D0/ C C BEGIN PROGRAM. BRANCH IF NO ACCELERATION IS TO BE TRIED. C IT=0 IF (TYPFUN.LT.2) GO TO 500 IF (NDEP.LT.NTDEP) GO TO 500 IF (ITOP.NE.NEFFP1) GO TO 500 FOLD(NP1)=-SCALE CALL SOLVE(BAS,NP1,NP1,LPERM,R,NP1,FOLD,PS,U,SCALPS,.FALSE.) DO 10 J=1,NP1 IJ=IPOINT(J) IF (IJ.EQ.0) GO TO 10 VTWT(IJ)=-PS(J) 10 CONTINUE C C VTWT GIVES THE WEIGHTS FOR THE FIRST BROYDEN STEP TO XNEW. C CALL LSTPCE(N,NP1,0,VDLT,VERT,ITOP,PI,SIGN,DELTA,XORIG,0,0,0, * VTWT,XNEW,0,0,0) C C CALCULATE THE STEP S (SCALED BY 1/DELTA) AND ITS NORM NMSTEP C DO 20 J=1,N S(J)=(XNEW(J)-XOLD(J))/DELTA(J) 20 CONTINUE NMSTEP=VSNORM(N,S) C C IF THE STEP IS TOO LONG OR (DUE TO ROUNDOFF) TOO SHORT, C BRANCH OUT. C IF (NMSTEP.GT.STPDWN) GO TO 500 IF (NMSTEP.EQ.ZERO) GO TO 500 CALL FNEVL1(N,XNEW,FNEW,0,INFO,NP1,NEFFP1,ITOP,VDLT,MOSTV, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT,IFACT) NTOTFN=NTOTFN+1 NMFNEW=VSNORM(N,FNEW) C C BRANCH OUT IF TERMINATION CRITERIA MET OR IF FUNCTION C REDUCTION INSUFFICIENT. C IF (NMFNEW.GE.NMFUN) GO TO 500 FLAG=0 IF (NMFNEW.LT.TOLFUN) FLAG=1 IF (NTOTFN.GT.MAXIT) FLAG=2 IF (FLAG.GT.0) GO TO 600 IF (NMFNEW.GE.P5*NMFUN) GO TO 400 IT=1 C C STEP WAS SUCCESSFUL. UPDATE APPROXIMATE DERIVATIVE MATRIX. C NEFF=NEFFP1-1 IF (NEFF.EQ.0) GO TO 400 QADRIV=.TRUE. C C DETERMINE WHETHER GOOD FACTORIZATION OF ADERIV IS EASILY AVAILABLE C IF ((NCYCIT.GT.2*NP1).OR.(NEFF.LT.N)) GO TO 100 IJ=0 IF (IPOINT(1).EQ.1) IJ=1 IF (IPOINT(1).EQ.NEFFP1) IJ=-1 IF (IJ.EQ.0) GO TO 100 IF ((IJ.LT.0).AND.(VDLT.LT.0)) GO TO 100 I=2 IF (IJ.LT.0) I=NEFF DO 30 J=2,NEFFP1 IF (IPOINT(J).NE.I) GO TO 100 I=I+IJ 30 CONTINUE I=1 IF (IJ.LT.0) I=NEFF DO 50 J=1,NEFF JP1=J+1 SIG=ONE IF (IJ.LT.0) SIG=-ONE K=PI(I) L=LPERM(J) LL=LPERM(JP1) IF (SIGN(K).LT.0) SIG=-SIG DO 40 M=1,NEFF ADERIV(M,K)=(BAS(M,LL)-BAS(M,L))*SIG R(M,J)=(R(M,JP1)-R(M,J))*SIG 40 CONTINUE KPOINT(J)=K I=I+IJ 50 CONTINUE R(NEFFP1,NEFF)=R(NEFFP1,NEFFP1)*SIG IF (VDLT.GE.0) GO TO 70 K=PI(1) L=PI(2) SIG=ONE IF (SIGN(L).LT.0) SIG=-ONE DO 60 M=1,NEFF ADERIV(M,K)=P5*ADERIV(M,K) R(M,1)=P5*R(M,1) ADERIV(M,L)=ADERIV(M,L)+SIG*ADERIV(M,K) R(M,2)=R(M,2)+SIG*R(M,1) 60 CONTINUE 70 CONTINUE CALL QRHESS(R,NP1,NEFFP1,1) DO 80 J=1,NEFF LPERM(J)=KPOINT(J) 80 CONTINUE GO TO 200 100 CONTINUE DO 120 J=1,NEFF DO 110 I=1,NEFF ADERIV(I,J)=ZERO 110 CONTINUE 120 CONTINUE NACTP1=NP1-NEFF DO 160 J=NACTP1,NP1 I=IPOINT(J) L=LPERM(J) IF (I.GE.NEFFP1) GO TO 140 K=PI(I) SIG= DBLE(SIGN(K)) IF ((I.LE.2).AND.(VDLT.LT.0)) SIG=P5*SIG DO 130 M=1,NEFF ADERIV(M,K)=ADERIV(M,K)-SIG*BAS(M,L) 130 CONTINUE IF ((I.GT.1).OR.(VDLT.GE.0)) GO TO 140 K=PI(2) SIG=P5* DBLE(SIGN(K)) DO 135 M=1,NEFF ADERIV(M,K)=ADERIV(M,K)-SIG*BAS(M,L) 135 CONTINUE 140 CONTINUE IF (I.EQ.1) GO TO 160 K=PI(I-1) SIG= DBLE(SIGN(K)) IF ((I.EQ.2).AND.(VDLT.LT.0)) SIG=P5*SIG DO 150 M=1,NEFF ADERIV(M,K)=ADERIV(M,K)+SIG*BAS(M,L) 150 CONTINUE 160 CONTINUE DO 180 J=1,N DO 170 I=1,N R(I,J)=ADERIV(I,J) 170 CONTINUE 180 CONTINUE DO 190 J=1,N LPERM(J)=J 190 CONTINUE CALL QRGEN(NEFF,NEFF,R,NP1,.TRUE.,LPERM,WEIGHT,PS) 200 TYPJAC=4 JACFAC=3 C C ADERIV IS NOW AN APPROXIMATE DERIVATIVE MATRIX WITH ITS C COLUMNS SCALED BY DELTA. PERFORM BROYDEN STEPS. C DO 210 J=1,NEFF XOLD(J)=XNEW(J) FOLD(J)=FNEW(J) 210 CONTINUE NMFUN=NMFNEW IF (NTIT.EQ.1) GO TO 500 C C CALCULATE THE PERMUTED STEP PS. C DO 220 J=1,NEFF LJ=LPERM(J) PS(J)=S(LJ) 220 CONTINUE FLAG=1 NTITM1=NTIT-1 DO 300 IT=1,NTITM1 IF (PRTLEV.LT.4) GO TO 230 WRITE (6,800) IT WRITE (6,900) (J,XOLD(J),FOLD(J),J=1,N) 230 STS=V2NORM(NEFF,PS)**2 C C CALCULATE THE SCALED PERMUTED STEP SCALPS. C DO 240 J=1,NEFF SCALPS(J)=PS(J)/STS 240 CONTINUE CALL SOLVE(ADERIV,N,NEFF,LPERM,R,NP1,FOLD,PS,U,S,.FALSE.) TEMP=ONE+DOTPRD(NEFF,PS,SCALPS) C C IF UPDATE MAKES MATRIX NEARLY SINGULAR, BRANCH OUT. C IF (TEMP.LT.P05) GO TO 500 C C UPDATE SCALED APPROXIMATE DERIVATIVE MATRIX. C DO 260 J=1,NEFF LJ=LPERM(J) DO 250 I=1,NEFF ADERIV(I,LJ)=ADERIV(I,LJ)+FNEW(I)*SCALPS(J) 250 CONTINUE 260 CONTINUE C C UPDATE UPPER TRIANGULAR FACTOR R. C CALL QRUPDT(NP1,NEFF,R,U,SCALPS) C C GET NEW PERMUTED STEP PS. C DO 270 J=1,NEFF PS(J)=-PS(J)/TEMP 270 CONTINUE C C GET NEW POINT XNEW. C ITEST=1 DO 280 J=1,NEFF L=LPERM(J) XNEW(L)=XOLD(L)+DELTA(L)*PS(J) IF (XNEW(L).NE.XOLD(L)) ITEST=0 280 CONTINUE IF (ITEST.EQ.1) GO TO 500 TEMP=VSNORM(NEFF,PS) IF (TEMP.GT.NMSTEP) GO TO 500 CALL FNEVL1(N,XNEW,FNEW,0,INFO,NP1,NEFFP1,ITOP,VDLT,MOSTV, * XDAT,FX,FXORIG,FF,FM,FDIAG,FDELTA,FB,FA,IDAT,IFACT) NTOTFN=NTOTFN+1 NMFNEW=VSNORM(N,FNEW) IF (NMFNEW.LT.TOLFUN) GO TO 600 IF (NMFNEW.GT.NMFUN) GO TO 500 IF (NMFNEW.GT.P5*NMFUN) GO TO 400 NMSTEP=TEMP FLAG=2 IF (NTOTFN.GT.MAXIT) GO TO 600 FLAG=1 DO 290 J=1,NEFF XOLD(J)=XNEW(J) FOLD(J)=FNEW(J) 290 CONTINUE NMFUN=NMFNEW 300 CONTINUE IT=NTIT GO TO 500 C C IF NEW POINT BETTER THAN OLD, UPDATE. C 400 CONTINUE DO 410 J=1,N XOLD(J)=XNEW(J) FOLD(J)=FNEW(J) 410 CONTINUE NMFUN=NMFNEW C C END OF ACCELERATION. DETERMINE WHETHER STEPS TAKEN. C 500 CONTINUE FLAG=0 IF (IT.GT.0) GO TO 550 C C NO STEPS TAKEN. C ITEST=0 DO 520 J=1,N DELTA(J)=DELTA(J)*STPDWN IF (ITEST.EQ.1) GO TO 510 IF (DABS(XOLD(J))+DELTA(J).NE.DABS(XOLD(J))) GO TO 510 WRITE (6,700) FLAG=2 GO TO 600 510 CONTINUE ITEST=1 XORIG(J)=XOLD(J) 520 CONTINUE IF (TYPFUN.GT.0) SCALE=SCALE*STPDWN IF (TYPJAC.EQ.4) TYPJAC=3 RETURN C C SOME STEPS TAKEN. C 550 CONTINUE N2=N*NEFF SCALE=VSNORM(N2,ADERIV)*NMSTEP C C RESCALE THE APPROXIMATE DERIVATIVE MATRIX AND ITS FACTOR. C DO 570 J=1,NEFF L=LPERM(J) TEMP=DELTA(L) DO 560 I=1,NEFF ADERIV(I,L)=ADERIV(I,L)/TEMP R(I,J)=R(I,J)/TEMP 560 CONTINUE 570 CONTINUE DO 580 J=1,N DELTA(J)=DELTA(J)*NMSTEP XORIG(J)=XOLD(J) 580 CONTINUE IF (PRTLEV.LE.0) RETURN WRITE (6,800) IT WRITE (6,900) (J,XOLD(J),FOLD(J),J=1,N) WRITE (6,1000) NMFUN WRITE (6,1100) NMSTEP RETURN C C TERMINATION. C 600 CONTINUE IF (FLAG.EQ.2) RETURN IT=IT+1 DO 610 J=1,N XOLD(J)=XNEW(J) FOLD(J)=FNEW(J) 610 CONTINUE NMFUN=NMFNEW WRITE (6,1200) IT IF (IT.EQ.1) RETURN C C RESCALE THE APPROXIMATE DERIVATIVE MATRIX AND ITS FACTOR. C DO 630 J=1,N L=LPERM(J) TEMP=DELTA(L) DO 620 I=1,N ADERIV(I,L)=ADERIV(I,L)/TEMP R(I,J)=R(I,J)/TEMP 620 CONTINUE 630 CONTINUE RETURN C C FORMATS C 700 FORMAT(50H0 STEEP FUNCTION -- STEP SIZE IS TOO SMALL ) 800 FORMAT(13H0 AFTER ,I3,34H BROYDEN STEP(S) THE SOLUTION AND * ,25H ITS FUNCTION VALUE ARE /) 900 FORMAT(6X,I2,5X,2D20.10) 1000 FORMAT(35H0 THE FUNCTION VALUE NORM IS ,D10.3) 1100 FORMAT(25H0 THE GRID RATIO IS ,D10.3) 1200 FORMAT (40H0 NUMBER OF BROYDEN STEPS TAKEN ,I3) C C LAST CARD OF ACCEL FOLLOWS. C END