* SUBROUTINE PA1MX2 ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION * WHICH IS DEFINED AS A MAXIMUM OF THE APPROXIMATED FUNCTIONS. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II NA NUMBER OF APPROXIMATED FUNCTIONS. * RI X(NF) VECTOR OF VARIABLES. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RA FA VALUE OF THE APPROXIMATED FUNCTION. * RO AF(NA) VECTOR WHOSE ELEMENTS ARE VALUES OF THE * APPROXIMATED FUNCTIONS. * RA GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. * RO AG(NF*NA) MATRIX WHOSE COLUMNS ARE GRADIENTS OF THE * APPROXIMATED FUNCTIONS. * RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * II KD DEGREE OF REQUIRED DERVATIVES. * IU LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II IEXT TYPE OF OBJECTIVE FUNCTION. IEXT<0-MAXIMUM OF POSITIVE * VALUES. IEXT=0-MAXIMUM OF ABSOLUTE VALUES. IEXT>0-MAXIMUM * OF NEGATIVE VALUES. * * SUBPROGRAMS USED : * SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION. * SE DER COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION. * S MXVCOP COPYING OF A VECTOR. * S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. * SUBROUTINE PA1MX2(NF,NA,X,F,FA,AF,GA,AG,G,KD,LD,IEXT) DOUBLE PRECISION F,FA INTEGER IEXT,KD,LD,NA,NF DOUBLE PRECISION AF(*),AG(*),G(*),GA(*),X(*) INTEGER NADD,NDECF,NFG,NFH,NFV,NIT,NRED,NREM,NRES DOUBLE PRECISION FTEMP,FVAL INTEGER K,KA,KAP,L,NAG,NAV COMMON /STAT/NDECF,NRES,NRED,NREM,NADD,NIT,NFV,NFG,NFH SAVE NAV,NAG,KAP IF (NIT.EQ.0) THEN NAV = 0 NAG = 0 END IF IF (KD.LE.LD) RETURN DO 20 KA = 1,NA IF (KD.LT.0) GO TO 20 IF (LD.GE.0) THEN FA = AF(KA) ELSE NAV = NAV + 1 CALL FUN(NF,KA,X,FA) AF(KA) = FA END IF IF (IEXT.EQ.0 .AND. FA.GE.0.0D0 .OR. IEXT.LT.0) THEN FTEMP = FA K = 1 ELSE FTEMP = -FA K =-1 END IF IF (KA.EQ.1) THEN FVAL = FTEMP KAP = KA L = K ELSE IF (FVAL.LT.FTEMP) THEN FVAL = FTEMP KAP = KA L = K END IF 10 IF (KD.LT.1) GO TO 20 NAG = NAG + 1 CALL DER(NF,KA,X,GA) CALL MXVCOP(NF,GA,AG((KA-1)*NF+1)) 20 CONTINUE IF (KD.GE.0 .AND. LD.LT.0) F = FVAL IF (KD.GE.1 .AND. LD.LT.1) THEN IF (L.GE.0) THEN CALL MXVCOP(NF,AG((KAP-1)*NF+1),G) ELSE CALL MXVNEG(NF,AG((KAP-1)*NF+1),G) END IF END IF NFV = NAV/NA NFG = NAG/NA LD = KD RETURN END * SUBROUTINE PF1HS1 ALL SYSTEMS 90/12/01 * PORTABILITY : ALL SYSTEMS * 90/12/01 LU : ORIGINAL VERSION * * PURPOSE : * NUMERICAL COMPUTATION OF THE HESSIAN MATRIX OF THE MODEL FUNCTION * USING ITS GRADIENTS. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * RI X(NF) VECTOR OF VARIABLES. * RO HF(M) HESSIAN MATRIX OF THE MODEL FUNCTION. * RI GF(NF) GRADIENT OF THE MODEL FUNCTION. * RA GO(NF) AUXILIARY VECTOR. * RI ETA1 PRECISION OF COMPUTED GRADIENTS. * * * SUBPROGRAMS USED : * SE FUNDER OBJECTIVE FUNCTION AND SUBGRADIENT EVALUATION. * S MXVCOP COPYING OF A VECTOR. * SUBROUTINE PF1HS1(NF,X,HF,GF,GO,ETA1) DOUBLE PRECISION ETA1 INTEGER NF DOUBLE PRECISION GF(*),GO(*),HF(*),X(*) INTEGER NADD,NDECF,NFG,NFH,NFV,NIT,NRED,NREM,NRES DOUBLE PRECISION ETA,FTEMP,XSTEP,XTEMP INTEGER IJ,IVAR,J COMMON /STAT/NDECF,NRES,NRED,NREM,NADD,NIT,NFV,NFG,NFH CALL MXVCOP(NF,GF,GO) ETA = SQRT(ETA1) DO 20 IVAR = 1,NF * * STEP SELECTION * XTEMP = X(IVAR) IF (XTEMP.GE.0.0D0) THEN XSTEP = ETA*MAX(ABS(XTEMP),1.0D0) ELSE XSTEP = -ETA*MAX(ABS(XTEMP),1.0D0) END IF X(IVAR) = XTEMP + XSTEP XSTEP = X(IVAR) - XTEMP CALL FUNDER(NF,X,FTEMP,GF) NFG = NFG + 1 * * NUMERICAL DIFFERENTIATION * DO 10 J = 1,NF IJ = MAX(IVAR,J)* (MAX(IVAR,J)-1)/2 + MIN(IVAR,J) C IJ = IND(IVAR,J) IF (J.GE.IVAR) THEN HF(IJ) = (GF(J)-GO(J))/XSTEP ELSE HF(IJ) = 0.5D0* (HF(IJ)+ (GF(J)-GO(J))/XSTEP) END IF 10 CONTINUE X(IVAR) = XTEMP 20 CONTINUE CALL MXVCOP(NF,GO,GF) RETURN END * SUBROUTINE PDDXQ1 ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD FOR MINIMAX * APPROXIMATION. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NA NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * II NC NUMBER OF LINEAR CONSTRAINTS. * 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 AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED * FUNCTIONS. * RO AFD(NA) VECTOR CONTAINING INCREMENTS OF THE APPROXIMATED * FUNCTIONS. * II IA(NA) VECTOR CONTAINING TYPES OF DEVIATIONS. * IO IAA(NF+1) VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS. * RI AG(NF*NA) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * APPROXIMATED FUNCTIONS. * RO AR((NF+1)*(NF+2)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RO AZ(NF+1) VECTOR OF LAGRANGE MULTIPLIERS. * 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. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RO G(NF+1) GRADIENT OF THE LAGRANGIAN FUNCTION. * RU H(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OR INVERSION OF THE * HESSIAN MATRIX APPROXIMATION. * RO S(NF+1) DIRECTION VECTOR. * RI F VALUE OF THE OBJECTIVE FUNCTION. * 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. * IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION. * IDECF=10-DIAGONAL MATRIX. * RI ETA0 MACHINE PRECISION. * 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. * RI TOLG TOLERANCE FOR THE GRADIENT OF THE LAGRANGIAN FUNCTION. * RO UMAX MAXIMUM LAGRANGE MULTIPLIER. * RO GMAX MAXIMUM PARTIAL DERIVATIVE. * RO GNORM NORM OF THE GRADIENT OF THE LAGRANGIAN FUNCTION. * RO SNORM NORM OF THE DIRECTION VECTOR. * RO XNORM VALUE OF LINEARIZED MINIMAX FUNCTION. * IO N DIMENSION OF MANIFOLD DEFINED BY ACTIVE CONSTRAINTS. * 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. * IO ITERD TYPE OF DIRECTION VECTOR. ITERD=1-CORRECT DIRECTION * VECTOR. ITERD=<0-FAILURE IN QUADRATIC PROGRAMMING. * IO ITERM CAUSE OF TERMINATION. * * SUBPROGRAMS USED : * S PLQDF1 DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD FOR * MINIMAX APPROXIMATION WITH LINEAR CONSTRAINTS. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * RF MXVMAX L-INFINITY NORM OF A VECTOR. * SUBROUTINE PDDXQ1(NF,NA,NC,X,IX,XL,XU,AF,AFD,IA,IAA,AG,AR,AZ,CF, + IC,CL,CU,CG,G,H,S,F,KBF,KBC,IDECF,ETA0,ETA2, + ETA9,EPS7,EPS9,TOLG,UMAX,GMAX,GNORM,SNORM,XNORM, + N,ITERQ,ITERD,ITERM) * * SPECIAL QUADRATIC PROGRAMMING SUBROUTINE * DOUBLE PRECISION EPS7,EPS9,ETA0,ETA2,ETA9,F,GMAX,GNORM,SNORM,TOLG, + UMAX,XNORM INTEGER IDECF,ITERD,ITERM,ITERQ,KBC,KBF,N,NA,NC,NF DOUBLE PRECISION AF(*),AFD(*),AG(*),AR(*),AZ(*),CF(*),CG(*),CL(*), + CU(*),G(*),H(*),S(*),X(*),XL(*),XU(*) INTEGER IA(*),IAA(*),IC(*),IX(*) INTEGER MFP DOUBLE PRECISION MXVDOT,MXVMAX MFP = 2 CALL PLQDF1(NF,NA,NC,X,IX,XL,XU,AF,AFD,IA,IAA,AG,AR,AZ,CF,IC,CL, + CU,CG,G,H,S,MFP,KBF,KBC,IDECF,ETA0,ETA2,ETA9,EPS7, + EPS9,XNORM,UMAX,GMAX,N,ITERQ) IF (ITERQ.LT.0) THEN ITERD = ITERQ - 10 RETURN END IF ITERD = 1 * * COMPUTATION OF VALUES FOR TERMINATION CRITERIA * GMAX = MXVMAX(NF,G) GNORM = SQRT(MXVDOT(NF,G,G)) SNORM = SQRT(MXVDOT(NF,S,S)) IF (GMAX.LE.1.0D-2*TOLG* (MIN(1.0D0,ABS(F)))) THEN ITERM = 4 END IF RETURN END * SUBROUTINE PDDBQ1 ALL SYSTEMS 96/12/01 * PORTABILITY : ALL SYSTEMS * 96/12/01 VL : ORIGINAL VERSION * * PURPOSE : * DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD FOR MINIMAX * APPROXIMATION AND BUNDLE UPDATING. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NA MAXIMUM BUNDLE DIMENSION. * II NC NUMBER OF LINEAR CONSTRAINTS. * 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. * RU F COMPUTED VALUE OF THE OBJECTIVE FUNCTION. * RI FO PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RU FP CURRENT MINIMUM VALUE OF THE OBJECTIVE FUNCTION. * RO FUB COMPARED VALUE OF THE OBJECTIVE FUNCTION. * RU AF(4*NA) VECTOR OF BUNDLE VALUES. * RO AFD(NA) VECTOR CONTAINING INCREMENTS OF BUNDLE FUNCTIONS. * IU IA(NA) VECTOR CONTAINING TYPES OF DEVIATIONS. * IO IAA(NF+1) VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS. * RU AG(NF*NA) MATRIX WHOSE COLUMNS ARE BUNDLE GRADIENTS. * RO AR((NF+1)*(NF+2)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RO AZ(NF+1) VECTOR OF LAGRANGE MULTIPLIERS. * 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. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI G(NF) SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RU H(NF) DIAGONAL MATRIX OF WEIGHT PARAMETERS. * RU S(NF+1) DIRECTION VECTOR. * RI XO(NF) INCREMENT VECTOR. * RU GO(NF+1) GRADIENT OF THE LAGRANGIAN FUNCTION. * RA XS(NF) AUXILIARY VECTOR. * RA GS(NF) AUXILIARY VECTOR. * RO P VALUE OF THE DIRECTIONAL DERIVATIVE. * RO R VALUE OF THE STEPSIZE PARAMETER. * RO RP VALUE OF THE STEPSIZE PARAMETER CORRESPONDING TO THE * CURRENT MINIMUM VALUE OF THE OBJECTIVE FUNCTION. * RU TO WEIGHT 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. * IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION. * IDECF=10-DIAGONAL MATRIX. * RI ETA0 MACHINE PRECISION. * 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. * RI TOLF TOLERANCE FOR CHANGE OF FUNCTION VALUES. * RI TOLG TOLERANCE FOR THE GRADIENT OF THE LAGRANGIAN FUNCTION. * RI ETA DISTANCE MEASURE PARAMETER. * RO UMAX MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER. * RO GMAX MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE. * RO GNORM NORM OF THE GRADIENT OF THE LAGRANGIAN FUNCTION. * RO SNORM NORM OF THE DIRECTION VECTOR. * RA XNORM AUXILIARY VARIABLE. * IO N DIMENSION OF MANIFOLD DEFINED BY ACTIVE CONSTRAINTS. * IU MAL CURRENT BUNDLE DIMENSION. * II NIT ACTUAL NUMBER OF ITERATIONS. * II MOS WEIGHT UPDATING METHOD SPECIFICATION. MOS=1-QUADRATIC * INTERPOLATION (MES2=1) OR LOCAL MINIMUM LOCALIZATION (MES2=2). * MOS=2-QUASI-NEWTON CONDITION. * IU NTESF ACTUAL NUMBER OF TESTS ON FUNCTION DECREASE. * IU NTESX ACTUAL NUMBER OF TESTS ON STEPLENGTH. * 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. * IO 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 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=9-SHORT STEP. ITERS=10-ZERO STEP. * ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL * DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION. * IO ITERM CAUSE OF TERMINATION. * * SUBPROGRAMS USED : * S PLQDF1 DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD FOR MINIMAX * APPROXIMATION WITH LINEAR CONSTRAINTS. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * S MXVCOP COPYING OF A VECTOR. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE PDDBQ1(NF,NA,NC,X,IX,XL,XU,F,FO,FP,FUB,AF,AFD,IA,IAA, + AG,AR,AZ,CF,IC,CL,CU,CG,G,H,S,XO,GO,XS,GS,P,R, + RP,TO,KBF,KBC,IDECF,ETA0,ETA2,ETA9,EPS7,EPS9, + TOLF,TOLG,ETA,UMAX,GMAX,GNORM,SNORM,XNORM,N,MAL, + NIT,MOS,NTESF,NTESX,ITERQ,ITERD,ITERS,ITERM) * * INITIALIZATION * DOUBLE PRECISION HALF,ONE,ZERO PARAMETER (HALF=5.0D-1,ONE=1.0D0,ZERO=0.0D0) DOUBLE PRECISION EPS7,EPS9,ETA,ETA0,ETA2,ETA9,F,FO,FP,FUB,GMAX, + GNORM,P,R,RP,SNORM,TO,TOLF,TOLG,UMAX,XNORM INTEGER IDECF,ITERD,ITERM,ITERQ,ITERS,KBC,KBF,MAL,MOS,N,NA,NC,NF, + NIT,NTESF,NTESX DOUBLE PRECISION AF(*),AFD(*),AG(*),AR(*),AZ(*),CF(*),CG(*),CL(*), + CU(*),G(*),GO(*),GS(*),H(*),S(*),X(*),XL(*), + XO(*),XS(*),XU(*) INTEGER IA(*),IAA(*),IC(*),IX(*) DOUBLE PRECISION DELF,FS,FU,PU,TOS,WK INTEGER I,IND,K,KA,L,MFP DOUBLE PRECISION MXVDOT,MXVMAX,MXVNM2 SAVE DELF,FU,PU,TOS,WK,IND IF (NIT.LE.1) THEN CALL MXVCOP(NF,G,GO) CALL MXVCOP(NF,GO,AG(1)) CALL MXVSET(NF,ZERO,S) SNORM = ZERO RP = ZERO R = ZERO FP = F FU = ZERO PU = ZERO IA(1) = 2 MAL = 0 WK = ZERO IND = 0 TO = ONE TOS = ONE IDECF = 10 END IF * MATRIX UPDATE * IF (MOS.EQ.2) THEN IF (IND.GE.1 .AND. ITERS.NE.9 .AND. NIT.GT.3) THEN WK = MXVDOT(NF,XS,XS) FS = (MXVDOT(NF,GO,XS)-MXVDOT(NF,GS,XS))/WK WK = MXVNM2(NF,GO,GS)/SQRT(WK) TOS = MIN(MAX(ABS(FS),1D-3),1D3) IF (FS.LT.WK*1D-3) TOS = MIN(MAX(WK,1D-3),1D4) IND = IND + 1 END IF TO = TOS IF (ITERS.EQ.5) THEN CALL MXVCOP(NF,XO,XS) CALL MXVCOP(NF,GO,GS) IND = 1 + MIN(1,IND) END IF END IF CALL MXVSET(NF,TO,H) * * BUNDLE VALUES UPDATE * FS = F IF (ITERS.GE.9) FS = FS - (R-RP)*P K = 1 DO 10 KA = 1,MAL AF(NA+NA+KA) = AF(NA+NA+KA) + RP*SNORM IF (IA(KA).GT.0) THEN ELSE IF (IA(KA).LT.0) THEN AFD(KA) = MXVDOT(NF,AG(K),S) END IF K = K + NF 10 CONTINUE CALL MXVDIR(MAL,RP,AFD,AF(NA+1),AF(NA+1)) PU = PU + RP*SNORM FU = FU + RP*MXVDOT(NF,AG,S) * * BUNDLE REDUCTION * DO 30 L = MAL,1,-1 IF (IA(L).NE.0) GO TO 30 K = (L-1)*NF + 1 DO 20 KA = L,MAL - 1 AF(NA+KA) = AF(NA+KA+1) AF(NA+NA+KA) = AF(NA+NA+KA+1) AF(NA*3+KA) = AF(NA*3+KA+1) IA(KA) = IA(KA+1) CALL MXVCOP(NF,AG(K+NF),AG(K)) K = K + NF 20 CONTINUE MAL = MAL - 1 30 CONTINUE IF (MAL.GE.NA) THEN K = NF + 1 DO 40 KA = 2,MAL - 1 AF(NA+KA) = AF(NA+KA+1) AF(NA+NA+KA) = AF(NA+NA+KA+1) AF(NA*3+KA) = AF(NA*3+KA+1) IA(KA) = IA(KA+1) CALL MXVCOP(NF,AG(K+NF),AG(K)) K = K + NF 40 CONTINUE MAL = MAL - 1 END IF * BUNDLE COMPLETION * MAL = MAX(2,MAL+1) K = (MAL-1)*NF + 1 AF(NA+1) = FU AF(NA+MAL) = FS AF(NA+NA+1) = PU AF(NA+NA+MAL) = (R-RP)*SNORM AF(NA*3+MAL) = SQRT(MXVDOT(NF,G,G)) DO 50 KA = 1,MAL AF(KA) = -MAX(ABS(AF(NA+KA)-FP),ETA*AF(NA+NA+KA)**2) 50 CONTINUE CALL MXVCOP(NF,G,AG(K)) IA(MAL) = 2 * * MAIN STOPPING CRITERION * F = FP IF (ITERS.LE.0) F = FO IF (MOS.EQ.1) WK = HALF* (SNORM*TO)**2 + + 5D2*MAX(ABS(FU-FP),PU*PU)/ (ABS(FP)+1D-3) IF (MOS.EQ.2) WK = HALF*GNORM*SNORM + + 5D1*MAX(ABS(FU-FP),ETA*PU*PU)/ (ABS(FP)+1D-3) IF (NIT.LE.1) WK = ONE IF (WK.LE.TOLG) THEN ITERM = 4 NIT = NIT - 1 GO TO 100 END IF * * PREPARATION FOR QUADRATIC PROGRAMMING * IF (NTESX.GT.0 .AND. NTESX.LT.NTESF) NTESF = NTESX IF (NIT.LE.1) DELF = ONE IF (FO.NE.F) DELF = ABS(FO-F)/MAX(ABS(F),ONE) MFP = 2 60 CALL PLQDF1(NF,MAL,NC,X,IX,XL,XU,AF,AFD,IA,IAA,AG,AR,AZ,CF,IC,CL, + CU,CG,GO,H,S,MFP,KBF,KBC,IDECF,ETA0,ETA2,ETA9,EPS7, + EPS9,XNORM,UMAX,GMAX,N,ITERQ) IF (ITERQ.LT.0) THEN ITERD = -2 GO TO 80 END IF ITERD = 1 GMAX = MXVMAX(NF,GO) GNORM = SQRT(MXVDOT(NF,GO,GO)) SNORM = SQRT(MXVDOT(NF,S,S)) IF (DELF.LE.TOLF .OR. ABS(FO-F)/MAX(ABS(F),ONE).GT.ONE) THEN ELSE IF (GNORM.LT.TOLG*TOLG*MXVMAX(MAL-1,AF(NA*3+2))) THEN * * REDEFINE BUNDLE FOR TOO SMALL DIRECTION VECTOR * I = MAL DO 70 KA = MAL,1,-1 IF (IA(KA).LT.0) I = KA 70 CONTINUE IF (I.EQ.MAL) GO TO 80 IA(I) = 0 GO TO 60 END IF 80 CONTINUE * * AGGREGATION * FU = ZERO PU = ZERO DO 90 KA = 1,NF - N L = IAA(KA) IF (L.GT.NC) THEN L = L - NC IF (L.GT.0) FU = FU - AZ(KA)*AF(NA+L) IF (L.GT.0) PU = PU - AZ(KA)*AF(NA+NA+L) END IF 90 CONTINUE CALL MXVCOP(NF,GO,AG(1)) * * PREPARATION FOR UYTRXD.I - TEST VALUES DETERMINATION * 100 FUB = FP RETURN END * SUBROUTINE PDDBQ2 ALL SYSTEMS 96/12/01 * PORTABILITY : ALL SYSTEMS * 96/12/01 VL : ORIGINAL VERSION * * PURPOSE : * DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD FOR MINIMAX * APPROXIMATION AND BUNDLE UPDATING * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NA MAXIMUM BUNDLE DIMENSION. * II NC NUMBER OF LINEAR CONSTRAINTS. * 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. * RU F COMPUTED VALUE OF THE OBJECTIVE FUNCTION. * RI FO PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RU FP CURRENT MINIMUM VALUE OF THE OBJECTIVE FUNCTION. * RO FUB COMPARED VALUE OF THE OBJECTIVE FUNCTION. * RU AF(5*NA) VECTOR OF BUNDLE VALUES. * RO AFD(NA) VECTOR CONTAINING INCREMENTS OF BUNDLE FUNCTIONS. * IU IA(NA) VECTOR CONTAINING TYPES OF DEVIATIONS. * IO IAA(NF+1) VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS. * RU AG(NF*NA) MATRIX WHOSE COLUMNS ARE BUNDLE GRADIENTS. * RO AR((NF+1)*(NF+2)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RO AZ(NF+1) VECTOR OF LAGRANGE MULTIPLIERS. * 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. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI G(NF) SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RI H(NF*(NF+1)/2) AGGREGATE HESSIAN MATRIX. * RI HF(NF*(NF+1)/2) HESSIAN MATRIX OF THE OBJECTIVE FUNCTION. * RU AH(NF*(NF+1)/2*NA) FIELD CONTAINING BUNDLE HESSIAN MATRICES. * RU S(NF+1) DIRECTION VECTOR. * RU GO(NF+1) GRADIENT OF THE LAGRANGIAN FUNCTION. * RO P VALUE OF THE DIRECTIONAL DERIVATIVE. * RO R VALUE OF THE STEPSIZE PARAMETER. * RO RP VALUE OF THE STEPSIZE PARAMETER CORRESPONDING TO THE * CURRENT MINIMUM VALUE OF THE OBJECTIVE FUNCTION. * II MFP TYPE OF FEASIBLE POINT. MFP=1-ARBITRARY FEASIBLE POINT. * MFP=2-OPTIMUM FEASIBLE POINT. * 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. * IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION. * IDECF=10-DIAGONAL MATRIX. * RI ETA0 MACHINE PRECISION. * 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. * RI TOLF TOLERANCE FOR CHANGE OF FUNCTION VALUES. * RI TOLG TOLERANCE FOR THE GRADIENT OF THE LAGRANGIAN FUNCTION. * RI ETA DISTANCE MEASURE PARAMETER. * RO UMAX MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER. * RO GMAX MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE. * RO GNORM NORM OF THE GRADIENT OF THE LAGRANGIAN FUNCTION. * RO SNORM NORM OF THE DIRECTION VECTOR. * RA XNORM AUXILIARY VARIABLE. * IO N DIMENSION OF MANIFOLD DEFINED BY ACTIVE CONSTRAINTS. * IU MAL CURRENT BUNDLE DIMENSION. * II NIT ACTUAL NUMBER OF ITERATIONS. * II MOS EXPONENT FOR DISTANCE MEASURE. * IU NTESF ACTUAL NUMBER OF TESTS ON FUNCTION DECREASE. * IU NTESX ACTUAL NUMBER OF TESTS ON STEPLENGTH. * 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. * IO 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 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=9-SHORT STEP. ITERS=10-ZERO STEP. * ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL * DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION. * IO ITERM CAUSE OF TERMINATION. * * SUBPROGRAMS USED : * SE FUNDER OBJECTIVE FUNCTION AND SUBGRADIENT EVALUATION. * S PLQDF1 DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD FOR MINIMAX * APPROXIMATION WITH LINEAR CONSTRAINTS. * S MXPDGF GILL-MURRAY DECOMPOSITION OF A DENSE SYMMETRIC MATRIX. * S MXPDGB BACK SUBSTITUTION AFTER GILL-MURRAY DECOMPOSITION. * S MXDPRB BACK SUBSTITUTION AFTER CHOLESKI DECOMPOSITION. * S MXDSMM MATRIX VECTOR PRODUCT. * RF MXDSMQ VALUE OF A QUADRATIC FORM WITH A DENSE SYMMETRIC MATRIX. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * S MXVCOP COPYING OF A VECTOR. * S MXVINA ABSOLUTE VALUES OF ELEMENTS OF INTEGER VECTOR. * S MXVINV CHANGE OF INTEGER VECTOR AFTER CONSTRAINT ADDITION. * RF MXVMAX L-INFINITY NORM OF A VECTOR. * S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE PDDBQ2(NF,NA,NC,X,IX,XL,XU,F,FO,FP,FUB,AF,AFD,IA,IAA, + AG,AR,AZ,CF,IC,CL,CU,CG,G,H,HF,AH,S,GO,P,R,RP, + KBF,KBC,IDECF,ETA0,ETA2,ETA9,EPS7,EPS9,TOLF, + TOLG,ETA,UMAX,GMAX,GNORM,SNORM,XNORM,N,MAL,NIT, + MOS,NTESF,NTESX,ITERQ,ITERD,ITERS,ITERM) DOUBLE PRECISION HALF,ONE,ZERO PARAMETER (HALF=5D-1,ONE=1D0,ZERO=0D0) DOUBLE PRECISION EPS7,EPS9,ETA,ETA0,ETA2,ETA9,F,FO,FP,FUB,GMAX, + GNORM,P,R,RP,SNORM,TOLF,TOLG,UMAX,XNORM INTEGER IDECF,ITERD,ITERM,ITERQ,ITERS,KBC,KBF,MAL,MOS,N,NA,NC,NF, + NIT,NTESF,NTESX DOUBLE PRECISION AF(*),AFD(*),AG(*),AH(*),AR(*),AZ(*),CF(*),CG(*), + CL(*),CU(*),G(*),GO(*),H(*),HF(*),S(*),X(*), + XL(*),XU(*) INTEGER IA(*),IAA(*),IC(*),IX(*) DOUBLE PRECISION DELF,FS,FU,FV,PP,PU,PV,WK INTEGER I,IND,IPOC,J,K,KA,L,MFP,NFF LOGICAL LA,LQ DOUBLE PRECISION MXDSMQ,MXVDOT,MXVMAX SAVE DELF,FS,FU,PU,PP,WK,IND,IPOC NFF = NF* (NF+1)/2 I = NA*4 * * INITIALIZATION * IF (NIT.LE.1) THEN CALL MXVCOP(NF,G,GO) CALL MXVCOP(NF,GO,AG) CALL MXVCOP(NFF,HF,H) CALL MXVCOP(NFF,HF,AH) CALL MXVSET(NF,ZERO,S) AF(I+1) = ZERO GNORM = ZERO SNORM = ZERO RP = ZERO R = ZERO FP = F FU = ZERO PU = ZERO IA(1) = 2 MAL = 0 IPOC = 0 IND = 0 WK = ZERO END IF * NULL STEPS COUNTER * IPOC = IPOC + 1 IF (ITERS.NE.10) IPOC = 0 LQ = IPOC .LE. 3 * * BUNDLE VALUES s, f LINEAR UPDATE * K = 1 DO 10 KA = 1,MAL AF(NA+NA+KA) = AF(NA+NA+KA) + RP*SNORM IF (IA(KA).GT.0) THEN ELSE IF (IA(KA).LT.0) THEN AFD(KA) = MXVDOT(NF,AG(K),S) END IF K = K + NF 10 CONTINUE CALL MXVDIR(MAL,RP,AFD,AF(NA+1),AF(NA+1)) PU = PU + RP*SNORM FU = FU + RP*MXVDOT(NF,AG,S) * * DAMPING PARAMETER ro (PP) DETERMINATION * FS = F PP = ONE IF (RP.NE.R) THEN FS = FS + (RP-R)*P FV = FS + HALF* (RP-R)**2*MXDSMQ(NF,HF,S,S) LA = (FV.GT.FS) .AND. (FV.GT.FP) IF (LA) PP = MAX(ZERO, (FP-FS)/ (FV-FS)) IF (LQ) FS = FS + PP* (FV-FS) ELSE IF (IND.GT.1) THEN PV = AF(NA+IND) FV = PV + HALF*R*R*MXDSMQ(NF,AH((IND-1)*NFF+1),S,S) LA = (FV.GT.PV) .AND. (FV.GT.FP) IF (LA) PP = MAX(ZERO, (FP-PV)/ (FV-PV)) END IF * BUNDLE VALUES f QUADRATIC UPDATE * PV = MXDSMQ(NF,AH,S,S) CALL MXVCOP(NFF,H,AH) DO 20 KA = 1,MAL AF(NA+KA) = AF(NA+KA) + AF(I+KA)*HALF*RP*RP* + MXDSMQ(NF,AH(1+ (KA-1)*NFF),S,S) 20 CONTINUE IF (RP.EQ.R .AND. IND.GT.1) AF(I+IND) = PP * * FU CORRECTION - COMPUTATION WITH NEW AGGREGATE MATRIX * FU = FU + AF(I+1)*HALF*RP*RP* (MXDSMQ(NF,AH,S,S)* (ONE+WK)-WK*PV) * * BUNDLE REDUCTION * DO 40 L = MAL,1,-1 IF (IA(L).NE.0) GO TO 40 K = (L-1)*NF + 1 DO 30 KA = L,MAL - 1 AF(NA+KA) = AF(NA+KA+1) AF(NA+NA+KA) = AF(NA+NA+KA+1) AF(NA*3+KA) = AF(NA*3+KA+1) AF(I+KA) = AF(I+KA+1) IA(KA) = IA(KA+1) CALL MXVCOP(NF,AG(K+NF),AG(K)) CALL MXVCOP(NFF,AH(L*NFF+1),AH(L*NFF+1-NFF)) K = K + NF 30 CONTINUE IND = IND - 1 MAL = MAL - 1 40 CONTINUE IF (MAL.GE.NA) THEN K = 1 + NF L = 1 + NFF DO 50 KA = 2,MAL - 1 AF(NA+KA) = AF(NA+KA+1) AF(NA+NA+KA) = AF(NA+NA+KA+1) AF(NA*3+KA) = AF(NA*3+KA+1) AF(I+KA) = AF(I+KA+1) IA(KA) = IA(KA+1) CALL MXVCOP(NF,AG(K+NF),AG(K)) CALL MXVCOP(NFF,AH(L+NFF),AH(L)) K = K + NF L = L + NFF 50 CONTINUE IND = IND - 1 MAL = MAL - 1 END IF * BUNDLE COMPLETION * MAL = MAX(2,MAL+1) K = (MAL-1)*NF + 1 L = (MAL-1)*NFF + 1 AF(NA+1) = FU AF(NA+MAL) = FS AF(NA+NA+1) = PU AF(NA+NA+MAL) = (R-RP)*SNORM AF(NA*3+MAL) = SQRT(MXVDOT(NF,G,G)) AF(I+1) = ONE DO 60 KA = 1,MAL AF(KA) = -MAX(ABS(AF(NA+KA)-FP),ETA*ABS(AF(NA+NA+KA))**MOS) 60 CONTINUE DO 70 KA = 0,MAL - 2 CALL MXDSMM(NF,AH(1+KA*NFF),S,AH(L)) CALL MXVDIR(NF,RP*AF(I+KA+1),AH(L),AG(1+KA*NF),AG(1+KA*NF)) 70 CONTINUE CALL MXVCOP(NF,G,AG(K)) IF (R.EQ.RP) THEN AF(I+MAL) = ONE IND = MAL ELSE AF(I+MAL) = PP IF (.NOT.LQ) AF(I+MAL) = ZERO CALL MXDSMM(NF,HF,S,AH(L)) CALL MXVDIR(NF, (RP-R)*AF(I+MAL),AH(L),AG(K),AG(K)) END IF CALL MXVCOP(NFF,HF,AH(L)) IA(MAL) = 2 * * MAIN STOPPING CRITERION - SEE ALSO UYFUT4.I * F = FP IF (ITERS.LE.0) F = FO WK = HALF*GNORM*SNORM + 5D1*MAX(ABS(FU-FP),ETA*PU*PU)/ + (ABS(FP)+1D-3) IF (NIT.LE.1) WK = ONE IF (WK.LE.TOLG) THEN ITERM = 4 NIT = NIT - 1 RETURN END IF CALL MXVCOP(NFF,H,HF) * * PREPARATION FOR QUADRATIC PROGRAMMING * IF (NTESX.GT.0 .AND. NTESX.LT.NTESF) NTESF = NTESX IF (NIT.LE.1) DELF = ONE IF (FO.NE.F) DELF = ABS(FO-F)/MAX(ABS(F),ONE) MFP = 2 80 CALL PLQDF1(NF,MAL,NC,X,IX,XL,XU,AF,AFD,IA,IAA,AG,AR,AZ,CF,IC,CL, + CU,CG,GO,HF,S,MFP,KBF,KBC,IDECF,ETA0,ETA2,ETA9,EPS7, + EPS9,XNORM,UMAX,GMAX,N,ITERQ) IF (ITERQ.LT.0) THEN ITERD = -2 GO TO 100 END IF ITERD = 1 GMAX = MXVMAX(NF,GO) GNORM = SQRT(MXVDOT(NF,GO,GO)) SNORM = SQRT(MXVDOT(NF,S,S)) IF (DELF.LE.TOLF .OR. ABS(FO-F)/MAX(ABS(F),ONE).GT.ONE) THEN ELSE IF (GNORM.LT.TOLG*TOLG*MXVMAX(MAL-1,AF(NA*3+2))) THEN * * REDEFINE BUNDLE FOR TOO SMALL DIRECTION VECTOR * J = MAL DO 90 KA = MAL,1,-1 IF (IA(KA).LT.0) J = KA 90 CONTINUE IF (J.EQ.MAL) GO TO 100 IA(J) = 0 GO TO 80 END IF 100 CONTINUE * * AGGREGATION * FU = ZERO PU = ZERO IF (NF.GT.N) CALL MXVSET(NFF,ZERO,H) K = 0 DO 110 KA = 1,NF - N L = IAA(KA) IF (L.GT.NC) THEN L = L - NC FU = FU - AZ(KA)*AF(NA+L) PU = PU - AZ(KA)*AF(NA+NA+L) CALL MXVDIR(NFF,-AZ(KA)*AF(I+L),AH((L-1)*NFF+1),H,H) IF (L.EQ.1) K = KA END IF 110 CONTINUE IDECF = 0 CALL MXVCOP(NF,GO,AG) * * CORRECTION FACTOR FOR FU (FIRST LAGRANGE MULTIPLIER, IF ACTIVE) * WK = ZERO IF (K.GT.0) WK = -AZ(K) * * PREPARATION FOR UYTRXD.I - TEST VALUES DETERMINATION * FUB = FP RETURN END * SUBROUTINE PLQDF1 ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD FOR MINIMAX * APPROXIMATION WITH LINEAR CONSTRAINTS. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NA NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * II NC NUMBER OF LINEAR CONSTRAINTS. * 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 AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED * FUNCTIONS. * RO AFD(NA) VECTOR CONTAINING INCREMENTS OF THE APPROXIMATED * FUNCTIONS. * II IA(NA) VECTOR CONTAINING TYPES OF DEVIATIONS. * IO IAA(NF+1) VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS. * RI AG(NF*NA) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * APPROXIMATED FUNCTIONS. * RO AR((NF+1)*(NF+2)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RO AZ(NF+1) VECTOR OF LAGRANGE MULTIPLIERS. * 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. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RO G(NF+1) GRADIENT OF THE LAGRANGIAN FUNCTION. * RU H(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OR INVERSION OF THE * HESSIAN MATRIX APPROXIMATION. * RO S(NF+1) 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. * IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION. * IDECF=10-DIAGONAL MATRIX. * RI ETA0 MACHINE PRECISION. * 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. * RO XNORM VALUE OF LINEARIZED MINIMAX FUNCTION. * RO UMAX MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER. * RO GMAX MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE. * IO N DIMENSION OF MANIFOLD DEFINED BY ACTIVE CONSTRAINTS. * 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 PLMINA DETERMINATION OF THE NEW ACTIVE FUNCTION. * S PLMINL DETERMINATION OF THE NEW ACTIVE LINEAR CONSTRAINT. * S PLMINS DETERMINATION OF THE NEW ACTIVE SIMPLE BOUND. * S PLMINT DETERMINATION OF THE NEW ACTIVE TRUST REGION BOUND. * S PLADF1 CONSTRAINT ADDITION. * S PLRMF0 CONSTRAINT DELETION. * S MXDPGF GILL-MURRAY DECOMPOSITION OF A DENSE SYMMETRIC MATRIX. * S MXDPGB BACK SUBSTITUTION AFTER GILL-MURRAY DECOMPOSITION. * S MXDPRB BACK SUBSTITUTION AFTER CHOLESKI DECOMPOSITION. * S MXDSMM MATRIX VECTOR PRODUCT. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * S MXVCOP COPYING OF A VECTOR. * S MXVINA ABSOLUTE VALUES OF ELEMENTS OF INTEGER VECTOR. * S MXVINV CHANGE OF INTEGER VECTOR AFTER CONSTRAINT ADDITION. * S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. * S MXVSET INITIATION OF A VECTOR. * * L.LUKSAN: DUAL METHOD FOR SOLVING A SPECIAL PROBLEM OF QUADRATIC * PROGRAMMING AS A SUBPROBLEM AT LINEARLY CONSTRAINED NONLINEAR MINIMAX * APPROXIMATION. kYBERNETIKA 20 (1984) 445-457. * SUBROUTINE PLQDF1(NF,NA,NC,X,IX,XL,XU,AF,AFD,IA,IAA,AG,AR,AZ,CF, + IC,CL,CU,CG,G,H,S,MFP,KBF,KBC,IDECF,ETA0,ETA2, + ETA9,EPS7,EPS9,XNORM,UMAX,GMAX,N,ITERQ) DOUBLE PRECISION EPS7,EPS9,ETA0,ETA2,ETA9,GMAX,UMAX,XNORM INTEGER IDECF,ITERQ,KBC,KBF,MFP,N,NA,NC,NF DOUBLE PRECISION AF(*),AFD(*),AG(*),AR(*),AZ(*),CF(*),CG(*),CL(*), + CU(*),G(*),H(*),S(*),X(*),XL(*),XU(*) INTEGER IA(*),IAA(*),IC(*),IX(*) INTEGER NADD,NDECF,NFG,NFH,NFV,NIT,NRED,NREM,NRES DOUBLE PRECISION BET,CON,E,GAM,PAR,Q,QO,SNORM,STEP,STEP1,STEP2,T, + TEMP INTEGER I,IER,INEW,INF,IOLD,J,K,KA,KC,KNEW,KOLD,KREM,L,NAA,NAR DOUBLE PRECISION MXVDOT,MXVMAX COMMON /STAT/NDECF,NRES,NRED,NREM,NADD,NIT,NFV,NFG,NFH T = 1.0D0 CON = ETA9 IF (IDECF.LT.0) IDECF = 1 IF (IDECF.EQ.0) THEN * * GILL-MURRAY DECOMPOSITION * TEMP = ETA2 CALL MXDPGF(NF,H,INF,TEMP,STEP) NDECF = NDECF + 1 IDECF = 1 END IF IF (IDECF.GE.2 .AND. IDECF.LE.8) THEN ITERQ = -10 RETURN END IF * * INITIATION * NRED = 0 IF (MFP.EQ.3) GO TO 10 N = NF NAA = 0 NAR = 0 XNORM = -ETA9 Q = 2.0D0 * XNORM CALL MXVINA(NA,IA) IF (KBF.GT.0) CALL MXVINA(NF,IX) IF (KBC.GT.0) CALL MXVINA(NC,IC) * * DIRECTION DETERMINATION * 10 CALL MXVSET(NF,0.0D0,S) DO 20 J = 1,NAA L = IAA(J) IF (L.GT.NC) THEN L = L - NC CALL MXVDIR(NF,AZ(J),AG((L-1)*NF+1),S,S) ELSE IF (L.GT.0) THEN CALL MXVDIR(NF,AZ(J),CG((L-1)*NF+1),S,S) ELSE L = -L S(L) = S(L) + AZ(J) END IF 20 CONTINUE CALL MXVCOP(NF,S,G) IF (NAA.GT.0) THEN IF (IDECF.EQ.1) THEN CALL MXDPGB(NF,H,S,0) ELSE IF (IDECF.EQ.9) THEN CALL MXDSMM(NF,H,G,S) ELSE CALL MXVMUL(NF,H,S,S,-1) END IF END IF * * INITIAL MINIMAX VARIABLE * IF (NAA.EQ.1) THEN TEMP = AF(INEW-NC) + MXVDOT(NF,AG((INEW-NC-1)*NF+1),S) XNORM = -SIGN(1,KNEW)*TEMP END IF * * CHECK OF FEASIBILITY * INEW = 0 PAR = 0.0D0 CALL PLMINA(NF,NA,NC,AF,AFD,IA,AG,S,INEW,KNEW,EPS9,XNORM,PAR) IF (NAA.GT.0) THEN CALL PLMINL(NF,NC,CF,IC,CL,CU,CG,S,KBC,INEW,KNEW,EPS9,PAR) CALL PLMINS(NF,IX,X,XL,XU,S,KBF,INEW,KNEW,EPS9,PAR) END IF IF (INEW.EQ.0) THEN * * SOLUTION ACHIEVED * CALL MXVNEG(NF,G,G) ITERQ = 2 RETURN ELSE QO = Q Q=0.5D0 * MXVDOT(NF,G,S) + XNORM IF (Q.LE.QO) THEN CALL MXVNEG(NF,G,G) ITERQ = 3 RETURN END IF SNORM = 0.0D0 END IF 30 IER = 0 * * STEPSIZE DETERMINATION * CALL PLADF1(NF,NC,IA,IAA,AG,AR,CG,H,S,G,IDECF,N,INEW,KNEW,IER, + EPS7,GMAX,UMAX,E,T) CALL PLDLAG(NF,NC,IA,IAA,S,N,KOLD) IF (KOLD.EQ.0) THEN * * ZERO STEPSIZE * STEP1 = 0.0D0 STEP = STEP1 SNORM = SIGN(1,KNEW) XNORM = XNORM - PAR ELSE * * PRIMAL STEPSIZE * CALL MXDPRB(NAA,AR,S,1) BET = E - MXVDOT(NAA,S,G) GAM = BET/MXVDOT(NAA,S,S) UMAX = BET*GAM + UMAX IF (UMAX.LE.EPS7*GMAX) THEN STEP1 = CON ELSE STEP1 = -PAR/UMAX END IF * * DUAL STEPSIZE * CALL MXDPRB(NAA,AR,S,-1) CALL MXDPRB(NAA,AR,G,-1) CALL MXVDIR(NAA,GAM,S,G,G) IF (KNEW.LT.0) CALL MXVNEG(NAA,G,G) STEP = MXVMAX(NAA,G) IOLD = 0 STEP2 = CON DO 40 J = 1,NAA 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) END IF 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 IF (ABS(G(J)).LE.ETA0*STEP) THEN ELSE TEMP = AZ(J)/G(J) IF (STEP2.GT.TEMP) THEN IOLD = J STEP2 = TEMP END IF END IF 40 CONTINUE * * FINAL STEPSIZE * STEP = MIN(STEP1,STEP2) IF (STEP.GE.CON) THEN * * FEASIBLE SOLUTION DOES NOT EXIST * ITERQ = -1 RETURN END IF * NEW LAGRANGE MULTIPLIERS * CALL MXVDIR(NAA,-STEP,G,AZ,AZ) SNORM = SNORM + SIGN(1,KNEW)*STEP XNORM = XNORM + SIGN(1,KNEW)*STEP*GAM PAR = PAR - (STEP/STEP1)*PAR END IF IF (STEP.EQ.STEP1) THEN IF (N.LT.0) THEN * * IMPOSSIBLE SITUATION * ITERQ = -5 RETURN END IF * * CONSTRAINT ADDITION * IF (IER.EQ.0) THEN N = N - 1 NAA = NAA + 1 NAR = NAR + NAA AZ(NAA) = SNORM END IF IF (INEW.GT.NC) THEN KA = INEW - NC CALL MXVINV(IA,KA,KNEW) ELSE 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 END IF NRED = NRED + 1 NADD = NADD + 1 GO TO 10 ELSE * * CONSTRAINT DELETION * DO 50 J = IOLD,NAA - 1 AZ(J) = AZ(J+1) 50 CONTINUE CALL PLRMF0(NF,NC,IX,IA,IAA,AR,IC,G,N,IOLD,KREM,IER) NAR = NAR - NAA NAA = NAA - 1 CALL MXVINA(NA,IA) IF (KBC.GT.0) CALL MXVINA(NC,IC) IF (KBF.GT.0) CALL MXVINA(NF,IX) DO 60 J = 1,NAA L = IAA(J) 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 60 CONTINUE GO TO 30 END IF END * SUBROUTINE PLMINA ALL SYSTEMS 90/12/01 * PORTABILITY : ALL SYSTEMS * 90/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF THE NEW ACTIVE FUNCTION. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II NA NUMBER OF CURRENT LINEAR APPROXIMATED FUNCTIONS. * II NC NUMBER OF CURRENT LINEAR CONSTRAINTS. * RI AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED * FUNCTIONS. * RO AFD(NA) VECTOR CONTAINING INCREMENTS OF THE APPROXIMATED * FUNCTIONS. * II IA(NA) VECTOR CONTAINING TYPES OF DEVIATIONS. * RI AG(NF*NA) VECTOR CONTAINING SCALING PARAMETERS. * RI S(NF) DIRECTION VECTOR. * IO INEW INDEX OF THE NEW ACTIVE FUNCTION. * IO KNEW SIGNUM OF THE NEW ACTIVE GRADIENT. * RI EPS9 TOLERANCE FOR ACTIVE FUNCTIONS. * RO XNORM VALUE OF LINEARIZED MINIMAX FUNCTION. * RA PAR AUXILIARY VARIABLE. * * SUBPROGRAMS USED : * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * SUBROUTINE PLMINA(NF,NA,NC,AF,AFD,IA,AG,S,INEW,KNEW,EPS9,XNORM, + PAR) DOUBLE PRECISION EPS9,PAR,XNORM INTEGER INEW,KNEW,NA,NC,NF DOUBLE PRECISION AF(*),AFD(*),AG(*),S(*) INTEGER IA(*) DOUBLE PRECISION POM,TEMP INTEGER JCG,KA DOUBLE PRECISION MXVDOT JCG = 1 DO 10 KA = 1,NA IF (IA(KA).GT.0) THEN TEMP = MXVDOT(NF,AG(JCG),S) AFD(KA) = TEMP TEMP = AF(KA) + TEMP IF (IA(KA).EQ.1 .OR. IA(KA).GE.3) THEN POM = XNORM + TEMP IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(XNORM),1.0D0))) THEN INEW = KA + NC KNEW = 1 PAR = POM END IF END IF IF (IA(KA).EQ.2 .OR. IA(KA).GE.3) THEN POM = XNORM - TEMP IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(XNORM),1.0D0))) THEN INEW = KA + NC KNEW = -1 PAR = POM END IF END IF END IF JCG = JCG + NF 10 CONTINUE RETURN END * SUBROUTINE PLMINL ALL SYSTEMS 90/12/01 * PORTABILITY : ALL SYSTEMS * 90/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF THE NEW ACTIVE LINEAR CONSTRAINT. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * RI CF(NC) 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. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI S(NF) DIRECTION VECTOR. * 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. * RI EPS9 TOLERANCE FOR ACTIVE CONSTRAINTS. * RA PAR AUXILIARY VARIABLE. * * SUBPROGRAMS USED : * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * SUBROUTINE PLMINL(NF,NC,CF,IC,CL,CU,CG,S,KBC,INEW,KNEW,EPS9,PAR) DOUBLE PRECISION EPS9,PAR INTEGER INEW,KBC,KNEW,NC,NF DOUBLE PRECISION CF(*),CG(*),CL(*),CU(*),S(*) INTEGER IC(*) DOUBLE PRECISION POM,TEMP INTEGER JCG,KC DOUBLE PRECISION MXVDOT IF (KBC.GT.0) THEN JCG = 1 DO 10 KC = 1,NC IF (IC(KC).GT.0) THEN TEMP = CF(KC) + MXVDOT(NF,CG(JCG),S) 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.0D0))) THEN INEW = KC KNEW = 1 PAR = POM END IF END IF 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.0D0))) THEN INEW = KC KNEW = -1 PAR = POM END IF END IF END IF JCG = JCG + NF 10 CONTINUE END IF 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 PLDLAG ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VECTOR OF LAGRANGE MULTIPLIERS FOR DUAL QP METHOD IS DETERMINED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II NC NUMBER OF LINEARIZED CONSTRAINTS. * II IA(NA) VECTOR CONTAINING TYPES OF DEVIATIONS. * II IAA(NF+1) VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS. * RO AZ(NF+1) OUTPUT VECTOR. * II N ACTUAL NUMBER OF VARIABLES. * IA KOLD AUXILIARY VARIABLE. * SUBROUTINE PLDLAG(NF,NC,IA,IAA,AZ,N,KOLD) INTEGER KOLD,N,NC,NF DOUBLE PRECISION AZ(*) INTEGER IA(*),IAA(*) DOUBLE PRECISION TEMP INTEGER J,L,NAA NAA = NF - N KOLD = 0 DO 10 J = 1,NAA L = IAA(J) IF (L.GT.NC) THEN L = L - NC TEMP = 1.0D0 IF (IA(L).EQ.-2 .OR. IA(L).EQ.-4) TEMP = -TEMP AZ(J) = TEMP KOLD = 1 ELSE AZ(J) = 0.0D0 END IF 10 CONTINUE RETURN END * SUBROUTINE PLADF1 ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * TRIANGULAR DECOMPOSITION OF KERNEL OF THE GENERAL PROJECTION * IS UPDATED AFTER FUNCTION OR CONSTRAINT ADDITION. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II NC NUMBER OF LINEARIZED CONSTRAINTS. * II IA(NA) VECTOR CONTAINING TYPES OF DEVIATIONS. * IU IAA(NF+1) VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS. * RI AG(NF*NA) MATRIX WHOSE COLUMNS ARE GRADIENTS OF THE LINEAR * APPROXIMATED FUNCTIONS. * RU AR((NF+1)*(NF+2)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF * THE ORTHOGONAL PROJECTION. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RI H(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OR INVERSION OF THE * HESSIAN MATRIX APPROXIMATION. * RA S(NF+1) AUXILIARY VECTOR. * RO G(NF+1) VECTOR USED IN THE DUAL RANGE SPACE QUADRATIC * PROGRAMMING METHOD. * IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION. * IDECF=10-DIAGONAL MATRIX. * IU N ACTUAL NUMBER OF VARIABLES. * II INEW INDEX OF THE NEW ACTIVE CONSTRAINT. * II KNEW SIGNUM OF THE NEW ACTIVE GRADIENT. * IO IER ERROR INDICATOR. * 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. * * SUBPROGRAMS USED : * S MXPDGB BACK SUBSTITUTION AFTER GILL-MURRAY DECOMPOSITION. * S MXDPRB BACK SUBSTITUTION AFTER CHOLESKI DECOMPOSITION. * 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 PLADF1(NF,NC,IA,IAA,AG,AR,CG,H,S,G,IDECF,N,INEW,KNEW, + IER,EPS7,GMAX,UMAX,E,T) DOUBLE PRECISION E,EPS7,GMAX,T,UMAX INTEGER IDECF,IER,INEW,KNEW,N,NC,NF DOUBLE PRECISION AG(*),AR(*),CG(*),G(*),H(*),S(*) INTEGER IA(*),IAA(*) INTEGER NADD,NDECF,NFG,NFH,NFV,NIT,NRED,NREM,NRES DOUBLE PRECISION POM,TEMP INTEGER J,JAG,JOB,K,L,NAA,NAR DOUBLE PRECISION MXVDOT COMMON /STAT/NDECF,NRES,NRED,NREM,NADD,NIT,NFV,NFG,NFH JOB = 1 E = 0.0D0 IF (INEW.GT.NC) E = SIGN(1,KNEW) IER = 0 IF (JOB.EQ.0 .AND. N.LT.0) IER = 2 IF (INEW.EQ.0) IER = 3 IF (IDECF.GE.2 .AND. IDECF.LE.8) IER = -2 IF (IER.NE.0) RETURN NAA = NF - N NAR = NAA* (NAA+1)/2 IF (INEW.GT.NC) THEN JAG = (INEW-NC-1)*NF + 1 IF (IDECF.EQ.1) THEN CALL MXVCOP(NF,AG(JAG),S) CALL MXDPGB(NF,H,S,0) ELSE IF (IDECF.EQ.9) THEN CALL MXDSMM(NF,H,AG(JAG),S) ELSE CALL MXVCOP(NF,AG(JAG),S) CALL MXVMUL(NF,H,S,S,-1) END IF GMAX = MXVDOT(NF,AG(JAG),S) + T ELSE IF (INEW.GT.0) THEN JAG = (INEW-1)*NF + 1 IF (IDECF.EQ.1) THEN CALL MXVCOP(NF,CG(JAG),S) CALL MXDPGB(NF,H,S,0) ELSE IF (IDECF.EQ.9) THEN CALL MXDSMM(NF,H,CG(JAG),S) ELSE CALL MXVCOP(NF,CG(JAG),S) CALL MXVMUL(NF,H,S,S,-1) END IF GMAX = MXVDOT(NF,CG(JAG),S) ELSE K = -INEW IF (IDECF.EQ.1) THEN CALL MXVSET(NF,0.0D0,S) S(K) = 1.0D0 CALL MXDPGB(NF,H,S,0) ELSE IF (IDECF.EQ.9) THEN CALL MXDSMV(NF,H,S,K) ELSE CALL MXVSET(NF,0.0D0,S) S(K) = 1.0D0/H(K) END IF GMAX = S(K) END IF IF (NAA.GT.0) THEN POM = T*E DO 10 J = 1,NAA L = IAA(J) IF (L.GT.NC) THEN L = L - NC G(J) = MXVDOT(NF,AG((L-1)*NF+1),S) IF (INEW.GT.NC) THEN TEMP = POM IF (IA(L).EQ.-2 .OR. IA(L).EQ.-4) TEMP = -TEMP G(J) = G(J) + TEMP END IF ELSE IF (L.GT.0) THEN G(J) = MXVDOT(NF,CG((L-1)*NF+1),S) ELSE L = -L G(J) = S(L) END IF 10 CONTINUE END IF IF (N.LT.0) THEN CALL MXDPRB(NAA,AR,G,1) UMAX = 0.0D0 IER = 2 RETURN ELSE IF (NAA.EQ.0) THEN UMAX = GMAX ELSE CALL MXDPRB(NAA,AR,G,1) UMAX = GMAX - MXVDOT(NAA,G,G) CALL MXVCOP(NAA,G,AR(NAR+1)) END IF IF (UMAX.LE.EPS7*GMAX) THEN IER = 1 RETURN ELSE NAA = NAA + 1 NAR = NAR + NAA IAA(NAA) = INEW AR(NAR) = SQRT(UMAX) IF (JOB.EQ.0) THEN N = N - 1 NADD = NADD + 1 END IF END IF 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,NDECF,NFG,NFH,NFV,NIT,NRED,NREM,NRES INTEGER L COMMON /STAT/NDECF,NRES,NRED,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 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) DOUBLE PRECISION EPS7,GMAX,UMAX INTEGER IER,INEW,N,NADD,NF DOUBLE PRECISION CG(*),CR(*),CZ(*),S(*) INTEGER ICA(*) 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) END IF DO 10 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 10 CONTINUE END IF 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) DOUBLE PRECISION EPS7,GMAX,UMAX INTEGER IDECF,IER,INEW,N,NADD,NF DOUBLE PRECISION CG(*),CR(*),CZ(*),H(*),S(*) INTEGER ICA(*) DOUBLE PRECISION CK,CL INTEGER I,J,K,L,N1 IF (IDECF.NE.0 .AND. IDECF.NE.9) THEN IER = -2 RETURN END IF 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) END IF DO 10 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 10 CONTINUE IF (IDECF.EQ.9) THEN L = N* (N+1)/2 IF (H(L+N1).NE.0.0D0) THEN CL = 1.0D0/H(L+N1) K = 0 DO 30 I = 1,N CK = CL*H(L+I) DO 20 J = 1,I K = K + 1 H(K) = H(K) - CK*H(L+J) 20 CONTINUE 30 CONTINUE END IF END IF END IF 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) DOUBLE PRECISION EPS7,GMAX,UMAX INTEGER IER,INEW,N,NADD,NF DOUBLE PRECISION CG(*),CR(*),S(*) INTEGER ICA(*) INTEGER I,J,K,L,NCA,NCR DOUBLE PRECISION MXVDOT 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 10 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) END IF 10 CONTINUE ELSE K = -INEW GMAX = 1.0D0 DO 20 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.0D0 END IF 20 CONTINUE END IF 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)) END IF 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 END IF 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) DOUBLE PRECISION STEP INTEGER KBC,NC DOUBLE PRECISION CF(*),CFD(*) INTEGER IC(*) INTEGER KC IF (KBC.GT.0) THEN DO 10 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) END IF 10 CONTINUE END IF 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) DOUBLE PRECISION STEP INTEGER KBF,NF DOUBLE PRECISION S(*),X(*) INTEGER IX(*) INTEGER I DO 10 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) END IF 10 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) DOUBLE PRECISION EPS9 INTEGER IND,INEW,KBF,NF DOUBLE PRECISION X(*),XL(*),XU(*) INTEGER IX(*) INTEGER I IND = 0 IF (KBF.GT.0) THEN DO 10 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) END IF 10 CONTINUE END IF 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) DOUBLE PRECISION EPS7,EPS9,ETA9,GMAX,UMAX INTEGER ITERL,KBC,KBF,MFP,N,NC,NF DOUBLE PRECISION CF(*),CFD(*),CG(*),CL(*),CR(*),CU(*),CZ(*),G(*), + GO(*),S(*),X(*),XL(*),XO(*),XU(*) INTEGER IC(*),ICA(*),IX(*) INTEGER NADD,NDECF,NFG,NFH,NFV,NIT,NRED,NREM,NRES DOUBLE PRECISION CON,DMAX,POM INTEGER I,IER,INEW,IOLD,IPOM,K,KC,KREM,NCA,NCR,NCZ COMMON /STAT/NDECF,NRES,NRED,NREM,NADD,NIT,NFV,NFG,NFH CON = ETA9 * * INITIATION * CALL MXVCOP(NF,X,XO) IPOM = 0 NRED = 0 KREM = 0 ITERL = 1 DMAX = 0.0D0 IF (MFP.EQ.3) GO TO 40 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 60 END IF N = 0 NCA = 0 NCZ = 0 DO 10 I = 1,NF IF (KBF.GT.0) THEN IF (IX(I).LT.0) THEN NCA = NCA + 1 ICA(NCA) = -I ELSE N = N + 1 CALL MXVSET(NF,0.0D0,CZ(NCZ+1)) CZ(NCZ+I) = 1.0D0 NCZ = NCZ + NF ENDIF ELSE N = N + 1 CALL MXVSET(NF,0.0D0,CZ(NCZ+1)) CZ(NCZ+I) = 1.0D0 NCZ = NCZ + NF END IF 10 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) IF (NF.GT.N) CALL PLSETC(NF,NC,X,XO,CF,IC,CG,S) DO 20 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 END IF 20 CONTINUE END IF 30 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 END IF IF (IPOM.EQ.0 .AND. ITERL.EQ.0) THEN * * FEASIBILITY ACHIEVED * ITERL = 1 CALL MXVCOP(NF,G,GO) IF (MFP.EQ.1) GO TO 60 END IF * * LAGRANGE MULTIPLIERS AND REDUCED GRADIENT DETERMINATION * 40 CALL PLTRBG(NF,N,NC,IX,IC,ICA,CG,CR,CZ,GO,S,EPS7,GMAX,UMAX,IOLD) INEW = 0 IF (GMAX.EQ.0.0D0) THEN * * OPTIMUM ON A LINEAR MANIFOLD OBTAINED * IF (IOLD.EQ.0) THEN IF (IPOM.EQ.0) THEN * * OPTIMAL SOLUTION ACHIEVED * ITERL = 2 GO TO 60 ELSE IPOM = 0 DO 50 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 END IF 50 CONTINUE IF (IPOM.EQ.0) THEN * * OPTIMAL SOLUTION ACHIEVED * CALL MXVCOP(NF,GO,G) ITERL = 2 GO TO 60 ELSE * * FEASIBLE SOLUTION DOES NOT EXIST * CALL MXVCOP(NF,GO,G) ITERL = -1 GO TO 60 END IF END IF 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) END IF DMAX = 0.0D0 GO TO 40 END IF 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 END IF GO TO 60 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) END IF DMAX = POM IF (DMAX.GT.0.0D0) NRED = NRED + 1 GO TO 30 END IF END IF 60 CONTINUE RETURN END * SUBROUTINE PLLPB2 ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/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. * 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 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. * 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. * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RO GO(NF) SAVED GRADIENT OF THE OBJECTIVE FUNCTION. * RI 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 PLGLAG GRADIENT OF THE LAGRANGIAN FUNCTION IS DETERMINED. * S PLSLAG NEGATIVE PROJECTED GRADIENT IS DETERMINED. * S PLTLAG THE OPTIMUM LAGRANGE MULTIPLIER IS DETERMINED. * S PLVLAG AN AUXILIARY VECTOR IS DETERMINED. * S PLADR0 CONSTRAINT ADDITION. * S PLRMF0 CONSTRAINT DELETION. * S MXDPRB BACK SUBSTITUTION AFTER CHOLESKI DECOMPOSITION. * S MXDSMI DETERMINATION OF THE INITIAL UNIT DENSE SYMMETRIC * MATRIX. * S MXVCOP COPYING OF A VECTOR. * S MXVDIF DIFFERENCE OF TWO VECTORS. * 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 PLLPB2(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) DOUBLE PRECISION EPS7,EPS9,ETA9,GMAX,UMAX INTEGER ITERL,KBC,KBF,MFP,N,NC,NF DOUBLE PRECISION CF(*),CFD(*),CG(*),CL(*),CR(*),CU(*),CZ(*),G(*), + GO(*),S(*),X(*),XL(*),XO(*),XU(*) INTEGER IC(*),ICA(*),IX(*) INTEGER NADD,NDECF,NFG,NFH,NFV,NIT,NRED,NREM,NRES DOUBLE PRECISION CON,DMAX,POM INTEGER I,IER,INEW,IOLD,IPOM,KC,KREM,MODE COMMON /STAT/NDECF,NRES,NRED,NREM,NADD,NIT,NFV,NFG,NFH CON = ETA9 * * INITIATION * CALL MXVCOP(NF,X,XO) CALL MXVCOP(NF,G,GO) IPOM = 0 NRED = 0 KREM = 0 ITERL = 1 DMAX = 0.0D0 IF (MFP.EQ.3) GO TO 40 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 60 END IF N = NF DO 10 I = 1,NF IF (KBF.GT.0) THEN IF (IX(I).LT.0) THEN N = N - 1 ICA(NF-N) = -I END IF END IF 10 CONTINUE CALL MXDSMI(NF-N,CR) IF (NC.GT.0) THEN * * ADDITION OF ACTIVE CONSTRAINTS AND INITIAL CHECK OF FEASIBILITY * CALL MXVINA(NC,IC) IF (NF.GT.N) CALL PLSETC(NF,NC,X,XO,CF,IC,CG,S) DO 20 KC = 1,NC IF (IC(KC).NE.0) THEN INEW = 0 CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW) CALL PLADR0(NF,N,ICA,CG,CR,S,EPS7,GMAX,UMAX,INEW,NADD, + IER) CALL MXVIND(IC,KC,IER) IF (IC(KC).LT.-10) IPOM = 1 END IF 20 CONTINUE END IF 30 IF (IPOM.EQ.1) THEN * * CHECK OF FEASIBILITY AND UPDATE OF THE FIRST PHASE OBJECTIVE * FUNCTION * CALL PLSETG(NF,NC,IC,CG,G,INEW) IF (INEW.EQ.0) IPOM = 0 END IF IF (IPOM.EQ.0 .AND. ITERL.EQ.0) THEN * * FEASIBILITY ACHIEVED * ITERL = 1 CALL MXVCOP(NF,GO,G) IF (MFP.EQ.1) GO TO 60 END IF * * LAGRANGE MULTIPLIERS DETERMINATION * 40 IF (NF.GT.N) THEN CALL PLVLAG(NF,N,NC,ICA,CG,CG,G,CZ) CALL MXDPRB(NF-N,CR,CZ,0) CALL PLTLAG(NF,N,NC,IX,IC,ICA,CZ,IC,EPS7,UMAX,IOLD) ELSE IOLD = 0 UMAX = 0.0D0 END IF * * PROJECTED GRADIENT DETERMINATION * IF (N.GT.0) THEN CALL MXVNEG(NF,G,S) CALL PLSLAG(NF,N,NC,ICA,CG,CZ,CG,S,EPS7,GMAX) ELSE GMAX = 0.0D0 END IF MODE = 1 - IPOM INEW = 0 IF (GMAX.EQ.0.0D0) THEN * * OPTIMUM ON A LINEAR MANIFOLD OBTAINED * IF (IOLD.EQ.0) THEN IF (IPOM.EQ.0) THEN * * OPTIMAL SOLUTION ACHIEVED * ITERL = 2 GO TO 60 ELSE IPOM = 0 DO 50 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 END IF 50 CONTINUE IF (IPOM.EQ.0) THEN * * OPTIMAL SOLUTION ACHIEVED * CALL MXVCOP(NF,GO,G) ITERL = 2 GO TO 60 ELSE * * FEASIBLE SOLUTION DOES NOT EXIST * CALL MXVCOP(NF,GO,G) ITERL = -1 GO TO 60 END IF END IF ELSE * * CONSTRAINT DELETION * CALL PLRMF0(NF,NC,IX,IC,ICA,CR,IC,S,N,IOLD,KREM,IER) DMAX = 0.0D0 GO TO 40 END IF ELSE * * 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 END IF GO TO 60 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 PLADR0(NF,N,ICA,CG,CR,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 PLADR0(NF,N,ICA,CG,CR,S,EPS7,GMAX,UMAX,INEW,NADD, + IER) CALL MXVIND(IX,I,IER) END IF DMAX = POM NRED = NRED + 1 GO TO 30 END IF END IF 60 CONTINUE 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) DOUBLE PRECISION STEP INTEGER INEW,KBC,KREM,NC,NF DOUBLE PRECISION CF(*),CFD(*),CG(*),CL(*),CU(*),S(*) INTEGER IC(*) DOUBLE PRECISION TEMP INTEGER JCG,KC DOUBLE PRECISION MXVDOT IF (KBC.GT.0) THEN JCG = 1 DO 10 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.0D0) 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 END IF END IF ELSE IF (TEMP.GT.0.0D0) 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 END IF END IF END IF ELSE IF (IC(KC).LT.-10) THEN TEMP = MXVDOT(NF,CG(JCG),S) CFD(KC) = TEMP IF (TEMP.GT.0.0D0) 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 END IF END IF ELSE IF (TEMP.LT.0.0D0) 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 END IF END IF END IF END IF JCG = JCG + NF 10 CONTINUE END IF 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) DOUBLE PRECISION STEP INTEGER INEW,KBF,KREM,NF DOUBLE PRECISION S(*),X(*),XL(*),XU(*) INTEGER IX(*) DOUBLE PRECISION TEMP INTEGER I IF (KBF.GT.0) THEN DO 10 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.0D0) 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 END IF END IF ELSE IF (S(I).GT.0.0D0) 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 END IF END IF END IF END IF 10 CONTINUE END IF 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) DOUBLE PRECISION EPS9 INTEGER INEW,KC DOUBLE PRECISION CF(*),CL(*),CU(*) INTEGER IC(*) 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.0D0) 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 END IF ELSE IF (IC(KC).EQ.2) THEN TEMP = EPS9*MAX(ABS(CU(KC)),1.0D0) 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 END IF ELSE IF (IC(KC).EQ.3 .OR. IC(KC).EQ.4) THEN TEMP = EPS9*MAX(ABS(CL(KC)),1.0D0) IF (CF(KC).GT.CL(KC)+TEMP) THEN TEMP = EPS9*MAX(ABS(CU(KC)),1.0D0) 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 END IF ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN IC(KC) = 13 INEW = KC ELSE IC(KC) = -13 END IF ELSE IF (IC(KC).EQ.5 .OR. IC(KC).EQ.6) THEN TEMP = EPS9*MAX(ABS(CL(KC)),1.0D0) IF (CF(KC).GT.CL(KC)+TEMP) THEN TEMP = EPS9*MAX(ABS(CU(KC)),1.0D0) 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 END IF ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN IC(KC) = 15 INEW = KC ELSE IC(KC) = -15 END IF 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) DOUBLE PRECISION EPS9 INTEGER I,INEW DOUBLE PRECISION X(*),XL(*),XU(*) INTEGER IX(*) DOUBLE PRECISION TEMP TEMP = 1.0D0 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 END IF 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 END IF 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 END IF IF (X(I).GE.XU(I)-EPS9*MAX(ABS(XU(I)),TEMP)) THEN IX(I) = 14 INEW = -I END IF END IF 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 IER,IOLD,KREM,N,NF,NREM DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*) INTEGER ICA(*) INTEGER I,J,KC,NCA,NCR,NCZ DOUBLE PRECISION MXVDOT 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.0D0,CZ(NCZ+1)) CZ(NCZ+NCA) = 1.0D0 CALL MXDPRB(NCA,CR,CZ(NCZ+1),-1) CALL MXVCOP(NCA,CZ(NCZ+1),CR(NCR+1)) N = N + 1 CALL MXVSET(NF,0.0D0,CZ(NCZ+1)) DO 10 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) END IF 10 CONTINUE GN(N) = MXVDOT(NF,CZ(NCZ+1),G) NREM = NREM + 1 IER = 0 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 NC,NF DOUBLE PRECISION CF(*),CG(*),S(*),X(*),XO(*) INTEGER IC(*) INTEGER JCG,KC DOUBLE PRECISION MXVDOT CALL MXVDIF(NF,X,XO,S) JCG = 0 DO 10 KC = 1,NC IF (IC(KC).NE.0) CF(KC) = CF(KC) + MXVDOT(NF,S,CG(JCG+1)) JCG = JCG + NF 10 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 INEW,NC,NF DOUBLE PRECISION CG(*),G(*) INTEGER IC(*) INTEGER KC CALL MXVSET(NF,0.0D0,G) INEW = 0 DO 10 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.0D0,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.0D0,CG((KC-1)*NF+1),G,G) INEW = 1 END IF 10 CONTINUE RETURN END * SUBROUTINE PLGLAG ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * GRADIENT OF THE LAGRANGIAN FUNCTION IS DETERMINED. * * 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*NC) MATRIX WHOSE COLUMNS ARE GRADIENTS OF THE LINEAR * APPROXIMATED FUNCTIONS. * RO AZ(NF+1) VECTOR OF LAGRANGE MPLTIPLIERS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RU G(NF) GRADIENT OF THE LAGRANGIAN FUNCTION. * * SUBPROGRAMS USED : * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * SUBROUTINE PLGLAG(NF,N,NC,IAA,AG,AZ,CG,G) INTEGER N,NC,NF DOUBLE PRECISION AG(*),AZ(*),CG(*),G(*) INTEGER IAA(*) INTEGER J,L,NAA NAA = NF - N DO 10 J = 1,NAA 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 RETURN END * SUBROUTINE PLSLAG ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * NEGATIVE PROJECTED GRADIENT IS DETERMINED USING LAGRANGE MULTIPLIERS. * * 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*NC) MATRIX WHOSE COLUMNS ARE GRADIENTS OF THE LINEAR * APPROXIMATED FUNCTIONS. * RO AZ(NF+1) VECTOR OF LAGRANGE MULTIPLIERS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RO S(NF) NEGATIVE PROJECTED GRADIENT OF THE QUADRATIC FUNCTION. * RI EPS7 TOLERANCE FOR LINEAR AND QUADRATIC PROGRAMMING. * RO GMAX NORM OF THE TRANSFORMED GRADIENT. * * SUBPROGRAMS USED : * S UXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * RF UXVMAX L-INFINITY NORM OF A VECTOR. * SUBROUTINE PLSLAG(NF,N,NC,IAA,AG,AZ,CG,S,EPS7,GMAX) DOUBLE PRECISION EPS7,GMAX INTEGER N,NC,NF DOUBLE PRECISION AG(*),AZ(*),CG(*),S(*) INTEGER IAA(*) INTEGER J,L,NAA DOUBLE PRECISION MXVMAX NAA = NF - N DO 10 J = 1,NAA L = IAA(J) IF (L.GT.NC) THEN L = L - NC CALL MXVDIR(NF,AZ(J),AG((L-1)*NF+1),S,S) ELSE IF (L.GT.0) THEN CALL MXVDIR(NF,AZ(J),CG((L-1)*NF+1),S,S) ELSE L = -L S(L) = S(L) + AZ(J) END IF 10 CONTINUE GMAX = MXVMAX(NF,S) IF (GMAX.LE.EPS7) GMAX = 0.0D0 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) DOUBLE PRECISION EPS7,UMAX INTEGER IOLD,N,NC,NF DOUBLE PRECISION AZ(*) INTEGER IA(*),IAA(*),IC(*),IX(*) DOUBLE PRECISION TEMP INTEGER J,K,L,NAA IOLD = 0 UMAX = 0.0D0 NAA = NF - N DO 10 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) END IF IF (K.LE.-5) THEN ELSE IF ((K.EQ.-1.OR.K.EQ.-3) .AND. UMAX+TEMP.GE.0.0D0) THEN ELSE IF ((K.EQ.-2.OR.K.EQ.-4) .AND. UMAX-TEMP.GE.0.0D0) THEN ELSE IOLD = J UMAX = ABS(TEMP) END IF 10 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) DOUBLE PRECISION EPS7,GMAX,UMAX INTEGER IOLD,N,NC,NF DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*) INTEGER IC(*),ICA(*),IX(*) INTEGER NCA,NCZ DOUBLE PRECISION MXVMAX GMAX = 0.0D0 IF (N.GT.0) THEN CALL MXDRMM(NF,N,CZ,G,GN) GMAX = MXVMAX(N,GN) END IF 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.0D0,GN) GMAX = 0.0D0 ELSE IOLD = 0 UMAX = 0.0D0 END IF 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 N,NC,NF DOUBLE PRECISION AG(*),CG(*),G(*),GN(*) INTEGER IAA(*) INTEGER J,L,NAA DOUBLE PRECISION MXVDOT NAA = NF - N DO 10 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) END IF 10 CONTINUE 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 MET SELECTION OF SELF SCALING. MET=1-SELF SCALING SUPPRESSED. * MET=2 INITIAL SELF SCALING. * II MOD CORRECTION IF THE NEGATIVE CURVATURE OCCURS. * MOD=1-CORRECTION SUPPRESSED. MOD=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,MOD) DOUBLE PRECISION PO,R INTEGER ITERH,KIT,MET,MOD,N,NIT DOUBLE PRECISION G(*),GO(*),H(*),S(*),XO(*) DOUBLE PRECISION A,B,C,DIS,GAM,PAR LOGICAL L1,L3 DOUBLE PRECISION MXDPGP,MXVDOT L1 = MET .GE. 3 .OR. MET .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 (MOD.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) END IF ELSE IF (B.LE.1.0D-4*C) THEN ITERH = 2 RETURN END IF END IF IF (L1) THEN * * DETERMINATION OF THE PARAMETER GAM (SELF SCALING) * PAR = C/B GAM = PAR IF (MET.GT.1) THEN IF (NIT.NE.KIT) THEN L3 = GAM .LT. 0.5D0 .OR. GAM .GT. 4.0D0 END IF END IF END IF IF (L3) THEN GAM = 1.0D0 PAR = GAM END IF * * BFGS UPDATE * CALL MXDPGU(N,H,PAR/B,GO,XO) CALL MXDPGU(N,H,-1.0D0/C,S,XO) ITERH = 0 IF (GAM.EQ.1.0D0) RETURN CALL MXDPGS(N,H,1.0D0/GAM) RETURN END * SUBROUTINE PUDVI2 ALL SYSTEMS 99/12/01 * PORTABILITY : ALL SYSTEMS * 99/12/01 VL : ORIGINAL VERSION * * PURPOSE : * NONSMOOTH VARIABLE METRIC UPDATE OF THE INVERSE HESSIAN MATRIX. * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * RU H(N*(N+1)/2) POSITIVE DEFINITE APPROXIMATION OF THE INVERSE * HESSIAN MATRIX. * RI S(N) DIRECTION VECTOR. * RI GV(N) AGGREGATED SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RU U(N) DIFFERENCE OF CURRENT AND PREVIOUS GRADIENTS. * RA V(N) AUXILIARY VECTOR. * RI T VALUE OF THE STEPSIZE PARAMETER. * RI Z DOT PRODUCT OF VECTORS S, U. * RI RO CORRECTION PARAMETER. * II JC CORRECTION INDICATOR. * II JU UPDATING INDICATOR. * II NNK CONSECUTIVE NULL STEPS COUNTER. * II JOB TYPE OF MINIMIZATION. JOB=0 - CONVEX MINIMIZATION. * JOB=1 - NONCONVEX MINIMIZATION. * * COMMON DATA : * * SUBPROGRAMS USED : * S MXDSMM MULTIPLICATION OF A DENSE SYMMETRIC MATRIX BY A VECTOR. * S MXDSMU UPDATE OF A DENSE SYMMETRIC MATRIX. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * RF MXVDEL SQUARED NORM OF A SHIFTED VECTOR. * RF MXVDOT DOT PRODUCT OF VECTORS. * S MXVLIN LINEAR COMBINATION OF TWO VECTORS. * SUBROUTINE PUDVI2(N,H,S,GV,U,V,Z,RO,JC,JU,NNK,JOB,NIT) DOUBLE PRECISION RO,Z INTEGER JC,JOB,JU,N,NIT,NNK DOUBLE PRECISION GV(*),H(*),S(*),U(*),V(*) DOUBLE PRECISION GAM,P,P1,P2,P3,Q,W INTEGER I,J,L LOGICAL LB,LR DOUBLE PRECISION MXVDEL,MXVDOT IF (JOB.GT.0) JU = 0 IF (Z.LE.0.0D0) RETURN CALL MXDSMM(N,H,U,V) W = MXVDOT(N,U,V) GAM = 1.0D0 IF (NIT.EQ.1) THEN Q = 1.0D0 IF (W.NE.0.0D0) Q = Z/W IF ((Q-2.5D-1)* (Q-3.0D0).GT.0.0D0) GAM = MIN(3.0D0, + MAX(2.0D-2,Q)) END IF P1 = MXVDEL(N,-1.0D0,S,V) P2 = MXVDOT(N,GV,V) - MXVDOT(N,GV,S) P3 = MXVDOT(N,GV,GV) P = W - Z LB = NNK .EQ. 0 LR = NNK .NE. 0 .AND. P2 .LT. 0.0D0 IF (JC.EQ.1 .AND. (P1.LT.RO*P*N.OR.P2*P2.LT.RO*P*P3)) LR = .FALSE. IF (LB) THEN CALL MXVLIN(N,1.0D0/Z,V,-0.5D0* (1.0D0/GAM+W/Z)/Z,S,U) L = 1 IF (JOB.GT.0) JU = 1 DO 20 I = 1,N DO 10 J = 1,I H(L) = (H(L)-U(I)*S(J)-S(I)*U(J))*GAM L = L + 1 10 CONTINUE 20 CONTINUE ELSE IF (LR) THEN CALL MXVDIF(N,V,S,U) JU = 1 CALL MXDSMU(N,H,-1.0D0/P,U) END IF RETURN END * SUBROUTINE PS0LA2 ALL SYSTEMS 94/12/01 * PORTABILITY : ALL SYSTEMS * 94/12/01 LU : ORIGINAL VERSION * * PURPOSE : * EXTENDED LINE SEARCH WITHOUT DIRECTIONAL DERIVATIVES. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II NA NUMBER OF APPROXIMATED FUNCTIONS. * RU X(NF) VECTOR OF VARIABLES. * RI XO(NF) OLD VECTOR OF VARIABLES. * RI S(NF) DIRECTION VECTOR. * RO R VALUE OF THE STEPSIZE PARAMETER. * RO RO PREVIOUS 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 RMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER. * RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER. * RI FMIN LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. * RI FMAX UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. * RA FA VALUE OF THE APPROXIMATED FUNCTION. * RO AF(NA) VECTOR WHOSE ELEMENTS ARE VALUES OF THE * APPROXIMATED FUNCTIONS. * RA GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. * RO AG(NF*NA) MATRIX WHOSE COLUMNS ARE GRADIENTS OF THE * APPROXIMATED FUNCTIONS. * RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * II KD DEGREE OF REQUIRED DERVATIVES. * IU LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II IEXT TYPE OF OBJECTIVE FUNCTION. IEXT<0-MAXIMUM OF POSITIVE * VALUES. IEXT=0-MAXIMUM OF ABSOLUTE VALUES. IEXT>0-MAXIMUM * OF NEGATIVE VALUES. * II NIT ACTUAL NUMBER OF ITERATIONS. * II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. * RI TOLS TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE * CHANGE OF THE FUNCTION VALUE). * II MES METHOD SELECTION. MES=1-BISECTION. MES=2-TWO POINT * QUADRATIC INTERPOLATION. MES=3-THREE POINT QUADRATIC * INTERPOLATION. * IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * II MRED MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * 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. * * SUBPROGRAM USED : * S UA1MX2 COMPUTATION OF THE VALUE AND THE GRADIENT OF THE * OBJECTIVE FUNCTION WHICH IS DEFINED AS A MAXIMUM OF THE * APPROXIMATED FUNCTIONS. * S PNINT3 EXTRAPOLATION OR INTERPOLATION WITHOUT DIRECTIONAL * DERIVATIVES. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * * METHOD : * SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH EXTENDED TERMINATION * CRITERIA. * SUBROUTINE PS0LA2(NF,NA,X,XO,S,R,RO,F,FO,PO,RMIN,RMAX,FMIN,FMAX, + FA,AF,GA,AG,G,KD,LD,IEXT,NIT,KIT,TOLS,MES,NRED, + MRED,ITERS) DOUBLE PRECISION F,FA,FMAX,FMIN,FO,PO,R,RMAX,RMIN,RO,TOLS INTEGER IEXT,ITERS,KD,KIT,LD,MES,MRED,NA,NF,NIT,NRED DOUBLE PRECISION AF(*),AG(*),G(*),GA(*),S(*),X(*),XO(*) DOUBLE PRECISION FI,FL,FU,RI,RL,RU INTEGER MERR,MODE,MTYP LOGICAL L1,L2,L3 IF (PO.GE.0.0D0) THEN R = 0.0D0 ITERS = -2 RETURN END IF IF (RMAX.LE.0.0D0) THEN ITERS = 0 RETURN END IF * * INITIAL STEPSIZE SELECTION * R = 1.0D0 R = MAX(R,RMIN) R = MIN(R,RMAX) MODE = 0 RU = 0.0D0 FU = FO RI = 0.0D0 FI = FO * * NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION) * 10 CALL PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR) IF (MERR.GT.0) THEN ITERS = -MERR RETURN ELSE IF (MODE.EQ.1) THEN NRED = NRED - 1 R = MIN(R,RMAX) ELSE IF (MODE.EQ.2) THEN NRED = NRED + 1 END IF * * COMPUTATION OF THE NEW FUNCTION VALUE * KD = 0 LD = -1 CALL MXVDIR(NF,R,S,XO,X) CALL PA1MX2(NF,NA,X,F,FA,AF,GA,AG,G,KD,LD,IEXT) KD = 1 IF (F.LE.FMIN) THEN ITERS = 7 RETURN 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.0D1 END IF * * TEST ON TERMINATION * IF (L1 .AND. .NOT.L3) THEN ITERS = 0 RETURN ELSE IF (L2 .AND. .NOT.F.GE.FU) THEN ITERS = 7 RETURN ELSE IF (L3) THEN ITERS = 5 RETURN ELSE IF (ABS(NRED).GE.MRED) THEN ITERS = -1 RETURN ELSE MODE = MAX(MODE,1) MTYP = ABS(MES) IF (F.GE.FMAX) MTYP = 1 END IF IF (MODE.EQ.1) THEN * * INTERVAL CHANGE AFTER EXTRAPOLATION * RL = RI FL = FI RI = RU FI = FU RU = R FU = F IF (F.GE.FI) THEN NRED = 0 MODE = 2 END IF * * INTERVAL CHANGE AFTER INTERPOLATION ELSE IF (R.LE.RI) THEN IF (F.LE.FI) THEN RU = RI FU = FI RI = R FI = F ELSE RL = R FL = F END IF ELSE IF (F.LE.FI) THEN RL = RI FL = FI RI = R FI = F ELSE RU = R FU = F END IF END IF GO TO 10 END * SUBROUTINE PS1L05 ALL SYSTEMS 94/12/01 * PORTABILITY : ALL SYSTEMS * 94/12/01 VL : ORIGINAL VERSION * * PURPOSE : * SPECIAL LINE SEARCH WITH DIRECTIONAL DERIVATIVES FOR BUNDLE METHODS. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * RU X(NF) VECTOR OF VARIABLES. * RI XO(NF) OLD VECTOR OF VARIABLES. * RU S(NF+1) DIRECTION VECTOR. * RO R VALUE OF THE STEPSIZE PARAMETER. * RO RP PREVIOUS VALUE OF THE STEPSIZE PARAMETER. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RI FO PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RU FP CURRENT MINIMUM 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. * RU TO WEIGHT PARAMETER * RU G(NF+1) SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RO SNORM NORM OF THE DIRECTION VECTOR. * RI RMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER * RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER * RI FMIN LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. * RI FMAX UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. * RI TOLS LINE SEARCH PARAMETER (IN TEST ON THE DECREASE * OF THE FUNCTION VALUE). * RI TOLP LINE SEARCH PARAMETER (IN TERMINATION CONDITION * FOR NULL AND SHORT STEPS). * RI ETA DISTANCE MEASURE PARAMETER. * II MES INTERPOLATION 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. * II MES2 WEIGHT UPDATING METHOD SELECTION. MES2=1-QUADRATIC * INTERPOLATION. MES2=2-LOCAL MINIMUM LOCALIZATION. * 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=9-SHORT STEP. ITERS=10-ZERO STEP. * ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL * DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION. * * SUBPROGRAM USED : * SE FUNDER OBJECTIVE FUNCTION AND SUBGRADIENT EVALUATION. * S PNINT1 EXTRAPOLATION OR INTERPOLATION WITH DIRECTIONAL * DERIVATIVES. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * * METHOD : * SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH SPECIAL TERMINATION * CRITERIA. * SUBROUTINE PS1L05(NF,X,XO,S,R,RP,F,FO,FP,P,PO,PP,TO,G,SNORM,RMIN, + RMAX,FMIN,FMAX,TOLS,TOLP,ETA,MES,MES2,ITERS) DOUBLE PRECISION CON1 PARAMETER (CON1=2D-3) DOUBLE PRECISION ETA,F,FMAX,FMIN,FO,FP,P,PO,PP,R,RMAX,RMIN,RP, + SNORM,TO,TOLP,TOLS INTEGER ITERS,MES,MES2,NF DOUBLE PRECISION G(*),S(*),X(*),XO(*) INTEGER NADD,NDECF,NFG,NFH,NFV,NIT,NRED,NREM,NRES DOUBLE PRECISION FAUX,FL,FU,PL,PU,RL,RTEMP,RU,TOAUX,TOP,TOP1,TOP2, + TOPB,TOPC INTEGER IPOC,IPOCM,IPOCN,MERR,MODE,MTYP LOGICAL L3,L5,M0,M3,M4 DOUBLE PRECISION MXVDOT COMMON /STAT/NDECF,NRES,NRED,NREM,NADD,NIT,NFV,NFG,NFH SAVE IPOC,M4,TOP,TOP1,TOP2,TOAUX,FAUX,IPOCN,IPOCM IF (NIT.LE.1) THEN IPOC = 0 IPOCN = 0 TOAUX = 1.0D0 END IF IF (NIT.EQ.NIT/5*5+1) THEN FAUX = 0.0D0 IPOCM = IPOCN END IF M4 = RP .EQ. R IF (.NOT.M4) PP = P IF (RMAX.LE.0.0D0) THEN ITERS = 0 GO TO 40 END IF NRED = 0 * * INITIAL STEPSIZE SELECTION * RTEMP = FMIN - F RP = 0.0D0 FP = FO R = 1.0D0 R = MIN(MAX(R,RMIN),RMAX) MODE = 0 MTYP = 10 RL = 0.0D0 FL = FO PL = PO RU = 0.0D0 FU = FO PU = PO * * NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION) * 10 IF (RL.EQ.0.0D0) THEN MTYP = MIN(MTYP,2) PL = PP CALL PNINT1(RL,RU,FL,FU,PL,PU,R,MODE,MTYP,MERR) ELSE MERR = 0 R = 0.5D0* (RU+RL) END IF IF (MERR.GT.0) THEN ITERS = -MERR GO TO 40 ELSE IF (MODE.EQ.1) THEN NRED = NRED - 1 R = MIN(R,RMAX) ELSE IF (MODE.EQ.2) THEN NRED = NRED + 1 END IF * * NEW FUNCTION VALUE AND NEW DIRECTIONAL DERIVATIVE * CALL MXVDIR(NF,R,S,XO,X) CALL FUNDER(NF,X,F,G) NFV = NFV + 1 NFG = NFG + 1 P = MXVDOT(NF,G,S) IF (F.LE.FMIN) THEN ITERS = 7 GO TO 40 END IF L3 = (F-FO.LE.TOLS*R*PO) .OR. (F.LE.FO-1.0D0) MODE = MAX(MODE,1) IF (MODE.EQ.1) THEN * * INTERVAL CHANGE AFTER EXTRAPOLATION * TOP2 = 2.0D0* (P*R+FO-F)/ (R*SNORM)**2 TOP = TOP2 TOP1 = -1.0D0 IF (M4) TOP1 = 2.0D0* (F-FO-PP*R)/ (R*SNORM)**2 IF (TOP.LT.-CON1) THEN TOP = TOP1 ELSE IF (TOP1.GE.-CON1) TOP = MIN(TOP1,TOP) END IF IF (MES2.EQ.2 .AND. PP.LT.0.0D0) THEN IF (P.GE.0.0D0) TOP = MIN(TOP,TO/R* (1.0D0-P/PP)) IF (P.LT.0.0D0) TOP = MIN(TOP,TO/ (1.0D0+1.0D1*P/PP)) END IF M0 = M4 .OR. IPOC .GE. 6 .OR. L3 RL = RU FL = FU PL = PU RU = R FU = F PU = P IF (.NOT.L3) THEN NRED = 0 MODE = 2 END IF ELSE * * INTERVAL CHANGE AFTER INTERPOLATION * TOPB = 2.0D0* (P*R+FO-F)/ (R*SNORM)**2 TOPC = 2.0D0* (FU-F-P* (RU-R))/ ((RU-R)*SNORM)**2 IF (TOPB.LT.-CON1) THEN TOP = TOPC ELSE IF (TOPC.GE.-CON1) TOP = MIN(TOPB,TOPC) IF (TOPC.LT.-CON1) TOP = TOPB END IF IF (L3) THEN RL = R FL = F PL = P ELSE RU = R FU = F PU = P END IF END IF * * TEST ON TERMINATION * IF (F.LT.FP) THEN RP = R FP = F END IF RTEMP = MAX(ABS(FP-F+P* (R-RP)),ETA* ((R-RP)*SNORM)**2) M3 = P - RTEMP .GE. TOLP*PO L5 = RL .EQ. 0.0D0 MTYP = 1 IF (L5) MTYP = MES IF (F.GE.FMAX) MTYP = 1 IF (MODE.EQ.1 .AND. L3) L5 = .FALSE. IF (L3 .AND. RP.GT.RMIN) GO TO 30 IF (M3 .AND. (R.LT.1D0.OR. (IPOC.LE.6.AND.NIT.GT.1))) GO TO 20 IF (M3 .AND. .NOT.M0) GO TO 20 IF (RU-RL.GE.RMIN) GO TO 10 20 IF (RP.EQ.0.0D0) THEN ITERS = 10 IPOC = IPOC + 1 ELSE ITERS = 9 IPOC = 0 END IF GO TO 40 30 ITERS = 5 IPOC = 0 IF (MODE.EQ.1) IPOCN = IPOCN + 1 R = RP F = FP 40 CONTINUE * * WEIGHT UPDATING * FAUX = MAX(FAUX,F-FP) IF (NIT.EQ.NIT/5*5 .AND. MES2.EQ.2) THEN IF (IPOCN-IPOCM.GE.4) THEN TOAUX = MAX(TOAUX/1.75D0,CON1) ELSE IF (IPOCN.EQ.IPOCM .AND. FAUX.GE.-R*PO .AND. + NIT.GT.10) THEN TOAUX = MIN(TOAUX*1.25D0,1.0D0/CON1) END IF END IF TOP = TOP*TOAUX RU = (TOP+CON1*1D-1)* (TOP*CON1*1D-2-1.0D0) IF (M0 .AND. RU.LE.0.0D0 .AND. MES2.EQ. + 1) TO = MIN(MAX(TOP,CON1,TO/1.0D1),1.0D0/CON1,TO*1.0D1) IF (M0 .AND. RU.LE.0.0D0 .AND. MES2.EQ. + 2) TO = MIN(MAX(TOP,CON1/1.0D1,TO/1.0D1),1.0D0/CON1,TO*1.0D1) RETURN END * SUBROUTINE PS1L07 ALL SYSTEMS 99/12/01 * PORTABILITY : ALL SYSTEMS * 99/12/01 VL : ORIGINAL VERSION * * PURPOSE : * SPECIAL LINE SEARCH FOR NONSMOOTH CONVEX VARIABLE METRIC METHOD. * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * II MA DECLARED NUMBER OF LINEAR APPROXIMATED FUNCTIONS * II MAL CURRENT NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * RU X(N) VECTOR OF VARIABLES. * RO G(N) SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RI S(N) DIRECTION VECTOR. * RU U(N) PREVIOUS VECTOR OF VARIABLES. * RI V(N) PREVIOUS SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RI AF(4*MA) VECTOR OF BUNDLE FUNCTIONS VALUES. * RI AG(N*MA) MATRIX WHOSE COLUMNS ARE BUNDLE SUBGRADIENTS. * RI AY(N*MA) MATRIX WHOSE COLUMNS ARE VARIABLE VECTORS. * RO T VALUE OF THE STEPSIZE PARAMETER. * RO TB BUNDLE PARAMETER FOR MATRIX SCALING. * RO FO PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RU PO PREVIOUS DIRECTIONAL DERIVATIVE. * RI TMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER. * RI TMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER. * RI D1 AUXILIARY VALUE FOR NULL/DESCENT STEP TEST. * RU DF AUXILIARY VALUE FOR TEST ON FUNCTION DECREASE. * RI ETA9 MAXIMUM FOR REAL NUMBERS. * RI TOLF LOWER BOUND FOR FUNCTION DECREASE. * IO JL TERMINATION INDICATOR. * IO NE EXTRAPOLATION COUNTER. * IO NK NULL STEP INDICATOR. NK=0-DESCENT STEP. NK=2-NULL STEP. * IU NV AUXILIARY NUMBER OF FUNCTION EVALUATIONS. * IU NTESF ACTUAL NUMBER OF TESTS ON FUNCTION DECREASE. * II MTESF MAXIMUM NUMBER OF TESTS ON FUNCTION DECREASE. * IO ITERS NULL STEP INDICATOR. ITERS=0-NULL STEP. ITERS=1-DESCENT * STEP. * * SUBPROGRAMS USED : * S MXVCOP COPYING OF A VECTOR. * RF MXVMX2 L-INFINITY NORM OF VECTOR DIFFERENCE. * S PNSTP2 STEPSIZE DETERMINATION FOR DESCENT STEPS. * S PNSTP3 STEPSIZE DETERMINATION FOR NULL STEPS. * * METHOD : * SPECIAL METHOD OF STEP LENGTH DETERMINATION. * SUBROUTINE PS1L07(N,MA,MAL,X,G,S,U,V,AF,AG,AY,T,TB,FO,F,PO,TMIN, + TMAX,D1,DF,ETA9,TOLF,JL,NE,NV,NTESF,MTESF,ITERS) DOUBLE PRECISION D1,DF,ETA9,F,FO,PO,T,TB,TMAX,TMIN,TOLF INTEGER ITERS,JL,MA,MAL,MTESF,N,NE,NTESF,NV DOUBLE PRECISION AF(*),AG(*),AY(*),G(*),S(*),U(*),V(*),X(*) INTEGER NADD,NDECF,NFG,NFH,NFV,NIT,NRED,NREM,NRES DOUBLE PRECISION DT,W DOUBLE PRECISION MXVMX2 COMMON /STAT/NDECF,NRES,NRED,NREM,NADD,NIT,NFV,NFG,NFH T = MIN(1.0D0,TMAX) IF (PO.NE.0.0D0) THEN IF (ITERS.EQ.1) THEN CALL PNSTP2(N,MA,MAL,X,AF,AG,AY,S,F,PO,T,TB,ETA9) ELSE CALL PNSTP3(N,MA,MAL,X,AF,AG,AY,S,F,PO,T,TB,ETA9) END IF END IF T = MIN(MAX(T,TMIN),TMAX) DT = T NE = 0 * * FUNCTION AND GRADIENT EVALUATION AT A NEW POINT * 10 CALL MXVDIR(N,T,S,U,X) CALL FUNDER(N,X,F,G) NFV = NFV + 1 NFG = NFG + 1 NV = NV + 1 * * NULL/DESCENT STEP TEST (ITERS=0/1) * ITERS = 1 IF (F.GT.FO-T*D1) ITERS = 0 W = DF IF (ABS(FO-F).GE.DF*1.0D-5) W = ABS(FO-F) IF (ITERS.EQ.1) DF = W IF (W/MAX(ABS(F),1.0D0).LE.TOLF .OR. FO.EQ.F) THEN NTESF = NTESF + 1 JL = 1 IF (NTESF.GE.MTESF) RETURN ELSE NTESF = 0 END IF W = MXVMX2(N,G,V) T = DT DT = DT + T IF (W.EQ.0.0D0 .AND. ITERS.EQ.1 .AND. NE.LE.9 .AND. + DT.LE.TMAX) THEN * * EXTRAPOLATION * NE = NE + 1 FO = F CALL MXVCOP(N,X,U) GO TO 10 END IF JL = 0 RETURN END * SUBROUTINE PS1L08 ALL SYSTEMS 99/12/01 * PORTABILITY : ALL SYSTEMS * 99/12/01 VL : ORIGINAL VERSION * * PURPOSE : * SPECIAL LINE SEARCH FOR NONSMOOTH NONCONVEX VARIABLE METRIC METHOD. * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * II MA DECLARED NUMBER OF LINEAR APPROXIMATED FUNCTIONS * II MAL CURRENT NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * RU X(N) VECTOR OF VARIABLES. * RO G(N) SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RI S(N) DIRECTION VECTOR. * RU U(N) PREVIOUS VECTOR OF VARIABLES. * RI AF(4*MA) VECTOR OF BUNDLE FUNCTIONS VALUES. * RI AG(N*MA) MATRIX WHOSE COLUMNS ARE BUNDLE SUBGRADIENTS. * RI AY(N*MA) MATRIX WHOSE COLUMNS ARE VARIABLE VECTORS. * RO T VALUE OF THE STEPSIZE PARAMETER. * RO TB BUNDLE PARAMETER FOR MATRIX SCALING. * RO FO PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RU PO PREVIOUS DIRECTIONAL DERIVATIVE. * RU P DIRECTIONAL DERIVATIVE. * RI TMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER. * RI TMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER. * RI SNORM EUCLIDEAN NORM OF THE DIRECTION VECTOR. * RI WK STOPPING PARAMETER. * RI EPS1 TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE * CHANGE OF THE FUNCTION VALUE). * RI EPS2 TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE * DIRECTIONAL DERIVATIVE). * RI ETA5 DISTANCE MEASURE PARAMETER. * RI ETA9 MAXIMUM FOR REAL NUMBERS. * II JE EXTRAPOLATION INDICATOR. * RI MOS3 LOCALITY MEASURE PARAMETER. * IO ITERS NULL STEP INDICATOR. ITERS=0-NULL STEP. ITERS=1-DESCENT * STEP. * SUBPROGRAMS USED : * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * S PNSTP4 STEPSIZE DETERMINATION FOR DESCENT STEPS. * S PNSTP5 STEPSIZE DETERMINATION FOR NULL STEPS. * S PNINT1 EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH * WITH DIRECTIONAL DERIVATIVES. * * METHOD : * SPECIAL METHOD OF STEP LENGTH DETERMINATION. * SUBROUTINE PS1L08(N,MA,MAL,X,G,S,U,AF,AG,AY,T,TB,FO,F,PO,P,TMIN, + TMAX,SNORM,WK,EPS1,EPS2,ETA5,ETA9,JE,MOS3,ITERS) DOUBLE PRECISION EPS1,EPS2,ETA5,ETA9,F,FO,P,PO,SNORM,T,TB,TMAX, + TMIN,WK INTEGER ITERS,JE,MA,MAL,MOS3,N DOUBLE PRECISION AF(*),AG(*),AY(*),G(*),S(*),U(*),X(*) INTEGER NADD,NDECF,NFG,NFH,NFV,NIT,NRED,NREM,NRES DOUBLE PRECISION BET,FL,FU,PL,PU,TL,TU INTEGER IER DOUBLE PRECISION MXVDOT COMMON /STAT/NDECF,NRES,NRED,NREM,NADD,NIT,NFV,NFG,NFH IF (JE.GT.0) T = DBLE(2-JE/99)*T IF (JE.LE.0) T = MIN(1.0D0,TMAX) IF (PO.EQ.0.0D0 .OR. JE.GT.0) GO TO 10 IF (ITERS.EQ.1) THEN CALL PNSTP4(N,MA,MAL,X,AF,AG,AY,S,F,PO,T,TB,ETA5,ETA9,MOS3) ELSE CALL PNSTP5(N,MA,MAL,X,AF,AG,AY,S,F,PO,T,TB,ETA5,ETA9,MOS3) END IF 10 T = MIN(MAX(T,TMIN),TMAX) TL = 0.0D0 TU = T FL = FO PL = PO * * FUNCTION AND GRADIENT EVALUATION AT A NEW POINT * 20 CALL MXVDIR(N,T,S,U,X) CALL FUNDER(N,X,F,G) NFV = NFV + 1 NFG = NFG + 1 P = MXVDOT(N,G,S) * * NULL/DESCENT STEP TEST (ITERS=0/1) * ITERS = 1 IF (F.LE.FO-T* (EPS1+EPS1)*WK) THEN TL = T FL = F PL = P ELSE TU = T FU = F PU = P END IF BET = MAX(ABS(FO-F+P*T),ETA5* (SNORM*T)**MOS3) IF (F.LE.FO-T*EPS1*WK .AND. (T.GE.TMIN.OR. + BET.GT.EPS1*WK)) GO TO 40 IF (P-BET.GE.-EPS2*WK .OR. TU-TL.LT.TMIN*1.0D-1) GO TO 30 IF (TL.EQ.0.0D0 .AND. PL.LT.0.0D0) THEN CALL PNINT1(TL,TU,FL,FU,PL,PU,T,2,2,IER) ELSE T = 5.0D-1* (TU+TL) END IF GO TO 20 30 ITERS = 0 40 CONTINUE RETURN END * SUBROUTINE PNINT1 ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 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 ZERO,HALF,ONE,TWO,THREE,C1L,C1U,C2L,C2U,C3L,FOUR 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=0.1D0, + FOUR=4.0D0) DOUBLE PRECISION FL,FU,PL,PU,R,RL,RU INTEGER MERR,MODE,MTYP DOUBLE PRECISION A,B,C,D,DEN,DIS INTEGER NTYP MERR = 0 IF (MODE.LE.0) RETURN IF (PL.GE.ZERO) THEN MERR = 2 RETURN ELSE IF (RU.LE.RL) THEN MERR = 3 RETURN END IF DO 10 NTYP = MTYP,1,-1 IF (NTYP.EQ.1) THEN * * BISECTION * IF (MODE.EQ.1) THEN R = FOUR*RU RETURN ELSE R = HALF* (RL+RU) RETURN END IF ELSE IF (NTYP.EQ.MTYP) THEN A = (FU-FL)/ (PL* (RU-RL)) B = PU/PL END IF IF (NTYP.EQ.2) THEN * * QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH ONE DIRECTIONAL * DERIVATIVE * DEN = TWO* (ONE-A) ELSE IF (NTYP.EQ.3) THEN * * QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH TWO DIRECTIONAL * DERIVATIVES * DEN = ONE - B ELSE IF (NTYP.EQ.4) THEN * * CUBIC EXTRAPOLATION OR INTERPOLATION * C = B - TWO*A + ONE D = B - THREE*A + TWO DIS = D*D - THREE*C IF (DIS.LT.ZERO) GO TO 10 DEN = D + SQRT(DIS) ELSE IF (NTYP.EQ.5) THEN * * CONIC EXTRAPOLATION OR INTERPOLATION * DIS = A*A - B IF (DIS.LT.ZERO) GO TO 10 DEN = A + SQRT(DIS) IF (DEN.LT.ZERO) GO TO 10 DEN = ONE - B* (ONE/DEN)**3 END IF IF (MODE.EQ.1 .AND. DEN.GT.ZERO .AND. DEN.LT.ONE) THEN * * EXTRAPOLATION ACCEPTED * R = RL + (RU-RL)/DEN R = MAX(R,C1L*RU) R = MIN(R,C1U*RU) RETURN ELSE IF (MODE.EQ.2 .AND. DEN.GT.ONE) THEN * * INTERPOLATION ACCEPTED * R = RL + (RU-RL)/DEN IF (RL.EQ.ZERO) 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)) RETURN END IF 10 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 PNSTP2 ALL SYSTEMS 99/12/01 * PORTABILITY : ALL SYSTEMS * 99/12/01 VL : ORIGINAL VERSION * * PURPOSE : * STEPSIZE SELECTION USING POLYHEDRAL APPROXIMATION * FOR DESCENT STEP IN CONVEX VARIABLE METRIC METHOD. * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * II MA DECLARED NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * II MAL CURRENT NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * RU X(N) VECTOR OF VARIABLES. * RI AF(4*MA) VECTOR OF BUNDLE FUNCTIONS VALUES. * RI AG(N*MA) MATRIX WHOSE COLUMNS ARE BUNDLE SUBGRADIENTS. * RI AY(N*MA) MATRIX WHOSE COLUMNS ARE VARIABLE VECTORS. * RI S(N) DIRECTION VECTOR. * RI F VALUE OF THE OBJECTIVE FUNCTION. * RI DF DIRECTIONAL DERIVATIVE. * RO T VALUE OF THE STEPSIZE PARAMETER. * RO TB BUNDLE PARAMETER FOR MATRIX SCALING. * RI ETA9 MAXIMUM FOR REAL NUMBERS. * SUBROUTINE PNSTP2(N,MA,MAL,X,AF,AG,AY,S,F,DF,T,TB,ETA9) DOUBLE PRECISION DF,ETA9,F,T,TB INTEGER MA,MAL,N DOUBLE PRECISION AF(*),AG(*),AY(*),S(*),X(*) DOUBLE PRECISION ALF,ALFL,ALFR,BET,BETL,BETR,Q,R,W INTEGER I,J,JN,K,L,LQ W = DF*T* (1.0D0-T*0.5D0) * * INITIAL CHOICE OF POSSIBLY ACTIVE LINES * K = 0 L = -1 JN = 0 TB = SQRT(ETA9) BETR = -ETA9 DO 20 J = 1,MAL - 1 BET = 0.0D0 ALFL = AF(J) - F DO 10 I = 1,N Q = AG(JN+I) ALFL = ALFL + (X(I)-AY(JN+I))*Q BET = BET + S(I)*Q 10 CONTINUE ALF = ABS(ALFL) R = 1.0D0 - BET/DF IF (R*R+ (ALF+ALF)/DF.GT.1.0D-6) THEN K = K + 1 AF(MA+K) = ALF AF(MA+MA+K) = BET R = T*BET - ALF IF (R.GT.W) THEN W = R L = K END IF END IF IF (BET.GT.0.0D0) TB = MIN(TB,ALF/ (BET-DF)) BETR = MAX(BETR,BET-ALF) JN = JN + N 20 CONTINUE IF (L.LT.0 .OR. BETR.LE.DF*0.5D0) RETURN LQ = 1 BETR = AF(MA+MA+L) IF (BETR.LE.0.0D0) THEN IF (T.LT.1.0D0 .OR. BETR.EQ.0.0D0) RETURN LQ = 2 END IF ALFR = AF(MA+L) * * ITERATION LOOP * 30 IF (LQ.GE.1) THEN Q = 1.0D0 - BETR/DF R = Q + SQRT(Q*Q+ (ALFR+ALFR)/DF) IF (BETR.GE.0.0D0) R = - (ALFR+ALFR)/ (DF*R) R = MIN(1.95D0,MAX(0.0D0,R)) ELSE IF (ABS(BETR-BETL)+ABS(ALFR-ALFL).LT.-1.0D-4*DF) RETURN R = (ALFR-ALFL)/ (BETR-BETL) END IF IF (ABS(T-R).LT.1.0D-4) RETURN T = R AF(MA+L) = -1.0D0 W = T*BETR - ALFR L = -1 DO 40 J = 1,K ALF = AF(MA+J) IF (ALF.LT.0.0D0) GO TO 40 BET = AF(MA+MA+J) R = T*BET - ALF IF (R.GT.W) THEN W = R L = J END IF 40 CONTINUE IF (L.LT.0) RETURN BET = AF(MA+MA+L) IF (BET.EQ.0.0D0) RETURN * * NEW INTERVAL SELECTION * ALF = AF(MA+L) IF (BET.LT.0.0D0) THEN IF (LQ.EQ.2) THEN ALFR = ALF BETR = BET ELSE ALFL = ALF BETL = BET LQ = 0 END IF ELSE IF (LQ.EQ.2) THEN ALFL = ALFR BETL = BETR LQ = 0 END IF ALFR = ALF BETR = BET END IF GO TO 30 END * SUBROUTINE PNSTP3 ALL SYSTEMS 99/12/01 * PORTABILITY : ALL SYSTEMS * 99/12/01 VL : ORIGINAL VERSION * * PURPOSE : * STEPSIZE SELECTION USING POLYHEDRAL APPROXIMATION * FOR NULL STEP IN CONVEX VARIABLE METRIC METHOD. * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * II MA DECLARED NUMBER OF LINEAR APPROXIMATED FUNCTIONS * II MAL CURRENT NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * RU X(N) VECTOR OF VARIABLES. * RI AF(4*MA) VECTOR OF BUNDLE FUNCTIONS VALUES. * RI AG(N*MA) MATRIX WHOSE COLUMNS ARE BUNDLE SUBGRADIENTS. * RI AY(N*MA) MATRIX WHOSE COLUMNS ARE VARIABLE VECTORS. * RI S(N) DIRECTION VECTOR. * RI F VALUE OF THE OBJECTIVE FUNCTION. * RI DF DIRECTIONAL DERIVATIVE. * RO T VALUE OF THE STEPSIZE PARAMETER. * RO TB BUNDLE PARAMETER FOR MATRIX SCALING. * RI ETA9 MAXIMUM FOR REAL NUMBERS. * SUBROUTINE PNSTP3(N,MA,MAL,X,AF,AG,AY,S,F,DF,T,TB,ETA9) DOUBLE PRECISION DF,ETA9,F,T,TB INTEGER MA,MAL,N DOUBLE PRECISION AF(*),AG(*),AY(*),S(*),X(*) DOUBLE PRECISION ALF,ALFL,ALFR,BET,BETL,BETR,R,TP,W INTEGER I,J,JN,K,L W = DF*T TP = T * * INITIAL CHOICE OF POSSIBLY ACTIVE LINES * K = 0 L = -1 JN = 0 TB = SQRT(ETA9) DO 20 J = 1,MAL - 1 BET = 0.0D0 ALFL = AF(J) - F DO 10 I = 1,N R = AG(JN+I) ALFL = ALFL + (X(I)-AY(JN+I))*R BET = BET + S(I)*R 10 CONTINUE ALF = ABS(ALFL) R = T*BET - ALF IF (R.GT.DF*T) THEN K = K + 1 AF(MA+K) = ALF AF(MA+MA+K) = BET IF (R.GT.W) THEN W = R L = K END IF END IF IF (BET.GT.0.0D0) TB = MIN(TB,ALF/ (BET-DF)) JN = JN + N 20 CONTINUE IF (L.LT.0) RETURN BETR = AF(MA+MA+L) IF (BETR.LE.0.0D0) RETURN ALFR = AF(MA+L) ALF = ALFR BET = BETR ALFL = 0.0D0 BETL = DF * * ITERATION LOOP * 30 IF (ABS(BETR-BETL)+ABS(ALFR-ALFL).LT.-1.0D-4*DF) RETURN R = T IF (BETR-BETL.NE.0.0D0) T = MIN((ALFR-ALFL)/ (BETR-BETL),TP) IF (ABS(T-R).LT.1.0D-3) RETURN AF(MA+L) = -1.0D0 W = T*BET - ALF L = -1 DO 40 J = 1,K ALF = AF(MA+J) IF (ALF.LT.0.0D0) GO TO 40 BET = AF(MA+MA+J) R = T*BET - ALF IF (R.GT.W) THEN W = R L = J END IF 40 CONTINUE IF (L.LT.0) RETURN * * NEW INTERVAL SELECTION * BET = AF(MA+MA+L) ALF = AF(MA+L) IF (BET.LE.0.0D0) THEN ALFL = ALF BETL = BET ELSE ALFR = ALF BETR = BET END IF GO TO 30 END * SUBROUTINE PNSTP4 ALL SYSTEMS 99/12/01 * PORTABILITY : ALL SYSTEMS * 99/12/01 VL : ORIGINAL VERSION * * PURPOSE : * STEPSIZE SELECTION USING POLYHEDRAL APPROXIMATION * FOR DESCENT STEP IN NONCONVEX VARIABLE METRIC METHOD. * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * II MA DECLARED NUMBER OF LINEAR APPROXIMATED FUNCTIONS * II MAL CURRENT NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * RU X(N) VECTOR OF VARIABLES. * RI AF(4*MA) VECTOR OF BUNDLE FUNCTIONS VALUES. * RI AG(N*MA) MATRIX WHOSE COLUMNS ARE BUNDLE SUBGRADIENTS. * RI AY(N*MA) MATRIX WHOSE COLUMNS ARE VARIABLE VECTORS. * RI F VALUE OF THE OBJECTIVE FUNCTION. * RI DF DIRECTIONAL DERIVATIVE. * RO T VALUE OF THE STEPSIZE PARAMETER. * RO TB BUNDLE PARAMETER FOR MATRIX SCALING. * RI ETA5 DISTANCE MEASURE PARAMETER. * RI ETA9 MAXIMUM FOR REAL NUMBERS. * RI MOS3 LOCALITY MEASURE PARAMETER. * SUBROUTINE PNSTP4(N,MA,MAL,X,AF,AG,AY,S,F,DF,T,TB,ETA5,ETA9,MOS3) DOUBLE PRECISION DF,ETA5,ETA9,F,T,TB INTEGER MA,MAL,MOS3,N DOUBLE PRECISION AF(*),AG(*),AY(*),S(*),X(*) DOUBLE PRECISION ALF,ALFL,ALFR,BET,BETL,BETR,DX,Q,R,W INTEGER I,J,JN,K,L,LQ W = DF*T* (1.0D0-T*0.5D0) * * INITIAL CHOICE OF POSSIBLY ACTIVE LINES * K = 0 L = -1 JN = 0 TB = SQRT(ETA9) BETR = -ETA9 DO 20 J = 1,MAL - 1 R = 0.0D0 BET = 0.0D0 ALFL = AF(J) - F DO 10 I = 1,N DX = X(I) - AY(JN+I) Q = AG(JN+I) R = R + DX*DX ALFL = ALFL + DX*Q BET = BET + S(I)*Q 10 CONTINUE IF (MOS3.NE.2) R = R** (DBLE(MOS3)*0.5D0) ALF = MAX(ABS(ALFL),ETA5*R) R = 1.0D0 - BET/DF IF (R*R+ (ALF+ALF)/DF.GT.1.0D-6) THEN K = K + 1 AF(MA+K) = ALF AF(MA+MA+K) = BET R = T*BET - ALF IF (R.GT.W) THEN W = R L = K END IF END IF IF (BET.GT.0.0D0) TB = MIN(TB,ALF/ (BET-DF)) BETR = MAX(BETR,BET-ALF) JN = JN + N 20 CONTINUE LQ = -1 IF (BETR.LE.DF*0.5D0) RETURN LQ = 1 IF (L.LT.0) RETURN BETR = AF(MA+MA+L) IF (BETR.LE.0.0D0) THEN IF (T.LT.1.0D0 .OR. BETR.EQ.0.0D0) RETURN LQ = 2 END IF ALFR = AF(MA+L) * * ITERATION LOOP * 30 IF (LQ.GE.1) THEN Q = 1.0D0 - BETR/DF R = Q + SQRT(Q*Q+ (ALFR+ALFR)/DF) IF (BETR.GE.0.0D0) R = - (ALFR+ALFR)/ (DF*R) R = MIN(1.95D0,MAX(0.0D0,R)) ELSE IF (ABS(BETR-BETL)+ABS(ALFR-ALFL).LT.-1.0D-4*DF) RETURN R = (ALFR-ALFL)/ (BETR-BETL) END IF IF (ABS(T-R).LT.1.0D-4) RETURN T = R AF(MA+L) = -1.0D0 W = T*BETR - ALFR L = -1 DO 40 J = 1,K ALF = AF(MA+J) IF (ALF.LT.0.0D0) GO TO 40 BET = AF(MA+MA+J) R = T*BET - ALF IF (R.GT.W) THEN W = R L = J END IF 40 CONTINUE IF (L.LT.0) RETURN BET = AF(MA+MA+L) IF (BET.EQ.0.0D0) RETURN * * NEW INTERVAL SELECTION * ALF = AF(MA+L) IF (BET.LT.0.0D0) THEN IF (LQ.EQ.2) THEN ALFR = ALF BETR = BET ELSE ALFL = ALF BETL = BET LQ = 0 END IF ELSE IF (LQ.EQ.2) THEN ALFL = ALFR BETL = BETR LQ = 0 END IF ALFR = ALF BETR = BET END IF GO TO 30 END * SUBROUTINE PNSTP5 ALL SYSTEMS 99/12/01 * PORTABILITY : ALL SYSTEMS * 99/12/01 VL : ORIGINAL VERSION * * PURPOSE : * STEPSIZE SELECTION USING POLYHEDRAL APPROXIMATION * FOR NULL STEP IN NONCONVEX VARIABLE METRIC METHOD. * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * II MA DECLARED NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * II MAL CURRENT NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * RU X(N) VECTOR OF VARIABLES. * RI AF(4*MA) VECTOR OF BUNDLE FUNCTIONS VALUES. * RI AG(N*MA) MATRIX WHOSE COLUMNS ARE BUNDLE SUBGRADIENTS. * RI AY(N*MA) MATRIX WHOSE COLUMNS ARE VARIABLE VECTORS. * RI S(N) DIRECTION VECTOR. * RI F VALUE OF THE OBJECTIVE FUNCTION. * RI DF DIRECTIONAL DERIVATIVE. * RO T VALUE OF THE STEPSIZE PARAMETER. * RO TB BUNDLE PARAMETER FOR MATRIX SCALING. * RI ETA5 DISTANCE MEASURE PARAMETER. * RI ETA9 MAXIMUM FOR REAL NUMBERS. * RI MOS3 LOCALITY MEASURE PARAMETER. * SUBROUTINE PNSTP5(N,MA,MAL,X,AF,AG,AY,S,F,DF,T,TB,ETA5,ETA9,MOS3) DOUBLE PRECISION DF,ETA5,ETA9,F,T,TB INTEGER MA,MAL,MOS3,N DOUBLE PRECISION AF(*),AG(*),AY(*),S(*),X(*) DOUBLE PRECISION ALF,ALFL,ALFR,BET,BETL,BETR,DX,Q,R,W INTEGER I,J,JN,K,L W = DF*T * * INITIAL CHOICE OF POSSIBLY ACTIVE PARABOLAS * K = 0 L = -1 JN = 0 TB = SQRT(ETA9) BETR = -ETA9 DO 20 J = 1,MAL - 1 BET = 0.0D0 R = 0.0D0 ALFL = AF(J) - F DO 10 I = 1,N DX = X(I) - AY(JN+I) R = R + DX*DX Q = AG(JN+I) ALFL = ALFL + DX*Q BET = BET + S(I)*Q 10 CONTINUE IF (MOS3.NE.2) R = R** (DBLE(MOS3)*0.5D0) ALF = MAX(ABS(ALFL),ETA5*R) IF (BET+BET.GT.DF) TB = MIN(TB,ALF/ (BET-DF)) BETR = MAX(BETR,BET-ALF) IF (ALF.LT.BET-DF) THEN K = K + 1 R = T*BET - ALF AF(MA+K) = ALF AF(MA+MA+K) = BET IF (R.GT.W) THEN W = R L = K END IF END IF JN = JN + N 20 CONTINUE IF (L.LT.0) RETURN BETR = AF(MA+MA+L) ALFR = AF(MA+L) ALF = ALFR BET = BETR ALFL = 0.0D0 BETL = DF * * ITERATION LOOP * 30 W = BET/DF IF (ABS(BETR-BETL)+ABS(ALFR-ALFL).LT.-1.0D-4*DF) RETURN IF (BETR-BETL.EQ.0.0D0) STOP 11 R = (ALFR-ALFL)/ (BETR-BETL) IF (ABS(T-W).LT.ABS(T-R)) R = W Q = T T = R IF (ABS(T-Q).LT.1.0D-3) RETURN AF(MA+L) = -1.0D0 W = T*BET - ALF L = -1 DO 40 J = 1,K ALF = AF(MA+J) IF (ALF.LT.0.0D0) GO TO 40 BET = AF(MA+MA+J) R = T*BET - ALF IF (R.GT.W) THEN W = R L = J END IF 40 CONTINUE IF (L.LT.0) RETURN BET = AF(MA+MA+L) Q = BET - T*DF IF (Q.EQ.0.0D0) RETURN * * NEW INTERVAL SELECTION * ALF = AF(MA+L) IF (Q.LT.0.0D0) THEN ALFL = ALF BETL = BET ELSE ALFR = ALF BETR = BET END IF GO TO 30 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) DOUBLE PRECISION EPS7,EPS9,GMAX,R,UMAX INTEGER IER,INEW,ITERM,KBC,KBF,N,NC,NF DOUBLE PRECISION CF(*),CFD(*),CG(*),CL(*),CR(*),CU(*),CZ(*),H(*), + S(*),X(*),XL(*),XU(*) INTEGER IC(*),ICA(*),IX(*) INTEGER NADD,NDECF,NFG,NFH,NFV,NIT,NRED,NREM,NRES DOUBLE PRECISION DEN,TEMP INTEGER I,IJ,IK,J,K,KC,KJ,KK,L,LL COMMON /STAT/NDECF,NRES,NRED,NREM,NADD,NIT,NFV,NFG,NFH IF (KBC.GT.0) THEN IF (R.NE.0.0D0) CALL PLDIRL(NC,CF,CFD,IC,R,KBC) IF (INEW.NE.0) THEN IF (KBF.GT.0) THEN DO 10 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 END IF 10 CONTINUE END IF DO 20 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 END IF 20 CONTINUE END IF ELSE IF (KBF.GT.0) THEN K = 0 DO 70 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.0D0) THEN IJ = 0 KJ = KK DO 40 J = 1,N IF (J.LE.K) THEN KJ = KJ + 1 ELSE KJ = KJ + J - 1 END IF IF (J.NE.K) TEMP = H(KJ)/DEN IK = KK DO 30 I = 1,J IF (I.LE.K) THEN IK = IK + 1 ELSE IK = IK + I - 1 END IF IJ = IJ + 1 IF (I.NE.K .AND. J.NE.K) H(IJ) = H(IJ) + + TEMP*H(IK) 30 CONTINUE 40 CONTINUE END IF LL = KK + K DO 60 I = K + 1,N DO 50 J = 1,I LL = LL + 1 IF (J.NE.K) THEN KK = KK + 1 H(KK) = H(LL) END IF 50 CONTINUE 60 CONTINUE N = N - 1 END IF 70 CONTINUE END IF RETURN END * SUBROUTINE PYAGB1 ALL SYSTEMS 99/12/01 * PORTABILITY : ALL SYSTEMS * 99/12/01 VL : ORIGINAL VERSION * * PURPOSE : * SUBGRADIENT AGGREGATION FOR NONSMOOTH VARIABLE METRIC METHOD. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * RI H(N*(N+1)/2) POSITIVE DEFINITE APPROXIMATION OF THE INVERSE * HESSIAN MATRIX. * RI G(NF) SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RI GO(NF) PREVIOUS SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RU GV(NF) AGGREGATED SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RU GN(N) REDUCED AGGREGATED SUBGRADIENT. * RI SN(N) REDUCED DIRECTION VECTOR. * RI CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. * RU S(N) REDUCED SUBGRADIENT. * RU U(N) PREVIOUS REDUCED SUBGRADIENT. * RO ALF LINEARIZATION TERM. * RU ALFV AGGREGATED LINEARIZATION TERM. * 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. * SUBROUTINE PYAGB1(NF,N,IX,H,G,GO,GV,GN,SN,CZ,S,U,ALF,ALFV,KBF,KBC) DOUBLE PRECISION ALF,ALFV INTEGER KBC,KBF,N,NF DOUBLE PRECISION CZ(*),G(*),GN(*),GO(*),GV(*),H(*),S(*),SN(*),U(*) INTEGER IX(*) DOUBLE PRECISION A,ALFM,B,LAM1,LAM2,PQ,PR,PRQR,QQP,QR,RR,RRP,RRQ, + W,W1,W2 INTEGER I,J,K,L ALFM = 0.0D0 * * GENERAL ROUTINE - HERE ALWAYS INPUT PARAMETER ALFM=0 * RR = ALFV + ALFV RRP = ALFV - ALFM RRQ = ALFV - ALF PQ = 0.0D0 PR = 0.0D0 QR = 0.0D0 QQP = ALF - ALFM PRQR = 0.0D0 IF (KBC.GT.0) THEN CALL MXDRMM(NF,N,CZ,G,S) CALL MXDRMM(NF,N,CZ,GO,U) DO 10 I = 1,N RR = RR - SN(I)*GN(I) U(I) = U(I) - GN(I) GN(I) = S(I) - GN(I) RRP = RRP + SN(I)*U(I) RRQ = RRQ + SN(I)*GN(I) 10 CONTINUE DO 40 I = 1,N L = I* (I-1)/2 + 1 W1 = 0.0D0 W2 = 0.0D0 DO 20 J = 1,I - 1 W = H(L) W1 = W1 + W*GN(J) W2 = W2 + W*U(J) L = L + 1 20 CONTINUE DO 30 J = I,N W = H(L) W1 = W1 + W*GN(J) W2 = W2 + W*U(J) L = L + J 30 CONTINUE PR = PR + U(I)*W2 QR = QR + GN(I)*W1 PQ = PQ + (GN(I)-U(I))* (W1-W2) QQP = QQP + S(I)* (W1-W2) PRQR = PRQR + U(I)*W1 40 CONTINUE ELSE IF (KBF.GT.0) THEN K = 0 DO 50 I = 1,NF IF (IX(I).GE.0) THEN K = K + 1 RR = RR - S(I)*GV(I) U(K) = GO(I) - GV(I) W = G(I) - GV(I) RRP = RRP + S(I)*U(K) RRQ = RRQ + S(I)*W S(K) = W END IF 50 CONTINUE K = 0 DO 80 I = 1,NF IF (IX(I).GE.0) THEN K = K + 1 L = K* (K-1)/2 + 1 W1 = 0.0D0 W2 = 0.0D0 DO 60 J = 1,K - 1 W = H(L) W1 = W1 + W*S(J) W2 = W2 + W*U(J) L = L + 1 60 CONTINUE DO 70 J = K,N W = H(L) W1 = W1 + W*S(J) W2 = W2 + W*U(J) L = L + J 70 CONTINUE PR = PR + U(K)*W2 QR = QR + S(K)*W1 PQ = PQ + (S(K)-U(K))* (W1-W2) QQP = QQP + G(I)* (W1-W2) PRQR = PRQR + U(K)*W1 END IF 80 CONTINUE ELSE DO 90 I = 1,NF RR = RR - S(I)*GV(I) U(I) = GO(I) - GV(I) GV(I) = G(I) - GV(I) RRP = RRP + S(I)*U(I) RRQ = RRQ + S(I)*GV(I) 90 CONTINUE DO 120 I = 1,NF L = I* (I-1)/2 + 1 W1 = 0.0D0 W2 = 0.0D0 DO 100 J = 1,I - 1 W = H(L) W1 = W1 + W*GV(J) W2 = W2 + W*U(J) L = L + 1 100 CONTINUE DO 110 J = I,NF W = H(L) W1 = W1 + W*GV(J) W2 = W2 + W*U(J) L = L + J 110 CONTINUE PR = PR + U(I)*W2 QR = QR + GV(I)*W1 PQ = PQ + (GV(I)-U(I))* (W1-W2) QQP = QQP + G(I)* (W1-W2) PRQR = PRQR + U(I)*W1 120 CONTINUE END IF IF (PR.LE.0.0D0 .OR. QR.LE.0.0D0) GO TO 130 A = RRQ/QR B = PRQR/QR W = PRQR*B - PR IF (W.EQ.0.0D0) GO TO 130 LAM1 = (A*PRQR-RRP)/W LAM2 = A - LAM1*B IF (LAM1* (LAM1-1.0D0).LT.0.0D0 .AND. + LAM2* (LAM1+LAM2-1.0D0).LT.0.0D0) GO TO 140 * * MINIMUM ON THE BOUNDARY * 130 LAM1 = 0.0D0 LAM2 = 0.0D0 IF (ALF.LE.ALFV) LAM2 = 1.0D0 IF (QR.GT.0.0D0) LAM2 = MIN(1.0D0,MAX(0.0D0,RRQ/QR)) W = (LAM2*QR-RRQ-RRQ)*LAM2 A = 0.0D0 IF (ALFM.LE.ALFV) A = 1.0D0 IF (PR.GT.0.0D0) A = MIN(1.0D0,MAX(0.0D0,RRP/PR)) B = (A*PR-RRP-RRP)*A IF (B.LT.W) THEN W = B LAM1 = A LAM2 = 0.0D0 END IF IF (QQP* (QQP-PQ).GE.0.0D0) GO TO 140 IF (QR-RRQ-RRQ-QQP*QQP/PQ.GE.W) GO TO 140 LAM1 = QQP/PQ LAM2 = 1.0D0 - LAM1 140 IF (LAM1.EQ.0.0D0 .AND. LAM2* (LAM2-1.0D0).LT.0.0D0 .AND. + RRP-LAM2*PRQR.GT.0.0D0 .AND. PR.GT. + 0.0D0) LAM1 = MIN(1.0D0-LAM2, (RRP-LAM2*PRQR)/PR) A = 1.0D0 - LAM1 - LAM2 B = 1.0D0 - LAM2 IF (KBC.GT.0) THEN DO 150 I = 1,N GN(I) = S(I) + LAM1*U(I) - B*GN(I) 150 CONTINUE DO 160 I = 1,NF GV(I) = LAM1*GO(I) + LAM2*G(I) + A*GV(I) 160 CONTINUE ELSE IF (KBF.GT.0) THEN K = 0 DO 170 I = 1,NF GV(I) = LAM1*GO(I) + LAM2*G(I) + A*GV(I) IF (IX(I).GE.0) THEN K = K + 1 GN(K) = GV(I) END IF 170 CONTINUE ELSE DO 180 I = 1,NF GV(I) = LAM1*GO(I) + (1.0D0-LAM1)*G(I) - A*GV(I) 180 CONTINUE END IF ALFV = LAM1*ALFM + LAM2*ALF + A*ALFV RETURN END * SUBROUTINE PYAGB2 ALL SYSTEMS 99/12/01 * PORTABILITY : ALL SYSTEMS * 99/12/01 VL : ORIGINAL VERSION * * PURPOSE : * SIMPLIFIED AGGREGATION FOR NONSMOOTH VARIABLE METRIC METHOD. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * RI H(M) POSITIVE DEFINITE APPROXIMATION OF THE INVERSE HESSIAN * MATRIX. * RI G(N) SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RU GV(N) AGGREGATED SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RU GN(N) REDUCED AGGREGATED SUBGRADIENT. * RI SN(N) REDUCED DIRECTION VECTOR. * RI CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. * RU S(N) REDUCED SUBGRADIENT. * RO ALF LINEARIZATION TERM. * RU ALFV AGGREGATED LINEARIZATION TERM. * 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. * SUBROUTINE PYAGB2(NF,N,IX,H,G,GV,GN,SN,CZ,S,ALF,ALFV,KBF,KBC) DOUBLE PRECISION ALF,ALFV INTEGER KBC,KBF,N,NF DOUBLE PRECISION CZ(*),G(*),GN(*),GV(*),H(*),S(*),SN(*) INTEGER IX(*) DOUBLE PRECISION LAM,P,Q,W INTEGER I,J,K,L P = ALFV - ALF IF (KBC.GT.0) THEN CALL MXDRMM(NF,N,CZ,G,S) DO 10 I = 1,N GN(I) = S(I) - GN(I) P = P + SN(I)*GN(I) 10 CONTINUE Q = 0.0D0 DO 40 I = 1,N L = I* (I-1)/2 + 1 W = 0.0D0 DO 20 J = 1,I - 1 W = W + H(L)*GN(J) L = L + 1 20 CONTINUE DO 30 J = I,N W = W + H(L)*GN(J) L = L + J 30 CONTINUE Q = Q + W*GN(I) 40 CONTINUE LAM = 0.5D0 + SIGN(0.5D0,P) IF (Q.GT.0.0D0) LAM = MIN(1.0D0,MAX(0.0D0,P/Q)) P = 1.0D0 - LAM DO 50 I = 1,N GN(I) = S(I) - P*GN(I) 50 CONTINUE DO 60 I = 1,NF GV(I) = LAM*G(I) + P*GV(I) 60 CONTINUE ELSE IF (KBF.GT.0) THEN K = 0 DO 70 I = 1,NF IF (IX(I).GE.0) THEN K = K + 1 Q = G(I) - GV(I) P = P + S(I)*Q S(K) = Q END IF 70 CONTINUE Q = 0.0D0 K = 0 DO 100 I = 1,NF IF (IX(I).GE.0) THEN K = K + 1 L = K* (K-1)/2 + 1 W = 0.0D0 DO 80 J = 1,K - 1 W = W + H(L)*S(J) L = L + 1 80 CONTINUE DO 90 J = K,N W = W + H(L)*S(J) L = L + J 90 CONTINUE Q = Q + W*S(K) END IF 100 CONTINUE LAM = 0.5D0 + SIGN(0.5D0,P) IF (Q.GT.0.0D0) LAM = MIN(1.0D0,MAX(0.0D0,P/Q)) P = 1.0D0 - LAM K = 0 DO 110 I = 1,NF GV(I) = LAM*G(I) + P*GV(I) IF (IX(I).GE.0) THEN K = K + 1 GN(K) = GV(I) END IF 110 CONTINUE ELSE DO 120 I = 1,NF GV(I) = G(I) - GV(I) P = P + S(I)*GV(I) 120 CONTINUE Q = 0.0D0 DO 150 I = 1,NF L = I* (I-1)/2 + 1 W = 0.0D0 DO 130 J = 1,I - 1 W = W + H(L)*GV(J) L = L + 1 130 CONTINUE DO 140 J = I,NF W = W + H(L)*GV(J) L = L + J 140 CONTINUE Q = Q + W*GV(I) 150 CONTINUE LAM = 0.5D0 + SIGN(0.5D0,P) IF (Q.GT.0.0D0) LAM = MIN(1.0D0,MAX(0.0D0,P/Q)) P = 1.0D0 - LAM DO 160 I = 1,NF GV(I) = G(I) - P*GV(I) 160 CONTINUE END IF ALFV = LAM*ALF + P*ALFV RETURN END * SUBROUTINE PYBUN1 ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 VL : ORIGINAL VERSION * * PURPOSE : * VARIABLE METRIC UPDATE OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX * WITH THE POSSIBILITY OF MATRIX INNOVATION. * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * II MA DECLARED NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * II MAL CURRENT NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * RI X(N) VECTOR OF VARIABLES. * RI G(N) SUBGRADIENT OF THE OBJECTIVE FUNCTION. * RI F VALUE OF THE OBJECTIVE FUNCTION. * RU AY(N*MA) MATRIX WHOSE COLUMNS ARE VARIABLE VECTORS. * RU AG(N*MA) MATRIX WHOSE COLUMNS ARE BUNDLE SUBGRADIENTS. * RU AF(4*MA) VECTOR OF VALUES OF BUNDLE FUNCTIONS. * IO ITERS NULL STEP INDICATOR. ITERS=0-NULL STEP. ITERS=1-DESCENT * STEP. * * SUBPROGRAMS USED : * RF MXVDOT DOT PRODUCT OF VECTORS. * SUBROUTINE PYBUN1(N,MA,MAL,X,G,F,AY,AG,AF,ITERS) DOUBLE PRECISION F INTEGER ITERS,MA,MAL,N DOUBLE PRECISION AF(*),AG(*),AY(*),G(*),X(*) INTEGER I,IND,K,KN,L L = 0 IF (ITERS.EQ.0) L = 1 * * BUNDLE REDUCTION * KN = 0 IF (MAL.GE.MA) THEN DO 20 K = 1,MAL - 1 KN = K*N - N DO 10 I = 1,N IF (G(I).NE.AG(KN+I)) GO TO 20 10 CONTINUE IND = K GO TO 30 20 CONTINUE IND = 1 30 DO 40 K = IND,MAL - 1 AF(K) = AF(K+1) AF(K+MA*3) = AF(K+1+MA*3) KN = K*N + 1 CALL MXVCOP(N,AG(KN),AG(KN-N)) CALL MXVCOP(N,AY(KN),AY(KN-N)) 40 CONTINUE MAL = MAL - 1 END IF * * BUNDLE COMPLETION * IF (L.GT.0 .AND. KN.EQ.0) THEN AF(MAL+1) = AF(MAL) AF(3*MA+MAL+1) = AF(3*MA+MAL) KN = MAL*N + 1 CALL MXVCOP(N,AG(KN-N),AG(KN)) CALL MXVCOP(N,AY(KN-N),AY(KN)) END IF MAL = MAL + 1 KN = MAL - L AF(KN) = F AF(KN+MA*3) = L K = (KN-1)*N + 1 CALL MXVCOP(N,G,AG(K)) CALL MXVCOP(N,X,AY(K)) 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) DOUBLE PRECISION EPS8,GMAX,UMAX INTEGER IER,IOLD,ITERM,KBC,KBF,KOLD,KREM,N,NF DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),H(*) INTEGER IC(*),ICA(*),IX(*) INTEGER NADD,NDECF,NFG,NFH,NFV,NIT,NRED,NREM,NRES INTEGER I,J,K,KC,L COMMON /STAT/NDECF,NRES,NRED,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.0D0,H(K+1)) H(K+N) = 1.0D0 KC = ICA(NF-N+1) IF (KC.GT.0) THEN IC(KC) = -IC(KC) ELSE K = -KC IX(K) = -IX(K) END IF END IF ELSE IOLD = 0 END IF ELSE IF (KBF.GT.0) THEN IF (UMAX.GT.EPS8*GMAX) THEN IX(IOLD) = MIN(ABS(IX(IOLD)),3) DO 10 I = N,KOLD,-1 GN(I+1) = GN(I) 10 CONTINUE GN(KOLD) = G(IOLD) N = N + 1 K = N* (N-1)/2 L = K + N DO 30 I = N,KOLD,-1 DO 20 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.0D0 L = L - 1 ELSE H(L) = 0.0D0 L = L - 1 END IF 20 CONTINUE 30 CONTINUE ELSE IOLD = 0 KOLD = 0 END IF END IF 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) DOUBLE PRECISION DMAX,F,FO,P,PO,R INTEGER ITERS,KBC,KBF,N,NF DOUBLE PRECISION CZ(*),G(*),GO(*),SN(*),X(*),XO(*) INTEGER IX(*) 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) END IF DMAX = 0.0D0 IF (KBC.GT.0) THEN DO 10 I = 1,NF DMAX = MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D0)) 10 CONTINUE IF (N.GT.0) THEN CALL MXVSCL(N,R,SN,XO) CALL MXVCOP(NF,GO,SN) CALL MXDRMM(NF,N,CZ,SN,GO) END IF ELSE IF (KBF.GT.0) THEN K = 0 DO 20 I = 1,NF IF (IX(I).LT.0) GO TO 20 K = K + 1 DMAX = MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D0)) XO(K) = XO(I) GO(K) = GO(I) 20 CONTINUE ELSE DO 30 I = 1,NF DMAX = MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D0)) 30 CONTINUE END IF 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) DOUBLE PRECISION GMAX,UMAX INTEGER IOLD,KBC,KBF,KOLD,N,NF DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*) INTEGER IC(*),ICA(*),IX(*) DOUBLE PRECISION TEMP INTEGER I,J,K,KC,NCA,NCZ DOUBLE PRECISION MXVDOT,MXVMAX IOLD = 0 KOLD = 0 UMAX = 0.0D0 GMAX = 0.0D0 IF (KBC.GT.0) THEN IF (NF.GT.N) THEN NCA = NF - N NCZ = N*NF CALL MXVCOP(NF,G,GN) DO 10 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) END IF 10 CONTINUE CALL MXDPRB(NCA,CR,CZ(NCZ+1),0) DO 20 J = 1,NCA TEMP = CZ(NCZ+J) KC = ICA(J) IF (KC.GT.0) THEN K = IC(KC) ELSE I = -KC K = IX(I) END IF IF (K.LE.-5) THEN ELSE IF ((K.EQ.-1.OR.K.EQ.-3) .AND. + UMAX+TEMP.GE.0.0D0) THEN ELSE IF ((K.EQ.-2.OR.K.EQ.-4) .AND. + UMAX-TEMP.GE.0.0D0) THEN ELSE IOLD = J UMAX = ABS(TEMP) END IF 20 CONTINUE END IF IF (N.GT.0) THEN CALL MXDRMM(NF,N,CZ,G,GN) GMAX = MXVMAX(N,GN) END IF ELSE IF (KBF.GT.0) THEN J = 0 IOLD = 0 KOLD = 0 DO 30 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.0D0) THEN ELSE IF ((K.EQ.-2.OR.K.EQ.-4) .AND. + UMAX-TEMP.GE.0.0D0) THEN ELSE IOLD = I KOLD = J + 1 UMAX = ABS(TEMP) END IF 30 CONTINUE N = J ELSE DO 40 I = 1,NF TEMP = G(I) GMAX = MAX(GMAX,ABS(TEMP)) 40 CONTINUE N = NF END IF 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) DOUBLE PRECISION F,FO,FP,P,PO,RMAX,RO INTEGER INEW,KBC,KBF,KREM,N,NC,NF DOUBLE PRECISION CF(*),CFD(*),CG(*),CL(*),CU(*),CZ(*),G(*),GO(*), + S(*),SN(*),X(*),XL(*),XO(*),XU(*) INTEGER IC(*),IX(*) INTEGER I,K FP = FO RO = 0.0D0 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.0D0,S) END IF ELSE IF (KBF.GT.0) THEN K = N + 1 DO 10 I = NF,1,-1 IF (IX(I).LT.0) THEN S(I) = 0.0D0 ELSE K = K - 1 S(I) = SN(K) END IF 10 CONTINUE INEW = 0 CALL PLMAXS(NF,X,IX,XL,XU,S,RMAX,KBF,KREM,INEW) END IF 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