* SUBROUTINE PA0GS1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE: * NUMERICAL COMPUTATION OF THE GRADIENT OF THE MODEL FUNCTION. * * PARAMETERS : * II N NUMBER OF VARIABLES. * II KA INDEF OF THE APPROXIMATED FUNCTION. * RI X(N) VECTOR OF VARIABLES. * RO GA(N) GRADIENT OF THE APPROXIMATED FUNCTION. * RI FA VALUE OF THE APPROXIMATED FUNCTION. * RI ETA1 PRECISION OF THE COMPUTED VALUES. * IU NAV NUMBER OF APPROXIMATED FUNCTION EVALUATIONS. * SUBROUTINE PA0GS1(N,KA,X,GA,FA,ETA1,NAV) INTEGER N,KA,NAV DOUBLE PRECISION X(*),GA(*),FA,ETA1 DOUBLE PRECISION XSTEP,XTEMP,FTEMP,ETA INTEGER IVAR ETA=SQRT(ETA1) FTEMP=FA DO 4 IVAR=1,N C C STEP SELECTION C XSTEP=1.0D 0 XSTEP=ETA*MAX(ABS(X(IVAR)),XSTEP)*SIGN(1.0D 0,X(IVAR)) XTEMP=X(IVAR) X(IVAR)=X(IVAR)+XSTEP XSTEP=X(IVAR)-XTEMP NAV=NAV+1 CALL FUN(N,KA,X,FA) C C NUMERICAL DIFFERENTIATION C GA(IVAR)=(FA-FTEMP)/XSTEP X(IVAR)=XTEMP 4 CONTINUE FA=FTEMP RETURN END * SUBROUTINE PA1SQ1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION * WHICH IS DEFINED AS A SUM OF SQUARES. * * PARAMETERS: * II N NUMBER OF VARIABLES. * RI X(N) VECTOR OF VARIABLES. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RO AF(N) VALUES OF THE APPROXIMATED FUNCTIONS. * RI GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. * RI AG(N*N) RECTANGULAR MATRIX WHICH IS USED FOR THE DIRECTION * VECTOR DETERMINATION. * RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RI ETA1 PRECISION OF THE COMPUTES FUNCTION VALUES. * II KDA DEGREE OF COMPUTED DERIVATIVES. * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * IU NFV NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED. * IU NFG NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED. * * SUBPROGRAMS USED : * S MXVCOP COPYING OF A VECTOR. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE PA1SQ1(N,X,F,AF,GA,AG,G,ETA1,KDA,KD,LD,NFV,NFG) INTEGER N,KDA,KD,LD,NFV,NFG DOUBLE PRECISION X(*),F,AF(*),GA(*),AG(*),G(*),ETA1 INTEGER KA,NAV DOUBLE PRECISION FA IF (KD.LE.LD) RETURN IF (KD.GE.0.AND.LD.LT.0) THEN F=0.0D0 NFV=NFV+1 ENDIF IF (KD.GE.1.AND.LD.LT.1) THEN CALL MXVSET(N,0.0D0,G) IF (KDA.GT.0) NFG=NFG+1 ENDIF NAV=0 DO 30 KA=1,N IF (KD.LT.0) GO TO 30 IF (LD.GE.0) THEN FA = AF(KA) GO TO 10 ELSE CALL FUN(N,KA,X,FA) AF(KA) = FA ENDIF IF (LD.GE.0) GO TO 10 F=F+FA*FA 10 IF (KD.LT.1) GO TO 30 IF (KDA.GT.0) THEN CALL DFUN(N,KA,X,GA) ELSE CALL PA0GS1(N,KA,X,GA,FA,ETA1,NAV) ENDIF CALL MXVDIR(N,FA,GA,G,G) CALL MXVCOP(N,GA,AG((KA-1)*N+1)) 30 CONTINUE NFV=NFV+NAV/N IF (KD.GE.0.AND.LD.LT.0) F=0.5D0*F LD=KD RETURN END * SUBROUTINE PA2SQ1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF THE VALUE AND THE GRADIENT AND THE HESSIAN MATRIX * OF THE OBJECTIVE FUNCTION WHICH IS DEFINED AS A SUM OF SQUARES. * * PARAMETERS: * II NF NUMBER OF VARIABLES. * II NA NUMBER OF APPROXIMATED FUNCTIONS. * RI GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. * RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RI HA(NF*(NF+1)/2) HESSIAN MATRIX OF THE APPROXIMATED FUNCTION. * RO H(NF*(NF+1)/2) HESSIAN MATRIX OF THE OBJECTIVE FUNCTION. * RI FA VALUE OF THE SELECTED FUNCTION. * RO F VALUE OF THE OBJECTIVE FUNCTION. * * COMMON DATA : * II KDA DEGREE OF ANALYTICALLY COMPUTED DERIVATIVES. * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * IU NAV NUMBER OF APPROXIMATED FUNCTION VALUES COMPUTED. * IU NAG NUMBER OF APPROXIMATED FUNCTION GRADIENTS COMPUTED. * IU NAH NUMBER OF APPROXIMATED FUNCTION HESSIAN MATRICES COMPUTED. * IU NFV NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED. * IU NFG NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED. * IU NFH NUMBER OF OBJECTIVE FUNCTION HESSIAN MATRICES COMPUTED. * IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * * STATUS VARIABLES : * NS,ISP,TSS * * SUBPROGRAMS USED : * S UYPRO1 PROLOGUE. * S UYEPI1 EPILOGUE. * S UYSET0 STATUS DEFINITION. * S MXVSET INITIATION OF A VECTOR. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * S MXDSMO INITIATION OF A DENSE SYMMETRIC MATRIX. * S MXDSMA DENSE SYMMETRIC MATRIX AUGMENTED BY THE SCALED DENSE * SYMMETRIC MATRIX. * S MXDSMU CORRECTION OF A DENSE SYMMETRIC MATRIX. * SUBROUTINE PA2SQ1(NF,NA,X,F,AF,GA,G,H,ETA1,KDA,KD,LD,NFV,NFG) INTEGER NF,NA,KDA,KD,LD,NFV,NFG DOUBLE PRECISION X(*),F,AF(*),GA(*),G(*),H(*),ETA1 INTEGER KA,NAV DOUBLE PRECISION FA IF (KD.LE.LD) RETURN IF (KD.GE.0.AND.LD.LT.0) THEN F=0.0D0 NFV=NFV+1 ENDIF IF (KD.GE.1.AND.LD.LT.1) THEN CALL MXVSET(NF,0.0D0,G) IF (KDA.GT.0) NFG=NFG+1 ENDIF IF (KD.GE.2.AND.LD.LT.2) CALL MXVSET(NF*(NF+1)/2,0.0D0,H) NAV=0 DO 30 KA=1,NA IF (KD.LT.0) GO TO 30 IF (LD.GE.0) THEN FA = AF(KA) GO TO 10 ELSE CALL FUN(NF,KA,X,FA) AF(KA) = FA ENDIF IF (LD.GE.0) GO TO 10 F=F+FA*FA 10 IF (KD.LT.1) GO TO 30 IF (KDA.GT.0) THEN CALL DFUN(NF,KA,X,GA) ELSE CALL PA0GS1(NF,KA,X,GA,FA,ETA1,NAV) ENDIF CALL MXVDIR(NF,FA,GA,G,G) IF (KD.LT.2) GO TO 30 CALL MXDSMU(NF,H,1.0D0,GA) 30 CONTINUE NFV=NFV+NAV/NA IF (KD.GE.0.AND.LD.LT.0) F=0.5D0*F LD=KD RETURN END * SUBROUTINE PC1F01 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF THE VALUE AND THE GRADIENT OF THE CONSTRAINT FUNCTION. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * RI X(NF) VECTOR OF VARIABLES. * RI FC VALUE OF THE SELECTED CONSTRAINT FUNCTION. * RI CF(NC) VECTOR CONTAINING VALUES OF CONSTRAINT FUNCTIONS. * RI CL(NC) VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CU(NC) VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * RI GC(NF) GRADIENT OF THE SELECTED CONSTRAINT FUNCTION. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE GRADIENTS OF CONSTRAINT * FUNCTIONS. * RO CMAX MAXIMUM CONSTRAINT VIOLATION. * II KD DEGREE OF REQUIRED DERVATIVES. * II LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * * SUBPROGRAMS USED : * S MXVCOP COPYING OF A VECTOR. * SUBROUTINE PC1F01(NF,NC,X,FC,CF,CL,CU,IC,GC,CG,CMAX,KD,LD) DOUBLE PRECISION FC,CMAX INTEGER KD,LD,NC,NF DOUBLE PRECISION CF(*),CG(*),CL(*),CU(*),GC(*),X(*) INTEGER IC(*) DOUBLE PRECISION POM,TEMP INTEGER KC IF (KD.LE.LD) RETURN IF (LD.LT.0) CMAX = 0.0D0 DO 20 KC = 1,NC IF (KD.LT.0) GO TO 20 IF (LD.GE.0) THEN FC = CF(KC) GO TO 10 ELSE CALL CON(NF,KC,X,FC) CF(KC) = FC END IF IF (IC(KC).GT.0) THEN POM = 0.0D0 TEMP = CF(KC) IF (IC(KC).EQ.1.OR.IC(KC).GE.3) POM = MIN(POM, + TEMP - CL(KC)) IF (IC(KC).EQ.2.OR.IC(KC).GE.3) POM = MIN(POM, + CU(KC) - TEMP) IF (POM.LT.0.0D0) THEN CMAX = MAX(CMAX,-POM) END IF END IF 10 IF (KD.LT.1) GO TO 20 IF (LD.GE.1) THEN CALL MXVCOP(NF,CG((KC - 1) * NF + 1),GC) GO TO 20 ELSE CALL DCON(NF,KC,X,GC) CALL MXVCOP(NF,GC,CG((KC - 1) * NF + 1)) END IF 20 CONTINUE LD = KD RETURN END * SUBROUTINE PF1F01 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION. * * PARAMETERS: * II NF NUMBER OF VARIABLES. * RI X(NF) VECTOR OF VARIABLES. * RI GF(NF) GRADIENT OF THE MODEL FUNCTION. * RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RI FF VALUE OF THE MODEL FUNCTION. * RO F VALUE OF THE OBJECTIVE FUNCTION. * II KD DEGREE OF REQUIRED DERIVATIVES. * II LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II IEXT TYPE OF EXTREMUM. IEXT=0-MINIMUM. IEXT=1-MAXIMUM. * * SUBPROGRAMS USED : * S MXVCOP COPYING OF A VECTOR. * S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. * SUBROUTINE PF1F01(NF,X,GF,G,FF,F,KD,LD,IEXT) DOUBLE PRECISION F,FF INTEGER IEXT,KD,LD,NF DOUBLE PRECISION GF(*),G(*),X(*) INTEGER NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES COMMON /STAT/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH IF (KD.LE.LD) RETURN IF (LD.GE.0) GO TO 10 NFV = NFV + 1 CALL OBJ(NF,X,FF) IF (IEXT.LE.0) THEN F = FF ELSE F = -FF END IF 10 IF (KD.LT.1) GO TO 20 IF (LD.GE.1) GO TO 20 NFG = NFG + 1 CALL DOBJ(NF,X,GF) IF (IEXT.GT.0) THEN CALL MXVNEG(NF,GF,G) END IF 20 LD = KD RETURN END * SUBROUTINE PLADB0 ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * NEW LINEAR CONSTRAINT OR A NEW SIMPLE BOUND IS ADDED TO THE * ACTIVE SET. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * IU N ACTUAL NUMBER OF VARIABLES. * IU ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RU CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RU CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. * RA S(NF) AUXILIARY VECTOR. * RI EPS7 TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS. * RO GMAX MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE. * RO UMAX MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER. * II INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * IU NADD NUMBER OF CONSTRAINT ADDITIONS. * IO IER ERROR INDICATOR. * * SUBPROGRAMS USED : * S PLADR0 CORRECTION OF KERNEL OF THE ORTHOGONAL PROJECTION * AFTER CONSTRAINT ADDITION. * S MXDRMM PREMULTIPLICATION OF A VECTOR BY A ROWWISE STORED DENSE * RECTANGULAR MATRIX. * S MXDRMV COPY OF THE SELECTED COLUMN OF A ROWWISE STORED DENSE * RECTANGULAR MATRIX. * S MXDRGR PLANE ROTATION OF A TRANSPOSED DENSE RECTANGULAR MATRIX. * S MXVORT DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR * PLANE ROTATION. * SUBROUTINE PLADB0(NF,N,ICA,CG,CR,CZ,S,EPS7,GMAX,UMAX,INEW,NADD, & IER) INTEGER NF,N,ICA(*),INEW,NADD,IER DOUBLE PRECISION CG(*),CR(*),CZ(*),S(*),EPS7,GMAX,UMAX DOUBLE PRECISION CK,CL INTEGER K,L,N1 CALL PLADR0(NF,N,ICA,CG,CR,S,EPS7,GMAX,UMAX,INEW,NADD,IER) IF (IER.NE.0) RETURN IF (N.GT.0) THEN N1=N+1 IF (INEW.GT.0) THEN CALL MXDRMM(NF,N1,CZ,CG((INEW-1)*NF+1),S) ELSE CALL MXDRMV(NF,N1,CZ,S,-INEW) ENDIF DO 1 L=1,N K=L+1 CALL MXVORT(S(K),S(L),CK,CL,IER) CALL MXDRGR(NF,CZ,K,L,CK,CL,IER) IF (IER.LT.0) RETURN 1 CONTINUE ENDIF IER=0 RETURN END * SUBROUTINE PLADB4 ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * NEW LINEAR CONSTRAINT OR A NEW SIMPLE BOUND IS ADDED TO THE ACTIVE * SET. TRANSFORMED HESSIAN MATRIX APPROXIMATION OR ITS INVERSION * IS UPDATED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * IU N ACTUAL NUMBER OF VARIABLES. * IU ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RU CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RU CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. * RU H(NF*(NF+1)/2) TRANSFORMED HESSIAN MATRIX APPROXIMATION OR * ITS INVERSION. * RA S(NF) AUXILIARY VECTOR. * RI EPS7 TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS. * RO GMAX MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE. * RO UMAX MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER. * IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION. * IDECF=10-DIAGONAL MATRIX. * II INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * IU NADD NUMBER OF CONSTRAINT ADDITIONS. * IO IER ERROR INDICATOR. * * SUBPROGRAMS USED : * S PLADR0 CORRECTION OF KERNEL OF THE ORTHOGONAL PROJECTION * AFTER CONSTRAINT ADDITION. * S MXDRMM PREMULTIPLICATION OF A VECTOR BY A ROWWISE STORED DENSE * RECTANGULAR MATRIX. * S MXDRMV COPY OF THE SELECTED COLUMN OF A ROWWISE STORED DENSE * RECTANGULAR MATRIX. * S MXDRGR PLANE ROTATION OF A TRANSPOSED DENSE RECTANGULAR MATRIX. * RECTANGULAR MATRIX. * S MXDSMR PLANE ROTATION OF A DENSE SYMMETRIC MATRIX. * S MXVORT DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR * PLANE ROTATION. * SUBROUTINE PLADB4(NF,N,ICA,CG,CR,CZ,H,S,EPS7,GMAX,UMAX,IDECF, & INEW,NADD,IER) INTEGER NF,N,ICA(*),IDECF,INEW,NADD,IER DOUBLE PRECISION CG(*),CR(*),CZ(*),H(*),S(*),EPS7,GMAX,UMAX DOUBLE PRECISION CK,CL INTEGER I,J,K,L,N1 IF (IDECF.NE.0.AND.IDECF.NE.9) THEN IER=-2 RETURN ENDIF CALL PLADR0(NF,N,ICA,CG,CR,S,EPS7,GMAX,UMAX,INEW,NADD,IER) IF (IER.NE.0) RETURN IF (N.GT.0) THEN N1=N+1 IF (INEW.GT.0) THEN CALL MXDRMM(NF,N1,CZ,CG((INEW-1)*NF+1),S) ELSE CALL MXDRMV(NF,N1,CZ,S,-INEW) ENDIF DO 1 L=1,N K=L+1 CALL MXVORT(S(K),S(L),CK,CL,IER) CALL MXDRGR(NF,CZ,K,L,CK,CL,IER) CALL MXDSMR(N1,H,K,L,CK,CL,IER) IF (IER.LT.0) RETURN 1 CONTINUE IF (IDECF.EQ.9) THEN L=N*(N+1)/2 IF (H(L+N1).NE.0.0D 0) THEN CL=1.0D 0/H(L+N1) K=0 DO 3 I=1,N CK=CL*H(L+I) DO 2 J=1,I K=K+1 H(K)=H(K)-CK*H(L+J) 2 CONTINUE 3 CONTINUE ENDIF ENDIF ENDIF IER=0 RETURN END * SUBROUTINE PLADR0 ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * TRIANGULAR DECOMPOSITION OF KERNEL OF THE ORTHOGONAL PROJECTION * IS UPDATED AFTER CONSTRAINT ADDITION. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * IU N ACTUAL NUMBER OF VARIABLES. * IU ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RU CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RA S(NF) AUXILIARY VECTOR. * RI EPS7 TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS. * RO GMAX MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE. * RO UMAX MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER. * II INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * IU NADD NUMBER OF CONSTRAINT ADDITIONS. * IO IER ERROR INDICATOR. * * SUBPROGRAMS USED : * S MXSPRB SPARSE BACK SUBSTITUTION. * S MXVCOP COPYING OF A VECTOR. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * SUBROUTINE PLADR0(NF,N,ICA,CG,CR,S,EPS7,GMAX,UMAX,INEW,NADD,IER) INTEGER NF,N,ICA(*),INEW,NADD,IER DOUBLE PRECISION CG(*),CR(*),S(*),EPS7,GMAX,UMAX DOUBLE PRECISION MXVDOT INTEGER NCA,NCR,I,J,K,L IER=0 IF (N.LE.0) IER=2 IF (INEW.EQ.0) IER=3 IF (IER.NE.0) RETURN NCA=NF-N NCR=NCA*(NCA+1)/2 IF (INEW.GT.0) THEN CALL MXVCOP(NF,CG((INEW-1)*NF+1),S) GMAX=MXVDOT(NF,CG((INEW-1)*NF+1),S) DO 1 J=1,NCA L=ICA(J) IF (L.GT.0) THEN CR(NCR+J)=MXVDOT(NF,CG((L-1)*NF+1),S) ELSE I=-L CR(NCR+J)=S(I) ENDIF 1 CONTINUE ELSE K=-INEW GMAX=1.0D 0 DO 2 J=1,NCA L=ICA(J) IF (L.GT.0) THEN CR(NCR+J)=CG((L-1)*NF+K)*GMAX ELSE CR(NCR+J)=0.0D 0 ENDIF 2 CONTINUE ENDIF IF (NCA.EQ.0) THEN UMAX=GMAX ELSE CALL MXDPRB(NCA,CR,CR(NCR+1),1) UMAX=GMAX-MXVDOT(NCA,CR(NCR+1),CR(NCR+1)) ENDIF IF (UMAX.LE.EPS7*GMAX) THEN IER=1 RETURN ELSE N=N-1 NCA=NCA+1 NCR=NCR+NCA ICA(NCA)=INEW CR(NCR)=SQRT(UMAX) NADD=NADD+1 ENDIF RETURN END * SUBROUTINE PLLPB1 ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF THE INITIAL FEASIBLE POINT AND THE LINEAR PROGRAMMING * SUBROUTINE. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NC NUMBER OF LINEAR CONSTRAINTS. * RU X(NF) VECTOR OF VARIABLES. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * RO XO(NF) SAVED VECTOR OF VARIABLES. * RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. * RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. * RU CF(NF) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS. * RA CFD(NF) VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT * FUNCTIONS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * IO ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CL(NC) VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CU(NC) VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RO CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RO CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RO GO(NF) SAVED GRADIENT OF THE OBJECTIVE FUNCTION. * RA S(NF) DIRECTION VECTOR. * II MFP TYPE OF FEASIBLE POINT. MFP=1-ARBITRARY FEASIBLE POINT. * MFP=2-OPTIMUM FEASIBLE POINT. MFP=3-REPEATED SOLUTION. * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. * II KBC SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR * CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO * SIDED LINEAR CONSTRAINTS. * RI ETA9 MAXIMUM FOR REAL NUMBERS. * RI EPS7 TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS. * RI EPS9 TOLERANCE FOR ACTIVITY OF CONSTRAINTS. * RO UMAX MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER. * RO GMAX MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE. * IO N DIMENSION OF THE MANIFOLD DEFINED BY ACTIVE CONSTRAINTS. * IO ITERL TYPE OF FEASIBLE POINT. ITERL=1-ARBITRARY FEASIBLE POINT. * ITERL=2-OPTIMUM FEASIBLE POINT. ITERL=-1 FEASIBLE POINT DOES * NOT EXISTS. ITERL=-2 OPTIMUM FEASIBLE POINT DOES NOT EXISTS. * * SUBPROGRAMS USED : * S PLINIT DETERMINATION OF INITIAL POINT SATISFYING SIMPLE BOUNDS. * S PLMAXL MAXIMUM STEPSIZE USING LINEAR CONSTRAINTS. * S PLMAXS MAXIMUM STEPSIZE USING SIMPLE BOUNDS. * S PLMAXT MAXIMUM STEPSIZE USING TRUST REGION BOUNDS. * S PLNEWL IDENTIFICATION OF ACTIVE LINEAR CONSTRAINTS. * S PLNEWS IDENTIFICATION OF ACTIVE SIMPLE BOUNDS. * S PLNEWT IDENTIFICATION OF ACTIVE TRUST REGION BOUNDS. * S PLDIRL NEW VALUES OF CONSTRAINT FUNCTIONS. * S PLDIRS NEW VALUES OF VARIABLES. * S PLSETC INITIAL VALUES OF CONSTRAINT FUNCTIONS. * S PLSETG DETERMINATION OF THE FIRST PHASE GRADIENT VECTOR. * S PLTRBG DETERMINATION OF LAGRANGE MULTIPLIERS AND COMPUTATION * S PLADB0 CONSTRAINT ADDITION. * S PLRMB0 CONSTRAINT DELETION. * S MXDCMM PREMULTIPLICATION OF A VECTOR BY A DENSE RECTANGULAR * MATRIX STORED BY COLUMNS. * S MXDRMM PREMULTIPLICATION OF A VECTOR BY A DENSE RECTANGULAR * MATRIX STORED BY ROWS. * S MXDSMI DETERMINATION OF THE INITIAL UNIT DENSE SYMMETRIC * MATRIX. * S MXVCOP COPYING OF A VECTOR. * S MXVINA ABSOLUTE VALUES OF ELEMENTS OF AN INTEGER VECTOR. * S MXVINC UPDATE OF AN INTEGER VECTOR. * S MXVIND CHANGE OF THE INTEGER VECTOR FOR CONSTRAINT ADDITION. * S MXVINT CHANGE OF THE INTEGER VECTOR FOR TRUST REGION BOUND * ADDITION. * S MXVMUL DIAGONAL PREMULTIPLICATION OF A VECTOR. * S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE PLLPB1(NF,NC,X,IX,XO,XL,XU,CF,CFD,IC,ICA,CL,CU,CG, & CR,CZ,G,GO,S,MFP,KBF,KBC,ETA9,EPS7,EPS9,UMAX,GMAX,N,ITERL) INTEGER NF,NC,IX(*),IC(*),ICA(*),MFP,KBF,KBC,N,ITERL DOUBLE PRECISION X(*),XO(*),XL(*),XU(*),CF(*),CFD(*),CL(*), & CU(*),CG(*),CR(*),CZ(*),G(*),GO(*),S(*),ETA9,EPS7,EPS9, & UMAX,GMAX DOUBLE PRECISION POM,CON,DMAX INTEGER NCA,NCR,NCZ,IPOM,I,K,IOLD,INEW,IER,KREM,KC,NRED INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH COMMON /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH CON=ETA9 * * INITIATION * CALL MXVCOP(NF,X,XO) IPOM=0 NRED=0 KREM=0 ITERL=1 DMAX=0.0D 0 IF (MFP.EQ.3) GO TO 5 IF (KBF.GT.0) CALL MXVINA(NF,IX) * * SHIFT OF VARIABLES FOR SATISFYING SIMPLE BOUNDS * CALL PLINIT(NF,X,IX,XL,XU,EPS9,KBF,INEW,ITERL) IF (ITERL.LT.0) THEN GO TO 6 ENDIF N=0 NCA=0 NCZ=0 DO 1 I=1,NF IF (KBF.GT.0.AND.IX(I).LT.0) THEN NCA=NCA+1 ICA(NCA)=-I ELSE N=N+1 CALL MXVSET(NF,0.0D 0,CZ(NCZ+1)) CZ(NCZ+I)=1.0D 0 NCZ=NCZ+NF ENDIF 1 CONTINUE CALL MXDSMI(NCA,CR) IF (NC.GT.0) THEN CALL MXDRMM(NF,NC,CG,X,CF) * * ADDITION OF ACTIVE CONSTRAINTS AND INITIAL CHECK OF FEASIBILITY * CALL MXVINA(NC,IC) C IF (NF.GT.N) CALL PLSETC(NF,NC,X,XO,CF,IC,CG,S) DO 2 KC=1,NC IF (IC(KC).NE.0) THEN INEW=0 CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW) CALL PLADB0(NF,N,ICA,CG,CR,CZ,S,EPS7,GMAX,UMAX,INEW,NADD,IER) CALL MXVIND(IC,KC,IER) IF (IC(KC).LT.-10) IPOM=1 ENDIF 2 CONTINUE ENDIF 3 IF (IPOM.EQ.1) THEN * * CHECK OF FEASIBILITY AND UPDATE OF THE FIRST PHASE OBJECTIVE * FUNCTION * CALL PLSETG(NF,NC,IC,CG,GO,INEW) IF (INEW.EQ.0) IPOM=0 ENDIF IF (IPOM.EQ.0.AND.ITERL.EQ.0) THEN * * FEASIBILITY ACHIEVED * ITERL=1 CALL MXVCOP(NF,G,GO) IF (MFP.EQ.1) GO TO 6 ENDIF * * LAGRANGE MULTIPLIERS AND REDUCED GRADIENT DETERMINATION * 5 CALL PLTRBG(NF,N,NC,IX,IC,ICA,CG,CR,CZ,GO,S,EPS7,GMAX,UMAX,IOLD) INEW=0 IF (GMAX.EQ.0.0D 0) THEN * * OPTIMUM ON A LINEAR MANIFOLD OBTAINED * IF (IOLD.EQ.0) THEN IF (IPOM.EQ.0) THEN * * OPTIMAL SOLUTION ACHIEVED * ITERL= 2 GO TO 6 ELSE IPOM=0 DO 22 KC=1,NC IF (IC(KC).LT.-10) THEN INEW=0 CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW) IF (IC(KC).LT.-10) IPOM=1 ENDIF 22 CONTINUE IF (IPOM.EQ.0) THEN * * OPTIMAL SOLUTION ACHIEVED * CALL MXVCOP(NF,GO,G) ITERL= 2 GO TO 6 ELSE * * FEASIBLE SOLUTION DOES NOT EXIST * CALL MXVCOP(NF,GO,G) ITERL=-1 GO TO 6 ENDIF ENDIF ELSE * * CONSTRAINT DELETION * CALL PLRMB0(NF,N,ICA,CG,CR,CZ,GO,S,IOLD,KREM,NREM,IER) KC=ICA(NF-N+1) IF (KC.GT.0) THEN IC(KC)=-IC(KC) ELSE K=-KC IX(K)=-IX(K) ENDIF DMAX=0.0D 0 GO TO 5 ENDIF ELSE * * DIRECTION DETERMINATION * NCA=NF-N NCR=NCA*(NCA+1)/2 CALL MXDCMM(NF,N,CZ,S,CR(NCR+1)) CALL MXVNEG(NF,CR(NCR+1),S) * * STEPSIZE SELECTION * POM=CON CALL PLMAXL(NF,NC,CF,CFD,IC,CL,CU,CG,S,POM,KBC,KREM,INEW) CALL PLMAXS(NF,X,IX,XL,XU,S,POM,KBF,KREM,INEW) IF (INEW.EQ.0) THEN IF (IPOM.EQ.0) THEN * * BOUNDED SOLUTION DOES NOT EXIST * ITERL=-2 ELSE * * FEASIBLE SOLUTION DOES NOT EXIST * ITERL=-3 ENDIF GO TO 6 ELSE * * STEP REALIZATION * CALL PLDIRS(NF,X,IX,S,POM,KBF) CALL PLDIRL(NC,CF,CFD,IC,POM,KBC) * * CONSTRAINT ADDITION * IF (INEW.GT.0) THEN KC=INEW INEW=0 CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW) CALL PLADB0(NF,N,ICA,CG,CR,CZ,S,EPS7,GMAX,UMAX,INEW,NADD,IER) CALL MXVIND(IC,KC,IER) ELSE IF (INEW+NF.GE.0) THEN I=-INEW INEW=0 CALL PLNEWS(X,IX,XL,XU,EPS9,I,INEW) CALL PLADB0(NF,N,ICA,CG,CR,CZ,S,EPS7,GMAX,UMAX,INEW,NADD,IER) CALL MXVIND(IX,I,IER) ENDIF DMAX=POM IF (DMAX.GT.0.0D 0) NRED=NRED+1 GO TO 3 ENDIF ENDIF 6 CONTINUE RETURN END * SUBROUTINE PLRMB0 ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * OLD LINEAR CONSTRAINT OR AN OLD SIMPLE BOUND IS REMOVED FROM THE * ACTIVE SET. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * II ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RU CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RU CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RU GN(NF) TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION. * II IOLD INDEX OF THE OLD ACTIVE CONSTRAINT. * IO KREM AUXILIARY VARIABLE. * IU NREM NUMBER OF CONSTRAINT DELETION. * IO IER ERROR INDICATOR. * * SUBPROGRAMS USED : * S PLRMR0 CORRECTION OF KERNEL OF THE ORTHOGONAL PROJECTION * AFTER CONSTRAINT DELETION. * S MXDPRB BACK SUBSTITUTION. * S MXVCOP COPYING OF A VECTOR. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * S MXVMUL DIAGONAL PREMULTIPLICATION OF A VECTOR. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE PLRMB0(NF,N,ICA,CG,CR,CZ,G,GN,IOLD,KREM,NREM,IER) INTEGER NF,N,ICA(*),IOLD,KREM,NREM,IER DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*) DOUBLE PRECISION MXVDOT INTEGER NCA,NCR,NCZ,I,J,KC IER=0 IF (N.EQ.NF) IER=2 IF (IOLD.EQ.0) IER=3 IF (IER.NE.0) RETURN NCA=NF-N NCR=NCA*(NCA-1)/2 NCZ=N*NF CALL PLRMR0(NF,ICA,CR,CZ(NCZ+1),N,IOLD,KREM,IER) CALL MXVSET(NCA,0.0D 0,CZ(NCZ+1)) CZ(NCZ+NCA)=1.0D 0 CALL MXDPRB(NCA,CR,CZ(NCZ+1),-1) CALL MXVCOP(NCA,CZ(NCZ+1),CR(NCR+1)) N=N+1 CALL MXVSET(NF,0.0D 0,CZ(NCZ+1)) DO 4 J=1,NCA KC=ICA(J) IF (KC.GT.0) THEN CALL MXVDIR(NF,CR(NCR+J),CG((KC-1)*NF+1),CZ(NCZ+1),CZ(NCZ+1)) ELSE I=-KC CZ(NCZ+I)=CZ(NCZ+I)+CR(NCR+J) ENDIF 4 CONTINUE GN(N)=MXVDOT(NF,CZ(NCZ+1),G) NREM=NREM+1 IER=0 RETURN END * SUBROUTINE PLQDB1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * IO N DIMENSION OF THE MANIFOLD DEFINED BY ACTIVE CONSTRAINTS. * II NC NUMBER OF LINEAR CONSTRAINTS. * RU X(NF) VECTOR OF VARIABLES. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * II IXA(NF) VECTOR CONTAINING INFORMATION ON TRUST REGION ACTIVITY. * RI XN(NF) VECTOR OF SCALING FACTORS. * RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. * RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. * RI CF(NF) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS. * RO CFD(NC) VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT * FUNCTIONS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * II ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CL(NC) VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CU(NC) VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RO CZ(NF) VECTOR OF LAGRANGE MULTIPLIERS. * RO G(NF) GRADIENT OF THE LAGRANGIAN FUNCTION. * RO GO(NF) SAVED GRADIENT OF THE OBJECTIVE FUNCTION. * RU H(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OR INVERSION OF THE * HESSIAN MATRIX APPROXIMATION. * RI S(NF) DIRECTION VECTOR. * RI ETA2 TOLERANCE FOR POSITIVE DEFINITENESS OF THE HESSIAN MATRIX. * RI ETA9 MAXIMUM FOR REAL NUMBERS. * RI EPS7 TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS. * RI EPS9 TOLERANCE FOR ACTIVITY OF CONSTRAINTS. * RU XDEL TRUST REGION BOUND. * RO UMAX MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER. * RO GMAX MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE. * II MFP TYPE OF FEASIBLE POINT. MFP=1-ARBITRARY FEASIBLE POINT. * MFP=2-OPTIMUM FEASIBLE POINT. MFP=3-REPEATED SOLUTION. * * COMMON DATA : * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. * II KBC SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR * CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO * SIDED LINEAR CONSTRAINTS. * II NORMF SCALING SPECIFICATION. NORMF=0-NO SCALING PERFORMED. * NORMF=1-SCALING FACTORS ARE DETERMINED AUTOMATICALLY. * NORMF=2-SCALING FACTORS ARE SUPPLIED BY USER. * IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION. * IDECF=10-DIAGONAL MATRIX. * IU NDECF NUMBER OF DECOMPOSITIONS. * IO ITERQ TYPE OF FEASIBLE POINT. ITERQ=1-ARBITRARY FEASIBLE POINT. * ITERQ=2-OPTIMUM FEASIBLE POINT. ITERQ=-1 FEASIBLE POINT DOES * NOT EXISTS. ITERQ=-2 OPTIMUM FEASIBLE POINT DOES NOT EXISTS. * * SUBPROGRAMS USED : * S PLMINS DETERMINATION OF THE NEW ACTIVE SIMPLE BOUND. * S PLMINL DETERMINATION OF THE NEW ACTIVE LINEAR CONSTRAINT. * S PLMINT DETERMINATION OF THE NEW ACTIVE TRUST REGION BOUND. * S PLADR1 ADDITION OF A NEW ACTIVE CONSTRAINT. * S PLRMR0 CONSTRAIN DELETION. * S PLSOB1 TRANSFORMATION OF THE LOCAL SOLUTION TO THE SOLUTION * OF THE ORIGINAL QP PROBLEM. * S MXDPGF GILL-MURRAY DECOMPOSITION OF A DENSE SYMMETRIC MATRIX. * S MXDPGB BACK SUBSTITUTION AFTER GILL-MURRAY DECOMPOSITION. * S MXDPRB BACK SUBSTITUTION. * S MXDSMM MATRIX VECTOR PRODUCT. * S MXVCOP COPYING OF A VECTOR. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * S MXVINA ABSOLUTE VALUES OF ELEMENTS OF AN INTEGER VECTOR. * S MXVINV CHANGE OF AN INTEGER VECTOR AFTER CONSTRAINT ADDITION. * S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. * SUBROUTINE PLQDB1(NF,NC,X,IX,XL,XU,CF,CFD,IC,ICA,CL,CU,CG,CR,CZ, & G,GO,H,S,MFP,KBF,KBC,IDECF,ETA2,ETA9,EPS7,EPS9,UMAX,GMAX,N,ITERQ) INTEGER NF,NC,IX(*),IC(*),ICA(*),MFP,KBF,KBC,IDECF,N,ITERQ DOUBLE PRECISION X(*),XL(*),XU(*),CF(*),CFD(*),CL(*),CU(*),CG(*), & CR(*),CZ(*),G(*),GO(*),H(*),S(*),ETA2,ETA9,EPS7,EPS9,UMAX,GMAX DOUBLE PRECISION CON,TEMP,STEP,STEP1,STEP2,DMAX,PAR,SNORM INTEGER NCA,NCR,I,J,K,IOLD,JOLD,INEW,JNEW,KNEW,INF,IER,KREM,KC, & NRED INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH COMMON /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH CON=ETA9 IF (IDECF.LT.0) IDECF=1 IF (IDECF.EQ.0) THEN C C GILL-MURRAY DECOMPOSITION C TEMP=ETA2 CALL MXDPGF(NF,H,INF,TEMP,STEP) NDEC=NDEC+1 IDECF=1 ENDIF IF (IDECF.GE.2.AND.IDECF.LE.8) THEN ITERQ=-10 RETURN ENDIF C C INITIATION C NRED=0 JOLD=0 JNEW=0 ITERQ=0 DMAX=0.0D0 IF (MFP.EQ.3) GO TO 1 N=NF NCA=0 NCR=0 IF (KBF.GT.0) CALL MXVINA(NF,IX) IF (KBC.GT.0) CALL MXVINA(NC,IC) C C DIRECTION DETERMINATION C 1 CALL MXVNEG(NF,GO,S) DO 2 J=1,NCA KC=ICA(J) IF (KC.GT.0) THEN CALL MXVDIR(NF,CZ(J),CG((KC-1)*NF+1),S,S) ELSE K=-KC S(K)=S(K)+CZ(J) ENDIF 2 CONTINUE CALL MXVCOP(NF,S,G) IF (IDECF.EQ.1) THEN CALL MXDPGB(NF,H,S,0) ELSE CALL MXDSMM(NF,H,G,S) ENDIF IF (ITERQ.EQ.3) GO TO 7 C C CHECK OF FEASIBILITY C INEW=0 PAR=0.0D0 CALL PLMINN(NF,NC,CF,CFD,IC,CL,CU,CG,S,EPS9,PAR,KBC,INEW,KNEW) CALL PLMINS(NF,IX,X,XL,XU,S,KBF,INEW,KNEW,EPS9,PAR) IF (INEW.EQ.0) THEN C C SOLUTION ACHIEVED C CALL MXVNEG(NF,G,G) ITERQ=2 GO TO 7 ELSE SNORM=0.0D0 ENDIF 4 IER=0 C C STEPSIZE DETERMINATION C CALL PLADR1(NF,N,ICA,CG,CR,H,S,G,EPS7,GMAX,UMAX,IDECF,INEW, & NADD,IER,1) CALL MXDPRB(NCA,CR,G,-1) IF (KNEW.LT.0) CALL MXVNEG(NCA,G,G) C C PRIMAL STEPSIZE C IF (IER.NE.0) THEN STEP1=CON ELSE STEP1=-PAR/UMAX ENDIF C C DUAL STEPSIZE C IOLD=0 STEP2=CON DO 5 J=1,NCA KC=ICA(J) IF (KC.GE.0) THEN K=IC(KC) ELSE I=-KC K=IX(I) ENDIF IF (K.LE.-5) THEN ELSE IF ((K.EQ.-1.OR.K.EQ.-3.).AND.G(J).LE.0.0D0) THEN ELSE IF ((K.EQ.-2.OR.K.EQ.-4.).AND.G(J).GE.0.0D0) THEN ELSE TEMP=CZ(J)/G(J) IF (STEP2.GT.TEMP) THEN IOLD=J STEP2=TEMP ENDIF ENDIF 5 CONTINUE C C FINAL STEPSIZE C STEP=MIN(STEP1,STEP2) IF (STEP.GE.CON) THEN C C FEASIBLE SOLUTION DOES NOT EXIST C ITERQ=-1 GO TO 7 ENDIF C C NEW LAGRANGE MULTIPLIERS C DMAX=STEP CALL MXVDIR(NCA,-STEP,G,CZ,CZ) SNORM=SNORM+SIGN(1,KNEW)*STEP PAR=PAR-(STEP/STEP1)*PAR IF (STEP.EQ.STEP1) THEN IF (N.LE.0) THEN C C IMPOSSIBLE SITUATION C ITERQ=-5 GO TO 7 ENDIF C C CONSTRAINT ADDITION C IF (IER.EQ.0) THEN N=N-1 NCA=NCA+1 NCR=NCR+NCA CZ(NCA)=SNORM ENDIF IF (INEW.GT.0) THEN KC=INEW CALL MXVINV(IC,KC,KNEW) ELSE IF (ABS(KNEW).EQ.1) THEN I=-INEW CALL MXVINV(IX,I,KNEW) ELSE I=-INEW IF (KNEW.GT.0) IX(I)=-3 IF (KNEW.LT.0) IX(I)=-4 ENDIF NRED=NRED+1 NADD=NADD+1 JNEW=INEW JOLD=0 GO TO 1 ELSE C C CONSTRAINT DELETION C DO 6 J=IOLD,NCA-1 CZ(J)=CZ(J+1) 6 CONTINUE CALL PLRMF0(NF,NC,IX,IC,ICA,CR,IC,G,N,IOLD,KREM,IER) NCR=NCR-NCA NCA=NCA-1 JOLD=IOLD JNEW=0 IF (KBC.GT.0) CALL MXVINA(NC,IC) IF (KBF.GT.0) CALL MXVINA(NF,IX) DO 8 J=1,NCA KC=ICA(J) IF (KC.GT.0) THEN IC(KC)=-IC(KC) ELSE KC=-KC IX(KC)=-IX(KC) ENDIF 8 CONTINUE GO TO 4 ENDIF 7 CONTINUE RETURN END * SUBROUTINE PLADR1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * TRIANGULAR DECOMPOSITION OF KERNEL OF THE GENERAL PROJECTION * IS UPDATED AFTER CONSTRAINT ADDITION. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * IU N ACTUAL NUMBER OF VARIABLES. * IU ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RU H(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OR INVERSION OF THE * HESSIAN MATRIX APPROXIMATION. * RA S(NF) AUXILIARY VECTOR. * RO G(NF) VECTOR USED IN THE DUAL RANGE SPACE QUADRATIC PROGRAMMING * METHOD. * RI EPS7 TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS. * RO GMAX MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE. * RO UMAX MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER. * RO E AUXILIARY VARIABLE. * RI T AUXILIARY VARIABLE. * IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION. * IDECF=10-DIAGONAL MATRIX. * II INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * IU NADD NUMBER OF CONSTRAINT ADDITIONS. * IO IER ERROR INDICATOR. * II JOB SPECIFICATION OF COMPUTATION. OUTPUT VECTOR G IS NOT OR IS * COMPUTED IN CASE WHEN N.LE.0 IF JOB=0 OR JOB=1 RESPECTIVELY. * * SUBPROGRAMS USED : * S MXDPGB BACK SUBSTITUTION. * S MXDPRB BACK SUBSTITUTION. * S MXDSMM MATRIX-VECTOR PRODUCT. * S MXDSMV COPYING OF A ROW OF DENSE SYMMETRIC MATRIX. * S MXVCOP COPYING OF A VECTOR. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * SUBROUTINE PLADR1(NF,N,ICA,CG,CR,H,S,G,EPS7,GMAX,UMAX,IDECF, & INEW,NADD,IER,JOB) INTEGER NF,N,ICA(*),IDECF,INEW,NADD,IER,JOB DOUBLE PRECISION CG(*),CR(*),H(*),S(*),G(*),EPS7,GMAX,UMAX DOUBLE PRECISION MXVDOT INTEGER NCA,NCR,JCG,J,K,L IER=0 IF (JOB.EQ.0.AND.N.LE.0) IER=2 IF (INEW.EQ.0) IER=3 IF (IDECF.NE.1.AND.IDECF.NE.9) IER=-2 IF (IER.NE.0) RETURN NCA=NF-N NCR=NCA*(NCA+1)/2 IF (INEW.GT.0) THEN JCG=(INEW-1)*NF+1 IF (IDECF.EQ.1) THEN CALL MXVCOP(NF,CG(JCG),S) CALL MXDPGB(NF,H,S,0) ELSE CALL MXDSMM(NF,H,CG(JCG),S) ENDIF GMAX=MXVDOT(NF,CG(JCG),S) ELSE K=-INEW IF (IDECF.EQ.1) THEN CALL MXVSET(NF,0.0D0,S) S(K)=1.0D 0 CALL MXDPGB(NF,H,S,0) ELSE CALL MXDSMV(NF,H,S,K) ENDIF GMAX=S(K) ENDIF DO 1 J=1,NCA L=ICA(J) IF (L.GT.0) THEN G(J)=MXVDOT(NF,CG((L-1)*NF+1),S) ELSE L=-L G(J)=S(L) ENDIF 1 CONTINUE IF (N.EQ.0) THEN CALL MXDPRB(NCA,CR,G,1) UMAX=0.0D0 IER=2 RETURN ELSE IF (NCA.EQ.0) THEN UMAX=GMAX ELSE CALL MXDPRB(NCA,CR,G,1) UMAX=GMAX-MXVDOT(NCA,G,G) CALL MXVCOP(NCA,G,CR(NCR+1)) ENDIF IF (UMAX.LE.EPS7*GMAX) THEN IER=1 RETURN ELSE NCA=NCA+1 NCR=NCR+NCA ICA(NCA)=INEW CR(NCR)=SQRT(UMAX) IF (JOB.EQ.0) THEN N=N-1 NADD=NADD+1 ENDIF ENDIF RETURN END * SUBROUTINE PLDIRL ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF THE NEW VALUES OF THE CONSTRAINT FUNCTIONS. * * PARAMETERS : * II NC NUMBER OF CONSTRAINTS. * RU CF(NF) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS. * RI CFD(NF) VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT * FUNCTIONS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * RI STEP CURRENT STEPSIZE. * II KBC SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR * CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO * SIDED LINEAR CONSTRAINTS. * SUBROUTINE PLDIRL(NC,CF,CFD,IC,STEP,KBC) INTEGER NC,IC(*),KBC DOUBLE PRECISION CF(*),CFD(*),STEP INTEGER KC IF (KBC.GT.0) THEN DO 1 KC=1,NC IF (IC(KC).GE.0.AND.IC(KC).LE.10) THEN CF(KC)=CF(KC)+STEP*CFD(KC) ELSE IF (IC(KC).LT.-10) THEN CF(KC)=CF(KC)+STEP*CFD(KC) ENDIF 1 CONTINUE ENDIF RETURN END * SUBROUTINE PLDIRS ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF THE NEW VECTOR OF VARIABLES. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * RU X(NF) VECTOR OF VARIABLES. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * RI S(NF) DIRECTION VECTOR. * RI STEP CURRENT STEPSIZE. * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. * SUBROUTINE PLDIRS(NF,X,IX,S,STEP,KBF) INTEGER NF,IX(*),KBF DOUBLE PRECISION X(*),S(*),STEP INTEGER I DO 1 I=1,NF IF (KBF.LE.0) THEN X(I)=X(I)+STEP*S(I) ELSE IF (IX(I).GE.0.AND.IX(I).LE.10) THEN X(I)=X(I)+STEP*S(I) ELSE IF (IX(I).LT.-10) THEN X(I)=X(I)+STEP*S(I) ENDIF 1 CONTINUE RETURN END * SUBROUTINE PLINIT ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF THE INITIAL POINT WHICH SATISFIES SIMPLE BOUNDS. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * RU X(NF) VECTOR OF VARIABLES. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. * RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. * RI EPS9 TOLERANCE FOR ACTIVE CONSTRAINTS. * IO INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * IO IND INDICATOR. IF IND.NE.0 THEN TRUST REGION BOUNDS CANNOT * BE SATISFIED. * * SUBPROGRAMS USED : * S PLNEWS TEST ON ACTIVITY OF A GIVEN SIMPLE BOUND. * SUBROUTINE PLINIT(NF,X,IX,XL,XU,EPS9,KBF,INEW,IND) INTEGER NF,IX(*),KBF,INEW,IND DOUBLE PRECISION X(*),XL(*),XU(*),EPS9 INTEGER I IND=0 IF (KBF.GT.0) THEN DO 1 I=1,NF CALL PLNEWS(X,IX,XL,XU,EPS9,I,INEW) IF (IX(I).LT.5) THEN ELSE IF (IX(I).EQ.5) THEN IX(I)=-5 ELSE IF (IX(I).EQ.11.OR.IX(I).EQ.13) THEN X(I)=XL(I) IX(I)=10-IX(I) ELSE IF (IX(I).EQ.12.OR.IX(I).EQ.14) THEN X(I)=XU(I) IX(I)=10-IX(I) ENDIF 1 CONTINUE ENDIF RETURN END * SUBROUTINE PLMAXL ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF THE MAXIMUM STEPSIZE USING LINEAR CONSTRAINTS. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II NC NUMBER OF CURRENT LINEAR CONSTRAINTS. * RI CF(NF) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCYIONS. * RO CFD(NF) VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT * FUNCTIONS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * RI CL(NC) VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CU(NC) VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI S(NF) DIRECTION VECTOR. * RO STEP MAXIMUM STEPSIZE. * II KBC SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR * CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO * SIDED LINEAR CONSTRAINTS. * II KREM INDICATION OF LINEARLY DEPENDENT GRADIENTS. * IO INEW INDEX OF THE NEW ACTIVE FUNCTION. * * SUBPROGRAMS USED : * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * SUBROUTINE PLMAXL(NF,NC,CF,CFD,IC,CL,CU,CG,S,STEP,KBC,KREM,INEW) INTEGER NF,NC,IC(*),KBC,KREM,INEW DOUBLE PRECISION CF(*),CFD(*),CL(*),CU(*),CG(*),S(*),STEP DOUBLE PRECISION TEMP,MXVDOT INTEGER JCG,KC IF (KBC.GT.0) THEN JCG=1 DO 1 KC=1,NC IF (KREM.GT.0.AND.IC(KC).GT.10) IC(KC)=IC(KC)-10 IF (IC(KC).GT.0.AND.IC(KC).LE.10) THEN TEMP=MXVDOT(NF,CG(JCG),S) CFD(KC)=TEMP IF (TEMP.LT.0.0D 0) THEN IF (IC(KC).EQ.1.OR.IC(KC).GE.3) THEN TEMP=(CL(KC)-CF(KC))/TEMP IF (TEMP.LE.STEP) THEN INEW=KC STEP=TEMP ENDIF ENDIF ELSE IF (TEMP.GT.0.0D 0) THEN IF (IC(KC).EQ.2.OR.IC(KC).GE.3) THEN TEMP=(CU(KC)-CF(KC))/TEMP IF (TEMP.LE.STEP) THEN INEW=KC STEP=TEMP ENDIF ENDIF ENDIF ELSE IF (IC(KC).LT.-10) THEN TEMP=MXVDOT(NF,CG(JCG),S) CFD(KC)=TEMP IF (TEMP.GT.0.0D 0) THEN IF (IC(KC).EQ.-11.OR.IC(KC).EQ.-13.OR.IC(KC).EQ.-15) THEN TEMP=(CL(KC)-CF(KC))/TEMP IF (TEMP.LE.STEP) THEN INEW=KC STEP=TEMP ENDIF ENDIF ELSE IF (TEMP.LT.0.0D 0) THEN IF (IC(KC).EQ.-12.OR.IC(KC).EQ.-14.OR.IC(KC).EQ.-16) THEN TEMP=(CU(KC)-CF(KC))/TEMP IF (TEMP.LE.STEP) THEN INEW=KC STEP=TEMP ENDIF ENDIF ENDIF ENDIF JCG=JCG+NF 1 CONTINUE ENDIF RETURN END * SUBROUTINE PLMAXS ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF THE MAXIMUM STEPSIZE USING THE SIMPLE BOUNDS * FOR VARIABLES. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * RI X(NF) VECTOR OF VARIABLES. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. * RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. * RI S(NF) DIRECTION VECTOR. * RO STEP MAXIMUM STEPSIZE. * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. * IO KREM INDICATION OF LINEARLY DEPENDENT GRADIENTS. * IO INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * SUBROUTINE PLMAXS(NF,X,IX,XL,XU,S,STEP,KBF,KREM,INEW) INTEGER NF,IX(*),KBF,KREM,INEW DOUBLE PRECISION X(*),XL(*),XU(*),S(*),STEP DOUBLE PRECISION TEMP INTEGER I IF (KBF.GT.0) THEN DO 1 I=1,NF IF (KREM.GT.0.AND.IX(I).GT.10) IX(I)=IX(I)-10 IF (IX(I).GT.0.AND.IX(I).LE.10) THEN IF (S(I).LT.0.0D 0) THEN IF (IX(I).EQ.1.OR.IX(I).GE.3) THEN TEMP=(XL(I)-X(I))/S(I) IF (TEMP.LE.STEP) THEN INEW=-I STEP=TEMP ENDIF ENDIF ELSE IF (S(I).GT.0.0D 0) THEN IF (IX(I).EQ.2.OR.IX(I).GE.3) THEN TEMP=(XU(I)-X(I))/S(I) IF (TEMP.LE.STEP) THEN INEW=-I STEP=TEMP ENDIF ENDIF ENDIF ENDIF 1 CONTINUE ENDIF KREM=0 RETURN END * SUBROUTINE PLNEWL ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * TEST ON ACTIVITY OF A GIVEN LINEAR CONSTRAINT. * * PARAMETERS : * II KC INDEX OF A GIVEN CONSTRAINT. * RI CF(NC) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS. * IU IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * RI CL(NC) VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CU(NC) VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI EPS9 TOLERANCE FOR ACTIVE CONSTRAINTS. * IO INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * SUBROUTINE PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW) INTEGER KC,IC(*),INEW DOUBLE PRECISION CF(*),CL(*),CU(*),EPS9 DOUBLE PRECISION TEMP IF (IC(KC).LT.-10) IC(KC)=-IC(KC)-10 IF (IC(KC).LE.0) THEN ELSE IF (IC(KC).EQ.1) THEN TEMP=EPS9*MAX(ABS(CL(KC)),1.0D 0) IF (CF(KC).GT.CL(KC)+TEMP) THEN ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN IC(KC)=11 INEW=KC ELSE IC(KC)=-11 ENDIF ELSE IF (IC(KC).EQ.2) THEN TEMP=EPS9*MAX(ABS(CU(KC)),1.0D 0) IF (CF(KC).LT.CU(KC)-TEMP) THEN ELSE IF (CF(KC).LE.CU(KC)+TEMP) THEN IC(KC)=12 INEW=KC ELSE IC(KC)=-12 ENDIF ELSE IF (IC(KC).EQ.3.OR.IC(KC).EQ.4) THEN TEMP=EPS9*MAX(ABS(CL(KC)),1.0D 0) IF (CF(KC).GT.CL(KC)+TEMP) THEN TEMP=EPS9*MAX(ABS(CU(KC)),1.0D 0) IF (CF(KC).LT.CU(KC)-TEMP) THEN ELSE IF (CF(KC).LE.CU(KC)+TEMP) THEN IC(KC)=14 INEW=KC ELSE IC(KC)=-14 ENDIF ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN IC(KC)=13 INEW=KC ELSE IC(KC)=-13 ENDIF ELSE IF (IC(KC).EQ.5.OR.IC(KC).EQ.6) THEN TEMP=EPS9*MAX(ABS(CL(KC)),1.0D 0) IF (CF(KC).GT.CL(KC)+TEMP) THEN TEMP=EPS9*MAX(ABS(CU(KC)),1.0D 0) IF (CF(KC).LT.CU(KC)-TEMP) THEN ELSE IF (CF(KC).LE.CU(KC)+TEMP) THEN IC(KC)=16 INEW=KC ELSE IC(KC)=-16 ENDIF ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN IC(KC)=15 INEW=KC ELSE IC(KC)=-15 ENDIF ENDIF RETURN END * SUBROUTINE PLMINN ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF THE NEW ACTIVE LINEAR CONSTRAINT. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * RI CF(NC) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS. * RO CFD(NC) VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT * FUNCTIONS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * RI CL(NC) VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CU(NC) VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI S(NF) DIRECTION VECTOR. * RI EPS9 TOLERANCE FOR ACTIVE CONSTRAINTS. * RA PAR AUXILIARY VARIABLE. * II KBC SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR * CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO * SIDED LINEAR CONSTRAINTS. * IO INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * IO KNEW SIGNUM OF THE NEW ACTIVE NORMAL. * * SUBPROGRAMS USED : * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * SUBROUTINE PLMINN(NF,NC,CF,CFD,IC,CL,CU,CG,S,EPS9,PAR,KBC,INEW, & KNEW) INTEGER NF,NC,IC(*),KBC,INEW,KNEW DOUBLE PRECISION CF(*),CFD(*),CL(*),CU(*),CG(*),S(*),EPS9,PAR DOUBLE PRECISION TEMP,POM,MXVDOT INTEGER JCG,KC IF (KBC.GT.0) THEN JCG=1 DO 1 KC=1,NC IF (IC(KC).GT.0) THEN TEMP=MXVDOT(NF,CG(JCG),S) CFD(KC)=TEMP TEMP=CF(KC)+TEMP IF (IC(KC).EQ.1.OR.IC(KC).GE.3) THEN POM=TEMP-CL(KC) IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(CL(KC)),1.0D 0))) THEN INEW=KC KNEW= 1 PAR=POM ENDIF ENDIF IF (IC(KC).EQ.2.OR.IC(KC).GE.3) THEN POM=CU(KC)-TEMP IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(CU(KC)),1.0D 0))) THEN INEW=KC KNEW=-1 PAR=POM ENDIF ENDIF ENDIF JCG=JCG+NF 1 CONTINUE ENDIF RETURN END * SUBROUTINE PLMINS ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF THE NEW ACTIVE SIMPLE BOUND. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * RI XO(NF) SAVED VECTOR OF VARIABLES. * RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. * RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. * RI S(NF) DIRECTION VECTOR. * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. * IO INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * IO KNEW SIGNUM OF THE NEW NORMAL. * RI EPS9 TOLERANCE FOR ACTIVE CONSTRAINTS. * RA PAR AUXILIARY VARIABLE. * SUBROUTINE PLMINS(NF,IX,XO,XL,XU,S,KBF,INEW,KNEW,EPS9,PAR) DOUBLE PRECISION EPS9,PAR INTEGER INEW,KBF,KNEW,NF DOUBLE PRECISION S(*),XL(*),XO(*),XU(*) INTEGER IX(*) DOUBLE PRECISION POM,TEMP INTEGER I IF (KBF.GT.0) THEN DO 10 I = 1,NF IF (IX(I).GT.0) THEN TEMP = 1.0D0 IF (IX(I).EQ.1 .OR. IX(I).GE.3) THEN POM = XO(I) + S(I)*TEMP - XL(I) IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(XL(I)), + TEMP))) THEN INEW = -I KNEW = 1 PAR = POM END IF END IF IF (IX(I).EQ.2 .OR. IX(I).GE.3) THEN POM = XU(I) - S(I)*TEMP - XO(I) IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(XU(I)), + TEMP))) THEN INEW = -I KNEW = -1 PAR = POM END IF END IF END IF 10 CONTINUE END IF RETURN END * SUBROUTINE PLNEWS ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * TEST ON ACTIVITY OF A GIVEN SIMPLE BOUND. * * PARAMETERS : * RI X(NF) VECTOR OF VARIABLES. * IU IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. * RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. * RI EPS9 TOLERANCE FOR ACTIVE CONSTRAINTS. * II I INDEX OF TESTED SIMPLE BOUND. * IO INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * SUBROUTINE PLNEWS(X,IX,XL,XU,EPS9,I,INEW) INTEGER IX(*),I,INEW DOUBLE PRECISION X(*),XL(*),XU(*),EPS9 DOUBLE PRECISION TEMP TEMP=1.0D 0 IF (IX(I).LE.0) THEN ELSE IF (IX(I).EQ.1) THEN IF (X(I).LE.XL(I)+EPS9*MAX(ABS(XL(I)),TEMP)) THEN IX(I)=11 INEW=-I ENDIF ELSE IF (IX(I).EQ.2) THEN IF (X(I).GE.XU(I)-EPS9*MAX(ABS(XU(I)),TEMP)) THEN IX(I)=12 INEW=-I ENDIF ELSE IF (IX(I).EQ.3.OR.IX(I).EQ.4) THEN IF (X(I).LE.XL(I)+EPS9*MAX(ABS(XL(I)),TEMP)) THEN IX(I)=13 INEW=-I ENDIF IF (X(I).GE.XU(I)-EPS9*MAX(ABS(XU(I)),TEMP)) THEN IX(I)=14 INEW=-I ENDIF ENDIF RETURN END * SUBROUTINE PLREDL ALL SYSTEMS 98/12/01 C PORTABILITY : ALL SYSTEMS C 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * TRANSFORMATION OF THE INCOMPATIBLE QUADRATIC PROGRAMMING SUBPROBLEM. * * PARAMETERS : * II NC NUMBER OF CURRENT LINEAR CONSTRAINTS. * RI CF(NF) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * RI CL(NC) VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CU(NC) VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS. * II KBC SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR * CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO * SIDED LINEAR CONSTRAINTS. * SUBROUTINE PLREDL(NC,CF,IC,CL,CU,KBC) INTEGER NC,IC(NC),KBC,K DOUBLE PRECISION CF(*),CL(*),CU(*) DOUBLE PRECISION TEMP INTEGER KC IF (KBC.GT.0) THEN DO 1 KC=1,NC K=IC(KC) IF (ABS(K).EQ.1.OR.ABS(K).EQ.3.OR.ABS(K).EQ.4) THEN TEMP=(CF(KC)-CL(KC)) IF (TEMP.LT.0) CF(KC)=CL(KC)+0.1D 0*TEMP ENDIF IF (ABS(K).EQ.2.OR.ABS(K).EQ.3.OR.ABS(K).EQ.4) THEN TEMP=(CF(KC)-CU(KC)) IF (TEMP.GT.0) CF(KC)=CU(KC)+0.1D 0*TEMP ENDIF IF (ABS(K).EQ.5.OR.ABS(K).EQ.6) THEN TEMP=(CF(KC)-CL(KC)) CF(KC)=CL(KC)+0.1D 0*TEMP ENDIF 1 CONTINUE ENDIF RETURN END * SUBROUTINE PLRMF0 ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * OPERATIONS AFTER CONSTRAINT DELETION. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * II IA(NA) VECTOR CONTAINING TYPES OF DEVIATIONS. * IU IAA(NF+1) VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS. * RU AR((NF+1)*(NF+2)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * RA S(NF+1) AUXILIARY VECTOR. * II N ACTUAL NUMBER OF VARIABLES. * II IOLD INDEX OF THE OLD ACTIVE CONSTRAINT. * IO KREM AUXILIARY VARIABLE. * IO IER ERROR INDICATOR. * * SUBPROGRAMS USED : * S PLRMR0 CORRECTION OF KERNEL OF THE ORTHOGONAL PROJECTION * AFTER CONSTRAINT DELETION. * SUBROUTINE PLRMF0(NF,NC,IX,IA,IAA,AR,IC,S,N,IOLD,KREM,IER) INTEGER IER,IOLD,KREM,N,NC,NF DOUBLE PRECISION AR(*),S(*) INTEGER IA(*),IAA(*),IC(*),IX(*) INTEGER NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES INTEGER L COMMON /STAT/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH CALL PLRMR0(NF,IAA,AR,S,N,IOLD,KREM,IER) N = N + 1 NREM = NREM + 1 L = IAA(NF-N+1) IF (L.GT.NC) THEN L = L - NC IA(L) = -IA(L) ELSE IF (L.GT.0) THEN IC(L) = -IC(L) ELSE L = -L IX(L) = -IX(L) END IF RETURN END * SUBROUTINE PLRMR0 ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * TRIANGULAR DECOMPOSITION OF KERNEL OF THE ORTHOGONAL PROJECTION IS * UPDATED AFTER CONSTRAINT DELETION. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * IU ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RU CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RA G(NF) AUXILIARY VECTOR. * II N ACTUAL NUMBER OF VARIABLES. * II IOLD INDEX OF THE OLD ACTIVE CONSTRAINT. * IO KREM AUXILIARY VARIABLE. * IO IER ERROR INDICATOR. * * SUBPROGRAMS USED : * S MXVCOP COPYING OF A VECTOR. * S MXVORT DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR * PLANE ROTATION. * S MXVROT PLANE ROTATION OF A VECTOR. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE PLRMR0(NF,ICA,CR,G,N,IOLD,KREM,IER) INTEGER IER,IOLD,KREM,N,NF DOUBLE PRECISION CR(*),G(*) INTEGER ICA(*) DOUBLE PRECISION CK,CL INTEGER I,J,K,KC,L,NCA NCA = NF - N IF (IOLD.LT.NCA) THEN K = IOLD* (IOLD-1)/2 KC = ICA(IOLD) CALL MXVCOP(IOLD,CR(K+1),G) CALL MXVSET(NCA-IOLD,0.0D0,G(IOLD+1)) K = K + IOLD DO 20 I = IOLD + 1,NCA K = K + I CALL MXVORT(CR(K-1),CR(K),CK,CL,IER) CALL MXVROT(G(I-1),G(I),CK,CL,IER) L = K DO 10 J = I,NCA - 1 L = L + J CALL MXVROT(CR(L-1),CR(L),CK,CL,IER) 10 CONTINUE 20 CONTINUE K = IOLD* (IOLD-1)/2 DO 30 I = IOLD,NCA - 1 L = K + I ICA(I) = ICA(I+1) CALL MXVCOP(I,CR(L+1),CR(K+1)) K = L 30 CONTINUE ICA(NCA) = KC CALL MXVCOP(NCA,G,CR(K+1)) END IF KREM = 1 RETURN END * SUBROUTINE PLSETC ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF INITIAL VALUES OF THE CONSTRAINT FUNCTIONS. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NC NUMBER OF CURRENT LINEAR CONSTRAINTS. * RI X(NF) VECTOR OF VARIABLES. * RI XO(NF) SAVED VECTOR OF VARIABLES. * RU CF(NF) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCYIONS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * RI CG(NF*MCL) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RA S(NF) AUXILIARY VECTOR. * * SUBPROGRAMS USED : * S MXVDIF DIFFERENCE OF TWO VECTORS. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * S MXVMUL DIAGONAL PREMULTIPLICATION OF A VECTOR. * SUBROUTINE PLSETC(NF,NC,X,XO,CF,IC,CG,S) INTEGER NF,NC,IC(*) DOUBLE PRECISION X(*),XO(*),CF(*),CG(*),S(*) DOUBLE PRECISION MXVDOT INTEGER JCG,KC CALL MXVDIF(NF,X,XO,S) JCG=0 DO 1 KC=1,NC IF (IC(KC).NE.0) CF(KC)=CF(KC)+MXVDOT(NF,S,CG(JCG+1)) JCG=JCG+NF 1 CONTINUE RETURN END * SUBROUTINE PLSETG ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * GRADIENT DETERMINATION IN THE FIRST PHASE OF LP SUBROUTINE. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * IO INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * * SUBPROGRAMS USED : * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE PLSETG(NF,NC,IC,CG,G,INEW) INTEGER NF,NC,IC(*),INEW DOUBLE PRECISION CG(*),G(*) INTEGER KC CALL MXVSET(NF,0.0D 0,G) INEW=0 DO 4 KC=1,NC IF (IC(KC).GE.-10) THEN ELSE IF (IC(KC).EQ.-11.OR.IC(KC).EQ.-13.OR.IC(KC).EQ.-15) THEN CALL MXVDIR(NF,-1.0D 0,CG((KC-1)*NF+1),G,G) INEW=1 ELSE IF (IC(KC).EQ.-12.OR.IC(KC).EQ.-14.OR.IC(KC).EQ.-16) THEN CALL MXVDIR(NF, 1.0D 0,CG((KC-1)*NF+1),G,G) INEW=1 ENDIF 4 CONTINUE RETURN END * SUBROUTINE PLTLAG ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER IS * COMPUTED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * II NC NUMBER OF LINEARIZED CONSTRAINTS. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * II IA(NA) VECTOR CONTAINING TYPES OF DEVIATIONS. * II IAA(NF+1) VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS. * RI AZ(NF+1) VECTOR OF LAGRANGE MULTIPLIERS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * RI EPS7 TOLERANCE FOR LINEAR AND QUADRATIC PROGRAMMING. * RO UMAX MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER. * IO IOLD INDEX OF THE REMOVED CONSTRAINT. * SUBROUTINE PLTLAG(NF,N,NC,IX,IA,IAA,AZ,IC,EPS7,UMAX,IOLD) INTEGER NF,N,NC,IX(*),IA(*),IAA(*),IC(*),IOLD DOUBLE PRECISION AZ(*),EPS7,UMAX DOUBLE PRECISION TEMP INTEGER NAA,J,K,L IOLD=0 UMAX=0.0D 0 NAA=NF-N DO 2 J=1,NAA TEMP=AZ(J) L=IAA(J) IF (L.GT.NC) THEN L=L-NC K=IA(L) ELSE IF (L.GT.0) THEN K=IC(L) ELSE L=-L K=IX(L) ENDIF IF (K.LE.-5) THEN ELSE IF ((K.EQ.-1.OR.K.EQ.-3).AND.UMAX+TEMP.GE.0.0D 0) THEN ELSE IF ((K.EQ.-2.OR.K.EQ.-4).AND.UMAX-TEMP.GE.0.0D 0) THEN ELSE IOLD=J UMAX=ABS(TEMP) ENDIF 2 CONTINUE IF (UMAX.LE.EPS7) IOLD=0 RETURN END * SUBROUTINE PLTRBG ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * GRADIENT OF THE OBJECTIVE FUNCTION IS SCALED AND REDUCED. LAGRANGE * MULTIPLIERS ARE DETERMINED. TEST VALUES GMAX AND UMAX ARE COMPUTED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * II NC NUMBER OF CURRENT LINEAR CONSTRAINTS. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * II ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RU CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. VECTOR CZ(1,NF) CONTAINS LAGRANGE * MULTIPLIERS BEING DETERMINED. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RO GN(NF) TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION IF IT IS * NONZERO. * RI EPS7 TOLERANCE FOR LINEAR AND QUADRATIC PROGRAMMING. * RO GMAX NORM OF THE TRANSFORMED GRADIENT. * RO UMAX MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER. * IO IOLD INDEX OF THE REMOVED CONSTRAINT. * * SUBPROGRAMS USED : * S PLVLAG GRADIENT IS PREMULTIPLIED BY THE MATRIX WHOSE COLUMNS * ARE NORMALS OF THE ACTIVE CONSTRAINTS. * S PLTLAG COMPUTATION OF THE MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE * LAGRANGE MULTIPLIER. * S MXDRMM PREMULTIPLICATION OF A VECTOR BY A ROWWISE STORED DENSE * RECTANGULAR MATRIX. * S MXDPRB BACK SUBSTITUTION AFTER A CHOLESKI DECOMPOSITION. * RF MXVMAX L-INFINITY NORM OF A VECTOR. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE PLTRBG(NF,N,NC,IX,IC,ICA,CG,CR,CZ,G,GN,EPS7,GMAX,UMAX, & IOLD) INTEGER NF,N,NC,IX(*),IC(*),ICA(*),IOLD DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),EPS7,GMAX,UMAX DOUBLE PRECISION MXVMAX INTEGER NCA,NCZ GMAX=0.0D 0 IF (N.GT.0) THEN CALL MXDRMM(NF,N,CZ,G,GN) GMAX=MXVMAX(N,GN) ENDIF IF (NF.GT.N.AND.GMAX.LE.EPS7) THEN NCA=NF-N NCZ=N*NF CALL PLVLAG(NF,N,NC,ICA,CG,CG,G,CZ(NCZ+1)) CALL MXDPRB(NCA,CR,CZ(NCZ+1),0) CALL PLTLAG(NF,N,NC,IX,IC,ICA,CZ(NCZ+1),IC,EPS7,UMAX,IOLD) IF (UMAX.LE.EPS7) IOLD=0 CALL MXVSET(N,0.0D 0,GN) GMAX=0.0D 0 ELSE IOLD=0 UMAX=0.0D 0 ENDIF RETURN END * SUBROUTINE PLVLAG ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * GRADIENT OF THE OBJECTIVE FUNCTION IS PREMULTIPLIED BY TRANSPOSE * OF THE MATRIX WHOSE COLUMNS ARE NORMALS OF CURRENT ACTIVE CONSTRAINTS * AND GRADIENTS OF CURRENT ACTIVE FUNCTIONS. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * II NC NUMBER OF LINEARIZED CONSTRAINTS. * II IAA(NF+1) VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS. * RI AG(NF*NA) VECTOR CONTAINING SCALING PARAMETERS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RO GN(NF+1) OUTPUT VECTOR. * * SUBPROGRAMS USED : * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * SUBROUTINE PLVLAG(NF,N,NC,IAA,AG,CG,G,GN) INTEGER NF,N,NC,IAA(*) DOUBLE PRECISION AG(*),CG(*),G(*),GN(*) DOUBLE PRECISION MXVDOT INTEGER NAA,J,L NAA=NF-N DO 1 J=1,NAA L=IAA(J) IF (L.GT.NC) THEN L=L-NC GN(J)=MXVDOT(NF,AG((L-1)*NF+1),G) ELSE IF (L.GT.0) THEN GN(J)=MXVDOT(NF,CG((L-1)*NF+1),G) ELSE L=-L GN(J)=G(L) ENDIF 1 CONTINUE RETURN END * SUBROUTINE PNINT1 ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITH DIRECTIONAL * DERIVATIVES. * * PARAMETERS : * RI RL LOWER VALUE OF THE STEPSIZE PARAMETER. * RI RU UPPER VALUE OF THE STEPSIZE PARAMETER. * RI FL VALUE OF THE OBJECTIVE FUNCTION FOR R=RL. * RI FU VALUE OF THE OBJECTIVE FUNCTION FOR R=RU. * RI PL DIRECTIONAL DERIVATIVE FOR R=RL. * RI PU DIRECTIONAL DERIVATIVE FOR R=RU. * RO R VALUE OF THE STEPSIZE PARAMETER OBTAINED. * II MODE MODE OF LINE SEARCH. * II MTYP METHOD SELECTION. MTYP=1-BISECTION. MTYP=2-QUADRATIC * INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE). * MTYP=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL * DERIVATIVES). MTYP=4-CUBIC INTERPOLATION. MTYP=5-CONIC * INTERPOLATION. * IO MERR ERROR INDICATOR. MERR=0 FOR NORMAL RETURN. * * METHOD : * EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS. * SUBROUTINE PNINT1(RL,RU,FL,FU,PL,PU,R,MODE,MTYP,MERR) DOUBLE PRECISION RL, RU, FL, FU, PL, PU, R INTEGER MODE,MTYP,MERR,NTYP DOUBLE PRECISION A,B,C,D,DIS,DEN DOUBLE PRECISION C1L,C1U,C2L,C2U,C3L PARAMETER (C1L=1.1D 0,C1U=1.0D 3,C2L=1.0D-2,C2U=0.9D 0, & C3L=0.1D 0) MERR=0 IF (MODE.LE.0) RETURN IF (PL.GE.0.0D 0) THEN MERR=2 RETURN ELSE IF (RU.LE.RL) THEN MERR=3 RETURN ENDIF DO 1 NTYP=MTYP,1,-1 IF (NTYP.EQ.1) THEN C C BISECTION C IF (MODE.EQ.1) THEN R=4.0D 0*RU RETURN ELSE R=0.5D 0*(RL+RU) RETURN ENDIF ELSE IF (NTYP.EQ.MTYP) THEN A = (FU-FL)/(PL*(RU-RL)) B = PU/PL ENDIF IF (NTYP.EQ.2) THEN C C QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH ONE DIRECTIONAL C DERIVATIVE C DEN = 2.0D 0*(1.0D 0-A) ELSE IF (NTYP.EQ.3) THEN C C QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH TWO DIRECTIONAL C DERIVATIVES C DEN = 1.0D 0 - B ELSE IF (NTYP.EQ.4) THEN C C CUBIC EXTRAPOLATION OR INTERPOLATION C C = B - 2.0D 0*A + 1.0D 0 D = B - 3.0D 0*A + 2.0D 0 DIS = D*D - 3.0D 0*C IF (DIS.LT.0.0D 0) GO TO 1 DEN = D + SQRT(DIS) ELSE IF (NTYP.EQ.5) THEN C C CONIC EXTRAPOLATION OR INTERPOLATION C DIS = A*A - B IF (DIS.LT.0.0D 0) GO TO 1 DEN = A + SQRT(DIS) IF (DEN.LE.0.0D 0) GO TO 1 DEN = 1.0D 0 - B*(1.0D 0/DEN)**3 ENDIF IF (MODE.EQ.1.AND.DEN.GT.0.0D 0.AND.DEN.LT.1.0D 0) THEN C C EXTRAPOLATION ACCEPTED C R = RL + (RU-RL)/DEN R = MAX(R,C1L*RU) R = MIN(R,C1U*RU) RETURN ELSE IF (MODE.EQ.2.AND.DEN.GT.1.0D 0) THEN C C INTERPOLATION ACCEPTED C R = RL + (RU-RL)/DEN IF (RL.EQ.0.0D 0) THEN R = MAX(R,RL+C2L*(RU-RL)) ELSE R = MAX(R,RL+C3L*(RU-RL)) ENDIF R = MIN(R,RL+C2U*(RU-RL)) RETURN ENDIF 1 CONTINUE END * SUBROUTINE PNINT3 ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITHOUT DIRECTIONAL * DERIVATIVES. * * PARAMETERS : * RI RO INITIAL VALUE OF THE STEPSIZE PARAMETER. * RI RL LOWER VALUE OF THE STEPSIZE PARAMETER. * RI RU UPPER VALUE OF THE STEPSIZE PARAMETER. * RI RI INNER VALUE OF THE STEPSIZE PARAMETER. * RI FO VALUE OF THE OBJECTIVE FUNCTION FOR R=RO. * RI FL VALUE OF THE OBJECTIVE FUNCTION FOR R=RL. * RI FU VALUE OF THE OBJECTIVE FUNCTION FOR R=RU. * RI FI VALUE OF THE OBJECTIVE FUNCTION FOR R=RI. * RO PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE. * RO R VALUE OF THE STEPSIZE PARAMETER OBTAINED. * II MODE MODE OF LINE SEARCH. * II MTYP METHOD SELECTION. MTYP=1-BISECTION. MTYP=2-TWO POINT * QUADRATIC INTERPOLATION. MTYP=2-THREE POINT QUADRATIC * INTERPOLATION. * IO MERR ERROR INDICATOR. MERR=0 FOR NORMAL RETURN. * * METHOD : * EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS. * SUBROUTINE PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR) DOUBLE PRECISION ZERO,HALF,ONE,TWO,THREE,C1L,C1U,C2L,C2U,C3L PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0, + THREE=3.0D0,C1L=1.1D0,C1U=1.0D3,C2L=1.0D-2, + C2U=0.9D0,C3L=1.0D-1) DOUBLE PRECISION FI,FL,FO,FU,PO,R,RI,RL,RO,RU INTEGER MERR,MODE,MTYP DOUBLE PRECISION AI,AL,AU,DEN,DIS INTEGER NTYP LOGICAL L1,L2 MERR = 0 IF (MODE.LE.0) RETURN IF (PO.GE.ZERO) THEN MERR = 2 RETURN ELSE IF (RU.LE.RL) THEN MERR = 3 RETURN END IF L1 = RL .LE. RO L2 = RI .LE. RL DO 10 NTYP = MTYP,1,-1 IF (NTYP.EQ.1) THEN * * BISECTION * IF (MODE.EQ.1) THEN R = TWO*RU RETURN ELSE IF (RI-RL.LE.RU-RI) THEN R = HALF* (RI+RU) RETURN ELSE R = HALF* (RL+RI) RETURN END IF ELSE IF (NTYP.EQ.MTYP .AND. L1) THEN IF (.NOT.L2) AI = (FI-FO)/ (RI*PO) AU = (FU-FO)/ (RU*PO) END IF IF (L1 .AND. (NTYP.EQ.2.OR.L2)) THEN * * TWO POINT QUADRATIC EXTRAPOLATION OR INTERPOLATION * IF (AU.GE.ONE) GO TO 10 R = HALF*RU/ (ONE-AU) ELSE IF (.NOT.L1 .OR. .NOT.L2 .AND. NTYP.EQ.3) THEN * * THREE POINT QUADRATIC EXTRAPOLATION OR INTERPOLATION * AL = (FI-FL)/ (RI-RL) AU = (FU-FI)/ (RU-RI) DEN = AU - AL IF (DEN.LE.ZERO) GO TO 10 R = RI - HALF* (AU* (RI-RL)+AL* (RU-RI))/DEN ELSE IF (L1 .AND. .NOT.L2 .AND. NTYP.EQ.4) THEN * * THREE POINT CUBIC EXTRAPOLATION OR INTERPOLATION * DIS = (AI-ONE)* (RU/RI) DEN = (AU-ONE)* (RI/RU) - DIS DIS = AU + AI - DEN - TWO* (ONE+DIS) DIS = DEN*DEN - THREE*DIS IF (DIS.LT.ZERO) GO TO 10 DEN = DEN + SQRT(DIS) IF (DEN.EQ.ZERO) GO TO 10 R = (RU-RI)/DEN ELSE GO TO 10 END IF IF (MODE.EQ.1 .AND. R.GT.RU) THEN * * EXTRAPOLATION ACCEPTED * R = MAX(R,C1L*RU) R = MIN(R,C1U*RU) RETURN ELSE IF (MODE.EQ.2 .AND. R.GT.RL .AND. R.LT.RU) THEN * * INTERPOLATION ACCEPTED * IF (RI.EQ.ZERO .AND. NTYP.NE.4) THEN R = MAX(R,RL+C2L* (RU-RL)) ELSE R = MAX(R,RL+C3L* (RU-RL)) END IF R = MIN(R,RL+C2U* (RU-RL)) IF (R.EQ.RI) GO TO 10 RETURN END IF 10 CONTINUE END * SUBROUTINE PNSTEP ALL SYSTEMS 89/12/01 * PORTABILITY : ALL SYSTEMS * 89/01/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF A SCALING FACTOR FOR THE BOUNDARY STEP. * * PARAMETERS : * RI DEL MAXIMUM STEPSIZE. * RI A INPUT PARAMETER. * RI B INPUT PARAMETER. * RI C INPUT PARAMETER. * RO ALF SCALING FACTOR FOR THE BOUNDARY STEP SUCH THAT * A**2+2*B*ALF+C*ALF**2=DEL**2. * SUBROUTINE PNSTEP(DEL,A,B,C,ALF) DOUBLE PRECISION DEL, A, B, C, ALF DOUBLE PRECISION DEN, DIS DOUBLE PRECISION ZERO PARAMETER (ZERO = 0.0D 0) ALF = ZERO DEN = (DEL+A) * (DEL-A) IF (DEN .LE. ZERO) RETURN DIS = B*B + C*DEN IF (B .GE. ZERO) THEN ALF = DEN / (SQRT(DIS) + B) ELSE ALF = (SQRT(DIS) - B) / C ENDIF RETURN END * SUBROUTINE PP0AF8 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF VALUE OF THE AUGMENTED LAGRANGIAN FUNCTION. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II N DIMENSION OF THE CONSTRAINT NULL SPACE. * II NC NUMBER OF CONSTRAINTS. * RI CF(NC+1) VECTOR CONTAINING VALUES OF THE CONSTRAINTS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * II ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CL(NC) VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CU(NC) VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CZ(NC) VECTOR OF LAGRANGE MULTIPLIERS. * RI RPF PENALTY COEFFICIENT. * RO FC VALUE OF THE PENALTY TERM. * RO F VALUE OF THE PENALTY FUNCTION. * SUBROUTINE PP0AF8(NF,N,NC,CF,IC,ICA,CL,CU,CZ,RPF,FC,F) INTEGER NF,N,NC,IC(*),ICA(*) DOUBLE PRECISION CF(*),CL(*),CU(*),CZ(*),RPF,FC,F DOUBLE PRECISION POM,TEMP INTEGER J,KC FC=0.0D0 DO 1 KC=1,NC IF (IC(KC).GT.0) THEN POM=0.0D0 TEMP=CF(KC) IF (IC(KC).EQ.1.OR.IC(KC).GE.3) POM=MIN(POM,TEMP-CL(KC)) IF (IC(KC).EQ.2.OR.IC(KC).GE.3) POM=MIN(POM,CU(KC)-TEMP) FC=FC+RPF*ABS(POM) ENDIF 1 CONTINUE DO 2 J=1,NF-N KC=ICA(J) IF (KC.GT.0) THEN POM=0.0D0 TEMP=CF(KC) IF (IC(KC).EQ.1.OR.IC(KC).EQ.3.OR.IC(KC).EQ.5) & POM=MIN(POM,TEMP-CL(KC)) IF (IC(KC).EQ.2.OR.IC(KC).EQ.4.OR.IC(KC).EQ.6) & POM=MAX(POM,TEMP-CU(KC)) FC=FC-CZ(J)*POM ENDIF 2 CONTINUE F=CF(NC+1)+FC RETURN END * SUBROUTINE PPSET2 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF THE NEW PENALTY PARAMETERS. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * II ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CZ(NF) VECTOR OF LAGRANGE MULTIPLIERS. * RI CP(NC) VECTOR CONTAINING PENALTY PARAMETERS. * SUBROUTINE PPSET2(NF,N,NC,ICA,CZ,CP) INTEGER NF,N,NC,ICA(*) DOUBLE PRECISION CZ(*),CP(*) DOUBLE PRECISION TEMP INTEGER J,L,KC DO 1 KC=1,NC CP(KC)=0.5D 0*CP(KC) 1 CONTINUE DO 2 J=1,NF-N L=ICA(J) IF (L.GT.0) THEN TEMP=ABS(CZ(J)) CP(L)=MAX(TEMP,CP(L)+0.5D 0*TEMP) ENDIF 2 CONTINUE RETURN END * SUBROUTINE PS0G01 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SIMPLE SEARCH WITH TRUST REGION UPDATE. * * PARAMETERS : * RO R VALUE OF THE STEPSIZE PARAMETER. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RI FO INITIAL VALUE OF THE OBJECTIVE FUNCTION. * RI PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE. * RI PP QUADRATIC PART OF THE PREDICTED FUNCTION VALUE. * RU XDEL TRUST REGION BOUND. * RO XDELO PREVIOUS TRUST REGION BOUND. * RI XMAX MAXIMUM STEPSIZE. * RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER. * RI SNORM EUCLIDEAN NORM OF THE DIRECTION VECTOR. * RI BET1 LOWER BOUND FOR STEPSIZE REDUCTION. * RI BET2 UPPER BOUND FOR STEPSIZE REDUCTION. * RI GAM1 LOWER BOUND FOR STEPSIZE EXPANSION. * RI GAM2 UPPER BOUND FOR STEPSIZE EXPANSION. * RI EPS4 FIRST TOLERANCE FOR RATIO DF/DFPRED. STEP BOUND IS * DECREASED IF DF/DFPREDEPS5. * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * IU IDIR INDICATOR FOR DIRECTION DETERMINATION. * IDIR=0-BASIC DETERMINATION. IDIR=1-DETERMINATION * AFTER STEPSIZE REDUCTION. IDIR=2-DETERMINATION AFTER * STEPSIZE EXPANSION. * IO ITERS TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-STEP * BOUND WAS DECREASED. ITERS=2-STEP BOUND WAS UNCHANGED. * ITERS=3-STEP BOUND WAS INCREASED. ITERS=6-FIRST STEPSIZE. * II ITERD TERMINATION INDICATOR. ITERD<0-BAD DECOMPOSITION. * ITERD=0-DESCENT DIRECTION. ITERD=1-NEWTON LIKE STEP. * ITERD=2-INEXACT NEWTON LIKE STEP. ITERD=3-BOUNDARY STEP. * ITERD=4-DIRECTION WITH THE NEGATIVE CURVATURE. * ITERD=5-MARQUARDT STEP. * IO MAXST MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM * STEPSIZE WAS NOT OR WAS REACHED. * IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * II MRED MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * II KTERS TERMINATION SELECTION. KTERS=1-NORMAL TERMINATION. * KTERS=6-FIRST STEPSIZE. * II MES1 SWITCH FOR EXTRAPOLATION. MES1=1-CONSTANT INCREASING OF * THE INTERVAL. MES1=2-EXTRAPOLATION SPECIFIED BY THE PARAMETER * MES. MES1=3 SUPPRESSED EXTRAPOLATION. * II MES2 SWITCH FOR TERMINATION. MES2=1-NORMAL TERMINATION. * MES2=2-TERMINATION AFTER AT LEAST TWO STEPS (ASYMPTOTICALLY * PERFECT LINE SEARCH). * II MES3 SAFEGUARD AGAINST ROUNDING ERRORS. MES3=0-SAFEGUARD * SUPPRESSED. MES3=1-FIRST LEVEL OF SAFEGUARD. MES3=2-SECOND * LEVEL OF SAFEGUARD. * IU ISYS CONTROL PARAMETER. * * COMMON DATA : * * METHOD : * G.A.SCHULTZ, R.B.SCHNABEL, R.H.BYRD: A FAMILY OF TRUST-REGION-BASED * ALGORITHMS FOR UNCONSTRAINED MINIMIZATION WITH STRONG GLOBAL * CONVERGENCE PROPERTIES, SIAM J. NUMER.ANAL. 22 (1985) PP. 47-67. * SUBROUTINE PS0G01(R,F,FO,PO,PP,XDEL,XDELO,XMAX,RMAX,SNORM, & BET1,BET2,GAM1,GAM2,EPS4,EPS5,KD,LD,IDIR,ITERS,ITERD, & MAXST,NRED,MRED,KTERS,MES1,MES2,MES3,ISYS) INTEGER KD,LD,IDIR,ITERS,ITERD,MAXST,NRED,MRED,KTERS, & MES1,MES2,MES3,ISYS DOUBLE PRECISION R,F,FO,PO,PP,XDEL,XDELO,XMAX,RMAX,SNORM,BET1, & BET2,GAM1,GAM2,EPS4,EPS5 DOUBLE PRECISION DF,DFPRED INTEGER NRED1,NRED2 SAVE NRED1,NRED2 IF (ISYS.EQ.1) GO TO 2 C GO TO (1,2) ISYS+1 1 CONTINUE IF (IDIR.EQ.0) THEN NRED1=0 NRED2=0 ENDIF IDIR=0 XDELO=XDEL C C COMPUTATION OF THE NEW FUNCTION VALUE C R=MIN(1.0D 0,RMAX) KD= 0 LD=-1 ISYS=1 RETURN 2 CONTINUE IF(KTERS.LT.0.OR.KTERS.GT.5) THEN ITERS=6 ELSE DF=FO-F DFPRED=-R*(PO+R*PP) IF (DF.LT.EPS4*DFPRED) THEN C C STEP IS TOO LARGE, IT HAS TO BE REDUCED C IF (MES1.EQ.1) THEN XDEL=BET2*SNORM ELSE IF (MES1.EQ.2) THEN XDEL=BET2*MIN(0.5D 0*XDEL,SNORM) ELSE XDEL=0.5D 0*PO*SNORM/(PO+DF) XDEL=MAX(XDEL,BET1*SNORM) XDEL=MIN(XDEL,BET2*SNORM) ENDIF ITERS=1 IF (MES3.LE.1) THEN NRED2=NRED2+1 ELSE IF (ITERD.GT.2) NRED2=NRED2+1 ENDIF ELSE IF (DF.LE.EPS5*DFPRED) THEN C C STEP IS SUITABLE C ITERS=2 ELSE C C STEP IS TOO SMALL, IT HAS TO BE ENLARGED C IF (MES2.EQ.2) THEN XDEL=MAX(XDEL,GAM1*SNORM) ELSE IF (ITERD.GT.2) THEN XDEL=GAM1*XDEL ENDIF ITERS=3 ENDIF XDEL=MIN(XDEL,XMAX,GAM2*SNORM) IF (FO.LE.F) THEN IF (NRED1.GE.MRED) THEN ITERS=-1 ELSE IDIR=1 ITERS=0 NRED1=NRED1+1 ENDIF ENDIF ENDIF MAXST=0 IF (XDEL.GE.XMAX) MAXST=1 IF (MES3.EQ.0) THEN NRED=NRED1 ELSE NRED=NRED2 ENDIF ISYS=0 RETURN END * SUBROUTINE PS0L02 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * EXTENDED LINE SEARCH WITHOUT DIRECTIONAL DERIVATIVES. * * PARAMETERS : * RO R VALUE OF THE STEPSIZE PARAMETER. * RO RO INITIAL VALUE OF THE STEPSIZE PARAMETER. * RO RP PREVIOUS VALUE OF THE STEPSIZE PARAMETER. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RI FO INITIAL VALUE OF THE OBJECTIVE FUNCTION. * RO FP PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RI PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE. * RO PP PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE. * RI FMIN LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. * RI FMAX UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. * RI RMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER * RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER * RI TOLS TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE * CHANGE OF THE FUNCTION VALUE). * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II NIT ACTUAL NUMBER OF ITERATIONS. * II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. * IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * II MRED MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * IO MAXST MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM * STEPSIZE WAS NOT OR WAS REACHED. * II IEST LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND * IS NOT OR IS GIVEN. * II INITS CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE * IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL * STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED * INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL * STEPSIZE. * IO ITERS TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT * LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY * STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE. * ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE. * ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION. * ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL * DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION. * II KTERS TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH. * KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE. * KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE. * KTERS=6-FIRST STEPSIZE. * II MES METHOD SELECTION. MES=1-BISECTION. MES=2-QUADRATIC * INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE). * MES=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL * DERIVATIVES). MES=4-CUBIC INTERPOLATION. MES=5-CONIC * INTERPOLATION. * IU ISYS CONTROL PARAMETER. * * SUBPROGRAM USED : * S PNINT3 EXTRAPOLATION OR INTERPOLATION WITHOUT DIRECTIONAL * DERIVATIVES. * * METHOD : * SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH EXTENDED TERMINATION * CRITERIA. * SUBROUTINE PS0L02(R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX, & TOLS,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS, & MES,ISYS) INTEGER KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS, & MES,ISYS DOUBLE PRECISION R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,TOLS DOUBLE PRECISION RL,FL,RU,FU,RI,FI,RTEMP,TOL INTEGER MTYP,MERR,MODE,INIT1,MES1,MES2 LOGICAL L1,L2,L3,L4,L6,L7 PARAMETER(TOL=1.0D-4) SAVE MTYP,MODE,MES1,MES2 SAVE RL,FL,RU,FU,RI,FI IF (ISYS.EQ.1) GO TO 3 C GO TO (1,3) ISYS+1 1 CONTINUE MES1=2 MES2=2 ITERS=0 IF (PO.GE.0.0D 0) THEN R=0.0D 0 ITERS=-2 GO TO 4 ENDIF IF(RMAX.LE.0.0D 0) THEN ITERS= 0 GO TO 4 ENDIF C C INITIAL STEPSIZE SELECTION C IF (INITS.GT.0) THEN RTEMP=FMIN-F ELSE IF (IEST.EQ.0) THEN RTEMP=F-FP ELSE RTEMP=MAX(F-FP,1.0D 1*(FMIN-F)) ENDIF INIT1=ABS(INITS) RP=0.0D 0 FP=FO PP=PO IF (INIT1.EQ.0) THEN ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN R=1.0D 0 ELSE IF (INIT1.EQ.2) THEN R=MIN(1.0D 0,4.0D 0*RTEMP/PO) ELSE IF (INIT1.EQ.3) THEN R=MIN(1.0D 0, 2.0D 0*RTEMP/PO) ELSE IF (INIT1.EQ.4) THEN R=2.0D 0*RTEMP/PO ENDIF RTEMP=R R=MAX(R,RMIN) R=MIN(R,RMAX) MODE=0 RL=0.0D 0 FL=FO RU=0.0D 0 FU=FO RI=0.0D 0 FI=FO C C NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION) C 2 CALL PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR) IF (MERR.GT.0) THEN ITERS=-MERR GO TO 4 ELSE IF (MODE.EQ.1) THEN NRED=NRED-1 R=MIN(R,RMAX) ELSE IF (MODE.EQ.2) THEN NRED=NRED+1 ENDIF C C COMPUTATION OF THE NEW FUNCTION VALUE C KD= 0 LD=-1 ISYS=1 RETURN 3 CONTINUE IF (ITERS.NE.0) GO TO 4 IF (F.LE.FMIN) THEN ITERS=7 GO TO 4 ELSE L1=R.LE.RMIN.AND.NIT.NE.KIT L2=R.GE.RMAX L3=F-FO.LE.TOLS*R*PO.OR.F-FMIN.LE.(FO-FMIN)/1.0D 1 L4=F-FO.GE.(1.0D 0-TOLS)*R*PO.OR.MES2.EQ.2.AND.MODE.EQ.2 L6=RU-RL.LE.TOL*RU.AND.MODE.EQ.2 L7=MES2.LE.2.OR.MODE.NE.0 MAXST=0 IF (L2) MAXST=1 ENDIF C C TEST ON TERMINATION C IF (L1.AND..NOT.L3) THEN ITERS=0 GO TO 4 ELSE IF (L2.AND..NOT.F.GE.FU) THEN ITERS=7 GO TO 4 ELSE IF (L6) THEN ITERS=1 GO TO 4 ELSE IF (L3.AND.L7.AND.KTERS.EQ.5) THEN ITERS=5 GO TO 4 ELSE IF (L3.AND.L4.AND.L7.AND.(KTERS.EQ.2.OR.KTERS.EQ.3.OR. * KTERS.EQ.4)) THEN ITERS=2 GO TO 4 ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN ITERS=6 GO TO 4 ELSE IF(ABS(NRED).GE.MRED) THEN ITERS=-1 GO TO 4 ELSE RP=R FP=F MODE=MAX(MODE,1) MTYP=ABS(MES) IF(F.GE.FMAX) MTYP=1 ENDIF IF (MODE.EQ.1) THEN C C INTERVAL CHANGE AFTER EXTRAPOLATION C RL=RI FL=FI RI=RU FI=FU RU=R FU=F IF (F.GE.FI) THEN NRED=0 MODE=2 ELSE IF ( MES1 .EQ. 1) THEN MTYP=1 ENDIF C C INTERVAL CHANGE AFTER INTERPOLATION C ELSE IF (R.LE.RI) THEN IF (F.LE.FI) THEN RU=RI FU=FI RI=R FI=F ELSE RL=R FL=F ENDIF ELSE IF (F.LE.FI) THEN RL=RI FL=FI RI=R FI=F ELSE RU=R FU=F ENDIF ENDIF GO TO 2 4 ISYS=0 RETURN END * SUBROUTINE PS1L01 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * STANDARD LINE SEARCH WITH DIRECTIONAL DERIVATIVES. * * PARAMETERS : * RO R VALUE OF THE STEPSIZE PARAMETER. * RO RP PREVIOUS VALUE OF THE STEPSIZE PARAMETER. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RI FO INITIAL VALUE OF THE OBJECTIVE FUNCTION. * RO FP PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RO P VALUE OF THE DIRECTIONAL DERIVATIVE. * RI PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE. * RO PP PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE. * RI FMIN LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. * RI FMAX UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. * RI RMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER * RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER * RI TOLS TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE * CHANGE OF THE FUNCTION VALUE). * RI TOLP TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE * CHANGE OF THE DIRECTIONAL DERIVATIVE). * RO PAR1 PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC * UPDATES. * RO PAR2 PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC * UPDATES. * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II NIT ACTUAL NUMBER OF ITERATIONS. * II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. * IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * II MRED MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * IO MAXST MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM * STEPSIZE WAS NOT OR WAS REACHED. * II IEST LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND * IS NOT OR IS GIVEN. * II INITS CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE * IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL * STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED * INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL * STEPSIZE. * IO ITERS TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT * LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY * STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE. * ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE. * ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION. * ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL * DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION. * II KTERS TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH. * KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE. * KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE. * KTERS=6-FIRST STEPSIZE. * II MES METHOD SELECTION. MES=1-BISECTION. MES=2-QUADRATIC * INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE). * MES=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL * DERIVATIVES). MES=4-CUBIC INTERPOLATION. MES=5-CONIC * INTERPOLATION. * IU ISYS CONTROL PARAMETER. * * SUBPROGRAM USED : * S PNINT1 EXTRAPOLATION OR INTERPOLATION WITH DIRECTIONAL * DERIVATIVES. * * METHOD : * SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH STANDARD TERMINATION * CRITERIA. * SUBROUTINE PS1L01(R,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX, & TOLS,TOLP,PAR1,PAR2,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS, & ITERS,KTERS,MES,ISYS) INTEGER KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS, & MES,ISYS DOUBLE PRECISION R,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX, & TOLS,TOLP,PAR1,PAR2 DOUBLE PRECISION RL,FL,PL,RU,FU,PU,RTEMP INTEGER MTYP,MERR,MODE,INIT1,MES1,MES2,MES3 LOGICAL L1,L2,L3,L5,L7,M1,M2,M3 DOUBLE PRECISION CON,CON1 PARAMETER (CON=1.0D-2,CON1=1.0D-13) SAVE MTYP,MODE,MES1,MES2,MES3 SAVE RL,FL,PL,RU,FU,PU IF (ISYS.EQ.1) GO TO 3 C GO TO (1,3) ISYS+1 1 CONTINUE MES1=2 MES2=2 MES3=2 ITERS=0 IF (PO.GE.0.0D 0) THEN R=0.0D 0 ITERS=-2 GO TO 4 ENDIF IF(RMAX.LE.0.0D 0) THEN ITERS=0 GO TO 4 ENDIF C C INITIAL STEPSIZE SELECTION C IF (INITS.GT.0) THEN RTEMP=FMIN-F ELSE IF (IEST.EQ.0) THEN RTEMP=F-FP ELSE RTEMP=MAX(F-FP,1.0D 1*(FMIN-F)) ENDIF INIT1=ABS(INITS) RP=0.0D 0 FP=FO PP=PO IF (INIT1.EQ.0) THEN ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN R=1.0D 0 ELSE IF (INIT1.EQ.2) THEN R=MIN(1.0D 0,4.0D 0*RTEMP/PO) ELSE IF (INIT1.EQ.3) THEN R=MIN(1.0D 0, 2.0D 0*RTEMP/PO) ELSE IF (INIT1.EQ.4) THEN R=2.0D 0*RTEMP/PO ENDIF R=MAX(R,RMIN) R=MIN(R,RMAX) MODE=0 RU=0.0D 0 FU=FO PU=PO C C NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION) C 2 CALL PNINT1(RL,RU,FL,FU,PL,PU,R,MODE,MTYP,MERR) IF (MERR.GT.0) THEN ITERS=-MERR GO TO 4 ELSE IF (MODE.EQ.1) THEN NRED=NRED-1 R=MIN(R,RMAX) ELSE IF (MODE.EQ.2) THEN NRED=NRED+1 ENDIF C C COMPUTATION OF THE NEW FUNCTION VALUE AND THE NEW DIRECTIONAL C DERIVATIVE C KD= 1 LD=-1 ISYS=1 RETURN 3 CONTINUE IF (MODE.EQ.0) THEN PAR1=P/PO PAR2=F-FO ENDIF IF (ITERS.NE.0) GO TO 4 IF (F.LE.FMIN) THEN ITERS=7 GO TO 4 ELSE L1=R.LE.RMIN.AND.NIT.NE.KIT L2=R.GE.RMAX L3=F-FO.LE.TOLS*R*PO L5=P.GE.TOLP*PO.OR.MES2.EQ.2.AND.MODE.EQ.2 L7=MES2.LE.2.OR.MODE.NE.0 M1=.FALSE. M2=.FALSE. M3=L3 IF(MES3.GE.1) THEN M1=ABS(P).LE.CON*ABS(PO).AND.FO-F.GE.(CON1/CON)*ABS(FO) L3=L3.OR.M1 ENDIF IF(MES3.GE.2) THEN M2=ABS(P).LE.0.5D 0*ABS(PO).AND.ABS(FO-F).LE.2.0D 0*CON1*ABS(FO) L3=L3.OR.M2 ENDIF MAXST=0 IF (L2) MAXST=1 ENDIF C C TEST ON TERMINATION C IF (L1.AND..NOT.L3) THEN ITERS=0 GO TO 4 ELSE IF (L2.AND.L3.AND..NOT.L5) THEN ITERS=7 GO TO 4 ELSE IF (M3.AND.MES1.EQ.3) THEN ITERS=5 GO TO 4 ELSE IF (L3.AND.L5.AND.L7) THEN ITERS=4 GO TO 4 ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN ITERS=6 GO TO 4 ELSE IF (ABS(NRED).GE.MRED) THEN ITERS=-1 GO TO 4 ELSE RP=R FP=F PP=P MODE=MAX(MODE,1) MTYP=ABS(MES) IF(F.GE.FMAX) MTYP=1 ENDIF IF (MODE.EQ.1) THEN C C INTERVAL CHANGE AFTER EXTRAPOLATION C RL=RU FL=FU PL=PU RU=R FU=F PU=P IF (.NOT.L3) THEN NRED=0 MODE=2 ELSE IF ( MES1 .EQ. 1) THEN MTYP=1 ENDIF ELSE C C INTERVAL CHANGE AFTER INTERPOLATION C IF (.NOT.L3) THEN RU=R FU=F PU=P ELSE RL=R FL=F PL=P ENDIF ENDIF GO TO 2 4 ISYS=0 RETURN END * SUBROUTINE PUDBG1 ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VARIABLE METRIC UPDATE OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX * USING THE FACTORIZATION B=L*D*TRANS(L). * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * RU H(M) FACTORIZATION B=L*D*TRANS(L) OF A POSITIVE * DEFINITE APPROXIMATION OF THE HESSIAN MATRIX. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RA S(NF) AUXILIARY VECTOR. * RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. * RI GO(NF) GRADIENTS DIFFERENCE. * RI R VALUE OF THE STEPSIZE PARAMETER. * RI PO OLD VALUE OF THE DIRECTIONAL DERIVATIVE. * II NIT ACTUAL NUMBER OF ITERATIONS. * II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. * IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION. * ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS. * II MET1 SELECTION OF SELF SCALING. MET1=1-SELF SCALING SUPPRESSED. * MET1=2 INITIAL SELF SCALING. * II MEC CORRECTION IF THE NEGATIVE CURVATURE OCCURS. * MEC=1-CORRECTION SUPPRESSED. MEC=2-POWELL'S CORRECTION. * * SUBPROGRAMS USED : * S MXDPGU CORRECTION OF A DENSE SYMMETRIC POSITIVE DEFINITE * MATRIX IN THE FACTORED FORM B=L*D*TRANS(L). * S MXDPGS SCALING OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX * IN THE FACTORED FORM B=L*D*TRANS(L). * S MXVDIF DIFFERENCE OF TWO VECTORS. * RF MXVDOT DOT PRODUCT OF VECTORS. * S MXVSCL SCALING OF A VECTOR. * * METHOD : * BFGS VARIABLE METRIC METHOD. * SUBROUTINE PUDBG1(N,H,G,S,XO,GO,R,PO,NIT,KIT,ITERH,MET,MET1,MEC) DOUBLE PRECISION PO,R INTEGER ITERH,KIT,MET,MET1,MEC,N,NIT DOUBLE PRECISION G(*),GO(*),H(*),S(*),XO(*) DOUBLE PRECISION A,B,C,GAM,PAR,DEN,DIS LOGICAL L1,L3 DOUBLE PRECISION MXVDOT,MXDPGP L1 = MET1 .GE. 3 .OR. MET1 .EQ. 2 .AND. NIT .EQ. KIT L3 = .NOT. L1 * * DETERMINATION OF THE PARAMETERS B, C * B = MXVDOT(N,XO,GO) A = 0.0D0 IF (L1) THEN CALL MXVCOP(N,GO,S) CALL MXDPGB(N,H,S,1) A=MXDPGP(N,H,S,S) IF (A.LE.0.0D0) THEN ITERH=1 RETURN END IF END IF CALL MXVDIF(N,GO,G,S) CALL MXVSCL(N,R,S,S) C = -R*PO IF (C.LE.0.0D0) THEN ITERH = 3 RETURN END IF IF (MEC.GT.1) THEN IF (B.LE.1.0D-4*C) THEN * * POWELL'S CORRECTION * DIS=(1.0D0-0.1D0)*C/(C-B) CALL MXVDIF(N,GO,S,GO) CALL MXVDIR(N,DIS,GO,S,GO) B=C+DIS*(B-C) IF (L1) A=C+2.0D0*(1.0D0-DIS)*(B-C)+DIS*DIS*(A-C) ENDIF ELSE IF (B.LE.1.0D-4*C) THEN ITERH = 2 RETURN ENDIF ENDIF IF (L1) THEN * * DETERMINATION OF THE PARAMETER GAM (SELF SCALING) * IF (MET.EQ.1) THEN PAR = C/B ELSE IF (A.LE.0.0D 0) THEN PAR = C/B ELSE PAR=SQRT(C/A) ENDIF GAM = PAR IF (MET1.GT.1) THEN IF (NIT.NE.KIT) THEN L3=GAM.LT.0.5D0.OR.GAM.GT.4.0D0 ENDIF ENDIF ENDIF IF (L3) THEN GAM = 1.0D0 PAR = GAM END IF IF (MET.EQ.1) THEN * * BFGS UPDATE * CALL MXDPGU(N,H,PAR/B,GO,XO) CALL MXDPGU(N,H,-1.0D0/C,S,XO) ELSE * * HOSHINO UPDATE * DEN=PAR*B+C DIS=0.5D0*B CALL MXVDIR(N,PAR,GO,S,S) CALL MXDPGU(N,H,PAR/DIS,GO,XO) CALL MXDPGU(N,H,-1.0D0/DEN,S,XO) ENDIF ITERH = 0 IF (GAM.EQ.1.0D0) RETURN CALL MXDPGS(N,H,1.0D0/GAM) RETURN END * SUBROUTINE PUDBI1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VARIABLE METRIC UPDATES. * * PARAMETERS : * II N NUMBER OF VARIABLES. * RU H(N*(N+1)/2) UPDATED APPROXIMATION OF THE INVERSE HESSIAN * MATRIX. * RA S(N) AUXILIARY VECTOR. * RI XO(N) VECTOR OF VARIABLES DIFFERENCE. * RI GO(N) GRADIENT DIFFERENCE. * RO R VALUE OF THE STEPSIZE PARAMETER. * RI PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE. * RO PAR1 PARAMETER FOR CONTROL SCALING. * RO PAR2 PARAMETER FOR CONTROL SCALING. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RI FO INITIAL VALUE OF THE OBJECTIVE FUNCTION. * RI P CURRENT VALUE OF THE DIRECTIONAL DERIVATIVE. * II NIT NUMBER OF ITERATIONS. * II KIT INDEX OF THE ITERATION WITH THE LAST RESTART. * II MET VARIABLE METRIC UPDATE. * II MET1 SCALING STRATEGY. * II MET2 CORRECTION RULE. * IU IDECF DECOMPOSITION INDICATOR. * II ITERD TERMINATION INDICATOR. ITERD<0-BAD DECOMPOSITION. * ITERD=0-DESCENT DIRECTION. ITERD=1-NEWTON LIKE STEP. * ITERD=2-INEXACT NEWTON LIKE STEP. ITERD=3-BOUNDARY STEP. * ITERD=4-DIRECTION WITH THE NEGATIVE CURVATURE. * ITERD=5-MARQUARDT STEP. * IO ITERH UPDATE INDICATOR. ITERH=0-SUCCESSFUL UPDATE. * ITERH>0-UNSUCCESSFUL UPDATE. * * METHOD : * VARIOUS VARIABLE METRIC UPDATES INCLUDING BFGS UPDATE. * SUBROUTINE PUDBI1(N,H,S,XO,GO,R,PO,PAR1,PAR2,F,FO,P,NIT,KIT, & MET,MET1,MET2,IDECF,ITERD,ITERH) INTEGER N,NIT,KIT,MET,MET1,MET2,IDECF,ITERD,ITERH DOUBLE PRECISION H(*),S(*),XO(*),GO(*),R,PO DOUBLE PRECISION PAR1,PAR2 DOUBLE PRECISION F,FO,P DOUBLE PRECISION AA,CC DOUBLE PRECISION MXVDOT DOUBLE PRECISION DIS,POM,POM3,POM4,A,B,C,GAM,RHO,PAR DOUBLE PRECISION DEN LOGICAL L1,L2,L3 IF (MET.LE.0) GO TO 22 IF (IDECF.NE.9) THEN ITERH=-1 GO TO 22 ENDIF L1=ABS(MET1).GE.3.OR.ABS(MET1).EQ.2.AND.NIT.EQ.KIT L3=.NOT.L1 * * DETERMINATION OF THE PARAMETERS A, B, C * B=MXVDOT(N,XO,GO) IF (B.LE.0.0D 0) THEN ITERH=2 GO TO 22 ENDIF CALL MXDSMM(N,H,GO,S) A=MXVDOT(N,GO,S) IF (A.LE.0.0D 0) THEN ITERH=1 GO TO 22 ENDIF IF(MET.NE.1.OR.L1) THEN IF (ITERD.NE.1) THEN C=0.0D 0 ELSE C=-R*PO IF (C.LE.0.0D 0) THEN ITERH=3 GO TO 22 ENDIF ENDIF ELSE C=0.0D 0 ENDIF * * DETERMINATION OF THE PARAMETER RHO (NONQUADRATIC PROPERTIES) * IF (MET2.EQ.1) THEN RHO=1.0D 0 ELSE IF (FO-F+P.EQ.0) THEN RHO=1.0D 0 ELSE RHO=0.5D 0*B/(FO-F+P) ENDIF IF(RHO.LE.1.0D-2) RHO=1.0D 0 IF(RHO.GE.1.0D 2) RHO=1.0D 0 AA=A/B CC=C/B IF (L1) THEN * * DETERMINATION OF THE PARAMETER GAM (SELF SCALING) * PAR=A/B POM3=0.7D 0 POM4=6.0D 0 GAM=RHO/PAR IF (NIT.NE.KIT) THEN IF (MET1.EQ.3) THEN L2=PAR2.LE.0.0D 0 L3=L2.AND.ABS(PAR1).LE.0.2D 0 L3=L3.OR.(.NOT.L2.AND.GAM.GT.1.0D 0) L3=L3.OR.(L2.AND.PAR1.LT.0.0D 0.AND.GAM.GT.1.0D 0) L3=L3.OR.(L2.AND.PAR1.GT.0.0D 0.AND.GAM.LT.1.0D 0) L3=L3.OR.GAM.LT.POM3 L3=L3.OR.GAM.GT.POM4 ELSE IF (MET1.EQ.4) THEN L3=GAM.LT.POM3.OR.GAM.GT.POM4 ENDIF ENDIF ENDIF IF (L3) THEN GAM=1.0D 0 PAR=RHO/GAM ENDIF IF (MET.EQ.1) GO TO 17 * * NEW UPDATE * POM=1.0D 0/(AA*CC) IF (POM.LT.1.0D 0) THEN POM=MAX(1.0D-15,(SQRT(C/A)-POM)/(1.0D 0-POM)) GO TO 20 ENDIF 17 CONTINUE * * BFGS UPDATE * POM=1.0D 0 DIS=PAR+AA CALL MXVDIR(N,-DIS,XO,S,XO) DIS=1.0D 0/(B*DIS) CALL MXDSMU(N,H,DIS,XO) CALL MXDSMU(N,H,-DIS,S) GO TO 21 20 CONTINUE * * GENERAL UPDATE * DEN=PAR+POM*AA DIS=POM/DEN CALL MXDSMU(N,H,(PAR*DIS-1.0D 0)/A,S) CALL MXVDIR(N,-DIS,S,XO,S) CALL MXDSMU(N,H,DEN/B,S) 21 CONTINUE IF (GAM.EQ.1.0D 0) GO TO 22 * * SCALING * CALL MXDSMS(N,H,GAM) 22 CONTINUE ITERH=0 RETURN END * SUBROUTINE PUDBM2 ALL SYSTEMS 92/12/01 C PORTABILITY : ALL SYSTEMS C 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VARIABLE METRIC UPDATE OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX. * * PARAMETERS : * RU H(M) POSITIVE DEFINITE APPROXIMATION OF THE HESSIAN * MATRIX. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RA S(NF) AUXILIARY VECTOR. * RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. * RI GO(NF) GRADIENTS DIFFERENCE. * * COMMON DATA : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * II M NUMBER OF NONZERO ELEMENTS OF THE MATRIX. * II MET METHOD SELECTION. MET=1-BFGS UPDATE. MET=2-DFP UPDATE. * MET=3-HOSHINO UPDATE. * II MET1 SELECTION OF SELF SCALING. MET1=1-SELF SCALING SUPPRESSED. * MET1=2-INITIAL SELF SCALING. * II MET2 SELECTION OF THE LINE SEARCH MODEL. MET2=1-QUADRATIC MODEL. * MET2=2 USE OF TAYLOR EXPANSION. * II MET3 METHOD CORRECTION. MET3=1-NO CORRECTION. * MET3=2-POWELL'S CORRECTION. * II ITERD TERMINATION INDICATOR. ITERD<0-BAD DECOMPOSITION. * ITERD=0-DESCENT DIRECTION. ITERD=1-NEWTON LIKE STEP. * ITERD=2-INEXACT NEWTON LIKE STEP. ITERD=3-BOUNDARY STEP. * ITERD=4-DIRECTION WITH THE NEGATIVE CURVATURE. * ITERD=5-MARQUARDT STEP. * IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION. * ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS. * II IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=2-BUNCH-PARLETT * DECOMPOSITION. IDECF=3-INVERSION. * II ITRAN TRANSFORMATION INDICATOR. ITRAN=0 OR ITRAN=1 IF * TRANSFORMATION IS NOT OR IS USED. * II NIT ACTUAL NUMBER OF ITERATIONS. * II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. * RI R VALUE OF THE STEPSIZE PARAMETER. * RI F NEW VALUE OF THE OBJECTIVE FUNCTION. * RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. * RI P NEW VALUE OF THE DIRECTIONAL DERIVATIVE. * RI PO OLD VALUE OF THE DIRECTIONAL DERIVATIVE. * TO TUXX TEXT INFORMATION ON THE CORRECTION USED. * * SUBPROGRAMS USED : * S MXDSMM MATRIX-VECTOR PRODUCT. * S MXDSMU CORRECTION OF A DENSE SYMMETRIC MATRIX. * S MXDSMS SCALING OF A DENSE SYMMETRIC MATRIX. * S MXVDIF DIFFERENCE OF TWO VECTORS. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * RF MXVDOT DOT PRODUCT OF VECTORS. * S MXVNEG COPYING OF A VECTOR WITH THE CHANGE OF THE SIGN. * S MXVSCL SCALING OF A VECTOR. * S UOU1D1 PRINT OF ENTRY TO VARIABLE METRIC UPDATE. * S UOU1D2 PRINT OF EXIT FROM VARIABLE METRIC UPDATE. * * METHOD : * BASIC VARIABLE METRIC METHODS. * SUBROUTINE PUDBM2(NF,N,H,HH,S,XO,GO,SO,FO,PAR,MET1,MET3,IDECF, & ITERH) INTEGER NF,N,MET1,MET3,IDECF,ITERH DOUBLE PRECISION H(NF*(NF+1)/2),HH(NF*(NF+1)/2),S(NF),XO(NF), & GO(NF),SO(NF),FO,PAR DOUBLE PRECISION DEN,A,B,C,GAM,POM,MXVDOT LOGICAL L1 DOUBLE PRECISION CON PARAMETER (CON=1.0D-8) IF (IDECF.NE.0) THEN ITERH=-1 RETURN ENDIF L1=MET1.GE.2 C C DETERMINATION OF THE PARAMETERS B, C C CALL MXDSMM(N,H,XO,S) CALL MXVDIF(N,GO,S,SO) IF (MET3.EQ.2) CALL MXVSCL(N,1.0D0/SQRT(FO),SO,SO) B=MXVDOT(N,XO,SO) IF (B.LE.0.0D0) L1=.FALSE. A=0.0D0 CALL MXDSMM(N,HH,XO,S) C=MXVDOT(N,XO,S) IF (C.LE.0.0D0) L1=.FALSE. IF (L1) THEN C C DETERMINATION OF THE PARAMETER GAM (SELF SCALING) C GAM=C/B ELSE GAM=1.0D0 ENDIF PAR=GAM C C RANK ONE UPDATE C DEN=PAR*B-C IF (ABS(DEN).LE.CON*MAX(CON,ABS(PAR*B),ABS(C))) THEN IF (B.GT.0.0D0.AND.C.GT.0.0D0) GO TO 1 ITERH=4 RETURN ENDIF POM=PAR*B/DEN CALL MXVDIR(N,-PAR,SO,S,S) CALL MXDSMU(N,HH,1.0D0/DEN,S) GO TO 5 C C BFGS UPDATE C 1 POM=0.0D0 CALL MXDSMU(N,HH,PAR/B,SO) IF (C.GT.0.0D0) CALL MXDSMU(N,HH,-1.0D0/C,S) GO TO 5 5 ITERH=0 IF (GAM.NE.1.0D0) THEN CALL MXDSMS(N,HH,1.0D0/GAM) ENDIF RETURN END * SUBROUTINE PUDBQ1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * BROYDEN GOOD UPDATE OF A RECTANGULAR MATRIX AFTER THE QR * DECOMPOSITION. * * PARAMETERS : * II N NUMBER OF VARIABLES. * II NA NUMBER OF EQUATIONS. * RU H(N*(N+1)/2) UPDATED UPPER TRIANGULAR MATRIX. * RI ETA2 PARAMETER WHICH CONTROLS A NONSINGULARITY * RU AG(N*NA) UPDATED RECTANGULAR MATRIX. * RA S(N) AUXILIARY VECTOR. * RI XO(N) VECTOR OF VARIABLES DIFFERENCE. * RI AFO(NA) RIGHT HAND SIDES DIFFERENCE. * II MET VARIABLE METRIC UPDATE. * IO ITERH UPDATE INDICATOR. ITERH=0-SUCCESSFUL UPDATE. * ITERH>0-UNSUCCESSFUL UPDATE. * IU IDECA DECOMPOSITION INDICATOR. * II NDECA NUMBER OF DECOMPOSITIONS. * * METHOD : * VARIOUS VARIABLE METRIC UPDATES INCLUDING BFGS UPDATE. * SUBROUTINE PUDBQ1(N,NA,H,ETA2,AG,S,XO,AFO,MET,ITERH,IDECA, & NDECA) INTEGER N,NA,MET,INF,ITERH,IDECA,NDECA DOUBLE PRECISION H(*),ETA2,AG(*),S(*),XO(*),AFO(*) DOUBLE PRECISION DEN,MXVDOT IF (MET.LE.0) RETURN IF (IDECA.EQ.0) THEN * * QR DECOMPOSITION * DEN=ETA2 CALL MXDRQF(N,NA,AG,H) CALL MXDPRC(N,H,INF,DEN) NDECA=NDECA+1 IDECA=1 ELSE IF (IDECA.NE.1) THEN ITERH=-1 GO TO 22 ENDIF * * THE GOOD BROYDEN UPDATE * DEN=MXVDOT(N,XO,XO) IF (DEN.LE.0.0D 0) THEN ITERH=2 GO TO 22 ENDIF CALL MXVCOP(N,XO,S) 21 CONTINUE CALL MXVNEG(N,XO,XO) CALL MXDPRM(N,H,XO,1) CALL MXDRMD(N,NA,AG,XO,1.0D 0,AFO,AFO) CALL MXDRQU(N,NA,AG,H,1.0D 0/DEN,AFO,S,XO,INF) ITERH=0 22 CONTINUE RETURN END * SUBROUTINE PUDFM1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VARIABLE METRIC UPDATE OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX. * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * RU B(M) POSITIVE DEFINITE APPROXIMATION OF THE HESSIAN MATRIX. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RA S(NF) AUXILIARY VECTOR. * RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. * RI GO(NF) GRADIENTS DIFFERENCE. * RI F CURRENT VALUE OF THE OBJECTIVE FUNCTION. * RI FO PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RI ETA5 TOLERANCE FOR A HYBRID METHOD. * RI ETA9 MAXIMUM FOR REAL NUMBERS. * II IPOM1 METHOD INDICATOR. * IO IPOM2 INDICATOR FOR SCALING. * II MET METHOD SELECTION. MET=0-NO UPDATE. MET=1-BFGS UPDATE. * II MET1 SELECTION OF SELF SCALING. MET1=1-SELF SCALING SUPPRESSED. * MET1=2 SELF SCALING IN THE FIRST ITERATION AFTER RESTART. * MET1=3-SELF SCALING IN EACH ITERATION. * * COMMON DATA : * IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION. * ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS. * II IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * TO TUXX TEXT INFORMATION ON THE CORRECTION USED. * * SUBPROGRAMS USED : * S MXDSMM MATRIX-VECTOR PRODUCT. * S MXDSMU CORRECTION OF A DENSE SYMMETRIC MATRIX. * S MXDSMS SCALING OF A DENSE SYMMETRIC MATRIX. * RF MXVDOT DOT PRODUCT OF VECTORS. * S UOERR1 ERROR MESAGES. * S UOU1D1 PRINT OF ENTRY TO VARIABLE METRIC UPDATE. * S UOU1D2 PRINT OF EXIT FROM VARIABLE METRIC UPDATE. * * METHOD : * FLETCHER'S COMBINATION OF THE GAUSS-NEWTON AND THE BFGS METHODS. * SUBROUTINE PUDFM1(N,B,S,XO,GO,F,FO,ETA5,IPOM1,IPOM2,MET1,IDECF, & ITERH) INTEGER N,IPOM1,IPOM2,MET1,IDECF,ITERH DOUBLE PRECISION B(N*(N+1)/2),S(N),XO(N),GO(N),F,FO,ETA5 DOUBLE PRECISION MXVDOT DOUBLE PRECISION AB,BB,CB,GAM,PAR LOGICAL L1 IF (IDECF.NE.0) THEN ITERH=-1 RETURN ENDIF PAR=0.0D0 C C DETERMINATION OF THE PARAMETERS A,B,C C BB=MXVDOT(N,XO,GO) IF (BB.LE.0.0D0) THEN ITERH=2 IPOM1=0 RETURN ENDIF AB=0.0D0 CALL MXDSMM(N,B,XO,S) CB=MXVDOT(N,XO,S) IF (CB.LE.0.0D0) THEN ITERH=3 RETURN ENDIF L1=MET1.EQ.4.OR.MET1.EQ.3.AND.IPOM2.GE.1.OR.MET1.EQ.2 & .AND.IPOM2.EQ.1 IF (FO-F.GE.ETA5*FO) THEN IPOM1=0 ELSE IPOM1=1 ENDIF IF (L1) THEN C C DETERMINATION OF THE PARAMETER GAM (SELF SCALING) C GAM=CB/BB ELSE GAM=1.0D0 ENDIF C C BFGS UPDATE C CALL MXDSMU(N,B,GAM/BB,GO) CALL MXDSMU(N,B,-1.0D0/CB,S) ITERH=0 IPOM2=0 IF (L1) THEN CALL MXDSMS(N,B,1.0D0/GAM) ENDIF RETURN END * SUBROUTINE PUDRV1 ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DRIVER FOR HYBRID QUASI-NEWTON UPDATES. * * PARAMETERS: * RI R VALUE OF THE STEPSIZE PARAMETER. * RI FO PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RI F CURRENT VALUE OF THE OBJECTIVE FUNCTION. * RI PO PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE. * II IPOM1 UPDATE SELECTION. * II IPOM2 METHOD SELECTION. * IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * II IREST RESTART SPECIFICATION. IF IREST=0 DOES NOT HOLD THEN A * RESTART IS PERFORMED. * SUBROUTINE PUDRV1(R,FO,F,PO,IPOM1,IPOM2,NRED,IREST) INTEGER IPOM1,IPOM2,NRED,IREST DOUBLE PRECISION R,FO,F,PO DOUBLE PRECISION POM DOUBLE PRECISION CON1,CON2 PARAMETER(CON1=1.0D-1,CON2=1.0D-2) POM=(FO-F)/FO GO TO (1,2,3,4) IPOM2 1 IREST=1 IF (NRED.LE.0) THEN IPOM1=2 IREST=0 ELSE IPOM1=0 ENDIF GO TO 5 2 IREST=1 IF(POM.GE.CON2) THEN IPOM1=0 ELSE IF (F-FO.LE.R*PO) THEN IPOM1=0 ELSE IPOM1=1 IREST=0 ENDIF GO TO 5 3 IREST=1 IF (NRED.LE.0) THEN IF (IPOM1.NE.1) THEN IPOM1=2 IREST=0 ELSE IPOM1=0 ENDIF ELSE IF (POM.GE.CON2) THEN IPOM1=0 ELSE IF (IPOM1.NE.2) THEN IPOM1=1 IREST=0 ELSE IPOM1=0 ENDIF GO TO 5 4 IREST=1 IPOM1=0 5 CONTINUE RETURN END * SUBROUTINE PUDSD2 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * INITIATION OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * RU H(N*(N+1)/2) POSITIVE DEFINITE APPROXIMATION OF THE HESSIAN * MATRIX * RU B(N*(N+1)/2) POSITIVE DEFINITE APPROXIMATION OF THE HESSIAN * MATRIX * RI F CURRENT VALUE OF THE OBJECTIVE FUNCTION. * RI FO PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RI ETA5 TOLERANCE FOR A HYBRID METHOD. * II MET3 TYPE OF STRUCTURED UPDATE. MET3=1-STANDARD STRUCTURED * UPDATE. MET3=2-TOTALLY STRUCTURED UPDATE. * * COMMON DATA : * RU RAN RANDOM NUMBER. * II IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * II IREST RESTART SPECIFICATION. IF IREST=0 DOES NOT HOLD THEN A * RESTART IS PERFORMED. * II IDIR INDICATOR OF DIRECTION DETERMINATION. IDIR=0-BASIC * DETERMINATION. IDIR=1-DETERMINATION AFTER STEPSIZE * REDUCTION. IDIR=2-DETERMINATION AFTER STEPSIZE EXPANSION. * IU LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * TO TUXX TEXT INFORMATION ON THE RESTART USED. * * SUBPROGRAMS USED : * S MXDSMI GENERATION OF THE UNIT MATRIX. * S MXDSDO INITIATION OF A DIAGONAL MATRIX. * S MXDSMA DENSE SYMMETRIC MATRIX AUGMENTED * BY THE SCALED DENSE SYMMETRIC MATRIX. * S MXDSMO A SCALAR IS SET TO ALL ELEMENTS * OF A DENSE SYMMETRIC MATRIX. * S MXDSMS SCALING OF A DENSE SYMMETRIC MATRIX. * S UYSET3 DEFINITION OF THE RESTART VARIABLES. * SUBROUTINE PUDSD2(N,H,B,F,FO,ETA5,MET3,LD,IDIR,IDECF,IREST,IND) INTEGER N,MET3,LD,IDIR,IDECF,IREST,IND DOUBLE PRECISION H(N*(N+1)/2),B(N*(N+1)/2),F,FO,ETA5 INTEGER IUDSD SAVE IUDSD IND=0 IF (IREST.LT.0) THEN CALL MXDSMI(N,B) IF (F.LT.1.0D0) CALL MXDSMS(N,B,SQRT(F)) IUDSD=1 ELSE IF (IREST.EQ.0) THEN IF (IDIR.LE.0) THEN IF (FO-F.LE.ETA5*FO) THEN IF (MET3.EQ.2) THEN CALL MXDSMA(N,SQRT(F),B,H,H) ELSE CALL MXDSMA(N,1.0D0,B,H,H) ENDIF LD=MIN(LD,1) IUDSD=0 ENDIF ENDIF ELSEIF(IUDSD.EQ.0) THEN IF (MET3.EQ.2) THEN CALL MXDSMA(N,-SQRT(F),B,H,H) ELSE CALL MXDSMA(N,-1.0D0,B,H,H) ENDIF CALL MXDSMI(N,B) IF (F.LT.1.0D0) CALL MXDSMS(N,B,SQRT(F)) IUDSD=1 ELSE CALL MXDSMI(N,H) LD=MIN(LD,1) IUDSD=1 IND=1 ENDIF IDECF=0 RETURN END * SUBROUTINE PUDSD3 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * INITIATION OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * RU H(N*(N+1)/2) FACTORIZATION H=L*D*TRANS(L) OF A POSITIVE * SEMIDEFINITE HESSIAN MATRIX. * RU B(N*(N+1)/2) FACTORIZATION B=L*D*TRANS(L) OF A POSITIVE * DEFINITE APPROXIMATION OF THE HESSIAN MATRIX. * IU IPOM1 METHOD INDICATOR. * IU IPOM2 INDICATOR FOR SCALING. * * COMMON DATA : * RU RAN RANDOM NUMBER. * II IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * II IREST RESTART SPECIFICATION. IF IREST=0 DOES NOT HOLD THEN A * RESTART IS PERFORMED. * II ITERS TERMINATION INDICATOR. ITERS=0-ZERO STEP. * IU LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * * SUBPROGRAMS USED : * S MXDSMI GENERATION OF THE UNIT MATRIX. * S MXDSDO INITIATION OF A DIAGONAL MATRIX. * S UUDSMC COPYING OF A DENSE SYMMETRIC MATRIX. * RF UNRAN1 RANDOM NUMBER GENERATOR. * S UYSET3 DEFINITION OF THE RESTART VARIABLES. * SUBROUTINE PUDSD3(N,H,B,IPOM1,IPOM2,LD,IDECF,ITERS,IREST,IND) INTEGER N,IPOM1,IPOM2,LD,IDECF,ITERS,IREST,IND DOUBLE PRECISION H(N*(N+1)/2),B(N*(N+1)/2) INTEGER KDECF SAVE KDECF IND=0 IF (IREST.EQ.0.OR.IREST.GT.0.AND.IPOM1.EQ.0) THEN ELSE CALL MXDSMI(N,B) IPOM2= 1 KDECF=-1 ENDIF IF (IPOM1.EQ.1) THEN IF (ITERS.GT.0.OR.IREST.GT.0) THEN CALL MXDSMC(N,B,H) LD=MIN(LD,1) IF (IPOM2.EQ.1) IDECF=KDECF IF (IREST.GT.0) IND=1 ENDIF ELSE IF (IREST.LE.0) THEN ELSE IPOM1= 1 CALL MXDSMC(N,B,H) LD=MIN(LD,1) IF (IPOM2.EQ.1) IDECF=KDECF ENDIF ENDIF RETURN END * SUBROUTINE PYADB4 ALL SYSTEMS 98/12/01 * PORTABILITY : ALL SYSTEMS * 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * NEW LINEAR CONSTRAINTS OR NEW SIMPLE BOUNDS ARE ADDED TO THE ACTIVE * SET. GILL-MURRAY FACTORIZATION OF THE TRANSFORMED HESSIAN MATRIX * APPROXIMATION IS UPDATED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * IU N ACTUAL NUMBER OF VARIABLES. * II NC NUMBER OF LINEARIZED CONSTRAINTS. * RI X(NF) VECTOR OF VARIABLES. * IU IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. * RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. * RI CF(NC) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS. * RI CFD(NC) VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT FUNCTIONS. * IU IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * IU ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CL(NC) VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CU(NC) VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RU CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RU CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. * RU H(NF*(NF+1)/2) GILL-MURRAY FACTORIZATION OF THE TRANSFORMED * HESSIAN MATRIX APPROXIMATION. * RA S(NF) AUXILIARY VECTOR. * RI R VALUE OF THE STEPSIZE PARAMETER. * RI EPS7 TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS. * RI EPS9 TOLERANCE FOR ACTIVE CONSTRAINTS. * RO GMAX MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE. * RO UMAX MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER. * II KBF TYPE OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. KBF=1-ONE * SIDED SIMPLE BOUNDS. KBF=2-TWO SIDED SIMPLE BOUNDS. * II KBC TYPE OF CONSTRAINTS. KBC=0-NO CONSTRAINTS. KBC=1-CONSTRAINTS * WITH ONE SIDED BOUNDS. KBC=2-CONSTRAINTS WITH TWO SIDED * BOUNDS. * IU INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * IO IER ERROR INDICATOR. * IO ITERM TERMINATION INDICATOR. * * COMMON DATA : * IU NADD NUMBER OF CONSTRAINT ADDITIONS. * * SUBPROGRAMS USED : * S PLADB4 ADDITION OF A NEW ACTIVE CONSTRAINT. * S PLNEWS IDENTIFICATION OF ACTIVE UPPER BOUNDS. * S PLNEWL IDENTIFICATION OF ACTIVE LINEAR CONSTRAINRS. * S PLDIRL NEW VALUES OF CONSTRAINT FUNCTIONS. * S MXVIND CHANGE OF THE INTEGER VECTOR FOR CONSTRAINT ADDITION. * SUBROUTINE PYADB4(NF,N,NC,X,IX,XL,XU,CF,CFD,IC,ICA,CL,CU,CG, & CR,CZ,H,S,R,EPS7,EPS9,GMAX,UMAX,KBF,KBC,INEW,IER,ITERM) INTEGER NF,N,NC,IX(*),IC(*),ICA(*),KBF,KBC,INEW,IER,ITERM DOUBLE PRECISION X(*),XL(*),XU(*),CF(*),CFD(*),CL(*),CU(*), & CG(*),CR(*),CZ(*),H(*),S(*),R,EPS7,EPS9,GMAX,UMAX INTEGER I,J,K,L,IJ,IK,KC,KJ,KK,LL DOUBLE PRECISION DEN,TEMP INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH COMMON /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH IF (KBC.GT.0) THEN IF (R.NE.0.0D 0) CALL PLDIRL(NC,CF,CFD,IC,R,KBC) IF (INEW.NE.0) THEN IF (KBF.GT.0) THEN DO 1 I=1,NF INEW=0 CALL PLNEWS(X,IX,XL,XU,EPS9,I,INEW) CALL PLADB4(NF,N,ICA,CG,CR,CZ,H,S,EPS7,GMAX,UMAX,9,INEW, & NADD,IER) CALL MXVIND(IX,I,IER) IF (IER.LT.0) THEN ITERM=-15 RETURN ENDIF 1 CONTINUE ENDIF DO 2 KC=1,NC INEW=0 CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW) CALL PLADB4(NF,N,ICA,CG,CR,CZ,H,S,EPS7,GMAX,UMAX,9,INEW, & NADD,IER) CALL MXVIND(IC,KC,IER) IF (IER.LT.0) THEN ITERM=-15 RETURN ENDIF 2 CONTINUE ENDIF ELSE IF (KBF.GT.0) THEN K=0 DO 7 L=1,NF IF (IX(L).GE.0) K=K+1 INEW=0 CALL PLNEWS(X,IX,XL,XU,EPS9,L,INEW) IF (INEW.NE.0) THEN IX(L)=10-IX(L) KK=K*(K-1)/2 DEN=H(KK+K) IF (DEN.NE.0.0D 0) THEN IJ=0 KJ=KK DO 4 J=1,N IF (J.LE.K) THEN KJ=KJ+1 ELSE KJ=KJ+J-1 ENDIF IF (J.NE.K) TEMP=H(KJ)/DEN IK=KK DO 3 I=1,J IF (I.LE.K) THEN IK=IK+1 ELSE IK=IK+I-1 ENDIF IJ=IJ+1 IF (I.NE.K.AND.J.NE.K) H(IJ)=H(IJ)+TEMP*H(IK) 3 CONTINUE 4 CONTINUE ENDIF LL=KK+K DO 6 I=K+1,N DO 5 J=1,I LL=LL+1 IF (J.NE.K) THEN KK=KK+1 H(KK)=H(LL) ENDIF 5 CONTINUE 6 CONTINUE N=N-1 ENDIF 7 CONTINUE ENDIF RETURN END * SUBROUTINE PYFUT1 ALL SYSTEMS 98/12/01 C PORTABILITY : ALL SYSTEMS C 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * TERMINATION CRITERIA AND TEST ON RESTART. * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * RI F NEW VALUE OF THE OBJECTIVE FUNCTION. * RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. * RI UMAX MAXIMUN ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER. * RO GMAX NORM OF THE TRANSFORMED GRADIENT. * RI DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES. * RI TOLX LOWER BOUND FOR STEPLENGTH. * RI TOLF LOWER BOUND FOR FUNCTION DECREASE. * RI TOLB LOWER BOUND FOR FUNCTION VALUE. * RI TOLG LOWER BOUND FOR GRADIENT. * II KD DEGREE OF REQUIRED DERIVATIVES. * IU NIT ACTUAL NUMBER OF ITERATIONS. * II KIT NUMBER OF THE ITERATION AFTER RESTART. * II MIT MAXIMUM NUMBER OF ITERATIONS. * IU NFV ACTUAL NUMBER OF COMPUTED FUNCTION VALUES. * II MFV MAXIMUM NUMBER OF COMPUTED FUNCTION VALUES. * IU NFG ACTUAL NUMBER OF COMPUTED GRADIENT VALUES. * II MFG MAXIMUM NUMBER OF COMPUTED GRADIENT VALUES. * IU NTESX ACTUAL NUMBER OF TESTS ON STEPLENGTH. * II MTESX MAXIMUM NUMBER OF TESTS ON STEPLENGTH. * IU NTESF ACTUAL NUMBER OF TESTS ON FUNCTION DECREASE. * II MTESF MAXIMUM NUMBER OF TESTS ON FUNCTION DECREASE. * II ITES SYSTEM VARIBLE WHICH SPECIFIES TERMINATION. IF ITES=0 * THEN TERMINATION IS SUPPRESSED. * II IRES1 RESTART SPECIFICATION. RESTART IS PERFORMED AFTER * IRES1*N+IRES2 ITERATIONS. * II IRES2 RESTART SPECIFICATION. RESTART IS PERFORMED AFTER * IRES1*N+IRES2 ITERATIONS. * IU IREST RESTART INDICATOR. RESTART IS PERFORMED IF IREST>0. * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. * ITERS=0 FOR ZERO STEP. * IO ITERM TERMINATION INDICATOR. ITERM=1-TERMINATION AFTER MTESX * UNSUFFICIENT STEPLENGTHS. ITERM=2-TERMINATION AFTER MTESF * UNSUFFICIENT FUNCTION DECREASES. ITERM=3-TERMINATION ON LOWER * BOUND FOR FUNCTION VALUE. ITERM=4-TERMINATION ON LOWER BOUND * FOR GRADIENT. ITERM=11-TERMINATION AFTER MAXIMUM NUMBER OF * ITERATIONS. ITERM=12-TERMINATION AFTER MAXIMUM NUMBER OF * COMPUTED FUNCTION VALUES. * SUBROUTINE PYFUT1(N,F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLB,TOLG,KD, & NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,ITES,IRES1, & IRES2,IREST,ITERS,ITERM) INTEGER N,KD,NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF, & ITES,IRES1,IRES2,IREST,ITERS,ITERM DOUBLE PRECISION F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLG,TOLB DOUBLE PRECISION TEMP IF (ITERM.LT.0) RETURN IF (ITES .LE.0) GO TO 2 IF (ITERS.EQ.0) GO TO 1 IF (NIT.LE.0) FO=F+MIN(SQRT(ABS(F)),ABS(F)/1.0D 1) IF (F.LE.TOLB) THEN ITERM = 3 RETURN ENDIF IF (KD.GT.0) THEN IF (GMAX.LE.TOLG.AND.UMAX.LE.TOLG) THEN ITERM = 4 RETURN ENDIF ENDIF IF (NIT.LE.0) THEN NTESX = 0 NTESF = 0 ENDIF IF (DMAX.LE.TOLX) THEN ITERM = 1 NTESX = NTESX + 1 IF (NTESX .GE. MTESX) RETURN ELSE NTESX = 0 ENDIF TEMP=ABS(FO-F)/MAX(ABS(F),1.0D 0) IF (TEMP.LE.TOLF) THEN ITERM = 2 NTESF = NTESF+1 IF (NTESF.GE.MTESF) RETURN ELSE NTESF = 0 ENDIF 1 IF (NIT.GE.MIT) THEN ITERM = 11 RETURN ENDIF IF (NFV.GE.MFV) THEN ITERM = 12 RETURN ENDIF IF (NFG.GE.MFG) THEN ITERM = 13 RETURN ENDIF 2 ITERM = 0 IF (N.GT.0.AND.NIT-KIT.GE.IRES1*N+IRES2) THEN IREST=MAX(IREST,1) ENDIF NIT = NIT + 1 RETURN END * SUBROUTINE PYRMB1 ALL SYSTEMS 98/12/01 * PORTABILITY : ALL SYSTEMS * 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * OLD LINEAR CONSTRAINT OR AN OLD SIMPLE BOUND IS REMOVED FROM THE * ACTIVE SET. TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION AND * TRANSFORMED HESSIAN MATRIX APPROXIMATION ARE UPDATED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * IU IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * IU IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * IU ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RU CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RU CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RU GN(NF) TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION. * RU H(NF*(NF+1)/2) TRANSFORMED HESSIAN MATRIX APPROXIMATION. * RI EPS8 TOLERANCE FOR CONSTRAINT TO BE REMOVED. * RI UMAX MAXIMUN ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER. * RI GMAX NORM OF THE TRANSFORMED GRADIENT. * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. * II KBC SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR * CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO * SIDED LINEAR CONSTRAINTS. * II IOLD INDEX OF THE REMOVED CONSTRAINT. * IA KOLD AUXILIARY VARIABLE. * IA KREM AUXILIARY VARIABLE. * IO IER ERROR INDICATOR. * IO ITERM TERMINATION INDICATOR. * * COMMON DATA : * IU NREM NUMBER OF CONSTRAINT DELETIONS. * * SUBPROGRAMS USED : * S PLRMB0 CONSTRAINT DELETION. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE PYRMB1(NF,N,IX,IC,ICA,CG,CR,CZ,G,GN,H,EPS8,UMAX, & GMAX,KBF,KBC,IOLD,KOLD,KREM,IER,ITERM) INTEGER NF,N,IX(*),IC(*),ICA(*),KBF,KBC,IOLD,KOLD,KREM,IER, $ ITERM DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),H(*),EPS8,UMAX, & GMAX INTEGER I,J,K,KC,L INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH COMMON /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH IF (KBC.GT.0) THEN IF (UMAX.GT.EPS8*GMAX) THEN CALL PLRMB0(NF,N,ICA,CG,CR,CZ,G,GN,IOLD,KREM,NREM,IER) IF (IER.LT.0) THEN ITERM=-16 ELSE IF (IER.GT.0) THEN IOLD=0 ELSE K=N*(N-1)/2 CALL MXVSET(N,0.0D 0,H(K+1)) H(K+N)=1.0D 0 KC=ICA(NF-N+1) IF (KC.GT.0) THEN IC(KC)=-IC(KC) ELSE K=-KC IX(K)=-IX(K) ENDIF ENDIF ELSE IOLD=0 ENDIF ELSE IF (KBF.GT.0) THEN IF (UMAX.GT.EPS8*GMAX) THEN IX(IOLD)=MIN(ABS(IX(IOLD)),3) DO 1 I=N,KOLD,-1 GN(I+1)=GN(I) 1 CONTINUE GN(KOLD)=G(IOLD) N=N+1 K=N*(N-1)/2 L=K+N DO 3 I=N,KOLD,-1 DO 2 J=I,1,-1 IF (I.NE.KOLD.AND.J.NE.KOLD) THEN H(L)=H(K) K=K-1 L=L-1 ELSE IF (I.EQ.KOLD.AND.J.EQ.KOLD) THEN H(L)=1.0D 0 L=L-1 ELSE H(L)=0.0D 0 L=L-1 ENDIF 2 CONTINUE 3 CONTINUE ELSE IOLD=0 KOLD=0 ENDIF ENDIF RETURN END * SUBROUTINE PYTRBD ALL SYSTEMS 98/12/01 * PORTABILITY : ALL SYSTEMS * 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED * AND TRANSFORMED. TEST VALUE DMAX IS DETERMINED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * RI X(NF) VECTOR OF VARIABLES. * RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RU GO(NF) GRADIENTS DIFFERENCE. * RI CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM CURRENT * REDUCED SUBSPACE. * RU SN(NF) TRANSFORMED DIRECTION VECTOR. * RI R VALUE OF THE STEPSIZE PARAMETER. * RU F NEW VALUE OF THE OBJECTIVE FUNCTION. * RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. * RU P NEW VALUE OF THE DIRECTIONAL DERIVATIVE. * RU PO OLD VALUE OF THE DIRECTIONAL DERIVATIVE. * RO DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES. * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. * ITERS=0 FOR ZERO STEP. * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. * II KBC SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR * CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO * SIDED LINEAR CONSTRAINTS. * * SUBPROGRAMS USED : * S MXDRMM PREMULTIPLICATION OF A VECTOR BY TRANSPOSE OF A DENSE * RECTANGULAR MATRIX. * S MXVCOP COPYING OF A VECTOR. * S MXVDIF DIFFERENCE OF TWO VECTORS. * S MXVMUL DIAGONAL PREMULTIPLICATION OF A VECTOR. * S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE * SUBSTRACTED ONE. * S MXVSCL SCALING OF A VECTOR. * SUBROUTINE PYTRBD(NF,N,X,IX,XO,G,GO,CZ,SN,R,F,FO,P,PO,DMAX,ITERS, & KBF,KBC) INTEGER NF,N,IX(*),ITERS,KBF,KBC DOUBLE PRECISION X(*),XO(*),G(*),GO(*),CZ(*),SN(*),R,F,FO,P,PO, & DMAX INTEGER I,K IF (ITERS.GT.0) THEN CALL MXVDIF(NF,X,XO,XO) CALL MXVDIF(NF,G,GO,GO) PO=R*PO P=R*P ELSE F = FO P = PO CALL MXVSAV(NF,X,XO) CALL MXVSAV(NF,G,GO) ENDIF DMAX = 0.0D 0 IF (KBC.GT.0) THEN DO 1 I=1,NF DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0)) 1 CONTINUE IF (N.GT.0) THEN CALL MXVSCL(N,R,SN,XO) CALL MXVCOP(NF,GO,SN) CALL MXDRMM(NF,N,CZ,SN,GO) ENDIF ELSE IF (KBF.GT.0) THEN K=0 DO 2 I=1,NF IF (IX(I).LT.0) GO TO 2 K=K+1 DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0)) XO(K)=XO(I) GO(K)=GO(I) 2 CONTINUE ELSE DO 3 I=1,NF DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0)) 3 CONTINUE ENDIF RETURN END * SUBROUTINE PYTRBG ALL SYSTEMS 98/12/01 * PORTABILITY : ALL SYSTEMS * 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * GRADIENT OF THE OBJECTIVE FUNCTION IS SCALED AND REDUCED. * TEST VALUES GMAX AND UMAX ARE COMPUTED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * II ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RU CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RO GN(NF) TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION. * RI EPS7 TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS. * RO UMAX MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER. * RO GMAX NORM OF THE TRANSFORMED GRADIENT. * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. * II KBC SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR * CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO * SIDED LINEAR CONSTRAINTS. * II IOLD INDEX OF THE REMOVED CONSTRAINT. * IA KOLD AUXILIARY VARIABLE. * * SUBPROGRAMS USED : * S MXDRMM PREMULTIPLICATION OF A VECTOR BY A ROWWISE STORED DENSE * RECTANGULAR MATRIX. * S MXDPRB BACK SUBSTITUTION. * S MXVCOP COPYING OF A VECTOR. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * RF MXVMAX L-INFINITY NORM OF A VECTOR. * S MXVMUL DIAGONAL PREMULTIPLICATION OF A VECTOR. * SUBROUTINE PYTRBG(NF,N,IX,IC,ICA,CG,CR,CZ,G,GN,UMAX,GMAX, & KBF,KBC,IOLD,KOLD) INTEGER NF,N,IX(*),IC(*),ICA(*),KBF,KBC,IOLD,KOLD DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),UMAX,GMAX DOUBLE PRECISION TEMP,MXVMAX,MXVDOT INTEGER NCA,NCZ,I,J,K,KC IOLD=0 KOLD=0 UMAX=0.0D 0 GMAX=0.0D 0 IF (KBC.GT.0) THEN IF (NF.GT.N) THEN NCA=NF-N NCZ=N*NF CALL MXVCOP(NF,G,GN) DO 1 J=1,NCA K=ICA(J) IF (K.GT.0) THEN CZ(NCZ+J)=MXVDOT(NF,CG((K-1)*NF+1),GN) ELSE I=-K CZ(NCZ+J)=GN(I) ENDIF 1 CONTINUE CALL MXDPRB(NCA,CR,CZ(NCZ+1),0) DO 2 J=1,NCA TEMP=CZ(NCZ+J) KC=ICA(J) IF (KC.GT.0) THEN K=IC(KC) ELSE I=-KC K=IX(I) ENDIF IF (K.LE.-5) THEN ELSE IF ((K.EQ.-1.OR.K.EQ.-3).AND.UMAX+TEMP.GE.0.0D 0) THEN ELSE IF ((K.EQ.-2.OR.K.EQ.-4).AND.UMAX-TEMP.GE.0.0D 0) THEN ELSE IOLD=J UMAX=ABS(TEMP) ENDIF 2 CONTINUE ENDIF IF (N.GT.0) THEN CALL MXDRMM(NF,N,CZ,G,GN) GMAX=MXVMAX(N,GN) ENDIF ELSE IF (KBF.GT.0) THEN J=0 IOLD=0 KOLD=0 DO 3 I=1,NF TEMP=G(I) K=IX(I) IF (K.GE.0) THEN J=J+1 GN(J)=TEMP GMAX=MAX(GMAX,ABS(TEMP)) ELSE IF (K.LE.-5) THEN ELSE IF ((K.EQ.-1.OR.K.EQ.-3).AND.UMAX+TEMP.GE.0.0D 0) THEN ELSE IF ((K.EQ.-2.OR.K.EQ.-4).AND.UMAX-TEMP.GE.0.0D 0) THEN ELSE IOLD=I KOLD=J+1 UMAX=ABS(TEMP) ENDIF 3 CONTINUE N=J ELSE DO 4 I=1,NF TEMP=G(I) GMAX=MAX(GMAX,ABS(TEMP)) 4 CONTINUE N=NF ENDIF RETURN END * SUBROUTINE PYTRBH ALL SYSTEMS 98/12/01 C PORTABILITY : ALL SYSTEMS C 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * HESSIAN MATRIX OF THE OBJECTIVE FUNCTION OR ITS APPROXIMATION IS * SCALED AND REDUCED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * RI CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RI CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. * RI H(NF*(NF+1)/2) HESSIAN MATRIX OR ITS APPROXIMATION. * RA S(NF) AUXILIARY VECTOR. * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. * II KBC SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR * CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO * SIDED LINEAR CONSTRAINTS. * II LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. * * SUBPROGRAMS USED : * S MXDSMM MATRIX VECTOR PRODUCT. * S MXVCOP COPYING OF A VECTOR. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * SUBROUTINE PYTRBH(NF,N,IX,CR,CZ,H,S,KBF,KBC,LD,ITERS) INTEGER NF,N,IX(*),KBF,KBC,LD,ITERS DOUBLE PRECISION CR(*),CZ(*),H(*),S(*) DOUBLE PRECISION MXVDOT INTEGER NCA,NCR,ICZ,JCZ,I,J,K,L IF (LD.NE.2.OR.ITERS.EQ.0) RETURN IF (KBC.GT.0) THEN IF (N.LE.0) RETURN NCA=NF-N NCR=NCA*(NCA+1)/2 K=NCR JCZ=1 DO 4 J=1,N CALL MXDSMM(NF,H,CZ(JCZ),S) ICZ=1 DO 3 I=1,J K=K+1 CR(K)=MXVDOT(NF,CZ(ICZ),S) ICZ=ICZ+NF 3 CONTINUE JCZ=JCZ+NF 4 CONTINUE CALL MXVCOP(N*(N+1)/2,CR(NCR+1),H) ELSE IF (KBF.GT.0) THEN K=0 L=0 DO 2 I=1,NF DO 1 J=1,I K=K+1 IF(IX(I).LT.0.OR.IX(J).LT.0) GO TO 1 L=L+1 H(L)=H(K) 1 CONTINUE 2 CONTINUE ENDIF RETURN END * SUBROUTINE PYTRBS ALL SYSTEMS 98/12/01 * PORTABILITY : ALL SYSTEMS * 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SCALED AND REDUCED DIRECTION VECTOR IS BACK TRANSFORMED. * VECTORS X,G AND VALUES F,P ARE SAVED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * IU N ACTUAL NUMBER OF VARIABLES. * II NC NUMBER OF LINEARIZED CONSTRAINTS. * RI X(NF) VECTOR OF VARIABLES. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * RO XO(NF) SAVED VECTOR OF VARIABLES. * RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. * RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RO GO(NF) SAVED GRADIENT OF THE OBJECTIVE FUNCTION. * RI CF(NF) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCYIONS. * RO CFD(NF) VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT * FUNCTIONS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * RI CL(NC) VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CU(NC) VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. * RI SN(NF) TRANSFORMED DIRECTION VECTOR. * RO S(NF) DIRECTION VECTOR. * RO RO SAVED VALUE OF THE STEPSIZE PARAMETER. * RO FP PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RU FO SAVED VALUE OF THE OBJECTIVE FUNCTION. * RI F VALUE OF THE OBJECTIVE FUNCTION. * RO PO SAVED VALUE OF THE DIRECTIONAL DERIVATIVE. * RI P VALUE OF THE DIRECTIONAL DERIVATIVE. * RU RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER. * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. * II KBC SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR * CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO * SIDED LINEAR CONSTRAINTS. * IO KREM INDICATION OF LINEARLY DEPENDENT GRADIENTS. * IO INEW INDEX OF THE NEW ACTIVE FUNCTION. * * SUBPROGRAMS USED : * S PLMAXS DETERMINATION OF THE MAXIMUM STEPSIZE USING SIMPLE * BOUNDS. * S PLMAXL DETERMINATION OF THE MAXIMUM STEPSIZE USING LINEAR * CONSTRAINTS. * S MXDCMM MATRIX VECTOR PRODUCT. * S MXVCOP COPYING OF A VECTOR. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE PYTRBS(NF,N,NC,X,IX,XO,XL,XU,G,GO,CF,CFD,IC,CL,CU,CG, & CZ,SN,S,RO,FP,FO,F,PO,P,RMAX,KBF,KBC,KREM,INEW) INTEGER NF,N,NC,IX(*),IC(*),KBF,KBC,KREM,INEW DOUBLE PRECISION X(*),XO(*),XL(*),XU(*),G(*),GO(*),CF(*),CFD(*), & CL(*),CU(*),CG(*),CZ(*),SN(*),S(*),RO,FP,FO,F,PO,P,RMAX INTEGER I,K FP = FO RO = 0.0D 0 FO = F PO = P CALL MXVCOP(NF,X,XO) CALL MXVCOP(NF,G,GO) IF (KBC.GT.0) THEN IF (N.GT.0) THEN CALL MXDCMM(NF,N,CZ,SN,S) INEW=0 CALL PLMAXL(NF,NC,CF,CFD,IC,CL,CU,CG,S,RMAX,KBC,KREM,INEW) CALL PLMAXS(NF,X,IX,XL,XU,S,RMAX,KBF,KREM,INEW) ELSE CALL MXVSET(NF,0.0D 0,S) ENDIF ELSE IF (KBF.GT.0) THEN K=N+1 DO 1 I=NF,1,-1 IF (IX(I).LT.0) THEN S(I)=0.0D 0 ELSE K=K-1 S(I)=SN(K) ENDIF 1 CONTINUE INEW=0 CALL PLMAXS(NF,X,IX,XL,XU,S,RMAX,KBF,KREM,INEW) ENDIF RETURN END * SUBROUTINE PYTRFD ALL SYSTEMS 90/12/01 * PORTABILITY : ALL SYSTEMS * 90/12/01 LU : ORIGINAL VERSION * * PURPOSE : * PREPARATION OF VARIABLE METRIC UPDATE. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * RI X(NF) VECTOR OF VARIABLES. * RU XO(NF) SAVED VECTOR OF VARIABLES. * II IAA(NF+1) VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS. * RI AG(NF*NA) MATRIX WHOSE COLUMNS ARE GRADIENTS OF THE LINEAR * APPROXIMATED FUNCTIONS. * RI AZ(NF+1) VECTOR OF LAGRANGE MULTIPLIERS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI G(NF) GRADIENT OF THE LAGRANGIAN FUNCTION. * RU GO(NF) SAVED GRADIENT OF THE LAGRANGIAN FUNCTION. * II N ACTUAL NUMBER OF VARIABLES. * II KD DEGREE OF REQUIRED DERVATIVES. * IU LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * RU R VALUE OF THE STEPSIZE PARAMETER. * RU F VALUE OF THE OBJECTIVE FUNCTION. * RI FO SAVED VALUE OF THE OBJECTIVE FUNCTION. * RU P VALUE OF THE DIRECTIONAL DERIVATIVE. * RU PO SAVED VALUE OF THE DIRECTIONAL DERIVATIVE. * RO DMAX RELATIVE STEPSIZE. * IO ITERS TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT * LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY * STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE. * ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE. * ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION. * ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL * DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION. * * SUBPROGRAMS USED : * S MXVCOP COPYING OF A VECTOR. * S MXVDIF DIFFERENCE OF TWO VECTORS. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * S MXVSET INITIATION OF A VECTOR. * S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE * SUBSTRACTED ONE. * SUBROUTINE PYTRFD(NF,NC,X,XO,IAA,AG,AZ,CG,G,GO,N,KD,LD,R,F,FO,P, + PO,DMAX,ITERS) DOUBLE PRECISION DMAX,F,FO,P,PO,R INTEGER ITERS,KD,LD,N,NC,NF DOUBLE PRECISION AG(*),AZ(*),CG(*),G(*),GO(*),X(*),XO(*) INTEGER IAA(*) INTEGER I,J,L CALL MXVSET(NF,0.0D0,G) DO 10 J = 1,NF - N L = IAA(J) IF (L.GT.NC) THEN L = L - NC CALL MXVDIR(NF,-AZ(J),AG((L-1)*NF+1),G,G) ELSE IF (L.GT.0) THEN CALL MXVDIR(NF,-AZ(J),CG((L-1)*NF+1),G,G) ELSE L = -L G(L) = G(L) - AZ(J) END IF 10 CONTINUE IF (ITERS.GT.0) THEN CALL MXVDIF(NF,X,XO,XO) CALL MXVDIF(NF,G,GO,GO) PO = R*PO P = R*P ELSE R = 0.0D0 F = FO P = PO CALL MXVSAV(NF,X,XO) CALL MXVSAV(NF,G,GO) LD = KD END IF DMAX = 0.0D0 DO 20 I = 1,NF DMAX = MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D0)) 20 CONTINUE N = NF RETURN END * SUBROUTINE PYTRND ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD FOR MINIMAX * APPROXIMATION. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * RI X(NF) VECTOR OF VARIABLES. * RI XN(NF) VECTOR OF SCALING FACTORS. * RO XO(NF) SAVED VECTOR OF VARIABLES. * II ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RO CZ(NF) VECTOR OF LAGRANGE MULTIPLIERS. * RO CZS(NF) SAVED VECTOR OF LAGRANGE MULTIPLIERS. * RI G(NF) GRADIENT OF THE LAGRANGIAN FUNCTION. * RI GO(NF) SAVED GRADIENT OF THE LAGRANGIAN FUNCTION. * RO R VALUE OF THE STEPSIZE PARAMETER. * RO F NEW VALUE OF THE OBJECTIVE FUNCTION. * RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. * RO P NEW VALUE OF THE DIRECTIONAL DERIVATIVE. * RI PO OLD VALUE OF THE DIRECTIONAL DERIVATIVE. * RI CMAX VALUE OF THE CONSTRAINT VIOLATION. * RO CMAXO SAVED VALUE OF THE CONSTRAINT VIOLATION. * RO DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES. * * COMMON DATA : * II NORMF SCALING SPECIFICATION. NORMF=0-NO SCALING PERFORMED. * NORMF=1-SCALING FACTORS ARE DETERMINED AUTOMATICALLY. * NORMF=2-SCALING FACTORS ARE SUPPLIED BY USER. * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. * ITERS=0 FOR ZERO STEP. * * SUBPROGRAMS USED : * S MXVCOP COPYING OF A VECTOR. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * S MXVSET INITIATION OF A VECTOR. * S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE * SUBSTRACTED ONE. * SUBROUTINE PYTRND(NF,N,X,XO,ICA,CG,CZ,G,GO,R,F,FO,P,PO,CMAX,CMAXO, & DMAX,KD,LD,ITERS) INTEGER NF,N,KD,LD,ITERS INTEGER ICA(*) DOUBLE PRECISION X(*),XO(*),CG(*),CZ(*),G(*),GO(*),R,F,FO,P,PO, & CMAX,CMAXO,DMAX INTEGER I,J,L DO 2 J=1,NF-N L=ICA(J) IF (L.GT.0) THEN CALL MXVDIR(NF,-CZ(J),CG((L-1)*NF+1),G,G) ELSE L=-L G(L)=G(L)-CZ(J) ENDIF 2 CONTINUE IF (ITERS.GT.0) THEN CALL MXVDIF(NF,X,XO,XO) CALL MXVDIF(NF,G,GO,GO) PO=R*PO P=R*P ELSE F = FO P = PO CMAX = CMAXO CALL MXVSAV(NF,X,XO) CALL MXVSAV(NF,G,GO) LD = KD ENDIF DMAX = 0.0D0 DO 1 I=1,NF DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D0)) 1 CONTINUE N=NF RETURN END * SUBROUTINE PYTRUD ALL SYSTEMS 98/12/01 * PORTABILITY : ALL SYSTEMS * 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED * AND SCALED. TEST VALUE DMAX IS DETERMINED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * RI X(NF) VECTOR OF VARIABLES. * RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RU GO(NF) GRADIENTS DIFFERENCE. * RO R VALUE OF THE STEPSIZE PARAMETER. * RO F NEW VALUE OF THE OBJECTIVE FUNCTION. * RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. * RO P NEW VALUE OF THE DIRECTIONAL DERIVATIVE. * RI PO OLD VALUE OF THE DIRECTIONAL DERIVATIVE. * RO DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES. * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. * ITERS=0 FOR ZERO STEP. * * SUBPROGRAMS USED : * S PYSET1 DEGREE DEFINITION OF THE COMPUTED DERIVATIVES. * S MXVDIF DIFFERENCE OF TWO VECTORS. * S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE * SUBSTRACTED ONE. * SUBROUTINE PYTRUD(NF,X,XO,G,GO,R,F,FO,P,PO,DMAX,KD,LD,ITERS) INTEGER NF,KD,LD,ITERS DOUBLE PRECISION X(*),XO(*),G(*),GO(*),R,F,FO,P,PO,DMAX INTEGER I IF (ITERS.GT.0) THEN CALL MXVDIF(NF,X,XO,XO) CALL MXVDIF(NF,G,GO,GO) PO=R*PO P=R*P ELSE F = FO P = PO CALL MXVSAV(NF,X,XO) CALL MXVSAV(NF,G,GO) LD=KD ENDIF DMAX = 0.0D 0 DO 1 I=1,NF DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0)) 1 CONTINUE RETURN END * SUBROUTINE PYTRUF ALL SYSTEMS 98/12/01 * PORTABILITY : ALL SYSTEMS * 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VECTORS OF VARIABLES DIFFERENCE AND RIGHT HAND SIDES DIFFERENCE ARE * COMPUTED AND SCALED. TEST VALUE DMAX IS DETERMINED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II NA NUMBER OF APPROXIMATED FUNCTIONS. * RI X(NF) VECTOR OF VARIABLES. * RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. * RI AF(NA) VECTOR OF RIGHT HAND SIDES. * RI AFO(NA) VECTOR OF RIGHT HAND SIDES DIFFERENCE. * RO R VALUE OF THE STEPSIZE PARAMETER. * RO F NEW VALUE OF THE OBJECTIVE FUNCTION. * RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. * RO P NEW VALUE OF THE DIRECTIONAL DERIVATIVE. * RI PO OLD VALUE OF THE DIRECTIONAL DERIVATIVE. * RO DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES. * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. * ITERS=0 FOR ZERO STEP. * * SUBPROGRAMS USED : * S PYSET1 DEGREE DEFINITION OF THE COMPUTED DERIVATIVES. * S MXVDIF DIFFERENCE OF TWO VECTORS. * S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE * SUBSTRACTED ONE. * SUBROUTINE PYTRUF(NF,NA,X,XO,AF,AFO,R,F,FO,P,PO,DMAX,KD,LD, & ITERS) INTEGER NF,NA,KD,LD,ITERS DOUBLE PRECISION X(*),XO(*),AF(*),AFO(*),R,F,FO,P,PO,DMAX INTEGER I IF (ITERS.GT.0) THEN CALL MXVDIF(NF,X,XO,XO) CALL MXVDIF(NA,AF,AFO,AFO) PO=R*PO P=R*P ELSE R = 0.0D 0 F = FO P = PO CALL MXVSAV(NF,X,XO) CALL MXVSAV(NA,AF,AFO) LD=KD ENDIF DMAX = 0.0D 0 DO 1 I=1,NF DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0)) 1 CONTINUE RETURN END