************************************************************************ * SUBROUTINE PSEDU ALL SYSTEMS 97/01/22 * PURPOSE : * EASY TO USE SUBROUTINE FOR LARGE-SCALE UNCONSTRAINED MINIMIZATION * OF PARTIALLY SEPARABLE FUNCTIONS. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NA NUMBER OF APPROXIMATED FUNCTIONS. * IU MA NUMBER OF NONZERO ELEMENTS IN THE JACOBIAN MATRIX. * RI X(NF) VECTOR OF VARIABLES. * RO AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED * 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) TYPE OF VARIABLE METRIC UPDATE. IPAR(5)=1-THE BFGS * UPDATE. IPAR(5)=2-A COMBINATION OF THE BFGS AND THE RANK-ONE * UPDATES. IPAR(5)=3-THE DISCRETE NEWTON METHOD. * IPAR(6) THIS PARAMETER IS NOT USED IN THE SUBROUTINE PSED. * 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 PSED. * RPAR(8) THIS PARAMETER IS NOT USED IN THE SUBROUTINE PSED. * RPAR(9) THIS PARAMETER IS NOT USED IN THE SUBROUTINE PSED. * 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 PSED VARIABLE METRIC METHOD FOR PARTIABLY SEPERABLE FUNCTIONS. * 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 PSEDU(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 NB,LGA,LG,LH,LHA,LAH,LS,LXO,LGO,LAG,LAGO,LIH,LJH,LPSL, & LPERM,LINVP,LWN11,LWN12,LWN13,LWN14,MH,ML,IFIL,IER INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH 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,MH,ML,IAG) ALLOCATE (IA(8*NF+6+(IFIL+3)*MH)) IF (IPAR(5).LE.2) THEN ALLOCATE (RA(5*NF+2*MA+(IFIL+4)*MH)) ELSE ALLOCATE (RA(5*NF+ML+(IFIL+3)*MH)) END IF NB=0 * * POINTERS FOR AUXILIARY ARRAYS * LGA=1 LG=LGA+NF LS=LG+NF LXO=LS+NF LGO=LXO+NF IF (IPAR(5).LE.2) THEN LHA=LGO LAG=LGO+NF LAGO=LAG+MA LAH=LAGO+MA LH=LAH+MH ELSE LAG=LGO LAGO=LGO LAH=LGO LHA=LGO+NF LH=LHA+ML END IF LIH=1 LPSL=LIH+NF+1 LPERM=LPSL+NF+1 LINVP=LPERM+NF LWN11=LINVP+NF LWN12=LWN11+NF+1 LWN13=LWN12+NF+1 LWN14=LWN13+NF+1 LJH=LWN14+NF+1 CALL PSED(NF,NA,NB,(IFIL+3)*MH,X,IA,RA,RA,AF,RA(LGA),RA(LG), & RA(LHA),RA(LAH),RA(LH),IA(LIH),IA(LJH),RA(LAG),IAG,JAG,RA(LS), & RA(LXO),RA(LGO),RA(LAGO),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),GMAX,F,IPAR(1),IPAR(2),IPAR(3),IPAR(4),IPAR(5), & IPRNT,ITERM) DEALLOCATE (IA,RA) RETURN END ************************************************************************ * SUBROUTINE PSEDS ALL SYSTEMS 97/01/22 * PURPOSE : * EASY TO USE SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION * OF PARTIALLY SEPARABLE FUNCTIONS. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NA NUMBER OF APPROXIMATED FUNCTIONS. * IU MA NUMBER OF NONZERO ELEMENTS IN THE JACOBIAN MATRIX. * RI X(NF) VECTOR OF VARIABLES. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE * X(I) IS UNBOUNDED. IX(I)=1-LOWER BOUND XL(I).LE.X(I). * IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND * XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. * RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. * RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. * RO AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED * 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) TYPE OF VARIABLE METRIC UPDATE. IPAR(5)=1-THE BFGS * UPDATE. IPAR(5)=2-A COMBINATION OF THE BFGS AND THE RANK-ONE * UPDATES. IPAR(5)=3-THE DISCRETE NEWTON METHOD. * IPAR(6) THIS PARAMETER IS NOT USED IN THE SUBROUTINE PSED. * 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 PSED. * RPAR(8) THIS PARAMETER IS NOT USED IN THE SUBROUTINE PSED. * RPAR(9) THIS PARAMETER IS NOT USED IN THE SUBROUTINE PSED. * 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 PSED VARIABLE METRIC METHOD FOR PARTIABLY SEPERABLE FUNCTIONS. * 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 PSEDS(NF,NA,MA,X,IX,XL,XU,AF,IAG,JAG,IPAR,RPAR,F,GMAX, & ISPAS,IPRNT,ITERM) INTEGER NF,NA,MA,IX(*),IAG(*),JAG(*),IPAR(7),ISPAS,IPRNT,ITERM DOUBLE PRECISION X(*),XL(*),XU(*),AF(*),RPAR(9),F,GMAX INTEGER NB,LGA,LG,LH,LHA,LAH,LS,LXO,LGO,LAG,LAGO,LIH,LJH,LPSL, & LPERM,LINVP,LWN11,LWN12,LWN13,LWN14,MH,ML,IFIL,IER INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH 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,MH,ML,IAG) ALLOCATE (IA(8*NF+6+(IFIL+3)*MH)) IF (IPAR(5).LE.2) THEN ALLOCATE (RA(5*NF+2*MA+(IFIL+4)*MH)) ELSE ALLOCATE (RA(5*NF+ML+(IFIL+3)*MH)) END IF NB=1 * * POINTERS FOR AUXILIARY ARRAYS * LGA=1 LG=LGA+NF LS=LG+NF LXO=LS+NF LGO=LXO+NF IF (IPAR(5).LE.2) THEN LHA=LGO LAG=LGO+NF LAGO=LAG+MA LAH=LAGO+MA LH=LAH+MH ELSE LAG=LGO LAGO=LGO LAH=LGO LHA=LGO+NF LH=LHA+ML END IF LIH=1 LPSL=LIH+NF+1 LPERM=LPSL+NF+1 LINVP=LPERM+NF LWN11=LINVP+NF LWN12=LWN11+NF+1 LWN13=LWN12+NF+1 LWN14=LWN13+NF+1 LJH=LWN14+NF+1 CALL PSED(NF,NA,NB,(IFIL+3)*MH,X,IX,XL,XU,AF,RA(LGA),RA(LG), & RA(LHA),RA(LAH),RA(LH),IA(LIH),IA(LJH),RA(LAG),IAG,JAG,RA(LS), & RA(LXO),RA(LGO),RA(LAGO),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),GMAX,F,IPAR(1),IPAR(2),IPAR(3),IPAR(4),IPAR(5), & IPRNT,ITERM) DEALLOCATE (IA,RA) RETURN END ************************************************************************ * SUBROUTINE PSED ALL SYSTEMS 01/09/22 * PURPOSE : * GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION * OF PARTIALLY SEPARABLE FUNCTIONS. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NA NUMBER OF LINEAR APPROXIMATED FUNCTIONS. * II NB CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED. * NB>0-SIMPLE BOUNDS ACCEPTED. * II MMAX MAXIMUM DIMENSION OF THE SPARSE TABLEAU. * RI X(NF) VECTOR OF VARIABLES. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE * X(I) IS UNBOUNDED. IX(I)=1-LOWER BOUND XL(I).LE.X(I). * IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND * XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. * 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. * RA GA(NF) GRADIENT OF THE SELECTED APPROXIMATED FUNCTION. * RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RA HA(ML) HESSIAN MATRIX OF THE SELECTED APPROXIMATING FUNCTION. * RA AH(MH) ELEMENTS OF THE PARTITIONED HESSIAN MATRIX. * RA H(MMAX) NONZERO ELEMENTS OF THE SPARSE HESSIAN MATRIX * TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL * DIFFERENTIATION. * II IH(NF+1) POINTERS OF DIAGONAL ELEMENTS OF THE MATRIX H. * IU JH(MMAX) INDICES OF NONZERO ELEMENTS OF THE MATRIX H * TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL * DIFFERENTIATION. * RA AG(MA) JACOBIAN MATRIX OF THE PARTITIONED FUNCTION. * RI IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. * RI JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. * RO S(NF) DIRECTION VECTOR. * RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. * RI GO(NF) GRADIENTS DIFFERENCE. * RA AGO(MA) OLD JACOBIAN MATRIX OF THE PARTITIONED FUNCTION, * II 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. * 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 MET TYPE OF VARIABLE METRIC UPDATE. MET=1-THE BFGS UPDATE. * MET=2-A COMBINATION OF THE BFGS AND THE RANK-ONE UPDATES. * MET=3-THE 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 PA1SF3 COMPUTATION OF THE VALUE AND THE GRADIENT OF A * PARTIALLY SEPARABLE OBJECTIVE FUNCTION. * S PA2SF4 COMPUTATION OF THE VALUE, GRADIENT AND THE HESSIAN * MATRIX OF A PARTIALLY SEPARABLE OBJECTIVE FUNCTION. * S PCBS04 ELIMINATION OF BOX CONSTRAINT VIOLATIONS. * S PDSLM1 DIRECTION DETERMINATION USING THE GILL-MURRAY DIRECT * DECOMPOSITION METHOD. * S PFSET3 PREPARATION OF THE SPARSE HESSIAN MATRIX * S PFSET4 PREPARATION OF THE PARTITIONED HESSIAN MATRIX * S PS1L01 STEPSIZE SELECTION USING LINE SEARCH. * S PUBBM1 VARIABLE METRIC UPDATE OF THE PARTITIONED HESSIAN MATRIX. * S PYADC0 ADDITION OF A BOX CONSTRAINT. * S PYFUT1 TEST ON TERMINATION. * S PYRMC0 DELETION OF A BOX CONSTRAINT. * 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 HESSIAN MATRIX. * S MXSPCC SPARSE MATRIX REORDERING, SYMBOLIC FACTORIZATION, DATA * STRUCTURES TRANSFORMATION. INITIATION OF THE DIRECT SPARSE * SOLVER. * S MXSSMI SPARSE SYMMETRIC MATRIX IS REPLACED BY THE UNIT MATRIX. * S MXUDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * RF MXUDOT DOT PRODUCT OF TWO VECTORS. * S MXVCOP COPYING OF A VECTOR. * S MXVDIF DIFFERENCE OF TWO VECTORS. * S MXVINB PROJECTION OF A SPARSE SYMMETRIC MATRIX TO SATISFY BOX * CONSTRAINTS. * S MXVINE RESTORATION OF A SPARSE SYMMETRIC MATRIX OBTAINED BY * MXVINB * S MXVINS INITIATION OF THE INTEGER VECTOR. * S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE * SUBSTRACTED ONE. * S MXVSET INITIATINON 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 : * VARIABLE METRIC AND DISCRETE NEWTON METHODS FOR MINIMIZATION OF * LARGE-SCALE PARTIALLY SEPARABLE FUNCTIONS. * SUBROUTINE PSED(NF,NA,NB,MMAX,X,IX,XL,XU,AF,GA,G,HA,AH,H,IH,JH, & AG,IAG,JAG,S,XO,GO,AGO,PSL,PERM,INVP,WN11,WN12,WN13,WN14,XMAX, & TOLX,TOLF,TOLB,TOLG,FMIN,GMAX,F,MIT,MFV,MFG,IEST,MET,IPRNT, & ITERM) INTEGER NF,NA,NB,MMAX,IX(*),IH(*),JH(*),IAG(*),JAG(*),PSL(*), & PERM(*),INVP(*),WN11(*),WN12(*),WN13(*),WN14(*),IPRNT,MIT,MFV, & MFG,IEST,MET,ITERM DOUBLE PRECISION X(*),XL(*),XU(*),AF(*),AG(*),GA(*),G(*),HA(*), & AH(*),H(*),S(*),XO(*),GO(*),AGO(*),TOLX,TOLF,TOLG,TOLB,FMIN, & XMAX,GMAX,UMAX,F INTEGER IDECF,ITERD,ITERS,ITERH,KD,LD,NTESX,NTESF,MTESX,MTESF, & MRED,KIT,IREST,KBF,MET1,MES,MES1,MES2,MES3,MAXST,ISYS,ITES, & INITS,KTERS,IRES1,IRES2,NRED,INEW,IOLD,I,M,MA,MH,N,ICOR,ISNA DOUBLE PRECISION R,RO,RP,FO,FP,P,PO,PP,GNORM,SNORM,RMIN,RMAX, & FMAX,DMAX,ETA0,ETA2,ETA9,EPS8,EPS9,ALF1,ALF2,PAR1,PAR2,TOLD, & TOLS,TOLP DOUBLE PRECISION MXUDOT 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 PSED :'')') * * INITIATION * KBF=0 IF (NB.GT.0) KBF=2 NRES=0 NDEC=0 NIN=0 NIT=0 NFV=0 NFG=0 NFH=0 ISYS=0 ITES=1 NTESX=0 NTESF=0 MTESX=2 MTESF=2 INITS=2 ITERM=0 ITERD=0 ITERS=2 ITERH=0 KTERS=3 IREST=0 IRES1=999 IRES2=0 IDECF=0 MRED=10 MET1=1 MES=4 MES1=2 MES2=2 MES3=2 ETA0=1.0D-15 ETA2=1.0D-18 ETA9=1.0D 120 EPS8=1.00D 0 EPS9=1.00D-8 ALF1=1.0D-10 ALF2=1.0D 10 RMAX=ETA9 DMAX=ETA9 FMAX=1.0D 20 IF (IEST.LE.0) FMIN=-1.0D 60 IF (IEST.GT.0) IEST=1 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 (TOLG.LE.0.0D 0) TOLG=1.0D-6 IF (TOLB.LE.0.0D 0) TOLB=FMIN+1.0D-16 TOLD=1.0D-4 TOLS=1.0D-4 TOLP=0.9D 0 IF (MIT.LE.0) MIT=9000 IF (MFV.LE.0) MFV=9000 IF (MFG.LE.0) MFG=9000 IF (MET.LE.0) MET=2 IF (MET.LE.2) MFG=MFV KD=MAX(1,MET-1) LD=-1 ISNA=3-KD KIT=-(IRES1*NF+IRES2) FO=FMIN * * INITIAL OPERATIONS WITH SIMPLE BOUNDS * IF (KBF.GT.0) THEN DO 2 I = 1,NF IF ((IX(I).EQ.3.OR.IX(I).EQ.4) .AND. XU(I).LE.XL(I)) THEN XU(I) = XL(I) IX(I) = 5 ELSE IF (IX(I).EQ.5 .OR. IX(I).EQ.6) THEN XL(I) = X(I) XU(I) = X(I) IX(I) = 5 END IF 2 CONTINUE END IF MA=IAG(NA+1)-1 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 11190 ICOR=0 IF (KBF.GT.0) THEN CALL PCBS04(NF,X,IX,XL,XU,EPS9,KBF) CALL PYADC0(NF,N,X,IX,XL,XU,INEW) CALL MXVINE(IH(NF+1)-1,JH) END IF MH=0 CALL MXSPCC(NF,M,MH,MMAX,H,IH,JH,PSL,PERM,INVP,WN11,WN12,WN13, & WN14,ITERM) IF (ITERM.NE.0) GO TO 11190 IF (MET.LE.2) THEN CALL PA1SF3(NF,NA,X,GA,G,AG,IAG,JAG,F,AF,KD,LD,ISNA,NFV,NFG) ELSE CALL PA2SF4(NF,NA,X,IX,GA,G,GO,HA,H,IH,JH,IAG,JAG,AF,F,ETA0,KBF, & KD,LD,NFV,NFG,IDECF) END IF * * START OF THE ITERATION WITH TESTS FOR TERMINATION. * 11120 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 PYFUT1(NF,F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLB,TOLG,KD, & NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,ITES, & IRES1,IRES2,IREST,ITERS,ITERM) IF (ITERM.NE.0) GO TO 11190 IF (KBF.GT.0) CALL PYRMC0(NF,N,IX,G,EPS8,UMAX,GMAX,RMAX,IOLD, & IREST) 11130 CONTINUE IF (IREST.GT.0) THEN IF (MET.LE.2) THEN CALL MXBSMI(NA,AH,IAG) ELSE CALL MXSSMI(NF,H,IH) END IF IDECF=0 LD=MIN(LD,1) 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 11190 IF (KBF.GT.0) THEN CALL MXVINB(M,IX,JH) CALL MXVINB(MA,IX,JAG) END IF IF (ITERS.NE.0) THEN IF (MET.LE.2) THEN CALL MXVSET(IH(NF+1)-1,0.0D 0,H) CALL PFSET4(NA,H,IH,JH,AH,IAG,JAG) END IF END IF IF (KBF.GT.0) CALL PYTSCH(NF,IX,H,IH,JH,KBF) * * DIRECTION DETERMINATION * CALL PDSLM1(NF,MMAX,MH,IX,G,H,IH,JH,S,XO,PSL,PERM,WN11,WN12, & GNORM,SNORM,ETA2,KBF,IDECF,NDEC,ITERD,ITERM) * * TEST ON DESCENT DIRECTION AND PREPARATION OF LINE SEARCH * IF (KD.GT.0) P=MXUDOT(NF,G,S,IX,KBF) 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+TOLD*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 11190 IF (IREST.NE.0) GO TO 11130 IF (NIT.EQ.1) KIT=NIT CALL PYTRCS(NF,X,IX,XO,XL,XU,G,GO,S,RO,FP,FO,F,PO,P,RMAX,ETA9, & KBF) IF (RMAX.EQ.0.0D 0) GO TO 11175 IF (MET.LE.2) CALL MXVCOP(MA,AG,AGO) 11170 CONTINUE CALL PS1L01(R,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX, & TOLS,TOLP,PAR1,PAR2,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST, & INITS,ITERS,KTERS,MES,ISYS) IF (ISYS.EQ.0) GO TO 11174 CALL MXUDIR(NF,R,S,XO,X,IX,KBF) CALL PCBS04(NF,X,IX,XL,XU,EPS9,KBF) IF (MET.LE.2) THEN CALL PA1SF3(NF,NA,X,GA,G,AG,IAG,JAG,F,AF,KD,LD,ISNA,NFV,NFG) ELSE CALL PA2SF4(NF,NA,X,IX,GA,G,GO,HA,H,IH,JH,IAG,JAG,AF,F,ETA0,KBF, & KD,LD,NFV,NFG,IDECF) END IF P=MXUDOT(NF,G,S,IX,KBF) GO TO 11170 11174 CONTINUE IF (ITERS.LE.0) THEN R=0.0D 0 F=FO P=PO CALL MXVCOP(NF,XO,X) CALL MXVCOP(NF,GO,G) IF (MET.LE.2) CALL MXVCOP(MA,AGO,AG) IREST=MAX(IREST,1) LD=KD GO TO 11130 END IF KD=MAX(1,MET-1) IF (KD.GT.LD) THEN IF (MET.LE.2) THEN CALL PA1SF3(NF,NA,X,GA,G,AG,IAG,JAG,F,AF,KD,LD,ISNA,NFV,NFG) ELSE CALL PA2SF4(NF,NA,X,IX,GA,G,GO,HA,H,IH,JH,IAG,JAG,AF,F,ETA0,KBF, & KD,LD,NFV,NFG,IDECF) END IF END IF CALL PYTRCD(NF,X,IX,XO,G,GO,R,F,FO,P,PO,DMAX,KBF,KD,LD,ITERS) IF (MET.LE.2) THEN IF (ITERS.GT.0) THEN CALL MXVDIF(MA,AG,AGO,AGO) ELSE CALL MXVSAV(MA,AG,AGO) END IF IDECF=0 CALL PUBBM1(NA,AH,IAG,JAG,S,XO,AGO,ETA0,ETA9,ICOR,NIT,KIT,ITERH, & MET,MET1) END IF 11175 CONTINUE IF (ITERH.NE.0) IREST=MAX(IREST,1) IF (KBF.GT.0) CALL PYADC0(NF,N,X,IX,XL,XU,INEW) GO TO 11120 11190 CONTINUE IF (IPRNT.GT.1.OR.IPRNT.LT.0) & WRITE(6,'(1X,''EXIT FROM PSED :'')') 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) RETURN END