************************************************************************ * SUBROUTINE PINDU ALL SYSTEMS 97/01/22 * PURPOSE : * EASY TO USE SUBROUTINE FOR EQUALITY CONSTRAINED OPTIMIZATION. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * RI X(NF) VECTOR OF VARIABLES. * II IH(NF+1) POINTERS OF DIAGONAL ELEMENTS OF THE SPARSE HESSIAN * MATRIX. * II JH(NH) INDICES OF NONZERO ELEMENTS OF THE SPARSE HESSIAN MATRIX * TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL * DIFFERENTIATION. * RO CF(NC+1) VECTOR OF THE CONSTRAINT FUNCTION VALUES. * II ICG(NC+1) POSITION OF THE FIRST RORWS ELEMENTS IN THE * CONSTRAINT JACOBIAN MATRIX. * II JCG(ICG(NC+1)-1) COLUMN INDICES OF ELEMENTS IN THE CONSTRAINT * 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) CHOICE OF THE PRECONDITIONER FOR THE KKT SYSTEM. * IPAR(6)=1-INDEFINITE PRECONDITIONER WITH THE COMPLETE * CHOLESKI DECOMPOSITION. IPAR(6)=-1-INDEFINITE PRECONDITIONER * WITH THE INCOMPLETE CHOLESKI DECOMPOSITION. * IPAR(5) MAXIMUM NUMBER OF DENSE ROWS. * IPAR(6) MAXIMUM NUMBER OF NONZERO ELEMENTS IN A SPARSE ROW. * 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(5) REAL PARAMETERS: * RPAR(1) MAXIMUM STEPSIZE. * RPAR(2) TOLERANCE FOR CHANGE OF VARIABLES. * RPAR(3) TOLERANCE FOR CONSTRAINTS. * RPAR(4) TOLERANCE FOR THE GRADIENT OF THE LAGRANGIAN FUNCTION. * RPAR(5) PENALTY PARAMETER. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RO GMAX MAXIMUM PARTIAL DERIVATIVE. * RO CMAX MAXIMUM CONSTRAINT VIOLATION. * 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) SUBSEQUEBT ITERATIONS. * ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN * MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS. * ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB. * ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG. * ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. * ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. * IF ITERM=-6, THEN THE TERMINATION CRITERION HAS NOT BEEN * SATISFIED, BUT THE POINT OBTAINED IF USUALLY ACCEPTABLE. * 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 DECOMPOSITION. * IO NIIT 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 PIND FULLY ITERATIVE ALGORITHM FOR SOLVING SPARSE EQUALITY * CONSTRAINED NONLINEAR PROGRAMMING PROBLEMS. * * EXTERNAL SUBROUTINES : * SE OBJ COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION. * CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS A NUMBER * OF VARIABLES, X(NF) IS A VECTOR OF VARIABLES AND FF IS THE * VALUE OF THE OBJECTIVE FUNCTION. * SE DOBJ COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION. * CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS A NUMBER * OF VARIABLES, X(NF) IS A VECTOR OF VARIABLES AND GF(NF) IS * THE GRADIENT OF THE OBJECTIVE FUNCTION. * SE CON COMPUTATION OF THE VALUE OF THE CONSTRAINT FUNCTION. * CALLING SEQUENCE: CALL CON(NF,KC,X,FC) WHERE NF IS A NUMBER * OF VARIABLES, KC IS THE INDEX OF THE CONSTRAINT FUNCTION, * X(NF) IS A VECTOR OF VARIABLES AND FC IS THE VALUE OF THE * CONSTRAINT FUNCTION. * SE DCON COMPUTATION OF THE GRADIENT OF THE CONSTRAINT FUNCTION. * CALLING SEQUENCE: CALL DCON(NF,KC,X,GC) WHERE NF IS A NUMBER * OF VARIABLES, KC IS THE INDEX OF THE CONSTRAINT FUNCTION, * X(NF) IS A VECTOR OF VARIABLES AND GC(NF) IS THE GRADIENT OF * THE CONSTRAINT FUNCTION. * * METHOD : * PRECONDITIONED SMOOTHED CONJUGATE GRADIENT METHOD WITH A SPECIAL * TERMINATION CRITERION GUARANTEEING SUFFICIENT DESCENT OF THE * DIRECTION VECTOR. * SUBROUTINE PINDU(NF,NC,X,IH,JH,CF,ICG,JCG,IPAR,RPAR,F,GMAX,CMAX, & IPRNT,ITERM) INTEGER NF,NC,IH(*),JH(*),ICG(*),JCG(*),IPAR(7),IPRNT,ITERM DOUBLE PRECISION X(*),CF(*),RPAR(5),F,GMAX,CMAX INTEGER LGF,LH,LCG,LCFO,LGC,LG,LS,LXO,LGO,LXS,LCZ,LCZO,LCP,LB, & LCGD,LDD,LBD,LJH,LKCG,LCOL,LPSR,LPERC,LINVP,LWN11,LWN12,LWN13, & LWN14,LIB,LJB,MB,MC,MD,MH,IFIL INTEGER NRES,NDEC,NIIT,NIT,NFV,NFG,NFH COMMON /STAT/ NRES,NDEC,NIIT,NIT,NFV,NFG,NFH INTEGER IA(:) DOUBLE PRECISION RA(:) ALLOCATABLE IA,RA IFIL=IPAR(7) IF (IFIL.LE.0) IFIL=1 IF (IPAR(5).LE.0) IPAR(5)=10 CALL PFSET3(NF,NC,IH,JH,ICG,JCG) CALL PFSET2(NC,MB,MC,ICG) MC=ICG(NC+1)-1 MD=IPAR(7) MH=IH(NF+1)-1 C MB=MAX(MB,MH) ALLOCATE (IA(4*NF+6*NC+6+(2*IFIL+3)*MB+3*MH), & RA(7*NF+8*NC+MC+NC*MD+MD*(MD+3)/2+1+(2*IFIL+3)*MB+3*MH)) * * POINTERS FOR AUXILIARY ARRAYS * LGF=1 LH=LGF+NF LCG=LH+3*MH LCFO=LCG+MC LGC=LCFO+NC+1 LG=LGC+NF LS=LG+NF LXO=LS+NF+NC LGO=LXO+NF+NC LXS=LGO+NF+NC LCZ=LXS+NF+NC LCZO=LCZ+NC LCP=LCZO+NC LCGD=LCP+NC LDD=LCGD+NC*MD LBD=LDD+MD LB=LBD+MD*(MD+1)/2 LJH=1 LKCG=3*MH LCOL=LKCG+NF LPSR=LCOL+NF LPERC=LPSR+NC+1 LINVP=LPERC+NC LWN11=LINVP+NC LWN12=LWN11+NF+1 LWN13=LWN12+NF+1 LWN14=LWN13+NC+1 LIB=LWN14+NC+1 LJB=LIB+NC+1 CALL MXVICP(IH(NF+1)-1,JH,IA(LJH)) CALL PIND(NF,NC,3*MH,(2*IFIL+3)*MB,X,RA(LGF),RA(LH),IH,IA(LJH), & RA(LCG),ICG,JCG,IA(LKCG),CF,RA(LCFO),RA(LGC),RA(LG),RA(LS), & RA(LXO),RA(LGO),RA(LXS),RA(LCZ),RA(LCZO),RA(LCP),RA(LB),IA(LIB), & IA(LJB),RA(LCGD),RA(LDD),RA(LBD),IA(LCOL),IA(LPSR),IA(LPERC), & IA(LINVP),IA(LWN11),IA(LWN12),IA(LWN13),IA(LWN14),RPAR(1), & RPAR(2),RPAR(3),RPAR(4),RPAR(5),F,GMAX,CMAX,IPAR(1),IPAR(2), & IPAR(3),IPAR(4),IPAR(5),IPAR(6),IPRNT,ITERM) DEALLOCATE (IA,RA) RETURN END ************************************************************************ * SUBROUTINE PIND ALL SYSTEMS 96/12/01 C PORTABILITY : ALL SYSTEMS C 96/12/01 LU : ORIGINAL VERSION * * PURPOSE : * FULLY ITERATIVE ALGORITHM FOR SOLVING SPARSE EQUALITY CONSTRAINED * NONLINEAR PROGRAMMING PROBLEMS. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * II MAXH DIMENSION OF THE ARRAYS H AND JH. * II MAXB DIMENSION OF THE ARRAYS B AND JB. * RU X(NF) VECTOR OF VARIABLES. * RA GF(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RA H(MAXH) NONZERO ELEMENTS OF THE SPARSE SYMMETRIC APPROXIMATION * OF THE LAGRANGIAN FUNCTION 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(MAXH) INDICES OF NONZERO ELEMENTS OF THE MATRIX H TOGETHER * WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL * DIFFERENTIATION. * RA CG(ICG(NC+1)-1) NONZERO ELEMENTS OF THE SPARSE CONSTRAINT * JACOBIAN MATRIX. * II ICG(NC+1) POSITION OF THE FIRST RORWS ELEMENTS IN THE FIELD CG. * II JCG(ICG(NC+1)-1) COLUMN INDICES OF ELEMENTS IN THE FIELD CG. * IA KCG(ICG(NC+1)-1) AUXILIARY ARRAY. * RA CF(NC+1) VECTOR OF THE CONSTRAINT FUNCTION VALUES. * RA CFO(NC+1) SAVED VECTOR OF THE CONSTRAINT FUNCTION VALUES. * RA GC(NF) GRADIENT OF THE CONSTRAINT FUNCTION. * RA G(NF) GRADIENT OF THE LAGRANGIAN FUNCTION. * RA S(NF+NC) DIRECTION VECTOR. * RA XO(NF+NC) AUXILIARY VECTOR. * RA GO(NF+NC) AUXILIARY VECTOR. * RA XS(NF+NC) AUXILIARY VECTOR. * RA CZ(NC) VECTOR OF LAGRANGE MULTIPLIERS. * RA CZO(NC) SAVED VECTOR OF LAGRANGE MULTIPLIERS. * RA CP(NF+NC) AUXILIARY VECTOR. * RA B(MAXB) NONZERO ELEMENTS OF THE SPARSE SYMMETRIC MATRIX * B=TRANS(CG)*D*CG AND ITS INCOMPLETE CHOLESKI DECOMPOSITION. * IA IB(NC+1) POINTERS OF DIAGONAL ELEMENTS OF THE MATRIX B. * IA JB(MAXB) INDICES OF NONZERO ELEMENTS OF THE MATRIX B. * RI CGD(NF*ND) DENSE ROWS. * RA DD(ND) AUXILIARY ARRAY. * RO BD(ND*(ND+1)/2)) DECOMPOSITION OF THE UPDATED MATRIX. * IA COL(NF) VECTOR DISCERNING GROUPS OF THE HESSIAN COLUMNS OF THE * SAME COLOUR. * IA PSR(NC+1) POINTER VECTOR OF THE COMPACT FORM OF THE TRIANGULAR * FACTOR OF THE PRECONDITIONER. * IA PERC(NC) PERMUTATION VECTOR. * IA INVP(NC) INVERSE PERMUTATION VECTOR. * IA WN11(NC+1) AUXILIARY VECTOR. * IA WN12(NC+1) AUXILIARY VECTOR. * IA WN13(NC+1) AUXILIARY VECTOR. * IA WN14(NC+1) AUXILIARY VECTOR. * RI XMAX MAXIMUM STEPSIZE. * RI TOLX TOLERANCE FOR CHANGE OF VARIABLES. * RI TOLC TOLERANCE FOR THE FUNCTION FALUE. * RI TOLG TOLERANCE FOR THE GRADIENT OF THE LAGRANGIAN FUNCTION. * RI RPF1 PENALTY PARAMETER. * RO F VALUE OF THE OBJECTIVE FUNCTION * RO GMAX MAXIMUM PARTIAL DERIVATIVE. * RO CMAX MAXIMUM CONSTRAINT VIOLATION. * II MIT MAXIMUN NUMBER OF ITERATIONS. * II MFV MAXIMUN NUMBER OF FUNCTION EVALUATIONS. * II MFG MAXIMUN NUMBER OF GRADIENT EVALUATIONS. * II MOS SELECTION OF THE PRECONDITIONER FOR THE KKT SYSTEM. * MOS=1-INDEFINITE PRECONDITIONER WITH THE COMPLETE CHOLESKI * DECOMPOSITION. MOS=-1-INDEFINITE PRECONDITIONER WITH THE * INCOMPLETE CHOLESKI DECOMPOSITION. * II MD MAXIMUM NUMBER OF DENSE ROWS. * II MDE MAXIMUM NUMBER OF NONZERO ELEMENTS IN A SPARSE ROW. * 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) SUBSEQUEBT ITERATIONS. * ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN * MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS. * ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB. * ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG. * ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. * ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. * IF ITERM=-6, THEN THE TERMINATION CRITERION HAS NOT BEEN * SATISFIED, BUT THE POINT OBTAINED IF USUALLY ACCEPTABLE. * VALUES ITERM<=-40 DETECT A LACK OF SPACE. * * VARIABLES IN COMMON /STAT/ (STATISTICS) : * IO NRES NUMBER OF RESTARTS. * IO NDEC NUMBER OF MATRIX DECOMPOSITION. * IO NIIT 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 PFSET3 MODIFICATION OF THE SPARSE HESSIAN MATRIX STRUCTURE * USING THE SPARSE CONSTRAINT JACOBIAN MATRIX STRUCTURE. * S PFSED7 DETERMINATION OF THE SPARSE NORMAL EQUATION MATRIX * STRUCTURE USING THE SPARSE CONSTRAINT JACOBIAN MATRIX * STRUCTURE CONSIDERING DENSE ROWS. * S PFSED9 DETERMINATION OF THE SPARSE NORMAL EQUATION MATRIX * USING THE SPARSE CONSTRAINT JACOBIAN MATRIX CONSIDERING * DENSE ROWS. * S PF1HS2 NUMERICAL COMPUTATION OF THE HESSIAN MATRIX USING * DIFFERENCES OF GRADIENTS. * S PP0AF1 COMPUTATION OF THE VALUE OF THE AUGMENTED LAGRANGIAN * FUNCTION * S PP1LF3 DETERMINATION OF THE VALUE AND THE GRADIENT OF THE * LAGRANGIAN FUNCTION. * S PS0L02 LINE SEARCH USING ONLY FUNCTION VALUES. * S PYPTSH DETERMINATION OF GROUPS FOR NUMERICAL DIFFERENTIATION. * S MXSCMD MATRIX-VECTOR PRODUCT FOLLOWED BY THE ADDITION OF A * SCALED VECTOR. SPARSE RECTANGULAR MATRIX IS STORED COLUMNWISE. * S MXSPCC SPARSE MATRIX REORDERING, SYMBOLIC FACTORIZATION, DATA * STRUCTURES TRANSFORMATION. INITIATION OF THE DIRECT SPARSE * SOLVER. * S MXSPCF GILL-MURRAY DECOMPOSITION OF A SPARSE SYMMETRIC MATRIX. * S MXSPCT COPYING A SPARSE SYMMETRIC MATRIX INTO THE PERMUTED * FACTORIZED COMPACT SCHEME. * S MXSPCX GILL-MURRAY DECOMPOSITION OF THE DENSE PART OF THE * PRECONDITIONER. * S MXSPCZ BACK SUBSTITUTION USING THE DECOMPOSITION OBTAINED BY * MXSPCF AND MXSPCX. * S MXSPIF INCOMPLETE GILL-MURRAY DECOMPOSITION OF A SPARSE * SYMMETRIC MATRIX. * S MXSPIX INCOMPLETE GILL-MURRAY DECOMPOSITION OF THE DENSE PART OF * THE PRECONDITIONER. * S MXSPIZ BACK SUBSTITUTION USING THE DECOMPOSITION OBTAINED BY * MXSPIF AND MXSPIX. * S MXSRMD MATRIX-VECTOR PRODUCT FOLLOWED BY THE ADDITION OF A * SCALED VECTOR. SPARSE RECTANGULAR MATRIX IS STORED ROWWISE. * S MXSRMM MATRIX-VECTOR PRODUCT. SPARSE RECTANGULAR MATRIX IS * STORED ROWWISE. * S MXSSMM MATRIX-VECTOR PRODUCT. * S MXVCOP COPYING OF A VECTOR. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * S MXVDIF DIFFERENCE OF TWO VECTORS. * S MXVDOT DOT PRODUCT OF TWO VECTORS. * S MXVMUL MULTIPLICATION OF A VECTOR BY A DIAGONAL MATRIX. * S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. * S MXVSET INITIATION OF A VECTOR. * S MXVSUM SUM OF TWO VECTORS, * S MXWDOT DOT PRODUCT OF TWO VECTORS IN THE PACKED CASE. * * EXTERNAL SUBROUTINES : * SE OBJ COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION. * CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS A NUMBER * OF VARIABLES, X(NF) IS A VECTOR OF VARIABLES AND FF IS THE * VALUE OF THE OBJECTIVE FUNCTION. * SE DOBJ COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION. * CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS A NUMBER * OF VARIABLES, X(NF) IS A VECTOR OF VARIABLES AND GF(NF) IS * THE GRADIENT OF THE OBJECTIVE FUNCTION. * SE CON COMPUTATION OF THE VALUE OF THE CONSTRAINT FUNCTION. * CALLING SEQUENCE: CALL CON(NF,KC,X,FC) WHERE NF IS A NUMBER * OF VARIABLES, KC IS THE INDEX OF THE CONSTRAINT FUNCTION, * X(NF) IS A VECTOR OF VARIABLES AND FC IS THE VALUE OF THE * CONSTRAINT FUNCTION. * SE DCON COMPUTATION OF THE GRADIENT OF THE CONSTRAINT FUNCTION. * CALLING SEQUENCE: CALL DCON(NF,KC,X,GC) WHERE NF IS A NUMBER * OF VARIABLES, KC IS THE INDEX OF THE CONSTRAINT FUNCTION, * X(NF) IS A VECTOR OF VARIABLES AND GC(NF) IS THE GRADIENT OF * THE CONSTRAINT FUNCTION. * * METHOD : * PRECONDITIONED SMOOTHED CONJUGATE GRADIENT METHOD WITH A SPECIAL * TERMINATION CRITERION GUARANTEEING SUFFICIENT DESCENT OF THE * DIRECTION VECTOR. * SUBROUTINE PIND(NF,NC,MAXH,MAXB,X,GF,H,IH,JH,CG,ICG,JCG,KCG,CF, & CFO,GC,G,S,XO,GO,XS,CZ,CZO,CP,B,IB,JB,CGD,DD,BD,COL,PSR,PERC, & INVP,WN11,WN12,WN13,WN14,XMAX,TOLX,TOLC,TOLG,RPF1,F,GMAX,CMAX, & MIT,MFV,MFG,MOS,MD,MDE,IPRNT,ITERM) INTEGER NF,NC,MAXH,MAXB,IH(*),JH(*),ICG(*),JCG(*),KCG(*),IB(*), & JB(*),COL(*),PSR(*),PERC(*),INVP(*),WN11(*),WN12(*),WN13(*), & WN14(*),IPRNT,MIT,MFV,MFG,MES,MOS,MD,MDE,ITERM DOUBLE PRECISION X(*),GF(*),H(*),CG(*),CF(*),CFO(*),GC(*),G(*), & S(*),XO(*),GO(*),XS(*),CZ(*),CZO(*),CP(*),B(*),CGD(*),DD(*), & BD(*),XMAX,TOLX,TOLC,TOLG,TOLD,TOLS,RPF1,F,GMAX,CMAX LOGICAL L1 INTEGER KD,LD,IEST,INITS,ITERD,ITERS,NTESX,MTESX,MTESF,MRED, & NRED,IREST,KIT,I,J,K,L,M,JSTRT,JSTOP,KC,ISNC,MS,MH,MM,MMX,ND, & INF,INITD,KTERS,MAXST,ISYS DOUBLE PRECISION ETA0,ETA1,ETA2,ETA9,EPS6,SNORM,GNORM,R,RO,RP,FO, & FP,P,PO,PP,RMIN,RMAX,FMIN,FMAX,DMAX,FF,FC,POM,ALF1,ALF2,RPF2, & TOLB,PAR DOUBLE PRECISION RHO,RHO1,RHO2,SIG,SIG1,SIG2,ALF,BTB,BTR,MXVDOT, & MXWDOT INTEGER NRES,NDEC,NIIT,NIT,NFV,NFG,NFH COMMON /STAT/ NRES,NDEC,NIIT,NIT,NFV,NFG,NFH IF (ABS(IPRNT).GT.1) WRITE(6,'(1X,''ENTRY TO PIND :'')') * * INITIATION * NRES=0 NDEC=0 NIIT=0 NIT=0 NFV=0 NFG=0 NFH=0 ISYS=0 IEST=0 MTESX=2 MTESF=2 INITS=1 INITD=1 ITERM=0 ITERD=0 ITERS=2 IREST=0 KTERS=5 MRED=20 ETA0=1.0D-15 ETA1=1.0D-15 ETA2=5.0D-5 ETA9=1.0D 60 EPS6=0.9D 0 ALF1=1.0D-10 ALF2=1.0D 10 FMAX=1.0D 60 FMIN=-FMAX TOLB=-FMAX DMAX=ETA9 IF (TOLX.LE.0.0D 0) TOLX=1.0D-16 IF (TOLC.LE.0.0D 0) TOLC=1.0D-6 IF (TOLG.LE.0.0D 0) TOLG=1.0D-6 TOLD=1.0D-5 TOLS=1.0D-4 IF (XMAX.LE.0.0D 0) XMAX=1.0D 16 IF (RPF1.LE.0.0D 0) RPF1=1.0D-4 IF (MOS.EQ.0) MOS=1 MES=1 IF (MD .LE.0) MD =10 IF (MDE.LE.0) MDE=50 IF (MIT.LE.0) MIT=1000 IF (MFV.LE.0) MFV=1000 IF (MFG.LE.0) MFG=10000 LD=-1 KIT=0 FO=FMIN CALL MXVSET(NC,0.0D 0,CZ) C CALL PFSET3(NF,NC,IH,JH,ICG,JCG) CALL PFSED7(NF,NC,ND,MD,MDE,IB,JB,ICG,JCG,KCG,ITERM) * * SELECTION OF COLUMN GROUPS FOR GROUPWISE DIFFERENTIATION * CALL PYPTSH(NF,MAXH,IH,JH,COL,S,XO,GO,WN11,WN12,GF,ITERM) MH=0 M=IH(NF+1)-1 CALL MXVINS(M,0,JH(M+1)) IF (ITERM.NE.0) GO TO 90 IF (MOS.EQ.1) THEN * * SYMBOLIC GILL-MURRAY DECOMPOSITION * MS=0 L=IB(NC+1)-1 CALL MXSPCC(NC,L,MS,MAXB,B,IB,JB,PSR,PERC,INVP,WN11, & WN12,WN13,WN14,ITERM) IF (ITERM.NE.0) GO TO 90 END IF 10 CONTINUE * * COMPUTATION OF THE VALUE AND THE GRADIENT AND THE HESSIAN MATRIX * OF THE LAGRANGIAN FUNCTION * KD=1 ISNC=1 CALL PP1LF3(NF,NC,X,GF,CF,CG,ICG,JCG,GC,G,CZ,FF,F,CMAX,KD,LD, & NFV,NFG,ISNC) CF(NC+1)=FF 11 CONTINUE CALL PF1HS2(NF,MH,MAXH,X,IH,XO,H,IH,JH,GO,G,COL,WN11,WN12,XS, & FF,ETA0,0,ITERM,ISYS) IF (ISYS.GT.0) THEN KD=1 LD=0 ISNC=-1 CALL PP1LF3(NF,NC,X,GF,CF,CG,ICG,JCG,GC,GO,CZ,FF,F,CMAX,KD,LD, & NFV,NFG,ISNC) GO TO 11 END IF GMAX = 0.0D 0 DO 20 I=1,NF GMAX=MAX(GMAX,ABS(G(I))) 20 CONTINUE IF (ABS(IPRNT).GT.1) & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X, & ''C='', G13.6,2X,''G='',G13.6)') NIT,NFV,NFG,CMAX,GMAX * * START OF THE ITERATION WITH TESTS FOR TERMINATION. * IF (ITERM.NE.0) GO TO 90 IF (ITERS.EQ.0) GO TO 30 IF (GMAX.LE.TOLG.AND.CMAX.LE.TOLC) THEN ITERM=4 GO TO 90 END IF IF (F.LE.TOLB) THEN ITERM = 3 GO TO 90 END IF IF ( NIT .GT. 0) THEN IF ( DMAX .LE. TOLX) THEN ITERM = 1 NTESX = NTESX + 1 IF ( NTESX .GE. MTESX) GO TO 90 ELSE NTESX = 0 END IF END IF IF ( NIT .GE. MIT) THEN ITERM = 11 GO TO 90 END IF IF ( NFV .GE. MFV) THEN ITERM = 12 GO TO 90 END IF IF ( NFG .GE. MFG) THEN ITERM = 13 GO TO 90 END IF 30 ITERM = 0 NIT = NIT + 1 IF (ITERM.NE.0) GO TO 90 40 CONTINUE * * RESTART * IF(IREST.GT.0) THEN POM=GMAX/1.0D 1 DO 60 I=1,NF JSTRT=IH(I) JSTOP=IH(I+1)-1 H(JSTRT)=MIN(MAX(POM*ABS(H(JSTRT)),5.0D-3),5.0D 2) DO 50 J=JSTRT+1,JSTOP H(J)=0.0D 0 50 CONTINUE 60 CONTINUE 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 GO TO 90 END IF END IF IF (ITERM.NE.0) GO TO 90 * * DIRECTION DETERMINATION USING INDEFINITELY PRECONDITIONED * CONJUGATE GRADIENT METHOD APPLIED TO THE KARUSH-KUHN-TUCKER * SYSTEM OF LINEAR EQUATIONS * RHO1=MXVDOT(NF,G,G) SIG1=MXVDOT(NC,CF,CF) CALL MXVNEG(NF,G,XS) RHO=RHO1 SIG=SIG1 * * CONSTRUCTION OF PRECONDITIONER * L=IB(NC+1)-1 DO 14240 I=1,NF ALF=ABS(H(IH(I))) IF (ABS(ALF).LE.1.0D-3) ALF=1.0D-3 IF (ABS(ALF).GE.1.0D 6) ALF=1.0D 6 GC(I)=1.0D 0/ALF 14240 CONTINUE CALL PFSED9(NC,ND,B,IB,JB,CG,ICG,JCG,KCG,CGD,DD,GC) ALF=1.0D-15 IF (MOS.EQ.1) THEN CALL MXSPCT(NC,L,MS,MAXB,B,JB,PSR,ITERM) IF (ITERM.NE.0) GO TO 14090 * * GILL-MURRAY DECOMPOSITION * CALL MXSPCF(NC,B(L+1),PSR,JB(L+1),WN11,WN12,CP,INF,ALF,SIG) IF (ND.GT.0) THEN CALL MXSPCX(NC,ND,B(L+1),PSR,JB(L+1),PERC,CGD,DD,BD,S,S(NF+1) & ) END IF ELSE CALL MXSPIF(NC,B,IB,JB,INF,ALF,SIG,GO) IF (ND.GT.0) THEN CALL MXSPIX(NC,ND,B,IB,JB,CGD,DD,BD,S(NF+1)) END IF END IF NDEC=NDEC+1 MM=NF+NC GNORM=SQRT(RHO1+SIG1) INITD=MAX(ABS(INITD),1) PAR=MIN(EPS6,GNORM) IF (PAR.GT.1.0D 1*1.0D-3) THEN PAR=MIN(PAR,1.0D 0/DBLE(NIT)) END IF PAR=PAR*PAR * * CG INITIATION * CALL MXVNEG(NC,CF,XS(NF+1)) CALL MXVSET(MM,0.0D 0,S) * * CG PRECONDITIONING * CALL MXVMUL(NF,GC,XS,XO,1) CALL MXVCOP(NC,XS(NF+1),XO(NF+1)) CALL MXSRMD(NC,CG,ICG,JCG,XO,-1.0D 0,XO(NF+1),XO(NF+1)) IF (MOS.EQ.1) THEN CALL MXSPCZ(NC,ND,B(L+1),PSR,JB(L+1),PERC,CGD,DD,BD,XO,XO(NF+1), & CP) ELSE CALL MXSPIZ(NC,ND,B,IB,JB,CGD,DD,BD,XO(NF+1),CP) END IF CALL MXSCMD(NF,NC,CG,ICG,JCG,XO(NF+1),-1.0D 0,XS,XO) CALL MXVNEG(NF,XO,XO) CALL MXVMUL(NF,GC,XO,XO,1) RHO=MXVDOT(NF,XS,XO) SIG=MXVDOT(NC,XS(NF+1),XO(NF+1)) MMX=MM+3 * * CG ITERATIONS * DO 14015 NRED=1,MMX IF (RHO+SIG.EQ.0.0D 0) THEN ITERD=-4 GO TO 14090 END IF CALL MXSSMM(NF,H,IH,JH,XO,GO) CALL MXSCMD(NF,NC,CG,ICG,JCG,XO(NF+1),1.0D 0,GO,GO) CALL MXSRMM(NC,CG,ICG,JCG,XO,GO(NF+1)) ALF=MXVDOT(MM,XO,GO) IF (ALF.EQ.0.0D 0) THEN ITERD=-5 GO TO 14090 END IF * * CG STEP * ALF=(RHO+SIG)/ALF CALL MXVDIR(MM, ALF,XO,S,S) CALL MXVDIR(MM,-ALF,GO,XS,XS) NIIT=NIIT+1 * * TERMINATION CRITERION * L1=NRED.GE.MMX SIG2=MXVDOT(NC,XS(NF+1),XS(NF+1)) RHO2=MXVDOT(NF,XS,XS) IF (SIG2.LE.PAR*MAX(SIG1,1.0D 0/1.0D 60)) THEN IF (RHO2.LE.PAR*RHO1) L1=.TRUE. END IF RPF2=RPF1 IF (L1) THEN GO TO 14016 END IF * * CG PRECONDITIONING * CALL MXVMUL(NF,GC,XS,GO,1) CALL MXVCOP(NC,XS(NF+1),GO(NF+1)) CALL MXSRMD(NC,CG,ICG,JCG,GO,-1.0D 0,GO(NF+1),GO(NF+1)) IF (MOS.EQ.1) THEN CALL MXSPCZ(NC,ND,B(L+1),PSR,JB(L+1),PERC,CGD,DD,BD,GO,GO(NF+1), & CP) ELSE CALL MXSPIZ(NC,ND,B,IB,JB,CGD,DD,BD,GO(NF+1),CP) END IF CALL MXSCMD(NF,NC,CG,ICG,JCG,GO(NF+1),-1.0D 0,XS,GO) CALL MXVNEG(NF,GO,GO) CALL MXVMUL(NF,GC,GO,GO,1) BTB=ABS(SIG2) BTR=ABS(RHO2+SIG2) IF (BTR.LE.1.0D 0/ETA9) GO TO 14016 RHO2=MXVDOT(NF,XS,GO) SIG2=MXVDOT(NC,XS(NF+1),GO(NF+1)) * * SPECIAL TEST FOR NULL SPACE DIRECTION * IF (BTB.LE.PAR*(RHO1+SIG1).AND.ABS(RHO2+SIG2).LE.1.0D-12*BTR) THEN CALL MXVMUL(NF,GC,XS,GO,1) CALL MXSRMM(NC,CG,ICG,JCG,GO,GO(NF+1)) IF (MOS.EQ.1) THEN CALL MXSPCZ(NC,ND,B(L+1),PSR,JB(L+1),PERC,CGD,DD,BD,GO,GO(NF+1), & CP) ELSE CALL MXSPIZ(NC,ND,B,IB,JB,CGD,DD,BD,GO(NF+1),CP) END IF CALL MXVSUM(NC,S(NF+1),GO(NF+1),S(NF+1)) CALL MXSSMM(NF,H,IH,JH,S,GO) CALL MXSCMD(NF,NC,CG,ICG,JCG,S(NF+1),1.0D 0,GO,GO) CALL MXVSUM(NF,GO,G,GO) CALL MXSRMD(NC,CG,ICG,JCG,S,1.0D 0,CF,GO(NF+1)) RHO2=MXVDOT(NF,GO,GO) SIG2=MXVDOT(NC,GO(NF+1),GO(NF+1)) GO TO 14016 END IF ALF=(RHO2+SIG2)/(RHO+SIG) CALL MXVDIR(MM,ALF,XO,GO,XO) RHO=RHO2 SIG=SIG2 14015 CONTINUE * * AN INEXACT SOLUTION IS OBTAINED * 14016 CONTINUE * * COMPUTATION OF THE DIRECTIONAL DERIVATIVE * ITERD=2 GNORM=SQRT(RHO1) SNORM=SQRT(MXVDOT(NF,S,S)) P=MXVDOT(NF,G,S) DO 14770 KC=1,NC K=ICG(KC) L=ICG(KC+1)-K ALF=MXWDOT(L,JCG(K),CG(K),S,2) P=P+ALF*(S(NF+KC)+RPF2*CF(KC)) 14770 CONTINUE 14090 CONTINUE IF (ITERD.LT.0) ITERM=ITERD IF (ITERM.NE.0) GO TO 90 * * TEST FOR SUFFICIENT DESCENT * IREST=1 IF(SNORM.LE.0.0D 0) THEN ELSEIF(P+TOLD*GNORM*SNORM.LE.0.0D 0) THEN IREST=0 END IF IF (IREST.EQ.0) THEN NRED=0 RMIN=MIN(1.0D 0/1.0D 1,ALF1*GNORM/SNORM) RMAX=MIN(MAX(1.0D 1,ALF2*GNORM/SNORM),XMAX/SNORM) ELSE GO TO 40 END IF CALL MXVCOP(NC,CZ,CZO) CALL PP0AF1(NF,NC,CF,CZ,S,RPF2,FC,F) IF (ITERM.NE.0) GO TO 90 IF (IREST.NE.0) GO TO 40 IF (NIT.EQ.1) KIT=NIT * * PREPARATION OF LINE SEARCH * FP=FO RO=0.0D 0 FO=F PO=P CALL MXVCOP(NF,X,XO) CALL MXVCOP(NF,G,GO) CALL MXVCOP(NC+1,CF,CFO) * * LINE SEARCH WITHOUT DIRECTIONAL DERIVATIVES * 70 CONTINUE CALL PS0L02(R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,TOLS,KD,LD, & NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,MES,ISYS) IF (ISYS.GT.0) THEN ISNC=0 CALL MXVDIR(NF,R,S,XO,X) CALL PP1LF3(NF,NC,X,GF,CF,CG,ICG,JCG,GC,G,CZ,FF,F,CMAX,KD,LD, & NFV,NFG,ISNC) CF(NC+1)=FF CALL PP0AF1(NF,NC,CF,CZ,S,RPF2,FC,F) GO TO 70 END IF CALL MXVDIR(NC,R,S(NF+1),CZO,CZ) * * DECISION AFTER UNSUCCESSFUL LINE SEARCH * IF (ITERS.LE.0) THEN R=0.0D 0 F=FO P=PO CALL MXVCOP(NF,XO,X) FF=F CALL MXVCOP(NC+1,CFO,CF) CALL MXVCOP(NC,CZO,CZ) IREST=1 GO TO 40 END IF CALL MXVDIF(NF,X,XO,XO) PO=R*PO P=R*P DMAX=0.0D 0 DO 80 I=1,NF DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0)) 80 CONTINUE GO TO 10 90 IF (IPRNT.GT.1.OR.IPRNT.LT.0) & WRITE(6,'(1X,''EXIT FROM PIND :'')') IF (IPRNT.NE.0) & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X, & ''C='', G13.6,2X,''G='',G13.6,2X,''ITERM='',I3)') NIT,NFV,NFG, & CMAX,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