************************************************************************ * SUBROUTINE PSUMU ALL SYSTEMS 97/01/22 * PURPOSE : * EASY TO USE SUBROUTINE FOR LARGE-SCALE UNCONSTRAINED MINIMIZATION OF * SUMS OF ABSOLUTE VALUES 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 PSUM. * 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) INITIAL TRUST-REGION RADIUS. * RPAR(8) THIS PARAMETER IS NOT USED IN THE SUBROUTINE PSUM. * RPAR(9) MINIMUM PERMITTED VALUE OF THE BARRIER PARAMETER. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RO GMAX MAXIMUM PARTIAL DERIVATIVE. * 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 PSUM PRIMAL TRUST-REGION INTERIOR-POINT METHOD FOR LARGE-SCALE * PARTIALLY SEPARABLE SUMS OF ABSOLUTE VALUES. * 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 PSUMU(NF,NA,MA,X,AF,IAG,JAG,IPAR,RPAR,F,GMAX,ISPAS, & IPRNT,ITERM) INTEGER NF,NA,MA,IAG(*),JAG(*),IPAR(7),ISPAS,IPRNT,ITERM DOUBLE PRECISION X(*),AF(*),RPAR(9),F,GMAX INTEGER LAFO,LAG,LAGO,LGA,LAH,LAS,LAZ,LG,LH,LS,LXO,LGO,LGS,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+6*NF+1+(IFIL+4)*MB)) ELSE ALLOCATE (RA(MA+3*NA+6*NF+1+(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 LAS=LGA+NF LAZ=LAS+NA LG=LAZ+NA LS=LG+NF LXO=LS+NF LGO=LXO+NF LGS=LGO+NF+1 LH=LGS+NF 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+1 CALL PSUM(NF,NA,(IFIL+3)*MB,X,IA,AF,RA(LAFO),RA(LAG),RA(LAGO), & RA(LGA),RA(LAH),RA(LAS),RA(LAZ),RA(LG),RA(LH),IA(LIH),IA(LJH), & IA,IAG,JAG,RA(LS),RA(LXO),RA(LGO),RA(LGS),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(7),RPAR(9), & GMAX,F,IPAR(1),IPAR(2),IPAR(3),IPAR(4),IPAR(5),IPRNT,ITERM) DEALLOCATE (IA,RA) RETURN END ************************************************************************ * SUBROUTINE PSUM ALL SYSTEMS 01/09/22 * PURPOSE : * GENERAL SUBROUTINE FOR LARGE-SCALE UNCONSTRAINED MINIMIZATION OF * SUMS OF ABSOLUTE VALUES 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 AS(NA) AUXILIARY VECTOR. * RA AZ(NA) VECTOR OF 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. * 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 XDEL TRUST REGION STEPSIZE. * 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 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 PNSTEP REALIZATION OF THE BOUNDARY STEP. * S PF1HS2 NUMERICAL COMPUTATION OF THE HESSIAN MATRIX USING * DIFFERENCES OF GRADIENTS. * S PFSEB5 COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE * PARTITIONED HESSIAN MATRIX IN THE L-ONE CASE. * S PFSET2 NUMBER OF NONZERO ELEMENTS IN THE PARTITIONED HESSIAN * MATRIX. * S PFSET3 PREPARATION OF THE SPARSE NORMAL EQUATION MATRIX * STRUCTURE. * S PP0BA1 COMPUTATION OF THE VALUE OF THE BARRIER FUNCTION. * S PP1SA3 COMPUTATION OF THE VALUE AND THE GRADIENT OF THE * LAGRANGIAN FUNCTION. * S PPLAG1 DETERMINATION OF THE LAGRANGE MULTIPLIERS. * S PS0G01 STEPSIZE SELECTION USING TRUST REGION. * S PUBBM2 VARIABLE METRIC UPDATES OF THE PARTITIONED MATRIX. * S PYFUT8 TEST ON TERMINATION. * S PYTCAB SCALED DIFFERENCE OF THE JACOBIAN MATRICES IN THE L-ONE * CASE. * S PYPTSH DETERMINATION OF GROUPS FOR NUMERICAL DIFFERENTIATION. * 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 MXSPCM MATRIX-VECTOR PRODUCT USING THE SPARSE DECOMPOSITION * OBTAINED BY MXSPCF. * RF MXSPCQ GENERALIZED DOT PRODUCT USING THE SPARSE DECOMPOSITION * OBTAINED BY MXSPCF. * S MXSPCT COPYING A SPARSE SYMMETRIC MATRIX INTO THE PERMUTED * FACTORIZED COMPACT SCHEME. * RF MXSSMQ COMPUTATION OF THE SPARSE QUADRATIC TERM. * S MXVCOP COPYING OF A VECTOR. * S MXVDIF DIFFERENCE OF TWO VECTORS. * 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 MXVSCL SCALING 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 : * TRUST-REGION INTERIOR-POINT METHOD FOR LARGE SPARSE SUMS OF ABSOLUTE * VALUES. * SUBROUTINE PSUM(NF,NA,MMAX,X,IX,AF,AFO,AG,AGO,GA,AH,AS,AZ,G,H,IH, & JH,IA,IAG,JAG,S,XO,GO,GS,COL,PSL,PERM,INVP,WN11,WN12,WN13,WN14, & XMAX,TOLX,TOLF,TOLB,TOLG,FMIN,XDEL,ETA5,GMAX,F,MIT,MFV,MFG,IEST, & MED,IPRNT,ITERM) INTEGER NF,NA,MMAX,IX(*),IH(*),JH(*),IA(*),IAG(*),JAG(*),COL(*), & PSL(*),PERM(*),INVP(*),WN11(*),WN12(*),WN13(*),WN14(*),MIT,MFV, & MFG,IEST,MED,IPRNT,ITERM DOUBLE PRECISION X(*),AF(*),AFO(*),AG(*),AGO(*),GA(*),AH(*), & AS(*),AZ(*),G(*),H(*),S(*),XO(*),GO(*),GS(*),XMAX,TOLX,TOLF, & TOLG,TOLB,XDEL,ETA5,FMIN,GMAX,F INTEGER IDECF,ITERD,ITERS,KD,LD,NTESX,NTESF,MTESX,MTESF,MRED,KIT, & IREST,KBF,MAXST,IDIR,IOLD,INF,INITD,ITERH,IER,ISYS,KTERS,IRES1, & IRES2,NRED,MET,MET1,MET3,MET5,I,J,MES1,MES2,MES3,M,MA,MB,MM,MH, & JSTRT,JSTOP,KA,ISNA DOUBLE PRECISION R,RO,FF,FO,FP,FA,P,PO,PP,GNORM,SNORM,RMAX,DMAX, & UMAX,XDELO,ETA0,ETA2,ETA6,EPS4,EPS5,ETA9,ALF,BET1,BET2,GAM1, & GAM2,TOLP,RHO,RPF3,FFO,B1,B2,B3 DOUBLE PRECISION MXVDOT,MXSSMQ,MXSPCQ 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 PSUM :'')') * * INITIATION OF PROBLEM * KBF=0 NRES=0 NDEC=0 NIN=0 NIT=0 NFV=0 NFG=0 NFH=0 IDIR=0 ISYS=0 NTESX=0 NTESF=0 MTESX=2 MTESF=2 INITD=1 ITERM=0 ITERD=0 ITERS=2 ITERH=0 KTERS=0 IREST=0 IRES1=999 IRES2=0 MRED=10 MES1=3 MES2=2 MES3=1 ETA0=1.0D-15 ETA2=1.0D-18 ETA6=1.0D 0 ETA9=1.0D 120 EPS4=0.10D 0 EPS5=0.90D 0 BET1=0.05D 0 BET2=0.75D 0 GAM1=2.0D 0 GAM2=1.0D 6 RPF3=1.0D 0 RMAX=ETA9 DMAX=ETA9 IF (IEST.LE.0) FMIN=0.0D 0 FMIN=MAX(FMIN,0.0D 0) IEST=1 IF (XMAX.LE.0.0D 0) XMAX=1.0D 16 XDEL=MIN(XDEL,XMAX) IF (TOLX.LE.0.0D 0) TOLX=1.0D-16 IF (TOLF.LE.0.0D 0) TOLF=1.0D-12 IF (TOLB.LE.0.0D 0) TOLB=FMIN+1.0D-12 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 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 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 (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 IF (ETA5.LE.0.0D 0) THEN C TOLP=SQRT(ETA0) TOLP=1.0D-8 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 PP1SA3(NF,NA,X,GA,AG,IAG,JAG,G,AZ,FA,AF,FF,KD,LD,NFV,NFG, & ISNA) LD=0 CALL PPLAG1(NA,AF,AS,AZ,RPF3) CALL PP0BA1(NA,AS,RPF3,F) * * COMPUTATION OF THE GRADIENT AND THE HESSIAN MATRIX OF THE * LAGRANGIAN FUNCTION * KD=1 CALL PP1SA3(NF,NA,X,GA,AG,IAG,JAG,G,AZ,FA,AF,FF,KD,LD,NFV,NFG, & ISNA) 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 PP1SA3(NF,NA,X,GA,AG,IAG,JAG,GO,AZ,FA,AF,FF,KD,LD,NFV,NFG, & ISNA) 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 3*TOLG.OR.FF.LE.1.0D-8) 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 PFSEB5(NA,H,IH,JH,AH,IAG,JAG,AZ,MET5) ELSE CALL PYTSCH(NF,IX,H,IH,JH,KBF) END IF * * DIRECTION DETERMINATION * B2=MXVDOT(NF,G,G) GNORM=SQRT(B2) DO 13415 KA=1,NA ALF=2.0D 0*RPF3/(AS(KA)**2+AF(KA)**2) CALL PALNG3(AG,IAG,JAG,GO,KA) CALL PASSH3(H,IH,JH,IAG,JAG,GO,KA,ALF) 13415 CONTINUE IF (IDECF.NE.0.AND.IDECF.NE.1) THEN ITERD=-1 GO TO 13390 END IF INITD=MAX(ABS(INITD),1) MM=IH(NF+1)-1 IF (IDECF.EQ.0) THEN B1=MXSSMQ(NF,H,IH,JH,G,G) ELSE CALL MXVCOP(NF,G,GO) CALL MXVSFP(NF,PERM,GO,XO) CALL MXSPCM(NF,H(MM+1),PSL,JH(MM+1),GO,XO,1) B1=MXSPCQ(NF,H(MM+1),PSL,GO) END IF IF (XDEL.LE.0.0D 0) THEN * * INITIAL TRUST REGION BOUND * IF (B1.LE.0.0D 0) THEN XDEL=GNORM ELSE XDEL=(B2/B1)*GNORM END IF XDEL=MIN(XDEL,XMAX) END IF IF (B1.LE.0.0D 0.OR.B2*GNORM.GE.B1*XDEL) THEN * * SCALED STEEPEST DESCENT DIRECTION IS ACCEPTED * CALL MXVSCL(NF,-XDEL/GNORM,G,S) SNORM=XDEL ITERD=3 GO TO 13390 END IF IF (IDECF.EQ.0) THEN CALL MXSPCT(NF,MM,MH,MMAX,H,JH,PSL,ITERM) IF (ITERM.NE.0) THEN GO TO 13390 END IF * * SPARSE GILL-MURRAY DECOMPOSITION * RHO=ETA2 CALL MXSPCF(NF,H(MM+1),PSL,JH(MM+1),WN11,WN12,XO,INF,RHO,ALF) NDEC=NDEC+1 IDECF=1 END IF * * COMPUTATION OF THE NEWTON DIRECTION * CALL MXVNEG(NF,G,GO) CALL MXVSFP(NF,PERM,GO,XO) CALL MXSPCB(NF,H(MM+1),PSL,JH(MM+1),GO,0) CALL MXVSBP(NF,PERM,GO,XO) RHO=SQRT(MXVDOT(NF,GO,GO)) * * COMPUTATION OF THE STEEPEST DESCENT DIRECTION * B2=B2/B1 SNORM=B2*GNORM CALL MXVSCL(NF,-B2,G,S) CALL MXVDIF(NF,GO,S,XO) B1=MXVDOT(NF,S,XO) B2=MXVDOT(NF,XO,XO) IF (B2.LE.1.0D-8*XDEL*XDEL) THEN * * NEWTON AND THE STEEPEST DESCENT DIRECTION ARE * APPROXIMATELY EQUAL * CALL MXVCOP(NF,GO,S) SNORM=RHO ITERD=1 ELSE IF (B1.LE.0.0D 0) THEN * * BOUNDARY STEP WITH NEGATIVE INCREMENT * CALL PNSTEP(XDEL,SNORM,-B1,B2,B3) CALL MXVDIR(NF,-B3,XO,S,S) SNORM=XDEL ITERD=3 ELSE IF (RHO.LE.XDEL) THEN * * NEWTON DIRECTION IS ACCEPTED * CALL MXVCOP(NF,GO,S) SNORM=RHO ITERD=1 ELSE * * DOUBLE DOGLEG STRATEGY * RHO=XDEL/RHO B3=MXVDOT(NF,S,GO) RHO=MAX(RHO,SNORM*SNORM/B3) CALL MXVDIR(NF,-RHO,GO,S,XO) B1=SNORM*SNORM-RHO*B3 B2=MXVDOT(NF,XO,XO) CALL PNSTEP(XDEL,SNORM,-B1,B2,B3) CALL MXVDIR(NF,-B3,XO,S,S) SNORM=XDEL ITERD=3 END IF 13390 CONTINUE IF (IDECF.EQ.0) THEN PP=MXSSMQ(NF,H,IH,JH,S,S)*0.5D 0 ELSE CALL MXVCOP(NF,S,GO) CALL MXVSFP(NF,PERM,GO,XO) CALL MXSPCM(NF,H(MM+1),PSL,JH(MM+1),GO,XO,1) PP=MXSPCQ(NF,H(MM+1),PSL,GO)*0.5D 0 IF (ITERD.EQ.1.AND.INF.NE.0) ITERD=2 END IF * * END OF DIRECTION DETERMINATION * IF (KD.GT.0) P=MXVDOT(NF,G,S) * * TEST ON LOCALLY CONSTRAINED STEP AND PREPARATION OF STEPSIZE * SELECTION * IF (ITERD.LT.0) THEN ITERM=ITERD ELSE IF (SNORM.LE.0.0D 0) THEN IREST=MAX(IREST,1) ELSE IREST=0 END IF IF (IREST.EQ.0) THEN RMAX=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 CALL MXVCOP(NA,AF,AFO) IF (MED.EQ.1) CALL MXVCOP(MA,AG,AGO) 11060 CONTINUE CALL PS0G01(R,F,FO,PO,PP,XDEL,XDELO,XMAX,RMAX,SNORM, & BET1,BET2,GAM1,GAM2,EPS4,EPS5,KD,LD,IDIR,ITERS,ITERD, & MAXST,NRED,MRED,KTERS,MES1,MES2,MES3,ISYS) IF (ISYS.EQ.0) GO TO 11064 CALL MXVDIR(NF,R,S,XO,X) CALL PP1SA3(NF,NA,X,GA,AG,IAG,JAG,G,AZ,FA,AF,FF,KD,LD,NFV,NFG, & ISNA) LD=KD CALL PPLAG1(NA,AF,AS,AZ,RPF3) CALL PP0BA1(NA,AS,RPF3,F) 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) IF (ITERS.LT.0) THEN RPF3=TOLP IREST=MAX(IREST,1) LD=KD GO TO 11040 END IF FF=FFO 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 (ITERS.GE.2) THEN IF (GNORM.GE.ETA6) THEN ELSE IF (RPF3.GE.1.0D 2*GNORM**2) RPF3=GNORM**2 END IF RPF3=MAX(RPF3,TOLP) END IF 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 PYTCAB(NA,MA,AG,AGO,IAG,AZ,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 PSUM :'')') 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