************************************************************************ * SUBROUTINE PMAXU ALL SYSTEMS 97/01/22 * PURPOSE : * EASY TO USE SUBROUTINE FOR LARGE-SCALE UNCONSTRAINED MINIMAX * OPTIMIZATION WITH SPARSE JACOBIAN MATRICES. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NA NUMBER OF PARTIAL FUNCTIONS. * IU MA NUMBER OF NONZERO ELEMENTS IN THE JACOBIAN MATRIX. * RI X(NF) VECTOR OF VARIABLES. * RO AF(NA) VECTOR CONTAINING VALUES OF PARTIAL FUNCTIONS. * RI IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE JACOBIAN * MATRIX. * RI JAG(MA) COLUMN INDICES OF ELEMENTS IN THE JACOBIAN MATRIX. * II IPAR(7) INTEGER PAREMETERS: * IPAR(1) MAXIMUM NUMBER OF ITERATIONS. * IPAR(2) MAXIMUM NUMBER OF FUNCTION EVALUATIONS. * IPAR(3) MAXIMUM NUMBER OF GRADIENT EVALUATIONS. * IPAR(4) ESTIMATION INDICATOR. IPAR(4)=0-MINIMUM IS NOT * ESTIMATED. IPAR(4)=1-MINIMUM IS ESTIMATED BY THE VALUE * RPAR(6). * IPAR(5) METHOD USED. IPAR(5)=1-PARTITIONED VARIABLE METRIC * METHOD. IPAR(5)=2-DISCRETE NEWTON METHOD. * IPAR(6) THIS PARAMETER IS NOT USED IN THE SUBROUTINE PMAX. * IPAR(7) NUMBER DEFINING THE SPACE FOR FILL-IN (THE SIZE OF THIS * SPACE IS IFIL TIMES THE STANDARD SIZE). THE DEFAULT VALUE IS * IFIL=1. THE DEFAULT VALUE HAS TO BE INCREASED IF ITERM IS * LESS OR EQUAL TO -40. * RI RPAR(9) REAL PARAMETERS: * RPAR(1) MAXIMUM STEPSIZE. * RPAR(2) TOLERANCE FOR THE CHANGE OF VARIABLES. * RPAR(3) TOLERANCE FOR THE CHANGE OF FUNCTION VALUES. * RPAR(4) TOLERANCE FOR THE FUNCTION FALUE. * RPAR(5) TOLERANCE FOR THE GRADIENT NORM. * RPAR(6) ESTIMATION OF THE MINIMUM FUNCTION VALUE. * RPAR(7) THIS PARAMETER IS NOT USED IN THE SUBROUTINE PMAX. * RPAR(8) COEFFICIENT FOR THE BARRIER PARAMETER DECREASE. * RPAR(9) MINIMUM PERMITTED VALUE OF THE BARRIER PARAMETER. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RO GMAX MAXIMUM PARTIAL DERIVATIVE. * 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 ISPAS INPUT SPARSE STRUCTURE. ISPAS=1-STANDARD COORDINATE * FORM. ISPAS=2-SPARSE STRUCTURE COMPRESSED BY ROWS. * II IPRNT PRINT SPECIFICATION. IPRNT=0-NO PRINT. * ABS(IPRNT)=1-PRINT OF FINAL RESULTS. * ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS. * IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL * RESULTS. * IO ITERM VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. * ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN * MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS. * ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN * MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS. * ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB. * ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG. * ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, * BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. * ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. * ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. * VALUES ITERM<=-40 DETECT A LACK OF SPACE. IN THIS CASE, * PARAMETER IPAR(7) HAS TO BE INCREASED. * * VARIABLES IN COMMON /STAT/ (STATISTICS) : * IO NRES NUMBER OF RESTARTS. * IO NDEC NUMBER OF MATRIX DECOMPOSITIONS. * IO NIN NUMBER OF INNER ITERATIONS. * IO NIT NUMBER OF ITERATIONS. * IO NFV NUMBER OF FUNCTION EVALUATIONS. * IO NFG NUMBER OF GRADIENT EVALUATIONS. * IO NFH NUMBER OF HESSIAN EVALUATIONS. * * SUBPROGRAMS USED : * S PMAX PRIMAL LINE-SEARCH INTERIOR-POINT METHOD FOR LARGE-SCALE * PARTIALLY SEPARABLE MINIMAX OPTIMIZATION PROBLEMS. * S PASED3 COMPRESSED SPARSE STRUCTURE OF THE JACOBIAN MATRIX IS * COMPUTED FROM THE COORDINATE FORM. * S PFSET2 NUMBER OF NONZERO ELEMENTS IN THE PARTITIONED HESSIAN * MATRIX. * * EXTERNAL SUBROUTINES : * SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION. * CALLING SEQUENCE: CALL FUN(NF,KA,X,FA) WHERE NF IS A NUMBER * OF VARIABLES, KA IS THE INDEX OF THE APPROXIMATED FUNCTION, * X(NF) IS A VECTOR OF VARIABLES AND FA IS THE VALUE OF THE * APPROXIMATED FUNCTION. * SE DFUN COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION. * CALLING SEQUENCE: CALL DFUN(NF,KA,X,GA) WHERE NF IS A NUMBER * OF VARIABLES, KA IS THE INDEX OF THE APPROXIMATED FUNCTION, * X(NF) IS A VECTOR OF VARIABLES AND GA(NF) IS THE GRADIENT OF * THE APPROXIMATED FUNCTION. * SUBROUTINE PMAXU(NF,NA,MA,X,AF,IAG,JAG,IPAR,RPAR,F,GMAX,IEXT, & ISPAS,IPRNT,ITERM) INTEGER NF,NA,MA,IAG(*),JAG(*),IPAR(7),IEXT,ISPAS,IPRNT,ITERM DOUBLE PRECISION X(*),AF(*),RPAR(9),F,GMAX INTEGER LAFO,LAG,LAGO,LGA,LAH,LAZL,LAZU,LG,LH,LS,LXO,LGO,LGS,LGP, & LIH,LJH,LCOL,LPSL,LPERM,LINVP,LWN11,LWN12,LWN13,LWN14,MB,MC INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH,IFIL,IER COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH INTEGER IA(:) DOUBLE PRECISION RA(:) ALLOCATABLE IA,RA IFIL=IPAR(7) IF (IFIL.LE.0) IFIL=1 IF (ISPAS.LE.1) THEN CALL PASED3(NF,NA,MA,IAG,JAG,IER) IF (IER.NE.0) THEN WRITE (6,'(''INPUT ERROR : IER = '',I3)') IER STOP END IF ELSE MA=IAG(NA+1)-1 END IF CALL PFSET2(NA,MB,MC,IAG) ALLOCATE (IA(NA+9*NF+7+(IFIL+3)*MB)) IF (IPAR(5).LE.1) THEN ALLOCATE (RA(2*MA+3*NA+7*NF+7+(IFIL+4)*MB)) ELSE ALLOCATE (RA(MA+3*NA+7*NF+7+(IFIL+3)*MB)) END IF C C POINTERS FOR AUXILIARY ARRAYS C LAFO=1 LAG=LAFO+NA IF (IPAR(5).LE.1) THEN LAGO=LAG+MA LAH=LAGO+MA LGA=LAH+MB ELSE LAGO=LAG LAH=LAG LGA=LAG+MA END IF LAZL=LGA+NF LAZU=LAZL+NA LG=LAZU+NA LS=LG+NF+1 LXO=LS+NF+1 LGO=LXO+NF+1 LGS=LGO+NF+1 LGP=LGS+NF+1 LH=LGP+NF+1 LCOL=NA+1 LPSL=LCOL+NF LPERM=LPSL+NF+1 LINVP=LPERM+NF LWN11=LINVP+NF LWN12=LWN11+NF+1 LWN13=LWN12+NF+1 LWN14=LWN13+NF+1 LIH=LWN14+NF+1 LJH=LIH+NF+2 CALL PMAX(NF,NA,(IFIL+3)*MB,X,IA,AF,RA(LAFO),RA(LAG),RA(LAGO), & RA(LGA),RA(LAH),RA(LAZL),RA(LAZU),RA(LG),RA(LH),IA(LIH),IA(LJH), & IA,IAG,JAG,RA(LS),RA(LXO),RA(LGO),RA(LGS),RA(LGP),IA(LCOL), & IA(LPSL),IA(LPERM),IA(LINVP),IA(LWN11),IA(LWN12),IA(LWN13), & IA(LWN14),RPAR(1),RPAR(2),RPAR(3),RPAR(4),RPAR(5),RPAR(6), & RPAR(8),RPAR(9),GMAX,F,IPAR(1),IPAR(2),IPAR(3),IPAR(4),IPAR(5), & IEXT,IPRNT,ITERM) DEALLOCATE (IA,RA) RETURN END ************************************************************************ * SUBROUTINE PMAX ALL SYSTEMS 01/09/22 * PURPOSE : * GENERAL SUBROUTINE FOR LARGE-SCALE UNCONSTRAINED MINIMAX PROBLEMS * WITH SPARSE JACOBIAN MATRICES. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NA NUMBER OF PARTIAL FUNCTIONS. * II MMAX MAXIMUM DIMENSION OF THE SPARSE TABLEAU. * RI X(NF) VECTOR OF VARIABLES. * IA IX(NF) AUXILIARY VECTOR. * RO AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED * FUNCTIONS. * RA AFO(NA) AUXILIARY VECTOR. * RA AG(MA) JACOBIAN MATRIX OF THE PARTIALLY SEPARABLE FUNCTION. * RA AGO(NA) AUXILIARY VECTOR. * RA GA(NF) GRADIENT OF THE SELECTED PARTIAL FUNCTION. * RA AH(MB) HESSIAN MATRIX OF THE PARTIALLY SEPARABLE FUNCTION. * RA AZL(NA) VECTOR OF LOWER LAGRANGE MULTIPLIERS. * RA AZU(NA) VECTOR OF UPPER LAGRANGE MULTIPLIERS. * RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RA H(MMAX) NONZERO ELEMENTS OF THE APPROXIMATION OF THE SPARSE * HESSIAN MATRIX TOGETHER WITH AN ADDITIONAL SPACE USED FOR * THE NUMERICAL DIFFERENTIATION. * IA IH(NF+1) POINTERS OF DIAGONAL ELEMENTS OF THE MATRIX H. * IA JH(MMAX) INDICES OF NONZERO ELEMENTS OF THE MATRIX H * TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL * DIFFERENTIATION. * IA IA(NA) AUXILIARY VECTOR. * II IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. * II JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. * RA S(NF) DIRECTION VECTOR. * RA XO(NF) VECTORS OF VARIABLES DIFFERENCE. * RA GO(NF) GRADIENTS DIFFERENCE. * RA GS(NF) AUXILIARY VECTOR. * RA GP(NF) AUXILIARY VECTOR. * IA COL(NF) AUXILIARY ARRAY. * IA PSL(NF+1) POINTER VECTOR OF THE COMPACT FORM OF THE TRIANGULAR * FACTOR OF THE HESSIAN APPROXIMATION. * IA PERM(NF) PERMUTATION VECTOR. * IA INVP(NF) INVERSE PERMUTATION VECTOR. * IA WN11(NF+1) AUXILIARY VECTOR. * IA WN12(NF+1) AUXILIARY VECTOR. * IA WN13(NF+1) AUXILIARY VECTOR. * IA WN14(NF+1) AUXILIARY VECTOR. * RI XMAX MAXIMUM STEPSIZE. * RI TOLX TOLERANCE FOR CHANGE OF VARIABLES. * RI TOLF TOLERANCE FOR CHANGE OF FUNCTION VALUES. * RI TOLB TOLERANCE FOR THE FUNCTION VALUE. * RI TOLG TOLERANCE FOR THE GRADIENT NORM. * RI FMIN ESTIMATION OF THE MINIMUM FUNCTION VALUE. * RI ETA4 COEFFICIENT FOR THE BARRIER PARAMETER DECREASE. * RI ETA5 MINIMUM PERMITTED VALUE OF THE BARRIER PARAMETER. * RO GMAX MAXIMUM PARTIAL DERIVATIVE. * RO F VALUE OF THE OBJECTIVE FUNCTION. * II MIT MAXIMUM NUMBER OF ITERATIONS. * II MFV MAXIMUM NUMBER OF FUNCTION EVALUATIONS. * II MFG MAXIMUM NUMBER OF GRADIENT EVALUATIONS. * II IEST ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED. * IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN. * II MED METHOD USED. MED=1-PARTITIONED VARIABLE METRIC METHOD. * MED=2-DISCRETE NEWTON METHOD. * 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 IPRNT PRINT SPECIFICATION. IPRNT=0-NO PRINT. * ABS(IPRNT)=1-PRINT OF FINAL RESULTS. * ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS. * IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL * RESULTS. * IO ITERM VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. * ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN * MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS. * ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN * MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS. * ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB. * ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG. * ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, * BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. * ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. * ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. * VALUES ITERM<=-40 DETECT A LACK OF SPACE. * * VARIABLES IN COMMON /STAT/ (STATISTICS) : * IO NRES NUMBER OF RESTARTS. * IO NDEC NUMBER OF MATRIX DECOMPOSITIONS. * IO NIN NUMBER OF INNER ITERATIONS. * IO NIT NUMBER OF ITERATIONS. * IO NFV NUMBER OF FUNCTION EVALUATIONS. * IO NFG NUMBER OF GRADIENT EVALUATIONS. * IO NFH NUMBER OF HESSIAN EVALUATIONS. * * SUBPROGRAMS USED : * S PALNG3 EXTRACTION OF THE PARTIAL GRADIENT. * S PASSH3 MODIFICATION OF THE HESSIAN MATRIX. * S PF1HS2 NUMERICAL COMPUTATION OF THE HESSIAN MATRIX USING * DIFFERENCES OF GRADIENTS. * S PFSEB2 COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE * PARTITIONED HESSIAN MATRIX IN THE MINIMAX CASE. * S PFSET2 NUMBER OF NONZERO ELEMENTS IN THE PARTITIONED HESSIAN * MATRIX. * S PFSET3 PREPARATION OF THE SPARSE NORMAL EQUATION MATRIX * STRUCTURE. * S PNFUZ1 DETERMINATION OF THE LAGRANGE MULTIPLIERS. * S PNNEQ1 SOLUTION OF THE BASIC NONLINEAR EQUATION. * S PP0BX1 COMPUTATION OF THE VALUE OF THE BARRIER FUNCTION. * S PP1MX3 COMPUTATION OF THE VALUE AND THE GRADIENT OF THE * LAGRANGIAN FUNCTION. * S PS0L02 LINE SEARCH USING ONLY FUNCTION VALUES. * S PUBBM2 VARIABLE METRIC UPDATES OF THE PARTITIONED MATRIX. * S PYFUT8 TEST ON TERMINATION. * S PYPTSH DETERMINATION OF GROUPS FOR NUMERICAL DIFFERENTIATION. * S PYTCUB SCALED DIFFERENCE OF THE JACOBIAN MATRICES IN THE * MINIMAX CASE. * S PYTRCD COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE * METRIC UPDATE. * S PYTRCG COMPUTATION OF THE PROJECTED GRADIENT. * S PYTRCS COMPUTATION OF THE PROJECTED DIRECTION VECTOR. * S PYTSCH CORRECTION OF THE HESSIAN MATRIX. * S MXBSMI INITIATION OF THE PARTITIONED MATRIX. * S MXSPCB BACK SUBSTITUTION USING THE SPARSE DECOMPOSITION * OBTAINED BY MXSPCF. * S MXSPCC SPARSE MATRIX REORDERING, SYMBOLIC FACTORIZATION, DATA * STRUCTURES TRANSFORMATION. INITIATION OF THE DIRECT SPARSE * SOLVER. * S MXSPCF GILL-MURRAY DECOMPOSITION OD A SPARSE SYMMETRIC MATRIX. * S MXSPCT COPYING A SPARSE SYMMETRIC MATRIX INTO THE PERMUTED * FACTORIZED COMPACT SCHEME. * S MXVCOP COPYING OF A VECTOR. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * S MXVINE RESTORATION OF A SPARSE SYMMETRIC MATRIX OBTAINED BY * MXVINB * S MXVINS INITIATION OF THE INTEGER VECTOR. * S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. * S MXVSBP INVERSE PERMUTATION OF A VECTOR * S MXVSET INITIATION OF A VECTOR. * S MXVSFP PERMUTATION OF A VECTOR. * * EXTERNAL SUBROUTINES : * SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION. * CALLING SEQUENCE: CALL FUN(NF,KA,X,FA) WHERE NF IS A NUMBER * OF VARIABLES, KA IS THE INDEX OF THE APPROXIMATED FUNCTION, * X(NF) IS A VECTOR OF VARIABLES AND FA IS THE VALUE OF THE * APPROXIMATED FUNCTION. * SE DFUN COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION. * CALLING SEQUENCE: CALL DFUN(NF,KA,X,GA) WHERE NF IS A NUMBER * OF VARIABLES, KA IS THE INDEX OF THE APPROXIMATED FUNCTION, * X(NF) IS A VECTOR OF VARIABLES AND GA(NF) IS THE GRADIENT OF * THE APPROXIMATED FUNCTION. * * METHOD : * PRIMAL DISCRETE NEWTON INTERIOR-POINT ALGORITHM FOR LARGE-SCALE * PARTIALLY SEPARABLE MINIMAX OPTIMIZATION PROBLEMS. * SUBROUTINE PMAX(NF,NA,MMAX,X,IX,AF,AFO,AG,AGO,GA,AH,AZL,AZU,G,H, & IH,JH,IA,IAG,JAG,S,XO,GO,GS,GP,COL,PSL,PERM,INVP,WN11,WN12,WN13, & WN14,XMAX,TOLX,TOLF,TOLB,TOLG,FMIN,ETA4,ETA5,GMAX,F,MIT,MFV,MFG, & IEST,MED,IEXT,IPRNT,ITERM) INTEGER NA,NF,MMAX,IX(*),IH(*),JH(*),IA(*),IAG(*),JAG(*),COL(*), & PSL(*),PERM(*),INVP(*),WN11(*),WN12(*),WN13(*),WN14(*),MIT,MFV, & MFG,IEST,MED,IEXT,IPRNT,ITERM DOUBLE PRECISION X(*),AF(*),AFO(*),AG(*),AGO(*),GA(*),AH(*), & AZL(*),AZU(*),G(*),H(*),S(*),XO(*),GO(*),GS(*),GP(*),XMAX,TOLX, & TOLF,TOLB,TOLG,FMIN,ETA4,ETA5,GMAX,F INTEGER IDECF,ITERD,ITERS,KD,LD,NTESX,NTESF,MTESX,MTESF,MRED,KIT, & IREST,KBF,MAXST,IDIR,IOLD,INF,INITD,MEP,MER,MET,MET1,MET3,MET5, & IER,ICON,ISYS,KTERS,ITERH,IRES1,IRES2,NRED,I,J,JP,K,L,INITS,MES, & M,MA,MB,MM,MH,NNIT,JSTRT,JSTOP,KA,ISNA DOUBLE PRECISION R,RO,RP,FF,FO,FP,FA,P,PO,PP,GNORM,SNORM,RMIN, & RMAX,FMAX,DMAX,UMAX,ETA0,ETA2,ETA3,ETA6,ETA9,EPS0,EPS1,ALF,ALF1, & ALF2,BET,RHO,TOLP,RPF3,PAR,FFO DOUBLE PRECISION MXVDOT,PNFUZ1 INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH IF (ABS(IPRNT).GT.1) WRITE(6,'(1X,''ENTRY TO PMAX :'')') * * INITIATION OF PROBLEM * KBF=0 NRES=0 NDEC=0 NIN=0 NIT=0 NFV=0 NFG=0 NFH=0 ICON=0 ISYS=0 NTESX=0 NTESF=0 MTESX=2 MTESF=2 INITS=1 INITD=1 ITERM=0 ITERD=0 ITERS=2 KTERS=5 IREST=0 IRES1=999 IRES2=0 MRED=20 IDIR=0 MEP=1 MES=1 ETA0=1.0D-15 ETA2=1.0D-18 ETA3=1.0D-6 IF (ETA4.LE.0.0D 0) ETA4=8.5D-1 ETA6=1.0D 0 ETA9=1.0D 120 EPS0=1.0D-8 EPS1=1.0D-4 ALF1=1.0D-10 ALF2=1.0D 10 RPF3=1.0D 0 RMAX=ETA9 DMAX=ETA9 FMAX=1.0D 20 IF (IEXT.EQ.0) THEN IF (IEST.LE.0) FMIN=0.0D 0 FMIN=MAX(FMIN,0.0D 0) IEST=1 ELSE IEST=0 END IF IF (XMAX.LE.0.0D 0) XMAX=1.0D 16 IF (TOLX.LE.0.0D 0) TOLX=1.0D-16 IF (TOLF.LE.0.0D 0) TOLF=1.0D-14 IF (TOLB.LE.0.0D 0) TOLB=FMIN+1.0D-16 IF (TOLG.LE.0.0D 0) TOLG=1.0D-6 IF (MIT.LE.0) MIT=10000 IF (MFV.LE.0) MFV=10000 IF (MFG.LE.0) MFG=20000 IF (MED.LE.0) MED=1 IF (MED.EQ.1) THEN MER=2 MET=1 MET1=3 MET3=1 MET5=1 KIT=-(IRES1*NF+IRES2) CALL PFSET2(NA,MB,MA,IAG) MA=IAG(NA+1)-1 ELSE MED=2 MER=2 KIT=0 END IF CALL MXVINP(NF+1,IH) CALL MXVINP(NF,JH) CALL PFSET3(NF,NA,M,MMAX,IH,JH,IAG,JAG,ITERM) IF (ITERM.NE.0) GO TO 11080 CALL MXVINS(NA,3,IA) IF (IEXT.GT.0) CALL MXVINS(NA,1,IA) IF (IEXT.LT.0) CALL MXVINS(NA,2,IA) IF (MED.EQ.2) CALL PYPTSH(NF,MMAX,IH,JH,COL,S,XO,GO,WN11,WN12,GA, & ITERM) MH=0 CALL MXVINE(IH(NF+1)-1,JH) CALL MXSPCC(NF,M,MH,MMAX,H,IH,JH,PSL,PERM,INVP,WN11,WN12,WN13, & WN14,IER) IF (IER.NE.0) THEN ITERM=IER END IF * * SPARSE NEWTON METHOD * ISNA=2 KD=MED LD=-1 R=0.0D 0 FO=FMIN PAR=0.0D 0 IF (ETA5.LE.0.0D 0) THEN TOLP=ETA0**(2.0D 0/3.0D 0) ELSE TOLP=ETA5 END IF IF (ITERM.NE.0) GO TO 11080 11010 CONTINUE * * COMPUTATION OF THE VALUE OF THE LAGRANGIAN FUNCTION * KD=0 CALL PP1MX3(NF,NA,X,GA,AG,IAG,JAG,G,AZL,AZU,FA,AF,FF,KD,LD,NFV, & NFG,ISNA,IEXT) LD=0 IF (FF+RPF3.GT.(1.0D 0+ETA0)*FF) THEN NNIT=20 11015 CONTINUE CALL PNNEQ1(FF+RPF3,FF+RPF3*DBLE(2*NA),X(NF+1),G(NF+1), & 0.0D 0,MIN(ETA3,TOLG),NNIT,IER,ICON) IF (ICON.GT.0) THEN G(NF+1)=PNFUZ1(X(NF+1),NA,RPF3,AF,AZL,AZU,IEXT) GO TO 11015 END IF ELSE X(NF+1)=(1.0D 0+ETA0)*FF G(NF+1)=PNFUZ1(X(NF+1),NA,RPF3,AF,AZL,AZU,IEXT) END IF CALL PP0BX1(NA,X(NF+1),AF,F,FF,PAR,RPF3,MEP,IEXT) * * COMPUTATION OF THE GRADIENT AND THE HESSIAN MATRIX OF THE * LAGRANGIAN FUNCTION * KD=1 CALL PP1MX3(NF,NA,X,GA,AG,IAG,JAG,G,AZL,AZU,FA,AF,FF,KD,LD,NFV, & NFG,ISNA,IEXT) LD=1 IF (MED.EQ.1) GO TO 11025 11020 CONTINUE CALL PF1HS2(NF,MH,MMAX,X,IH,S,H,IH,JH,GO,G,COL,WN11,WN12,GS, & FF,ETA0,0,ITERM,ISYS) IF (ISYS.GT.0) THEN LD=0 ISNA=0 CALL PP1MX3(NF,NA,X,GA,AG,IAG,JAG,GO,AZL,AZU,FA,AF,FF,KD,LD,NFV, & NFG,ISNA,IEXT) GO TO 11020 END IF KD=2 LD=2 ISNA=2 IDECF=0 11025 CONTINUE IF (NIT.NE.0) GO TO 11070 11030 CONTINUE CALL PYTRCG(NF,NF,IX,G,UMAX,GMAX,KBF,IOLD) IF (ABS(IPRNT).GT.1) & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X, & ''F='', G16.9,2X,''G='',E10.3)') NIT,NFV,NFG,F,GMAX CALL PYFUT8(NF,FF,FFO,GMAX,DMAX,RPF3,TOLX,TOLF,TOLB,TOLG,TOLP,KD, & NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,IRES1,IRES2, & IREST,ITERS,ITERM) IF (ITERM.NE.0) GO TO 11080 11040 CONTINUE IF (IREST.GT.0) THEN IF (MED.EQ.1) THEN CALL MXBSMI(NA,AH,IAG) ELSE RHO=GMAX/1.0D 1 DO 20 I=1,NF JSTRT=IH(I) JSTOP=IH(I+1)-1 H(JSTRT)=MIN(MAX(RHO*ABS(H(JSTRT)),5.0D-3),5.0D 2) DO 10 J=JSTRT+1,JSTOP H(J)=0.0D 0 10 CONTINUE 20 CONTINUE END IF IDECF=0 IF (KIT.LT.NIT) THEN NRES=NRES+1 KIT=NIT ELSE ITERM=-10 IF (ITERS.LT.0) ITERM=ITERS-5 IF (GMAX.LE.1.0D 2*TOLG) ITERM=-ITERM END IF END IF IF (ITERM.NE.0) GO TO 11080 IF (MED.EQ.1) THEN CALL MXVSET(IH(NF+1)-1,0.0D 0,H) CALL PFSEB4(NA,H,IH,JH,AH,IAG,JAG,IA,AZL,AZU,MET5) ELSE CALL PYTSCH(NF,IX,H,IH,JH,KBF) END IF * * DIRECTION DETERMINATION * GNORM=SQRT(MXVDOT(NF,G,G)) CALL MXVSET(NF+1,0.0D 0,GS) DO 13315 KA=1,NA ALF=0.0D 0 BET=0.0D 0 IF (IEXT.LE.0) THEN RHO=AZU(KA)**2 ALF=ALF+RHO BET=BET+RHO END IF IF (IEXT.GE.0) THEN RHO=AZL(KA)**2 ALF=ALF+RHO BET=BET-RHO END IF ALF=ALF/RPF3 BET=BET/RPF3 CALL PALNG3(AG,IAG,JAG,GO,KA) CALL PASSH3(H,IH,JH,IAG,JAG,GO,KA,ALF) K=IAG(KA) L=IAG(KA+1)-K DO 13310 J=1,L JP=ABS(JAG(K)) GS(JP)=GS(JP)+BET*AG(K) K=K+1 13310 CONTINUE GS(NF+1)=GS(NF+1)+ALF 13315 CONTINUE IF (IDECF.NE.0.AND.IDECF.NE.1) THEN ITERD=-1 GO TO 13290 END IF INITD=MAX(ABS(INITD),1) MM=IH(NF+1)-1 IF (IDECF.NE.1) THEN CALL MXSPCT(NF,MM,MH,MMAX,H,JH,PSL,ITERM) IF (ITERM.NE.0) THEN GO TO 13290 END IF * * GILL-MURRAY DECOMPOSITION * RHO=ETA2 CALL MXSPCF(NF,H(MM+1),PSL,JH(MM+1),WN11,WN12,GO,INF,RHO,ALF) NDEC=NDEC+1 IDECF=1 END IF * * BACK SUBSTITUTIONS * CALL MXVCOP(NF,GS,GP) CALL MXVSFP(NF,PERM,GP,GO) CALL MXSPCB(NF,H(MM+1),PSL,JH(MM+1),GP,0) CALL MXVSBP(NF,PERM,GP,GO) CALL MXVNEG(NF,G,S) CALL MXVSFP(NF,PERM,S,GO) CALL MXSPCB(NF,H(MM+1),PSL,JH(MM+1),S,0) CALL MXVSBP(NF,PERM,S,GO) IF (MXVDOT(NF,GS,GP)-GS(NF+1).EQ.0.0D 0) THEN S(NF+1)=0.0D 0 ELSE S(NF+1)=-(MXVDOT(NF,GS,S)-G(NF+1))/(MXVDOT(NF,GS,GP)-GS(NF+1)) CALL MXVDIR(NF,S(NF+1),GP,S,S) END IF SNORM=SQRT(MXVDOT(NF,S,S)) * * COMPUTATION OF THE DIRECTIONAL DERIVATIVE * P=MXVDOT(NF,S,G) 13290 CONTINUE R=0.0D 0 * * END OF DIRECTION DETERMINATION * IF (KD.GT.0) P=MXVDOT(NF,G,S) * * TEST ON DESCENT DIRECTION AND PREPARATION OF LINE SEARCH * IF (ITERD.LT.0) THEN ITERM=ITERD ELSE * * TEST ON DESCENT DIRECTION * IF (SNORM.LE.0.0D 0) THEN IREST=MAX(IREST,1) ELSE IF (P+EPS0*GNORM*SNORM.LE.0.0D 0) THEN IREST=0 ELSE * * UNIFORM DESCENT CRITERION * IREST=MAX(IREST,1) END IF IF (IREST.EQ.0) THEN * * PREPARATION OF LINE SEARCH * NRED = 0 RMIN=ALF1*GNORM/SNORM RMAX=MIN(ALF2*GNORM/SNORM,XMAX/SNORM) END IF END IF IF (ITERM.NE.0) GO TO 11080 IF (IREST.NE.0) GO TO 11040 IF (NIT.EQ.1) KIT=NIT CALL PYTRCS(NF,X,IX,XO,X,X,G,GO,S,RO,FP,FO,F,PO,P,RMAX,ETA9, & KBF) FFO=FF XO(NF+1)=X(NF+1) GO(NF+1)=G(NF+1) CALL MXVCOP(NA,AF,AFO) IF (MED.EQ.1) CALL MXVCOP(MA,AG,AGO) 11060 CONTINUE CALL PS0L02(R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,EPS1, & KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,MES,ISYS) IF (ISYS.EQ.0) GO TO 11064 CALL MXVDIR(NF,R,S,XO,X) CALL PP1MX3(NF,NA,X,GA,AG,IAG,JAG,G,AZL,AZU,FA,AF,FF,KD,LD,NFV, & NFG,ISNA,IEXT) LD=KD IF (FF+RPF3.GT.(1.0D 0+ETA0)*FF) THEN NNIT=20 11063 CONTINUE CALL PNNEQ1(FF+RPF3,FF+RPF3*DBLE(2*NA),X(NF+1),G(NF+1), & 0.0D 0,MIN(ETA3,TOLG),NNIT,IER,ICON) IF (ICON.GT.0) THEN G(NF+1)=PNFUZ1(X(NF+1),NA,RPF3,AF,AZL,AZU,IEXT) GO TO 11063 END IF ELSE X(NF+1)=(1.0D 0+ETA0)*FF G(NF+1)=PNFUZ1(X(NF+1),NA,RPF3,AF,AZL,AZU,IEXT) END IF CALL PP0BX1(NA,X(NF+1),AF,F,FF,PAR,RPF3,MEP,IEXT) GO TO 11060 11064 CONTINUE KD=MED IF (ITERS.LE.0) THEN R=0.0D 0 F=FO P=PO CALL MXVCOP(NF,XO,X) FF=FFO X(NF+1)=XO(NF+1) G(NF+1)=GO(NF+1) CALL MXVCOP(NA,AFO,AF) IF (MED.EQ.1) CALL MXVCOP(MA,AGO,AG) IF (IDIR.EQ.0) IREST=MAX(IREST,1) LD=KD GO TO 11040 END IF IF (MER.EQ.1) THEN RPF3=MIN(ETA4*RPF3,GNORM**2) ELSE IF (MER.EQ.2) THEN RPF3=MIN(MAX(ETA4*RPF3,RPF3/(1.0D 2*RPF3+1.0D 0)),MAX(GNORM**2, & 1.0D-2**NIT)) ELSE IF (MER.EQ.3) THEN IF (GNORM.GE.ETA6) THEN ELSE RPF3=MIN(MAX(ETA4*RPF3,RPF3/(1.0D 2*RPF3+1.0D 0)),MAX(GNORM**2, & 1.0D-2**NIT)) END IF ELSE IF (MER.EQ.4) THEN IF (GNORM.GE.ETA6) THEN ELSE IF (RPF3.GE.1.0D 1*GNORM**2) RPF3=GNORM**2 END IF END IF RPF3=MAX(RPF3,TOLP) GO TO 11010 11070 CONTINUE CALL PYTRCD(NF,X,IX,XO,G,GO,R,F,FO,P,PO,DMAX,KBF,KD,LD,ITERS) IF (MED.EQ.1) THEN CALL PYTCUB(NA,MA,AG,AGO,IAG,IA,AZL,AZU,ITERS,MET5) IDECF=0 CALL PUBBM2(NA,AH,IAG,JAG,S,XO,AGO,ETA0,ETA9,NIT,KIT,ITERH,MET, & MET1,MET3) IF (ITERH.NE.0) IREST=MAX(IREST,1) END IF GO TO 11030 11080 CONTINUE F=FF IF (IPRNT.GT.1.OR.IPRNT.LT.0) & WRITE(6,'(1X,''EXIT FROM PMAX :'')') IF (IPRNT.NE.0) & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X, & ''F='', G16.9,2X,''G='',E10.3,2X,''ITERM='',I3)') NIT,NFV,NFG, & F,GMAX,ITERM IF (IPRNT.LT.0) & WRITE (6,'(1X,''X='',5(G14.7,1X):/(3X,5(G14.7,1X)))') & (X(I),I=1,NF) END