************************************************************************ * R port by Yohan Chalabi / 2008-01-29 * - added IMPLICIT NONE statement in all subroutines and functions * - removed common block /STAT/ and added array STAT as an argument in * subroutines : * VARIABLES IN ARRAY STAT (STATISTICS) : * IO STAT(1) NRES NUMBER OF RESTARTS. * IO STAT(2) NDEC NUMBER OF MATRIX DECOMPOSITION. * IO STAT(3) NREM NUMBER OF CONSTRAINT DELETIONS. * IO STAT(4) NADD NUMBER OF CONSTRAINT ADDITIONS. * IO STAT(5) NIT NUMBER OF ITERATIONS. * IO STAT(6) NFV NUMBER OF FUNCTION EVALUATIONS. * IO STAT(7) NFG NUMBER OF GRADIENT EVALUATIONS. * IO STAT(8) NFH NUMBER OF HESSIAN EVALUATIONS. * - changed WRITE statement to DBLEPR for R compatibitly * - added RCHKUSR() to allow interrupts from R ************************************************************************ ************************************************************************ * Copyright: * * Subroutines PMIN, PBUN, PNEW, PVAR, Copyright ACM, 2001. The original * versions were published in Transactions on Mathematical Software, * Vol.27, 2001, pp.193-213. Here are the author's modifications. They * are posted here by permission of ACM for your personal use. Not for * redistribution. Subroutines PLIP, PSEN, Copyright Jan Vlcek, 2007. * The remaining subroutines, Copyright Ladislav Luksan, 2007. Many of * sparse matrix modules were prepared by Miroslav Tuma. * * License: * * This library (with exception of PMIN, PBUN, PNEW, PVAR) is a free * software; you can redistribute it and/or modify it under the terms * of the GNU Lesser General Public License as published by the Free * Software Foundation; either version 2.1 of the License, or (at your * option) any later version (see http://www.gnu.org/copyleft/gpl.html). * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Lesser General Public License for more details. * * Permission is hereby granted to use or copy this program under the * terms of the GNU LGPL, provided that the Copyright, this License, * and the Availability of the original version is retained on all * copies. User documentation of any code that uses this code or any * modified version of this code must cite the Copyright, this License, * the Availability note, and "Used by permission." Permission to modify * the code and to distribute modified code is granted, provided the * Copyright, this License, and the Availability note are retained, * and a notice that the code was modified is included. * * Availability: * * http://www.cs.cas.cz/~luksan/subroutines.html * * Acknowledgements: * * This work was supported by the Grant Agency of the Czech Academy of * Sciences, under grant IAA1030405. ************************************************************************ ************************************************************************ * SUBROUTINE PSQPN ALL SYSTEMS 97/01/22 * PURPOSE : * EASY TO USE SUBROUTINE FOR GENERAL NONLINEAR PROGRAMMING PROBLEMS. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NB CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED. * NB>0-SIMPLE BOUNDS ACCEPTED. * II NC NUMBER OF LINEAR CONSTRAINTS. * 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-LOVER 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 CF(NC+1) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * IC(KC)=0-CONSTRAINT CF(KC) IS NOT USED. IC(KC)=1-LOVER * CONSTRAINT CL(KC).LE.CF(KC). IC(KC)=2-UPPER CONSTRAINT * CF(KC).LE.CU(KC). IC(KC)=3-TWO SIDE CONSTRAINT * CL(KC).LE.CF(KC).LE.CU(KC). IC(KC)=5-EQUALITY CONSTRAINT * CF(KC).EQ.CL(KC). * RI CL(NC) VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CU(NC) VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS. * II IPAR(6) INTEGER PAREMETERS: * IPAR(1) MAXIMUM NUMBER OF ITERATIONS. * IPAR(2) MAXIMUM NUMBER OF FUNCTION EVALUATIONS. * IPAR(3) THIS PARAMETER IS NOT USED IN THE SUBROUTINE PSQP. * IPAR(4) THIS PARAMETER IS NOT USED IN THE SUBROUTINE PSQP. * IPAR(5) VARIABLE METRIC UPDATE USED. IPAR(5)=1-THE BFGS UPDATE. * IPAR(5)-THE HOSHINO UPDATE. * IPAR(6) CORRECTION OF THE VARIABLE METRIC UPDATE IF A NEGATIVE * CURVATURE OCCURS. IPAR(6)=1-NO CORRECTION. IPAR(6)=2-POWELL'S * CORRECTION. * RI RPAR(5) REAL PARAMETERS: * RPAR(1) MAXIMUM STEPSIZE. * RPAR(2) TOLERANCE FOR CHANGE OF VARIABLES. * RPAR(3) TOLERANCE FOR CONSTRAINT VIOLATIONS. * RPAR(4) TOLERANCE FOR THE GRADIENT OF THE LAGRANGIAN FUNCTION. * RPAR(5) PENALTY COEFFICIENT. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RO GMAX MAXIMUM PARTIAL DERIVATIVE OF THE LAGRANGIAN FUNCTION. * 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. * * VARIABLES IN ARRAY STAT (STATISTICS) : * IO STAT(1) NRES NUMBER OF RESTARTS. * IO STAT(2) NDEC NUMBER OF MATRIX DECOMPOSITION. * IO STAT(3) NREM NUMBER OF CONSTRAINT DELETIONS. * IO STAT(4) NADD NUMBER OF CONSTRAINT ADDITIONS. * IO STAT(5) NIT NUMBER OF ITERATIONS. * IO STAT(6) NFV NUMBER OF FUNCTION EVALUATIONS. * IO STAT(7) NFG NUMBER OF GRADIENT EVALUATIONS. * IO STAT(8) NFH NUMBER OF HESSIAN EVALUATIONS. * * SUBPROGRAMS USED : * S PSQP RECURSIVE QUADRATIC PROGRAMMING METHOD WITH THE BFGS * VARIABLE METRIC UPDATE. * * EXTERNAL SUBROUTINES : * SE OBJ COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION. * CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE 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 THE NUMBER * OF VARIABLES, X(NF) IS A VECTOR OF VARIABLES AND GC(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 THE * 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 THE * 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. * SUBROUTINE PSQPN(NF,NB,NC,X,IX,XL,XU,CF,IC,CL,CU,IPAR,RPAR, + F,GMAX,CMAX,IPRNT,ITERM, STAT) * * POINTERS FOR AUXILIARY ARRAYS * IMPLICIT NONE DOUBLE PRECISION F,CMAX,GMAX INTEGER IPRNT,ITERM,NB,NC,NF DOUBLE PRECISION CF(*),CL(*),CU(*),RPAR(5),X(*),XL(*),XU(*) INTEGER IC(*),IPAR(6),IX(*) c INTEGER NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES INTEGER LCFD,LCFO,LCG,LCP,LCR,LCZ,LG,LGC,LGF,LGO,LH,LIA, + LS,LXO INTEGER STAT(8) c COMMON /STAT/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH INTEGER IA(:) DOUBLE PRECISION RA(:) ALLOCATABLE IA,RA ALLOCATE (IA(NF),RA((NF+NC+8)*NF+3*NC+1)) LCG = 1 LCFO = LCG + NF*NC LCFD = LCFO + NC + 1 LGC = LCFD + NC LCR = LGC + NF LCZ = LCR + NF*(NF+1)/2 LCP = LCZ + NF LGF = LCP + NC LG = LGF + NF LH = LG + NF LS = LH + NF* (NF+1)/2 LXO = LS + NF LGO = LXO + NF LIA = 1 CALL PSQP(NF,NB,NC,X,IX,XL,XU,CF,IC,CL,CU,RA,RA(LCFO),RA(LCFD), + RA(LGC),IA,RA(LCR),RA(LCZ),RA(LCP),RA(LGF),RA(LG), + RA(LH),RA(LS),RA(LXO),RA(LGO),RPAR(1),RPAR(2),RPAR(3), + RPAR(4),RPAR(5),CMAX,GMAX,F,IPAR(1),IPAR(2),IPAR(5), + IPAR(6),IPRNT,ITERM, STAT) DEALLOCATE(IA,RA) RETURN END ************************************************************************ * SUBROUTINE PSQP ALL SYSTEMS 97/01/22 * PURPOSE : * RECURSIVE QUADRATIC PROGRAMMING METHOD WITH THE BFGS VARIABLE METRIC * UPDATE FOR GENERAL NONLINEAR PROGRAMMING PROBLEMS. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NB CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED. * NB>0-SIMPLE BOUNDS ACCEPTED. * II NC NUMBER OF LINEAR CONSTRAINTS. * 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-LOVER 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 CF(NC+1) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * IC(KC)=0-CONSTRAINT CF(KC) IS NOT USED. IC(KC)=1-LOVER * CONSTRAINT CL(KC).LE.CF(KC). IC(KC)=2-UPPER CONSTRAINT * CF(KC).LE.CU(KC). IC(KC)=3-TWO SIDE CONSTRAINT * CL(KC).LE.CF(KC).LE.CU(KC). IC(KC)=5-EQUALITY CONSTRAINT * CF(KC).EQ.CL(KC). * 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. * RA CFO(NC) VECTOR CONTAINING SAVED VALUES OF THE CONSTRAINT * FUNCTIONS. * RA CFD(NC) VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT * FUNCTIONS. * RA GC(NF) GRADIENT OF THE SELECTED CONSTRAINT FUNCTION. * IO ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RO CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RO CZ(NF) VECTOR OF LAGRANGE MULTIPLIERS. * RO GF(NF) GRADIENT OF THE MODEL FUNCTION. * RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RU H(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OR INVERSION OF THE * HESSIAN MATRIX APPROXIMATION. * RO S(NF) DIRECTION VECTOR. * RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. * RI GO(NF) GRADIENTS DIFFERENCE. * RI XMAX MAXIMUM STEPSIZE. * RI TOLX TOLERANCE FOR CHANGE OF VARIABLES. * RI TOLC TOLERANCE FOR CONSTRAINT VIOLATIONS. * RI TOLG TOLERANCE FOR THE GRADIENT OF THE LAGRANGIAN FUNCTION. * RI RPF PENALTY COEFFICIENT. * RO CMAX MAXIMUM CONSTRAINT VIOLATION. * RO GMAX MAXIMUM PARTIAL DERIVATIVE OF THE LAGRANGIAN FUNCTION. * RO F VALUE OF THE OBJECTIVE FUNCTION. * FUNCTIONS. * II MIT MAXIMUN NUMBER OF ITERATIONS. * II MFV MAXIMUN NUMBER OF FUNCTION EVALUATIONS. * II MET VARIABLE METRIC UPDATE USED. MET=1-THE BFGS UPDATE. * MET=2-THE HOSHINO UPDATE. * II MEC CORRECTION IF THE NEGATIVE CURVATURE OCCURS. * MEC=1-CORRECTION SUPPRESSED. MEC=2-POWELL'S CORRECTION. * 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. * * VARIABLES IN ARRAY STAT (STATISTICS) : * IO STAT(1) NRES NUMBER OF RESTARTS. * IO STAT(2) NDEC NUMBER OF MATRIX DECOMPOSITION. * IO STAT(3) NREM NUMBER OF CONSTRAINT DELETIONS. * IO STAT(4) NADD NUMBER OF CONSTRAINT ADDITIONS. * IO STAT(5) NIT NUMBER OF ITERATIONS. * IO STAT(6) NFV NUMBER OF FUNCTION EVALUATIONS. * IO STAT(7) NFG NUMBER OF GRADIENT EVALUATIONS. * IO STAT(8) NFH NUMBER OF HESSIAN EVALUATIONS. * * SUBPROGRAMS USED : * S PC1F01 COMPUTATION OF THE VALUE AND THE GRADIENT OF THE * CONSTRAINT FUNCTION. * S PF1F01 COMPUTATION OF THE VALUE AND THE GRADIENT OF THE * OBJECTIVE FUNCTION. * S PLQDB1 GENERAL QUADRATIC PROGRAMMING SUBROUTINE BASED ON THE * GOLDFARB-IDNANI DUAL METHOD. * S PLNEWS IDENTIFICATION OF ACTIVE SIMPLE BOUNDS. * S PLREDL TRANSFORMATION OF THE INCOMPATIBLE QUADRATIC PROGRAMMING * SUBPROBLEM. * S PP0AF8 COMPUTATION OF VALUE OF THE AUGMENTED LAGRANGIAN * FUNCTION. * S PPSET2 COMPUTATION OF THE NEW PENALTY PARAMETERS. * S PS0L02 LINE SEARCH USING ONLY FUNCTION VALUES. * S PYTRND DETERMINATION OF DIFFERENCES FOR VARIABLE METRIC * UPDATES. * S PUDBG1 VARIABLE METRIC UPDATE AFTER GILL-MURRAY DECOMPOSITION. * S MXDSMI SYMMETRIC MATRIX IS REPLACED BY THE UNIT 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 AN INTEGER VECTOR. * RF MXVMAX L-INFINITY NORM OF A VECTOR. * S MXVSET INITIATION OF A VECTOR. * * EXTERNAL SUBROUTINES : * SE OBJ COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION. * CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE 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 THE NUMBER * OF VARIABLES, X(NF) IS A VECTOR OF VARIABLES AND GC(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 THE * 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 THE * 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 : * RECURSIVE QUADRATIC PROGRAMMING METHOD WITH THE BFGS VARIABLE METRIC * UPDATE. * SUBROUTINE PSQP(NF,NB,NC,X,IX,XL,XU,CF,IC,CL,CU,CG,CFO,CFD,GC, + ICA,CR,CZ,CP,GF,G,H,S,XO,GO,XMAX,TOLX,TOLC,TOLG, + RPF,CMAX,GMAX,F,MIT,MFV,MET,MEC,IPRNT,ITERM, STAT) IMPLICIT NONE DOUBLE PRECISION F,CMAX,GMAX,RPF,TOLC,TOLD,TOLG,TOLS,TOLX,XMAX INTEGER IPRNT,ITERM,MET,MET1,MEC,MES,MFV,MIT,NB,NC,NF DOUBLE PRECISION CF(*),CFD(*),CFO(*),CG(*),CL(*),CP(*),CR(*), + CZ(*),CU(*),G(*),GC(*),GF(*),GO(*),H(*),S(*), + X(*),XL(*),XO(*),XU(*) INTEGER IC(*),ICA(*),IX(*) c INTEGER NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES DOUBLE PRECISION ALF1,ALF2,CMAXO,DMAX,EPS7,EPS9,ETA0,ETA2,ETA9, + FMAX,FMIN,FO,GNORM,P,PO,R,RMAX,RMIN,RO,SNORM, + TOLB,TOLF,UMAX,RP,FP,PP,FF,FC INTEGER I,IDECF,IEXT,IREST,ITERD,ITERL,ITERH,ITERQ,ITERS, + KBC,KBF,KC,KD,KIT,LD,MRED,MTESF,MTESX,N,K, + NTESX,IEST,INITS,KTERS,MAXST,ISYS,MFP,NRED,IPOM, + LDS DOUBLE PRECISION MXVDOT, MXVMAX INTEGER STAT(8) c COMMON /STAT/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH C IF (ABS(IPRNT).GT.1) WRITE (6,'(1X,''ENTRY TO PSQP :'')') * * INITIATION * KBF = 0 KBC = 0 IF (NB.GT.0) KBF = 2 IF (NC.GT.0) KBC = 2 STAT(1) = 0 STAT(2) = 0 STAT(3) = 0 STAT(4) = 0 STAT(5) = 0 STAT(6) = 0 STAT(7) = 0 STAT(8) = 0 ISYS = 0 IEST = 0 IEXT = 0 MTESX = 2 MTESF = 2 INITS = 1 ITERM = 0 ITERS = 0 ITERD = 0 ITERQ = 0 MRED = 20 IREST = 1 ITERS = 2 KTERS = 5 IDECF = 1 ETA0 = 1.0D-15 ETA2 = 1.0D-15 ETA9 = 1.0D60 EPS7 = 1.0D-15 EPS9 = 1.0D-8 ALF1 = 1.0D-10 ALF2 = 1.0D10 FMAX = 1.0D60 FMIN = -FMAX TOLB = -FMAX DMAX = ETA9 TOLF = 1.0D-16 IF (XMAX.LE.0.0D0) XMAX = 1.0D+16 IF (TOLX.LE.0.0D0) TOLX = 1.0D-16 IF (TOLG.LE.0.0D0) TOLG = 1.0D-6 IF (TOLC.LE.0.0D0) TOLC = 1.0D-6 TOLD = 1.0D-8 TOLS = 1.0D-4 IF (RPF.LE.0.0D0) RPF = 1.0D-4 IF (MET.LE.0) MET = 1 MET1 = 2 IF (MEC.LE.0) MEC = 2 MES = 1 IF (MIT.LE.0) MIT = 1000 IF (MFV.LE.0) MFV = 2000 KD = 1 LD = -1 KIT = 0 CALL MXVSET(NC,0.0D0,CP) * * INITIAL OPERATIONS WITH SIMPLE BOUNDS * IF (KBF.GT.0) THEN DO 20 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 IF (IX(I).EQ.1 .OR. IX(I).EQ.3) X(I) = MAX(X(I),XL(I)) IF (IX(I).EQ.2 .OR. IX(I).EQ.3) X(I) = MIN(X(I),XU(I)) 20 CONTINUE END IF * INITIAL OPERATIONS WITH GENERAL CONSTRAINTS * IF (KBC.GT.0) THEN K = 0 DO 30 KC = 1,NC IF ((IC(KC).EQ.3.OR.IC(KC).EQ.4) .AND. + CU(KC).LE.CL(KC)) THEN CU(KC) = CL(KC) IC(KC) = 5 ELSE IF (IC(KC).EQ.5 .OR. IC(KC).EQ.6) THEN CU(KC) = CL(KC) IC(KC) = 5 END IF K = K + NF 30 CONTINUE END IF IF (KBF.GT.0) THEN DO 31 I = 1,NF IF (IX(I).GE.5) IX(I) = -IX(I) IF (IX(I).LE.0) THEN ELSE IF ((IX(I).EQ.1.OR.IX(I).EQ.3).AND.X(I).LE.XL(I)) + THEN X(I) = XL(I) ELSE IF ((IX(I).EQ.2.OR.IX(I).EQ.3).AND.X(I).GE.XU(I)) + THEN X(I) = XU(I) END IF CALL PLNEWS(X,IX,XL,XU,EPS9,I,ITERL) IF (IX(I).GT.10) IX(I) = 10 - IX(I) 31 CONTINUE END IF FO = FMIN GMAX = ETA9 DMAX = ETA9 40 CONTINUE LDS=LD CALL PF1F01(NF,X,GF,GF,FF,F,KD,LD,IEXT, STAT) LD=LDS CALL PC1F01(NF,NC,X,FC,CF,CL,CU,IC,GC,CG,CMAX,KD,LD) CF(NC+1)=F c IF (ABS(IPRNT).GT.1) c + WRITE (6,'(1X,''STAT(5)='',I4,2X,''STAT(6)='',I4,2X,''STAT(7)='',I4,2X, c + ''F='',G12.6,2X,''C='',E7.1,2X,''G='',E7.1)') c + STAT(5),STAT(6),STAT(7),F,CMAX,GMAX IF (IPRNT.EQ.1) THEN CALL DBLEPR("LLH improved to:", -1, F, 1) CALL DBLEPR("With X:", -1, X, NF) END IF c Allowing interrupts from R CALL RCHKUSR() * * START OF THE ITERATION WITH TESTS FOR TERMINATION. * IF (ITERM.LT.0) GO TO 80 IF (ITERS.EQ.0) GO TO 50 IF (F.LE.TOLB) THEN ITERM = 3 GO TO 80 END IF IF (DMAX.LE.TOLX) THEN ITERM = 1 NTESX = NTESX + 1 IF (NTESX.GE.MTESX) GO TO 80 ELSE NTESX = 0 END IF 50 IF (STAT(5).GE.MIT) THEN ITERM = 11 GO TO 80 END IF IF (STAT(6).GE.MFV) THEN ITERM = 12 GO TO 80 END IF ITERM = 0 STAT(5) = STAT(5) + 1 60 CONTINUE * * RESTART * N = NF IF (IREST.GT.0) THEN CALL MXDSMI(N,H) LD = MIN(LD,1) IDECF = 1 IF (KIT.LT.STAT(5)) THEN STAT(1) = STAT(1) + 1 KIT = STAT(5) ELSE ITERM = -10 IF (ITERS.LT.0) ITERM = ITERS - 5 GO TO 80 END IF END IF * * DIRECTION DETERMINATION USING A QUADRATIC PROGRAMMING PROCEDURE * CALL MXVCOP(NC+1,CF,CFO) MFP=2 IPOM=0 61 CONTINUE CALL PLQDB1(NF,NC,X,IX,XL,XU,CF,CFD,IC,ICA,CL,CU,CG,CR,CZ,G,GF, & H,S,MFP,KBF,KBC,IDECF,ETA2,ETA9,EPS7,EPS9,UMAX,GMAX, & N,ITERQ, STAT) IF (ITERQ.LT.0) THEN IF (IPOM.LT.10) THEN IPOM=IPOM+1 CALL PLREDL(NC,CF,IC,CL,CU,KBC) GO TO 61 END IF ITERD=ITERQ-10 GO TO 62 END IF IPOM=0 ITERD=1 GMAX=MXVMAX(NF,G) GNORM=SQRT(MXVDOT(NF,G,G)) SNORM=SQRT(MXVDOT(NF,S,S)) 62 CONTINUE IF (ITERD.LT.0) ITERM = ITERD IF (ITERM.NE.0) GO TO 80 CALL MXVCOP(NC+1,CFO,CF) * * TEST FOR SUFFICIENT DESCENT * P = MXVDOT(NF,G,S) IREST = 1 IF (SNORM.LE.0.0D0) THEN ELSE IF (P+TOLD*GNORM*SNORM.LE.0.0D0) THEN IREST = 0 END IF IF (IREST.EQ.0) THEN NRED = 0 RMIN = ALF1*GNORM/SNORM RMAX = MIN(ALF2*GNORM/SNORM,XMAX/SNORM) ELSE GO TO 60 END IF IF (GMAX.LE.TOLG.AND.CMAX.LE.TOLC) THEN ITERM=4 GO TO 80 ENDIF CALL PPSET2(NF,N,NC,ICA,CZ,CP) CALL MXVINA(NC,IC) CALL PP0AF8(NF,N,NC,CF,IC,ICA,CL,CU,CZ,RPF,FC,F) * * PREPARATION OF LINE SEARCH * RO = 0.0D0 FO = F PO = P CMAXO = CMAX CALL MXVCOP(NF,X,XO) CALL MXVCOP(NF,G,GO) CALL MXVCOP(NF,GF,CR) CALL MXVCOP(NC+1,CF,CFO) * * LINE SEARCH WITHOUT DIRECTIONAL DERIVATIVES * 11170 CONTINUE CALL PS0L02(R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX, & TOLS,KD,LD,STAT(5),KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS, & MES,ISYS) IF (ISYS.EQ.0) GO TO 11174 C GO TO (11174,11172) ISYS+1 11172 CONTINUE CALL MXVDIR(NF,R,S,XO,X) LDS=LD CALL PF1F01(NF,X,GF,G,FF,F,KD,LD,IEXT, STAT) LD=LDS CALL PC1F01(NF,NC,X,FC,CF,CL,CU,IC,GC,CG,CMAX,KD,LD) CF(NC+1)=F CALL PP0AF8(NF,N,NC,CF,IC,ICA,CL,CU,CZ,RPF,FC,F) GO TO 11170 11174 CONTINUE KD=1 * * DECISION AFTER UNSUCCESSFUL LINE SEARCH * IF (ITERS.LE.0) THEN R = 0.0D0 F = FO P = PO CALL MXVCOP(NF,XO,X) CALL MXVCOP(NF,CR,GF) CALL MXVCOP(NC+1,CFO,CF) IREST = 1 LD = KD GO TO 60 END IF * * COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE * FUNCTION TOGETHER WITH THE VALUES AND THE GRADIENTS OF THE * APPROXIMATED FUNCTIONS * IF (KD.GT.LD) THEN LDS=LD CALL PF1F01(NF,X,GF,GF,FF,F,KD,LD,IEXT, STAT) LD=LDS CALL PC1F01(NF,NC,X,FC,CF,CL,CU,IC,GC,CG,CMAX,KD,LD) END IF * * PREPARATION OF VARIABLE METRIC UPDATE * CALL MXVCOP(NF,GF,G) CALL PYTRND(NF,N,X,XO,ICA,CG,CZ,G,GO,R,F,FO,P,PO,CMAX,CMAXO, + DMAX,KD,LD,ITERS) * * VARIABLE METRIC UPDATE * CALL PUDBG1(N,H,G,S,XO,GO,R,PO,STAT(5),KIT,ITERH,MET,MET1,MEC) C IF (MER.GT.0.AND.ITERH.GT.0) IREST=1 * * END OF THE ITERATION * GO TO 40 80 CONTINUE c 80 IF (IPRNT.GT.1 .OR. IPRNT.LT.0) c + WRITE (6,'(1X,''EXIT FROM PSQP :'')') c IF (IPRNT.NE.0) c + WRITE (6,'(1X,''STAT(5)='',I4,2X,''STAT(6)='',I4,2X,''STAT(7)='',I4,2X, c + ''F='',G12.6,2X,''C='',E7.1,2X,''G='',E7.1,2X,''ITERM='',I3)') c + STAT(5),STAT(6),STAT(7),F,CMAX,GMAX,ITERM c IF (IPRNT.LT.0) c + WRITE (6,'(1X,''X='',5(G14.7,1X):/(3X,5(G14.7,1X)))') c + (X(I),I=1,NF) RETURN END * SUBROUTINE PA0GS1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE: * NUMERICAL COMPUTATION OF THE GRADIENT OF THE MODEL FUNCTION. * * PARAMETERS : * II N NUMBER OF VARIABLES. * II KA INDEF OF THE APPROXIMATED FUNCTION. * RI X(N) VECTOR OF VARIABLES. * RO GA(N) GRADIENT OF THE APPROXIMATED FUNCTION. * RI FA VALUE OF THE APPROXIMATED FUNCTION. * RI ETA1 PRECISION OF THE COMPUTED VALUES. * IU NAV NUMBER OF APPROXIMATED FUNCTION EVALUATIONS. * SUBROUTINE PA0GS1(N,KA,X,GA,FA,ETA1,NAV) IMPLICIT NONE INTEGER N,KA,NAV DOUBLE PRECISION X(*),GA(*),FA,ETA1 DOUBLE PRECISION XSTEP,XTEMP,FTEMP,ETA INTEGER IVAR ETA=SQRT(ETA1) FTEMP=FA DO 4 IVAR=1,N C C STEP SELECTION C XSTEP=1.0D 0 XSTEP=ETA*MAX(ABS(X(IVAR)),XSTEP)*SIGN(1.0D 0,X(IVAR)) XTEMP=X(IVAR) X(IVAR)=X(IVAR)+XSTEP XSTEP=X(IVAR)-XTEMP NAV=NAV+1 CALL FUN(N,KA,X,FA) C C NUMERICAL DIFFERENTIATION C GA(IVAR)=(FA-FTEMP)/XSTEP X(IVAR)=XTEMP 4 CONTINUE FA=FTEMP RETURN END * SUBROUTINE PA1SQ1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION * WHICH IS DEFINED AS A SUM OF SQUARES. * * PARAMETERS: * II N NUMBER OF VARIABLES. * RI X(N) VECTOR OF VARIABLES. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RO AF(N) VALUES OF THE APPROXIMATED FUNCTIONS. * RI GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. * RI AG(N*N) RECTANGULAR MATRIX WHICH IS USED FOR THE DIRECTION * VECTOR DETERMINATION. * RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RI ETA1 PRECISION OF THE COMPUTES FUNCTION VALUES. * II KDA DEGREE OF COMPUTED DERIVATIVES. * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * IU NFV NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED. * IU NFG NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED. * * SUBPROGRAMS USED : * S MXVCOP COPYING OF A VECTOR. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE PA1SQ1(N,X,F,AF,GA,AG,G,ETA1,KDA,KD,LD,NFV,NFG) IMPLICIT NONE INTEGER N,KDA,KD,LD,NFV,NFG DOUBLE PRECISION X(*),F,AF(*),GA(*),AG(*),G(*),ETA1 INTEGER KA,NAV DOUBLE PRECISION FA IF (KD.LE.LD) RETURN IF (KD.GE.0.AND.LD.LT.0) THEN F=0.0D0 NFV=NFV+1 ENDIF IF (KD.GE.1.AND.LD.LT.1) THEN CALL MXVSET(N,0.0D0,G) IF (KDA.GT.0) NFG=NFG+1 ENDIF NAV=0 DO 30 KA=1,N IF (KD.LT.0) GO TO 30 IF (LD.GE.0) THEN FA = AF(KA) GO TO 10 ELSE CALL FUN(N,KA,X,FA) AF(KA) = FA ENDIF IF (LD.GE.0) GO TO 10 F=F+FA*FA 10 IF (KD.LT.1) GO TO 30 IF (KDA.GT.0) THEN CALL DFUN(N,KA,X,GA) ELSE CALL PA0GS1(N,KA,X,GA,FA,ETA1,NAV) ENDIF CALL MXVDIR(N,FA,GA,G,G) CALL MXVCOP(N,GA,AG((KA-1)*N+1)) 30 CONTINUE NFV=NFV+NAV/N IF (KD.GE.0.AND.LD.LT.0) F=0.5D0*F LD=KD RETURN END * SUBROUTINE PC1F01 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF THE VALUE AND THE GRADIENT OF THE CONSTRAINT FUNCTION. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * RI X(NF) VECTOR OF VARIABLES. * RI FC VALUE OF THE SELECTED CONSTRAINT FUNCTION. * RI CF(NC) VECTOR CONTAINING VALUES OF CONSTRAINT FUNCTIONS. * RI CL(NC) VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS. * RI CU(NC) VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS. * II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. * RI GC(NF) GRADIENT OF THE SELECTED CONSTRAINT FUNCTION. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE GRADIENTS OF CONSTRAINT * FUNCTIONS. * RO CMAX MAXIMUM CONSTRAINT VIOLATION. * II KD DEGREE OF REQUIRED DERVATIVES. * II LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * * SUBPROGRAMS USED : * S MXVCOP COPYING OF A VECTOR. * SUBROUTINE PC1F01(NF,NC,X,FC,CF,CL,CU,IC,GC,CG,CMAX,KD,LD) IMPLICIT NONE DOUBLE PRECISION FC,CMAX INTEGER KD,LD,NC,NF DOUBLE PRECISION CF(*),CG(*),CL(*),CU(*),GC(*),X(*) INTEGER IC(*) DOUBLE PRECISION POM,TEMP INTEGER KC IF (KD.LE.LD) RETURN IF (LD.LT.0) CMAX = 0.0D0 DO 20 KC = 1,NC IF (KD.LT.0) GO TO 20 IF (LD.GE.0) THEN FC = CF(KC) GO TO 10 ELSE CALL CON(NF,KC,X,FC) CF(KC) = FC END IF IF (IC(KC).GT.0) THEN POM = 0.0D0 TEMP = CF(KC) IF (IC(KC).EQ.1.OR.IC(KC).GE.3) POM = MIN(POM, + TEMP - CL(KC)) IF (IC(KC).EQ.2.OR.IC(KC).GE.3) POM = MIN(POM, + CU(KC) - TEMP) IF (POM.LT.0.0D0) THEN CMAX = MAX(CMAX,-POM) END IF END IF 10 IF (KD.LT.1) GO TO 20 IF (LD.GE.1) THEN CALL MXVCOP(NF,CG((KC - 1) * NF + 1),GC) GO TO 20 ELSE CALL DCON(NF,KC,X,GC) CALL MXVCOP(NF,GC,CG((KC - 1) * NF + 1)) END IF 20 CONTINUE LD = KD RETURN END * SUBROUTINE PF1F01 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION. * * PARAMETERS: * II NF NUMBER OF VARIABLES. * RI X(NF) VECTOR OF VARIABLES. * RI GF(NF) GRADIENT OF THE MODEL FUNCTION. * RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. * RI FF VALUE OF THE MODEL FUNCTION. * RO F VALUE OF THE OBJECTIVE FUNCTION. * II KD DEGREE OF REQUIRED DERIVATIVES. * II LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II IEXT TYPE OF EXTREMUM. IEXT=0-MINIMUM. IEXT=1-MAXIMUM. * * SUBPROGRAMS USED : * S MXVCOP COPYING OF A VECTOR. * S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. * SUBROUTINE PF1F01(NF,X,GF,G,FF,F,KD,LD,IEXT, STAT) IMPLICIT NONE DOUBLE PRECISION F,FF INTEGER IEXT,KD,LD,NF DOUBLE PRECISION GF(*),G(*),X(*) c INTEGER NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES INTEGER STAT(8) c COMMON /STAT/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH IF (KD.LE.LD) RETURN IF (LD.GE.0) GO TO 10 STAT(6) = STAT(6) + 1 CALL OBJ(NF,X,FF) IF (IEXT.LE.0) THEN F = FF ELSE F = -FF END IF 10 IF (KD.LT.1) GO TO 20 IF (LD.GE.1) GO TO 20 STAT(7) = STAT(7) + 1 CALL DOBJ(NF,X,GF) IF (IEXT.GT.0) THEN CALL MXVNEG(NF,GF,G) END IF 20 LD = KD 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) IMPLICIT NONE INTEGER NF,N,ICA(*),INEW,NADD,IER DOUBLE PRECISION CG(*),CR(*),CZ(*),S(*),EPS7,GMAX,UMAX 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) ENDIF DO 1 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 1 CONTINUE ENDIF 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) IMPLICIT NONE INTEGER NF,N,ICA(*),IDECF,INEW,NADD,IER DOUBLE PRECISION CG(*),CR(*),CZ(*),H(*),S(*),EPS7,GMAX,UMAX DOUBLE PRECISION CK,CL INTEGER I,J,K,L,N1 IF (IDECF.NE.0.AND.IDECF.NE.9) THEN IER=-2 RETURN ENDIF 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) ENDIF DO 1 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 1 CONTINUE IF (IDECF.EQ.9) THEN L=N*(N+1)/2 IF (H(L+N1).NE.0.0D 0) THEN CL=1.0D 0/H(L+N1) K=0 DO 3 I=1,N CK=CL*H(L+I) DO 2 J=1,I K=K+1 H(K)=H(K)-CK*H(L+J) 2 CONTINUE 3 CONTINUE ENDIF ENDIF ENDIF 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) IMPLICIT NONE INTEGER NF,N,ICA(*),INEW,NADD,IER DOUBLE PRECISION CG(*),CR(*),S(*),EPS7,GMAX,UMAX DOUBLE PRECISION MXVDOT INTEGER NCA,NCR,I,J,K,L 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 1 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) ENDIF 1 CONTINUE ELSE K=-INEW GMAX=1.0D 0 DO 2 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.0D 0 ENDIF 2 CONTINUE ENDIF 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)) ENDIF 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 ENDIF 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, STAT) IMPLICIT NONE INTEGER NF,NC,IX(*),IC(*),ICA(*),MFP,KBF,KBC,N,ITERL DOUBLE PRECISION X(*),XO(*),XL(*),XU(*),CF(*),CFD(*),CL(*), & CU(*),CG(*),CR(*),CZ(*),G(*),GO(*),S(*),ETA9,EPS7,EPS9, & UMAX,GMAX DOUBLE PRECISION POM,CON,DMAX INTEGER NCA,NCR,NCZ,IPOM,I,K,IOLD,INEW,IER,KREM,KC,NRED c INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH INTEGER STAT(8) c COMMON /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH CON=ETA9 * * INITIATION * CALL MXVCOP(NF,X,XO) IPOM=0 NRED=0 KREM=0 ITERL=1 DMAX=0.0D 0 IF (MFP.EQ.3) GO TO 5 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 6 ENDIF N=0 NCA=0 NCZ=0 DO 1 I=1,NF IF (KBF.GT.0.AND.IX(I).LT.0) THEN NCA=NCA+1 ICA(NCA)=-I ELSE N=N+1 CALL MXVSET(NF,0.0D 0,CZ(NCZ+1)) CZ(NCZ+I)=1.0D 0 NCZ=NCZ+NF ENDIF 1 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) C IF (NF.GT.N) CALL PLSETC(NF,NC,X,XO,CF,IC,CG,S) DO 2 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,STAT(4),IER) CALL MXVIND(IC,KC,IER) IF (IC(KC).LT.-10) IPOM=1 ENDIF 2 CONTINUE ENDIF 3 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 ENDIF IF (IPOM.EQ.0.AND.ITERL.EQ.0) THEN * * FEASIBILITY ACHIEVED * ITERL=1 CALL MXVCOP(NF,G,GO) IF (MFP.EQ.1) GO TO 6 ENDIF * * LAGRANGE MULTIPLIERS AND REDUCED GRADIENT DETERMINATION * 5 CALL PLTRBG(NF,N,NC,IX,IC,ICA,CG,CR,CZ,GO,S,EPS7,GMAX,UMAX,IOLD) INEW=0 IF (GMAX.EQ.0.0D 0) THEN * * OPTIMUM ON A LINEAR MANIFOLD OBTAINED * IF (IOLD.EQ.0) THEN IF (IPOM.EQ.0) THEN * * OPTIMAL SOLUTION ACHIEVED * ITERL= 2 GO TO 6 ELSE IPOM=0 DO 22 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 ENDIF 22 CONTINUE IF (IPOM.EQ.0) THEN * * OPTIMAL SOLUTION ACHIEVED * CALL MXVCOP(NF,GO,G) ITERL= 2 GO TO 6 ELSE * * FEASIBLE SOLUTION DOES NOT EXIST * CALL MXVCOP(NF,GO,G) ITERL=-1 GO TO 6 ENDIF ENDIF ELSE * * CONSTRAINT DELETION * CALL PLRMB0(NF,N,ICA,CG,CR,CZ,GO,S,IOLD,KREM,STAT(3),IER) KC=ICA(NF-N+1) IF (KC.GT.0) THEN IC(KC)=-IC(KC) ELSE K=-KC IX(K)=-IX(K) ENDIF DMAX=0.0D 0 GO TO 5 ENDIF 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 ENDIF GO TO 6 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,STAT(4),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,STAT(4),IER) CALL MXVIND(IX,I,IER) ENDIF DMAX=POM IF (DMAX.GT.0.0D 0) NRED=NRED+1 GO TO 3 ENDIF ENDIF 6 CONTINUE 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) IMPLICIT NONE INTEGER NF,N,ICA(*),IOLD,KREM,NREM,IER DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*) DOUBLE PRECISION MXVDOT INTEGER NCA,NCR,NCZ,I,J,KC 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.0D 0,CZ(NCZ+1)) CZ(NCZ+NCA)=1.0D 0 CALL MXDPRB(NCA,CR,CZ(NCZ+1),-1) CALL MXVCOP(NCA,CZ(NCZ+1),CR(NCR+1)) N=N+1 CALL MXVSET(NF,0.0D 0,CZ(NCZ+1)) DO 4 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) ENDIF 4 CONTINUE GN(N)=MXVDOT(NF,CZ(NCZ+1),G) NREM=NREM+1 IER=0 RETURN END * SUBROUTINE PLQDB1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * IO N DIMENSION OF THE MANIFOLD DEFINED BY ACTIVE CONSTRAINTS. * II NC NUMBER OF LINEAR CONSTRAINTS. * RU X(NF) VECTOR OF VARIABLES. * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. * II IXA(NF) VECTOR CONTAINING INFORMATION ON TRUST REGION ACTIVITY. * RI XN(NF) VECTOR OF SCALING FACTORS. * 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 FUNCTIONS. * RO CFD(NC) 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. * RO G(NF) GRADIENT OF THE LAGRANGIAN FUNCTION. * RO GO(NF) SAVED GRADIENT OF THE OBJECTIVE FUNCTION. * RU H(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OR INVERSION OF THE * HESSIAN MATRIX APPROXIMATION. * RI S(NF) DIRECTION VECTOR. * 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. * RU XDEL TRUST REGION BOUND. * RO UMAX MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER. * RO GMAX MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE. * II MFP TYPE OF FEASIBLE POINT. MFP=1-ARBITRARY FEASIBLE POINT. * MFP=2-OPTIMUM FEASIBLE POINT. MFP=3-REPEATED SOLUTION. * * COMMON DATA : * 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 NORMF SCALING SPECIFICATION. NORMF=0-NO SCALING PERFORMED. * NORMF=1-SCALING FACTORS ARE DETERMINED AUTOMATICALLY. * NORMF=2-SCALING FACTORS ARE SUPPLIED BY USER. * IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. * IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION. * IDECF=10-DIAGONAL MATRIX. * IU NDECF NUMBER OF DECOMPOSITIONS. * 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 PLMINS DETERMINATION OF THE NEW ACTIVE SIMPLE BOUND. * S PLMINL DETERMINATION OF THE NEW ACTIVE LINEAR CONSTRAINT. * S PLMINT DETERMINATION OF THE NEW ACTIVE TRUST REGION BOUND. * S PLADR1 ADDITION OF A NEW ACTIVE CONSTRAINT. * S PLRMR0 CONSTRAIN DELETION. * S PLSOB1 TRANSFORMATION OF THE LOCAL SOLUTION TO THE SOLUTION * OF THE ORIGINAL QP PROBLEM. * S MXDPGF GILL-MURRAY DECOMPOSITION OF A DENSE SYMMETRIC MATRIX. * S MXDPGB BACK SUBSTITUTION AFTER GILL-MURRAY DECOMPOSITION. * S MXDPRB BACK SUBSTITUTION. * S MXDSMM MATRIX VECTOR PRODUCT. * S MXVCOP COPYING OF A VECTOR. * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * S MXVINA ABSOLUTE VALUES OF ELEMENTS OF AN INTEGER VECTOR. * S MXVINV CHANGE OF AN INTEGER VECTOR AFTER CONSTRAINT ADDITION. * S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. * SUBROUTINE PLQDB1(NF,NC,X,IX,XL,XU,CF,CFD,IC,ICA,CL,CU,CG,CR,CZ, & G,GO,H,S,MFP,KBF,KBC,IDECF,ETA2,ETA9,EPS7,EPS9,UMAX,GMAX,N,ITERQ, & STAT) IMPLICIT NONE INTEGER NF,NC,IX(*),IC(*),ICA(*),MFP,KBF,KBC,IDECF,N,ITERQ DOUBLE PRECISION X(*),XL(*),XU(*),CF(*),CFD(*),CL(*),CU(*),CG(*), & CR(*),CZ(*),G(*),GO(*),H(*),S(*),ETA2,ETA9,EPS7,EPS9,UMAX,GMAX DOUBLE PRECISION CON,TEMP,STEP,STEP1,STEP2,DMAX,PAR,SNORM INTEGER NCA,NCR,I,J,K,IOLD,JOLD,INEW,JNEW,KNEW,INF,IER,KREM,KC, & NRED c INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH INTEGER STAT(8) c COMMON /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH CON=ETA9 IF (IDECF.LT.0) IDECF=1 IF (IDECF.EQ.0) THEN C C GILL-MURRAY DECOMPOSITION C TEMP=ETA2 CALL MXDPGF(NF,H,INF,TEMP,STEP) STAT(2)=STAT(2)+1 IDECF=1 ENDIF IF (IDECF.GE.2.AND.IDECF.LE.8) THEN ITERQ=-10 RETURN ENDIF C C INITIATION C NRED=0 JOLD=0 JNEW=0 ITERQ=0 DMAX=0.0D0 IF (MFP.EQ.3) GO TO 1 N=NF NCA=0 NCR=0 IF (KBF.GT.0) CALL MXVINA(NF,IX) IF (KBC.GT.0) CALL MXVINA(NC,IC) C C DIRECTION DETERMINATION C 1 CALL MXVNEG(NF,GO,S) DO 2 J=1,NCA KC=ICA(J) IF (KC.GT.0) THEN CALL MXVDIR(NF,CZ(J),CG((KC-1)*NF+1),S,S) ELSE K=-KC S(K)=S(K)+CZ(J) ENDIF 2 CONTINUE CALL MXVCOP(NF,S,G) IF (IDECF.EQ.1) THEN CALL MXDPGB(NF,H,S,0) ELSE CALL MXDSMM(NF,H,G,S) ENDIF IF (ITERQ.EQ.3) GO TO 7 C C CHECK OF FEASIBILITY C INEW=0 PAR=0.0D0 CALL PLMINN(NF,NC,CF,CFD,IC,CL,CU,CG,S,EPS9,PAR,KBC,INEW,KNEW) CALL PLMINS(NF,IX,X,XL,XU,S,KBF,INEW,KNEW,EPS9,PAR) IF (INEW.EQ.0) THEN C C SOLUTION ACHIEVED C CALL MXVNEG(NF,G,G) ITERQ=2 GO TO 7 ELSE SNORM=0.0D0 ENDIF 4 IER=0 C C STEPSIZE DETERMINATION C CALL PLADR1(NF,N,ICA,CG,CR,H,S,G,EPS7,GMAX,UMAX,IDECF,INEW, & STAT(4),IER,1) CALL MXDPRB(NCA,CR,G,-1) IF (KNEW.LT.0) CALL MXVNEG(NCA,G,G) C C PRIMAL STEPSIZE C IF (IER.NE.0) THEN STEP1=CON ELSE STEP1=-PAR/UMAX ENDIF C C DUAL STEPSIZE C IOLD=0 STEP2=CON DO 5 J=1,NCA KC=ICA(J) IF (KC.GE.0) THEN K=IC(KC) ELSE I=-KC K=IX(I) ENDIF 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 TEMP=CZ(J)/G(J) IF (STEP2.GT.TEMP) THEN IOLD=J STEP2=TEMP ENDIF ENDIF 5 CONTINUE C C FINAL STEPSIZE C STEP=MIN(STEP1,STEP2) IF (STEP.GE.CON) THEN C C FEASIBLE SOLUTION DOES NOT EXIST C ITERQ=-1 GO TO 7 ENDIF C C NEW LAGRANGE MULTIPLIERS C DMAX=STEP CALL MXVDIR(NCA,-STEP,G,CZ,CZ) SNORM=SNORM+SIGN(1,KNEW)*STEP PAR=PAR-(STEP/STEP1)*PAR IF (STEP.EQ.STEP1) THEN IF (N.LE.0) THEN C C IMPOSSIBLE SITUATION C ITERQ=-5 GO TO 7 ENDIF C C CONSTRAINT ADDITION C IF (IER.EQ.0) THEN N=N-1 NCA=NCA+1 NCR=NCR+NCA CZ(NCA)=SNORM ENDIF 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 ENDIF NRED=NRED+1 STAT(4)=STAT(4)+1 JNEW=INEW JOLD=0 GO TO 1 ELSE C C CONSTRAINT DELETION C DO 6 J=IOLD,NCA-1 CZ(J)=CZ(J+1) 6 CONTINUE CALL PLRMF0(NF,NC,IX,IC,ICA,CR,IC,G,N,IOLD,KREM,IER, STAT) NCR=NCR-NCA NCA=NCA-1 JOLD=IOLD JNEW=0 IF (KBC.GT.0) CALL MXVINA(NC,IC) IF (KBF.GT.0) CALL MXVINA(NF,IX) DO 8 J=1,NCA KC=ICA(J) IF (KC.GT.0) THEN IC(KC)=-IC(KC) ELSE KC=-KC IX(KC)=-IX(KC) ENDIF 8 CONTINUE GO TO 4 ENDIF 7 CONTINUE RETURN END * SUBROUTINE PLADR1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * TRIANGULAR DECOMPOSITION OF KERNEL OF THE GENERAL 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. * RI CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RU H(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OR INVERSION OF THE * HESSIAN MATRIX APPROXIMATION. * RA S(NF) AUXILIARY VECTOR. * RO G(NF) VECTOR USED IN THE DUAL RANGE SPACE QUADRATIC PROGRAMMING * METHOD. * 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. * 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. * II JOB SPECIFICATION OF COMPUTATION. OUTPUT VECTOR G IS NOT OR IS * COMPUTED IN CASE WHEN N.LE.0 IF JOB=0 OR JOB=1 RESPECTIVELY. * * SUBPROGRAMS USED : * S MXDPGB BACK SUBSTITUTION. * S MXDPRB BACK SUBSTITUTION. * 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 PLADR1(NF,N,ICA,CG,CR,H,S,G,EPS7,GMAX,UMAX,IDECF, & INEW,NADD,IER,JOB) IMPLICIT NONE INTEGER NF,N,ICA(*),IDECF,INEW,NADD,IER,JOB DOUBLE PRECISION CG(*),CR(*),H(*),S(*),G(*),EPS7,GMAX,UMAX DOUBLE PRECISION MXVDOT INTEGER NCA,NCR,JCG,J,K,L IER=0 IF (JOB.EQ.0.AND.N.LE.0) IER=2 IF (INEW.EQ.0) IER=3 IF (IDECF.NE.1.AND.IDECF.NE.9) IER=-2 IF (IER.NE.0) RETURN NCA=NF-N NCR=NCA*(NCA+1)/2 IF (INEW.GT.0) THEN JCG=(INEW-1)*NF+1 IF (IDECF.EQ.1) THEN CALL MXVCOP(NF,CG(JCG),S) CALL MXDPGB(NF,H,S,0) ELSE CALL MXDSMM(NF,H,CG(JCG),S) ENDIF GMAX=MXVDOT(NF,CG(JCG),S) ELSE K=-INEW IF (IDECF.EQ.1) THEN CALL MXVSET(NF,0.0D0,S) S(K)=1.0D 0 CALL MXDPGB(NF,H,S,0) ELSE CALL MXDSMV(NF,H,S,K) ENDIF GMAX=S(K) ENDIF DO 1 J=1,NCA L=ICA(J) IF (L.GT.0) THEN G(J)=MXVDOT(NF,CG((L-1)*NF+1),S) ELSE L=-L G(J)=S(L) ENDIF 1 CONTINUE IF (N.EQ.0) THEN CALL MXDPRB(NCA,CR,G,1) UMAX=0.0D0 IER=2 RETURN ELSE IF (NCA.EQ.0) THEN UMAX=GMAX ELSE CALL MXDPRB(NCA,CR,G,1) UMAX=GMAX-MXVDOT(NCA,G,G) CALL MXVCOP(NCA,G,CR(NCR+1)) ENDIF IF (UMAX.LE.EPS7*GMAX) THEN IER=1 RETURN ELSE NCA=NCA+1 NCR=NCR+NCA ICA(NCA)=INEW CR(NCR)=SQRT(UMAX) IF (JOB.EQ.0) THEN N=N-1 NADD=NADD+1 ENDIF ENDIF 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) IMPLICIT NONE INTEGER NC,IC(*),KBC DOUBLE PRECISION CF(*),CFD(*),STEP INTEGER KC IF (KBC.GT.0) THEN DO 1 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) ENDIF 1 CONTINUE ENDIF 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) IMPLICIT NONE INTEGER NF,IX(*),KBF DOUBLE PRECISION X(*),S(*),STEP INTEGER I DO 1 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) ENDIF 1 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) IMPLICIT NONE INTEGER NF,IX(*),KBF,INEW,IND DOUBLE PRECISION X(*),XL(*),XU(*),EPS9 INTEGER I IND=0 IF (KBF.GT.0) THEN DO 1 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) ENDIF 1 CONTINUE ENDIF 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) IMPLICIT NONE INTEGER NF,NC,IC(*),KBC,KREM,INEW DOUBLE PRECISION CF(*),CFD(*),CL(*),CU(*),CG(*),S(*),STEP DOUBLE PRECISION TEMP,MXVDOT INTEGER JCG,KC IF (KBC.GT.0) THEN JCG=1 DO 1 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.0D 0) 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 ENDIF ENDIF ELSE IF (TEMP.GT.0.0D 0) 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 ENDIF ENDIF ENDIF ELSE IF (IC(KC).LT.-10) THEN TEMP=MXVDOT(NF,CG(JCG),S) CFD(KC)=TEMP IF (TEMP.GT.0.0D 0) 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 ENDIF ENDIF ELSE IF (TEMP.LT.0.0D 0) 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 ENDIF ENDIF ENDIF ENDIF JCG=JCG+NF 1 CONTINUE ENDIF 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) IMPLICIT NONE INTEGER NF,IX(*),KBF,KREM,INEW DOUBLE PRECISION X(*),XL(*),XU(*),S(*),STEP DOUBLE PRECISION TEMP INTEGER I IF (KBF.GT.0) THEN DO 1 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.0D 0) 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 ENDIF ENDIF ELSE IF (S(I).GT.0.0D 0) 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 ENDIF ENDIF ENDIF ENDIF 1 CONTINUE ENDIF 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) IMPLICIT NONE INTEGER KC,IC(*),INEW DOUBLE PRECISION CF(*),CL(*),CU(*),EPS9 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.0D 0) 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 ENDIF ELSE IF (IC(KC).EQ.2) THEN TEMP=EPS9*MAX(ABS(CU(KC)),1.0D 0) 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 ENDIF ELSE IF (IC(KC).EQ.3.OR.IC(KC).EQ.4) THEN TEMP=EPS9*MAX(ABS(CL(KC)),1.0D 0) IF (CF(KC).GT.CL(KC)+TEMP) THEN TEMP=EPS9*MAX(ABS(CU(KC)),1.0D 0) 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 ENDIF ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN IC(KC)=13 INEW=KC ELSE IC(KC)=-13 ENDIF ELSE IF (IC(KC).EQ.5.OR.IC(KC).EQ.6) THEN TEMP=EPS9*MAX(ABS(CL(KC)),1.0D 0) IF (CF(KC).GT.CL(KC)+TEMP) THEN TEMP=EPS9*MAX(ABS(CU(KC)),1.0D 0) 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 ENDIF ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN IC(KC)=15 INEW=KC ELSE IC(KC)=-15 ENDIF ENDIF RETURN END * SUBROUTINE PLMINN ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF THE NEW ACTIVE LINEAR CONSTRAINT. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * RI CF(NC) VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS. * RO CFD(NC) 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. * RI EPS9 TOLERANCE FOR ACTIVE CONSTRAINTS. * RA PAR AUXILIARY VARIABLE. * 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. * * SUBPROGRAMS USED : * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * SUBROUTINE PLMINN(NF,NC,CF,CFD,IC,CL,CU,CG,S,EPS9,PAR,KBC,INEW, & KNEW) IMPLICIT NONE INTEGER NF,NC,IC(*),KBC,INEW,KNEW DOUBLE PRECISION CF(*),CFD(*),CL(*),CU(*),CG(*),S(*),EPS9,PAR DOUBLE PRECISION TEMP,POM,MXVDOT INTEGER JCG,KC IF (KBC.GT.0) THEN JCG=1 DO 1 KC=1,NC IF (IC(KC).GT.0) THEN TEMP=MXVDOT(NF,CG(JCG),S) CFD(KC)=TEMP TEMP=CF(KC)+TEMP 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.0D 0))) THEN INEW=KC KNEW= 1 PAR=POM ENDIF ENDIF 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.0D 0))) THEN INEW=KC KNEW=-1 PAR=POM ENDIF ENDIF ENDIF JCG=JCG+NF 1 CONTINUE ENDIF 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) IMPLICIT NONE 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 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) IMPLICIT NONE INTEGER IX(*),I,INEW DOUBLE PRECISION X(*),XL(*),XU(*),EPS9 DOUBLE PRECISION TEMP TEMP=1.0D 0 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 ENDIF 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 ENDIF 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 ENDIF IF (X(I).GE.XU(I)-EPS9*MAX(ABS(XU(I)),TEMP)) THEN IX(I)=14 INEW=-I ENDIF ENDIF RETURN END * SUBROUTINE PLREDL ALL SYSTEMS 98/12/01 C PORTABILITY : ALL SYSTEMS C 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * TRANSFORMATION OF THE INCOMPATIBLE QUADRATIC PROGRAMMING SUBPROBLEM. * * PARAMETERS : * II NC NUMBER OF CURRENT LINEAR CONSTRAINTS. * 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. * II KBC SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR * CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO * SIDED LINEAR CONSTRAINTS. * SUBROUTINE PLREDL(NC,CF,IC,CL,CU,KBC) IMPLICIT NONE INTEGER NC,IC(NC),KBC,K DOUBLE PRECISION CF(*),CL(*),CU(*) DOUBLE PRECISION TEMP INTEGER KC IF (KBC.GT.0) THEN DO 1 KC=1,NC K=IC(KC) IF (ABS(K).EQ.1.OR.ABS(K).EQ.3.OR.ABS(K).EQ.4) THEN TEMP=(CF(KC)-CL(KC)) IF (TEMP.LT.0) CF(KC)=CL(KC)+0.1D 0*TEMP ENDIF IF (ABS(K).EQ.2.OR.ABS(K).EQ.3.OR.ABS(K).EQ.4) THEN TEMP=(CF(KC)-CU(KC)) IF (TEMP.GT.0) CF(KC)=CU(KC)+0.1D 0*TEMP ENDIF IF (ABS(K).EQ.5.OR.ABS(K).EQ.6) THEN TEMP=(CF(KC)-CL(KC)) CF(KC)=CL(KC)+0.1D 0*TEMP ENDIF 1 CONTINUE ENDIF 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, STAT) IMPLICIT NONE INTEGER IER,IOLD,KREM,N,NC,NF DOUBLE PRECISION AR(*),S(*) INTEGER IA(*),IAA(*),IC(*),IX(*) c INTEGER NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES INTEGER L INTEGER STAT(8) c COMMON /STAT/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH CALL PLRMR0(NF,IAA,AR,S,N,IOLD,KREM,IER) N = N + 1 STAT(3) = STAT(3) + 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) IMPLICIT NONE 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 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) IMPLICIT NONE INTEGER NF,NC,IC(*) DOUBLE PRECISION X(*),XO(*),CF(*),CG(*),S(*) DOUBLE PRECISION MXVDOT INTEGER JCG,KC CALL MXVDIF(NF,X,XO,S) JCG=0 DO 1 KC=1,NC IF (IC(KC).NE.0) CF(KC)=CF(KC)+MXVDOT(NF,S,CG(JCG+1)) JCG=JCG+NF 1 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) IMPLICIT NONE INTEGER NF,NC,IC(*),INEW DOUBLE PRECISION CG(*),G(*) INTEGER KC CALL MXVSET(NF,0.0D 0,G) INEW=0 DO 4 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.0D 0,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.0D 0,CG((KC-1)*NF+1),G,G) INEW=1 ENDIF 4 CONTINUE 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) IMPLICIT NONE INTEGER NF,N,NC,IX(*),IA(*),IAA(*),IC(*),IOLD DOUBLE PRECISION AZ(*),EPS7,UMAX DOUBLE PRECISION TEMP INTEGER NAA,J,K,L IOLD=0 UMAX=0.0D 0 NAA=NF-N DO 2 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) ENDIF IF (K.LE.-5) THEN ELSE IF ((K.EQ.-1.OR.K.EQ.-3).AND.UMAX+TEMP.GE.0.0D 0) THEN ELSE IF ((K.EQ.-2.OR.K.EQ.-4).AND.UMAX-TEMP.GE.0.0D 0) THEN ELSE IOLD=J UMAX=ABS(TEMP) ENDIF 2 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) IMPLICIT NONE INTEGER NF,N,NC,IX(*),IC(*),ICA(*),IOLD DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),EPS7,GMAX,UMAX DOUBLE PRECISION MXVMAX INTEGER NCA,NCZ GMAX=0.0D 0 IF (N.GT.0) THEN CALL MXDRMM(NF,N,CZ,G,GN) GMAX=MXVMAX(N,GN) ENDIF 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.0D 0,GN) GMAX=0.0D 0 ELSE IOLD=0 UMAX=0.0D 0 ENDIF 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) IMPLICIT NONE INTEGER NF,N,NC,IAA(*) DOUBLE PRECISION AG(*),CG(*),G(*),GN(*) DOUBLE PRECISION MXVDOT INTEGER NAA,J,L NAA=NF-N DO 1 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) ENDIF 1 CONTINUE RETURN END * SUBROUTINE PNINT1 ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 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) IMPLICIT NONE DOUBLE PRECISION RL, RU, FL, FU, PL, PU, R INTEGER MODE,MTYP,MERR,NTYP DOUBLE PRECISION A,B,C,D,DIS,DEN DOUBLE PRECISION C1L,C1U,C2L,C2U,C3L PARAMETER (C1L=1.1D 0,C1U=1.0D 3,C2L=1.0D-2,C2U=0.9D 0, & C3L=0.1D 0) MERR=0 IF (MODE.LE.0) RETURN IF (PL.GE.0.0D 0) THEN MERR=2 RETURN ELSE IF (RU.LE.RL) THEN MERR=3 RETURN ENDIF DO 1 NTYP=MTYP,1,-1 IF (NTYP.EQ.1) THEN C C BISECTION C IF (MODE.EQ.1) THEN R=4.0D 0*RU RETURN ELSE R=0.5D 0*(RL+RU) RETURN ENDIF ELSE IF (NTYP.EQ.MTYP) THEN A = (FU-FL)/(PL*(RU-RL)) B = PU/PL ENDIF IF (NTYP.EQ.2) THEN C C QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH ONE DIRECTIONAL C DERIVATIVE C DEN = 2.0D 0*(1.0D 0-A) ELSE IF (NTYP.EQ.3) THEN C C QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH TWO DIRECTIONAL C DERIVATIVES C DEN = 1.0D 0 - B ELSE IF (NTYP.EQ.4) THEN C C CUBIC EXTRAPOLATION OR INTERPOLATION C C = B - 2.0D 0*A + 1.0D 0 D = B - 3.0D 0*A + 2.0D 0 DIS = D*D - 3.0D 0*C IF (DIS.LT.0.0D 0) GO TO 1 DEN = D + SQRT(DIS) ELSE IF (NTYP.EQ.5) THEN C C CONIC EXTRAPOLATION OR INTERPOLATION C DIS = A*A - B IF (DIS.LT.0.0D 0) GO TO 1 DEN = A + SQRT(DIS) IF (DEN.LE.0.0D 0) GO TO 1 DEN = 1.0D 0 - B*(1.0D 0/DEN)**3 ENDIF IF (MODE.EQ.1.AND.DEN.GT.0.0D 0.AND.DEN.LT.1.0D 0) THEN C C EXTRAPOLATION ACCEPTED C R = RL + (RU-RL)/DEN R = MAX(R,C1L*RU) R = MIN(R,C1U*RU) RETURN ELSE IF (MODE.EQ.2.AND.DEN.GT.1.0D 0) THEN C C INTERPOLATION ACCEPTED C R = RL + (RU-RL)/DEN IF (RL.EQ.0.0D 0) THEN R = MAX(R,RL+C2L*(RU-RL)) ELSE R = MAX(R,RL+C3L*(RU-RL)) ENDIF R = MIN(R,RL+C2U*(RU-RL)) RETURN ENDIF 1 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) IMPLICIT NONE 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 PNSTEP ALL SYSTEMS 89/12/01 * PORTABILITY : ALL SYSTEMS * 89/01/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF A SCALING FACTOR FOR THE BOUNDARY STEP. * * PARAMETERS : * RI DEL MAXIMUM STEPSIZE. * RI A INPUT PARAMETER. * RI B INPUT PARAMETER. * RI C INPUT PARAMETER. * RO ALF SCALING FACTOR FOR THE BOUNDARY STEP SUCH THAT * A**2+2*B*ALF+C*ALF**2=DEL**2. * SUBROUTINE PNSTEP(DEL,A,B,C,ALF) IMPLICIT NONE DOUBLE PRECISION DEL, A, B, C, ALF DOUBLE PRECISION DEN, DIS DOUBLE PRECISION ZERO PARAMETER (ZERO = 0.0D 0) ALF = ZERO DEN = (DEL+A) * (DEL-A) IF (DEN .LE. ZERO) RETURN DIS = B*B + C*DEN IF (B .GE. ZERO) THEN ALF = DEN / (SQRT(DIS) + B) ELSE ALF = (SQRT(DIS) - B) / C ENDIF RETURN END * SUBROUTINE PP0AF8 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF VALUE OF THE AUGMENTED LAGRANGIAN FUNCTION. * * PARAMETERS : * II NF NUMBER OF VARIABLES. * II N DIMENSION OF THE CONSTRAINT NULL SPACE. * II NC NUMBER OF CONSTRAINTS. * RI CF(NC+1) VECTOR CONTAINING VALUES OF THE CONSTRAINTS. * 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 CZ(NC) VECTOR OF LAGRANGE MULTIPLIERS. * RI RPF PENALTY COEFFICIENT. * RO FC VALUE OF THE PENALTY TERM. * RO F VALUE OF THE PENALTY FUNCTION. * SUBROUTINE PP0AF8(NF,N,NC,CF,IC,ICA,CL,CU,CZ,RPF,FC,F) IMPLICIT NONE INTEGER NF,N,NC,IC(*),ICA(*) DOUBLE PRECISION CF(*),CL(*),CU(*),CZ(*),RPF,FC,F DOUBLE PRECISION POM,TEMP INTEGER J,KC FC=0.0D0 DO 1 KC=1,NC IF (IC(KC).GT.0) THEN POM=0.0D0 TEMP=CF(KC) IF (IC(KC).EQ.1.OR.IC(KC).GE.3) POM=MIN(POM,TEMP-CL(KC)) IF (IC(KC).EQ.2.OR.IC(KC).GE.3) POM=MIN(POM,CU(KC)-TEMP) FC=FC+RPF*ABS(POM) ENDIF 1 CONTINUE DO 2 J=1,NF-N KC=ICA(J) IF (KC.GT.0) THEN POM=0.0D0 TEMP=CF(KC) IF (IC(KC).EQ.1.OR.IC(KC).EQ.3.OR.IC(KC).EQ.5) & POM=MIN(POM,TEMP-CL(KC)) IF (IC(KC).EQ.2.OR.IC(KC).EQ.4.OR.IC(KC).EQ.6) & POM=MAX(POM,TEMP-CU(KC)) FC=FC-CZ(J)*POM ENDIF 2 CONTINUE F=CF(NC+1)+FC RETURN END * SUBROUTINE PPSET2 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF THE NEW PENALTY PARAMETERS. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * II ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CZ(NF) VECTOR OF LAGRANGE MULTIPLIERS. * RI CP(NC) VECTOR CONTAINING PENALTY PARAMETERS. * SUBROUTINE PPSET2(NF,N,NC,ICA,CZ,CP) IMPLICIT NONE INTEGER NF,N,NC,ICA(*) DOUBLE PRECISION CZ(*),CP(*) DOUBLE PRECISION TEMP INTEGER J,L,KC DO 1 KC=1,NC CP(KC)=0.5D 0*CP(KC) 1 CONTINUE DO 2 J=1,NF-N L=ICA(J) IF (L.GT.0) THEN TEMP=ABS(CZ(J)) CP(L)=MAX(TEMP,CP(L)+0.5D 0*TEMP) ENDIF 2 CONTINUE RETURN END * SUBROUTINE PS0G01 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SIMPLE SEARCH WITH TRUST REGION UPDATE. * * PARAMETERS : * RO R 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 PP QUADRATIC PART OF THE PREDICTED FUNCTION VALUE. * RU XDEL TRUST REGION BOUND. * RO XDELO PREVIOUS TRUST REGION BOUND. * RI XMAX MAXIMUM STEPSIZE. * RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER. * RI SNORM EUCLIDEAN NORM OF THE DIRECTION VECTOR. * RI BET1 LOWER BOUND FOR STEPSIZE REDUCTION. * RI BET2 UPPER BOUND FOR STEPSIZE REDUCTION. * RI GAM1 LOWER BOUND FOR STEPSIZE EXPANSION. * RI GAM2 UPPER BOUND FOR STEPSIZE EXPANSION. * RI EPS4 FIRST TOLERANCE FOR RATIO DF/DFPRED. STEP BOUND IS * DECREASED IF DF/DFPREDEPS5. * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * IU IDIR INDICATOR FOR DIRECTION DETERMINATION. * IDIR=0-BASIC DETERMINATION. IDIR=1-DETERMINATION * AFTER STEPSIZE REDUCTION. IDIR=2-DETERMINATION AFTER * STEPSIZE EXPANSION. * IO ITERS TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-STEP * BOUND WAS DECREASED. ITERS=2-STEP BOUND WAS UNCHANGED. * ITERS=3-STEP BOUND WAS INCREASED. ITERS=6-FIRST STEPSIZE. * II 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 MAXST MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM * STEPSIZE WAS NOT OR WAS REACHED. * IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * II MRED MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * II KTERS TERMINATION SELECTION. KTERS=1-NORMAL TERMINATION. * KTERS=6-FIRST STEPSIZE. * II MES1 SWITCH FOR EXTRAPOLATION. MES1=1-CONSTANT INCREASING OF * THE INTERVAL. MES1=2-EXTRAPOLATION SPECIFIED BY THE PARAMETER * MES. MES1=3 SUPPRESSED EXTRAPOLATION. * II MES2 SWITCH FOR TERMINATION. MES2=1-NORMAL TERMINATION. * MES2=2-TERMINATION AFTER AT LEAST TWO STEPS (ASYMPTOTICALLY * PERFECT LINE SEARCH). * II MES3 SAFEGUARD AGAINST ROUNDING ERRORS. MES3=0-SAFEGUARD * SUPPRESSED. MES3=1-FIRST LEVEL OF SAFEGUARD. MES3=2-SECOND * LEVEL OF SAFEGUARD. * IU ISYS CONTROL PARAMETER. * * COMMON DATA : * * METHOD : * G.A.SCHULTZ, R.B.SCHNABEL, R.H.BYRD: A FAMILY OF TRUST-REGION-BASED * ALGORITHMS FOR UNCONSTRAINED MINIMIZATION WITH STRONG GLOBAL * CONVERGENCE PROPERTIES, SIAM J. NUMER.ANAL. 22 (1985) PP. 47-67. * SUBROUTINE 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) IMPLICIT NONE INTEGER KD,LD,IDIR,ITERS,ITERD,MAXST,NRED,MRED,KTERS, & MES1,MES2,MES3,ISYS DOUBLE PRECISION R,F,FO,PO,PP,XDEL,XDELO,XMAX,RMAX,SNORM,BET1, & BET2,GAM1,GAM2,EPS4,EPS5 DOUBLE PRECISION DF,DFPRED INTEGER NRED1,NRED2 SAVE NRED1,NRED2 IF (ISYS.EQ.1) GO TO 2 C GO TO (1,2) ISYS+1 1 CONTINUE IF (IDIR.EQ.0) THEN NRED1=0 NRED2=0 ENDIF IDIR=0 XDELO=XDEL C C COMPUTATION OF THE NEW FUNCTION VALUE C R=MIN(1.0D 0,RMAX) KD= 0 LD=-1 ISYS=1 RETURN 2 CONTINUE IF(KTERS.LT.0.OR.KTERS.GT.5) THEN ITERS=6 ELSE DF=FO-F DFPRED=-R*(PO+R*PP) IF (DF.LT.EPS4*DFPRED) THEN C C STEP IS TOO LARGE, IT HAS TO BE REDUCED C IF (MES1.EQ.1) THEN XDEL=BET2*SNORM ELSE IF (MES1.EQ.2) THEN XDEL=BET2*MIN(0.5D 0*XDEL,SNORM) ELSE XDEL=0.5D 0*PO*SNORM/(PO+DF) XDEL=MAX(XDEL,BET1*SNORM) XDEL=MIN(XDEL,BET2*SNORM) ENDIF ITERS=1 IF (MES3.LE.1) THEN NRED2=NRED2+1 ELSE IF (ITERD.GT.2) NRED2=NRED2+1 ENDIF ELSE IF (DF.LE.EPS5*DFPRED) THEN C C STEP IS SUITABLE C ITERS=2 ELSE C C STEP IS TOO SMALL, IT HAS TO BE ENLARGED C IF (MES2.EQ.2) THEN XDEL=MAX(XDEL,GAM1*SNORM) ELSE IF (ITERD.GT.2) THEN XDEL=GAM1*XDEL ENDIF ITERS=3 ENDIF XDEL=MIN(XDEL,XMAX,GAM2*SNORM) IF (FO.LE.F) THEN IF (NRED1.GE.MRED) THEN ITERS=-1 ELSE IDIR=1 ITERS=0 NRED1=NRED1+1 ENDIF ENDIF ENDIF MAXST=0 IF (XDEL.GE.XMAX) MAXST=1 IF (MES3.EQ.0) THEN NRED=NRED1 ELSE NRED=NRED2 ENDIF ISYS=0 RETURN END * SUBROUTINE PS0L02 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * EXTENDED LINE SEARCH WITHOUT DIRECTIONAL DERIVATIVES. * * PARAMETERS : * RO R VALUE OF THE STEPSIZE PARAMETER. * RO RO INITIAL VALUE OF THE STEPSIZE PARAMETER. * RO RP PREVIOUS VALUE OF THE STEPSIZE PARAMETER. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RI FO INITIAL VALUE OF THE OBJECTIVE FUNCTION. * RO FP PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RI PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE. * RO PP PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE. * RI FMIN LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. * RI FMAX UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. * RI RMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER * RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER * RI TOLS TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE * CHANGE OF THE FUNCTION VALUE). * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II NIT ACTUAL NUMBER OF ITERATIONS. * II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. * IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * II MRED MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * IO MAXST MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM * STEPSIZE WAS NOT OR WAS REACHED. * II IEST LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND * IS NOT OR IS GIVEN. * II INITS CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE * IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL * STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED * INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL * 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. * II KTERS TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH. * KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE. * KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE. * KTERS=6-FIRST STEPSIZE. * II MES 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. * IU ISYS CONTROL PARAMETER. * * SUBPROGRAM USED : * S PNINT3 EXTRAPOLATION OR INTERPOLATION WITHOUT DIRECTIONAL * DERIVATIVES. * * METHOD : * SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH EXTENDED TERMINATION * CRITERIA. * SUBROUTINE 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) IMPLICIT NONE INTEGER KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS, & MES,ISYS DOUBLE PRECISION R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,TOLS DOUBLE PRECISION RL,FL,RU,FU,RI,FI,RTEMP,TOL INTEGER MTYP,MERR,MODE,INIT1,MES1,MES2 LOGICAL L1,L2,L3,L4,L6,L7 PARAMETER(TOL=1.0D-4) SAVE MTYP,MODE,MES1,MES2 SAVE RL,FL,RU,FU,RI,FI IF (ISYS.EQ.1) GO TO 3 C GO TO (1,3) ISYS+1 1 CONTINUE MES1=2 MES2=2 ITERS=0 IF (PO.GE.0.0D 0) THEN R=0.0D 0 ITERS=-2 GO TO 4 ENDIF IF(RMAX.LE.0.0D 0) THEN ITERS= 0 GO TO 4 ENDIF C C INITIAL STEPSIZE SELECTION C IF (INITS.GT.0) THEN RTEMP=FMIN-F ELSE IF (IEST.EQ.0) THEN RTEMP=F-FP ELSE RTEMP=MAX(F-FP,1.0D 1*(FMIN-F)) ENDIF INIT1=ABS(INITS) RP=0.0D 0 FP=FO PP=PO IF (INIT1.EQ.0) THEN ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN R=1.0D 0 ELSE IF (INIT1.EQ.2) THEN R=MIN(1.0D 0,4.0D 0*RTEMP/PO) ELSE IF (INIT1.EQ.3) THEN R=MIN(1.0D 0, 2.0D 0*RTEMP/PO) ELSE IF (INIT1.EQ.4) THEN R=2.0D 0*RTEMP/PO ENDIF RTEMP=R R=MAX(R,RMIN) R=MIN(R,RMAX) MODE=0 RL=0.0D 0 FL=FO RU=0.0D 0 FU=FO RI=0.0D 0 FI=FO C C NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION) C 2 CALL PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR) IF (MERR.GT.0) THEN ITERS=-MERR GO TO 4 ELSE IF (MODE.EQ.1) THEN NRED=NRED-1 R=MIN(R,RMAX) ELSE IF (MODE.EQ.2) THEN NRED=NRED+1 ENDIF C C COMPUTATION OF THE NEW FUNCTION VALUE C KD= 0 LD=-1 ISYS=1 RETURN 3 CONTINUE IF (ITERS.NE.0) GO TO 4 IF (F.LE.FMIN) THEN ITERS=7 GO TO 4 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.0D 1 L4=F-FO.GE.(1.0D 0-TOLS)*R*PO.OR.MES2.EQ.2.AND.MODE.EQ.2 L6=RU-RL.LE.TOL*RU.AND.MODE.EQ.2 L7=MES2.LE.2.OR.MODE.NE.0 MAXST=0 IF (L2) MAXST=1 ENDIF C C TEST ON TERMINATION C IF (L1.AND..NOT.L3) THEN ITERS=0 GO TO 4 ELSE IF (L2.AND..NOT.F.GE.FU) THEN ITERS=7 GO TO 4 ELSE IF (L6) THEN ITERS=1 GO TO 4 ELSE IF (L3.AND.L7.AND.KTERS.EQ.5) THEN ITERS=5 GO TO 4 ELSE IF (L3.AND.L4.AND.L7.AND.(KTERS.EQ.2.OR.KTERS.EQ.3.OR. * KTERS.EQ.4)) THEN ITERS=2 GO TO 4 ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN ITERS=6 GO TO 4 ELSE IF(ABS(NRED).GE.MRED) THEN ITERS=-1 GO TO 4 ELSE RP=R FP=F MODE=MAX(MODE,1) MTYP=ABS(MES) IF(F.GE.FMAX) MTYP=1 ENDIF IF (MODE.EQ.1) THEN C C INTERVAL CHANGE AFTER EXTRAPOLATION C RL=RI FL=FI RI=RU FI=FU RU=R FU=F IF (F.GE.FI) THEN NRED=0 MODE=2 ELSE IF ( MES1 .EQ. 1) THEN MTYP=1 ENDIF C C INTERVAL CHANGE AFTER INTERPOLATION C ELSE IF (R.LE.RI) THEN IF (F.LE.FI) THEN RU=RI FU=FI RI=R FI=F ELSE RL=R FL=F ENDIF ELSE IF (F.LE.FI) THEN RL=RI FL=FI RI=R FI=F ELSE RU=R FU=F ENDIF ENDIF GO TO 2 4 ISYS=0 RETURN END * SUBROUTINE PS1L01 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * STANDARD LINE SEARCH WITH DIRECTIONAL DERIVATIVES. * * PARAMETERS : * RO R VALUE OF THE STEPSIZE PARAMETER. * RO RP PREVIOUS VALUE OF THE STEPSIZE PARAMETER. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RI FO INITIAL VALUE OF THE OBJECTIVE FUNCTION. * RO FP PREVIOUS 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. * RI FMIN LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. * RI FMAX UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. * RI RMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER * RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER * RI TOLS TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE * CHANGE OF THE FUNCTION VALUE). * RI TOLP TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE * CHANGE OF THE DIRECTIONAL DERIVATIVE). * RO PAR1 PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC * UPDATES. * RO PAR2 PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC * UPDATES. * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II NIT ACTUAL NUMBER OF ITERATIONS. * II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. * IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * II MRED MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * IO MAXST MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM * STEPSIZE WAS NOT OR WAS REACHED. * II IEST LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND * IS NOT OR IS GIVEN. * II INITS CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE * IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL * STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED * INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL * 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. * II KTERS TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH. * KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE. * KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE. * KTERS=6-FIRST STEPSIZE. * II MES 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. * IU ISYS CONTROL PARAMETER. * * SUBPROGRAM USED : * S PNINT1 EXTRAPOLATION OR INTERPOLATION WITH DIRECTIONAL * DERIVATIVES. * * METHOD : * SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH STANDARD TERMINATION * CRITERIA. * SUBROUTINE 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) IMPLICIT NONE INTEGER KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS, & MES,ISYS DOUBLE PRECISION R,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX, & TOLS,TOLP,PAR1,PAR2 DOUBLE PRECISION RL,FL,PL,RU,FU,PU,RTEMP INTEGER MTYP,MERR,MODE,INIT1,MES1,MES2,MES3 LOGICAL L1,L2,L3,L5,L7,M1,M2,M3 DOUBLE PRECISION CON,CON1 PARAMETER (CON=1.0D-2,CON1=1.0D-13) SAVE MTYP,MODE,MES1,MES2,MES3 SAVE RL,FL,PL,RU,FU,PU IF (ISYS.EQ.1) GO TO 3 C GO TO (1,3) ISYS+1 1 CONTINUE MES1=2 MES2=2 MES3=2 ITERS=0 IF (PO.GE.0.0D 0) THEN R=0.0D 0 ITERS=-2 GO TO 4 ENDIF IF(RMAX.LE.0.0D 0) THEN ITERS=0 GO TO 4 ENDIF C C INITIAL STEPSIZE SELECTION C IF (INITS.GT.0) THEN RTEMP=FMIN-F ELSE IF (IEST.EQ.0) THEN RTEMP=F-FP ELSE RTEMP=MAX(F-FP,1.0D 1*(FMIN-F)) ENDIF INIT1=ABS(INITS) RP=0.0D 0 FP=FO PP=PO IF (INIT1.EQ.0) THEN ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN R=1.0D 0 ELSE IF (INIT1.EQ.2) THEN R=MIN(1.0D 0,4.0D 0*RTEMP/PO) ELSE IF (INIT1.EQ.3) THEN R=MIN(1.0D 0, 2.0D 0*RTEMP/PO) ELSE IF (INIT1.EQ.4) THEN R=2.0D 0*RTEMP/PO ENDIF R=MAX(R,RMIN) R=MIN(R,RMAX) MODE=0 RU=0.0D 0 FU=FO PU=PO C C NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION) C 2 CALL PNINT1(RL,RU,FL,FU,PL,PU,R,MODE,MTYP,MERR) IF (MERR.GT.0) THEN ITERS=-MERR GO TO 4 ELSE IF (MODE.EQ.1) THEN NRED=NRED-1 R=MIN(R,RMAX) ELSE IF (MODE.EQ.2) THEN NRED=NRED+1 ENDIF C C COMPUTATION OF THE NEW FUNCTION VALUE AND THE NEW DIRECTIONAL C DERIVATIVE C KD= 1 LD=-1 ISYS=1 RETURN 3 CONTINUE IF (MODE.EQ.0) THEN PAR1=P/PO PAR2=F-FO ENDIF IF (ITERS.NE.0) GO TO 4 IF (F.LE.FMIN) THEN ITERS=7 GO TO 4 ELSE L1=R.LE.RMIN.AND.NIT.NE.KIT L2=R.GE.RMAX L3=F-FO.LE.TOLS*R*PO L5=P.GE.TOLP*PO.OR.MES2.EQ.2.AND.MODE.EQ.2 L7=MES2.LE.2.OR.MODE.NE.0 M1=.FALSE. M2=.FALSE. M3=L3 IF(MES3.GE.1) THEN M1=ABS(P).LE.CON*ABS(PO).AND.FO-F.GE.(CON1/CON)*ABS(FO) L3=L3.OR.M1 ENDIF IF(MES3.GE.2) THEN M2=ABS(P).LE.0.5D 0*ABS(PO).AND.ABS(FO-F).LE.2.0D 0*CON1*ABS(FO) L3=L3.OR.M2 ENDIF MAXST=0 IF (L2) MAXST=1 ENDIF C C TEST ON TERMINATION C IF (L1.AND..NOT.L3) THEN ITERS=0 GO TO 4 ELSE IF (L2.AND.L3.AND..NOT.L5) THEN ITERS=7 GO TO 4 ELSE IF (M3.AND.MES1.EQ.3) THEN ITERS=5 GO TO 4 ELSE IF (L3.AND.L5.AND.L7) THEN ITERS=4 GO TO 4 ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN ITERS=6 GO TO 4 ELSE IF (ABS(NRED).GE.MRED) THEN ITERS=-1 GO TO 4 ELSE RP=R FP=F PP=P MODE=MAX(MODE,1) MTYP=ABS(MES) IF(F.GE.FMAX) MTYP=1 ENDIF IF (MODE.EQ.1) THEN C C INTERVAL CHANGE AFTER EXTRAPOLATION C RL=RU FL=FU PL=PU RU=R FU=F PU=P IF (.NOT.L3) THEN NRED=0 MODE=2 ELSE IF ( MES1 .EQ. 1) THEN MTYP=1 ENDIF ELSE C C INTERVAL CHANGE AFTER INTERPOLATION C IF (.NOT.L3) THEN RU=R FU=F PU=P ELSE RL=R FL=F PL=P ENDIF ENDIF GO TO 2 4 ISYS=0 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 MET1 SELECTION OF SELF SCALING. MET1=1-SELF SCALING SUPPRESSED. * MET1=2 INITIAL SELF SCALING. * II MEC CORRECTION IF THE NEGATIVE CURVATURE OCCURS. * MEC=1-CORRECTION SUPPRESSED. MEC=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,MET1,MEC) IMPLICIT NONE DOUBLE PRECISION PO,R INTEGER ITERH,KIT,MET,MET1,MEC,N,NIT DOUBLE PRECISION G(*),GO(*),H(*),S(*),XO(*) DOUBLE PRECISION A,B,C,GAM,PAR,DEN,DIS LOGICAL L1,L3 DOUBLE PRECISION MXVDOT,MXDPGP L1 = MET1 .GE. 3 .OR. MET1 .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 (MEC.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) ENDIF ELSE IF (B.LE.1.0D-4*C) THEN ITERH = 2 RETURN ENDIF ENDIF IF (L1) THEN * * DETERMINATION OF THE PARAMETER GAM (SELF SCALING) * IF (MET.EQ.1) THEN PAR = C/B ELSE IF (A.LE.0.0D 0) THEN PAR = C/B ELSE PAR=SQRT(C/A) ENDIF GAM = PAR IF (MET1.GT.1) THEN IF (NIT.NE.KIT) THEN L3=GAM.LT.0.5D0.OR.GAM.GT.4.0D0 ENDIF ENDIF ENDIF IF (L3) THEN GAM = 1.0D0 PAR = GAM END IF IF (MET.EQ.1) THEN * * BFGS UPDATE * CALL MXDPGU(N,H,PAR/B,GO,XO) CALL MXDPGU(N,H,-1.0D0/C,S,XO) ELSE * * HOSHINO UPDATE * DEN=PAR*B+C DIS=0.5D0*B CALL MXVDIR(N,PAR,GO,S,S) CALL MXDPGU(N,H,PAR/DIS,GO,XO) CALL MXDPGU(N,H,-1.0D0/DEN,S,XO) ENDIF ITERH = 0 IF (GAM.EQ.1.0D0) RETURN CALL MXDPGS(N,H,1.0D0/GAM) RETURN END * SUBROUTINE PUDBI1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VARIABLE METRIC UPDATES. * * PARAMETERS : * II N NUMBER OF VARIABLES. * RU H(N*(N+1)/2) UPDATED APPROXIMATION OF THE INVERSE HESSIAN * MATRIX. * RA S(N) AUXILIARY VECTOR. * RI XO(N) VECTOR OF VARIABLES DIFFERENCE. * RI GO(N) GRADIENT DIFFERENCE. * RO R VALUE OF THE STEPSIZE PARAMETER. * RI PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE. * RO PAR1 PARAMETER FOR CONTROL SCALING. * RO PAR2 PARAMETER FOR CONTROL SCALING. * RO F VALUE OF THE OBJECTIVE FUNCTION. * RI FO INITIAL VALUE OF THE OBJECTIVE FUNCTION. * RI P CURRENT VALUE OF THE DIRECTIONAL DERIVATIVE. * II NIT NUMBER OF ITERATIONS. * II KIT INDEX OF THE ITERATION WITH THE LAST RESTART. * II MET VARIABLE METRIC UPDATE. * II MET1 SCALING STRATEGY. * II MET2 CORRECTION RULE. * IU IDECF DECOMPOSITION INDICATOR. * II 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 ITERH UPDATE INDICATOR. ITERH=0-SUCCESSFUL UPDATE. * ITERH>0-UNSUCCESSFUL UPDATE. * * METHOD : * VARIOUS VARIABLE METRIC UPDATES INCLUDING BFGS UPDATE. * SUBROUTINE PUDBI1(N,H,S,XO,GO,R,PO,PAR1,PAR2,F,FO,P,NIT,KIT, & MET,MET1,MET2,IDECF,ITERD,ITERH) IMPLICIT NONE INTEGER N,NIT,KIT,MET,MET1,MET2,IDECF,ITERD,ITERH DOUBLE PRECISION H(*),S(*),XO(*),GO(*),R,PO DOUBLE PRECISION PAR1,PAR2 DOUBLE PRECISION F,FO,P DOUBLE PRECISION AA,CC DOUBLE PRECISION MXVDOT DOUBLE PRECISION DIS,POM,POM3,POM4,A,B,C,GAM,RHO,PAR DOUBLE PRECISION DEN LOGICAL L1,L2,L3 IF (MET.LE.0) GO TO 22 IF (IDECF.NE.9) THEN ITERH=-1 GO TO 22 ENDIF L1=ABS(MET1).GE.3.OR.ABS(MET1).EQ.2.AND.NIT.EQ.KIT L3=.NOT.L1 * * DETERMINATION OF THE PARAMETERS A, B, C * B=MXVDOT(N,XO,GO) IF (B.LE.0.0D 0) THEN ITERH=2 GO TO 22 ENDIF CALL MXDSMM(N,H,GO,S) A=MXVDOT(N,GO,S) IF (A.LE.0.0D 0) THEN ITERH=1 GO TO 22 ENDIF IF(MET.NE.1.OR.L1) THEN IF (ITERD.NE.1) THEN C=0.0D 0 ELSE C=-R*PO IF (C.LE.0.0D 0) THEN ITERH=3 GO TO 22 ENDIF ENDIF ELSE C=0.0D 0 ENDIF * * DETERMINATION OF THE PARAMETER RHO (NONQUADRATIC PROPERTIES) * IF (MET2.EQ.1) THEN RHO=1.0D 0 ELSE IF (FO-F+P.EQ.0) THEN RHO=1.0D 0 ELSE RHO=0.5D 0*B/(FO-F+P) ENDIF IF(RHO.LE.1.0D-2) RHO=1.0D 0 IF(RHO.GE.1.0D 2) RHO=1.0D 0 AA=A/B CC=C/B IF (L1) THEN * * DETERMINATION OF THE PARAMETER GAM (SELF SCALING) * PAR=A/B POM3=0.7D 0 POM4=6.0D 0 GAM=RHO/PAR IF (NIT.NE.KIT) THEN IF (MET1.EQ.3) THEN L2=PAR2.LE.0.0D 0 L3=L2.AND.ABS(PAR1).LE.0.2D 0 L3=L3.OR.(.NOT.L2.AND.GAM.GT.1.0D 0) L3=L3.OR.(L2.AND.PAR1.LT.0.0D 0.AND.GAM.GT.1.0D 0) L3=L3.OR.(L2.AND.PAR1.GT.0.0D 0.AND.GAM.LT.1.0D 0) L3=L3.OR.GAM.LT.POM3 L3=L3.OR.GAM.GT.POM4 ELSE IF (MET1.EQ.4) THEN L3=GAM.LT.POM3.OR.GAM.GT.POM4 ENDIF ENDIF ENDIF IF (L3) THEN GAM=1.0D 0 PAR=RHO/GAM ENDIF IF (MET.EQ.1) GO TO 17 * * NEW UPDATE * POM=1.0D 0/(AA*CC) IF (POM.LT.1.0D 0) THEN POM=MAX(1.0D-15,(SQRT(C/A)-POM)/(1.0D 0-POM)) GO TO 20 ENDIF 17 CONTINUE * * BFGS UPDATE * POM=1.0D 0 DIS=PAR+AA CALL MXVDIR(N,-DIS,XO,S,XO) DIS=1.0D 0/(B*DIS) CALL MXDSMU(N,H,DIS,XO) CALL MXDSMU(N,H,-DIS,S) GO TO 21 20 CONTINUE * * GENERAL UPDATE * DEN=PAR+POM*AA DIS=POM/DEN CALL MXDSMU(N,H,(PAR*DIS-1.0D 0)/A,S) CALL MXVDIR(N,-DIS,S,XO,S) CALL MXDSMU(N,H,DEN/B,S) 21 CONTINUE IF (GAM.EQ.1.0D 0) GO TO 22 * * SCALING * CALL MXDSMS(N,H,GAM) 22 CONTINUE ITERH=0 RETURN END * SUBROUTINE PUDBQ1 ALL SYSTEMS 97/12/01 C PORTABILITY : ALL SYSTEMS C 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * BROYDEN GOOD UPDATE OF A RECTANGULAR MATRIX AFTER THE QR * DECOMPOSITION. * * PARAMETERS : * II N NUMBER OF VARIABLES. * II NA NUMBER OF EQUATIONS. * RU H(N*(N+1)/2) UPDATED UPPER TRIANGULAR MATRIX. * RI ETA2 PARAMETER WHICH CONTROLS A NONSINGULARITY * RU AG(N*NA) UPDATED RECTANGULAR MATRIX. * RA S(N) AUXILIARY VECTOR. * RI XO(N) VECTOR OF VARIABLES DIFFERENCE. * RI AFO(NA) RIGHT HAND SIDES DIFFERENCE. * II MET VARIABLE METRIC UPDATE. * IO ITERH UPDATE INDICATOR. ITERH=0-SUCCESSFUL UPDATE. * ITERH>0-UNSUCCESSFUL UPDATE. * IU IDECA DECOMPOSITION INDICATOR. * II NDECA NUMBER OF DECOMPOSITIONS. * * METHOD : * VARIOUS VARIABLE METRIC UPDATES INCLUDING BFGS UPDATE. * SUBROUTINE PUDBQ1(N,NA,H,ETA2,AG,S,XO,AFO,MET,ITERH,IDECA, & NDECA) IMPLICIT NONE INTEGER N,NA,MET,INF,ITERH,IDECA,NDECA DOUBLE PRECISION H(*),ETA2,AG(*),S(*),XO(*),AFO(*) DOUBLE PRECISION DEN,MXVDOT IF (MET.LE.0) RETURN IF (IDECA.EQ.0) THEN * * QR DECOMPOSITION * DEN=ETA2 CALL MXDRQF(N,NA,AG,H) CALL MXDPRC(N,H,INF,DEN) NDECA=NDECA+1 IDECA=1 ELSE IF (IDECA.NE.1) THEN ITERH=-1 GO TO 22 ENDIF * * THE GOOD BROYDEN UPDATE * DEN=MXVDOT(N,XO,XO) IF (DEN.LE.0.0D 0) THEN ITERH=2 GO TO 22 ENDIF CALL MXVCOP(N,XO,S) 21 CONTINUE CALL MXVNEG(N,XO,XO) CALL MXDPRM(N,H,XO,1) CALL MXDRMD(N,NA,AG,XO,1.0D 0,AFO,AFO) CALL MXDRQU(N,NA,AG,H,1.0D 0/DEN,AFO,S,XO,INF) ITERH=0 22 CONTINUE RETURN END * SUBROUTINE PUDRV1 ALL SYSTEMS 97/12/01 * PORTABILITY : ALL SYSTEMS * 97/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DRIVER FOR HYBRID QUASI-NEWTON UPDATES. * * PARAMETERS: * RI R VALUE OF THE STEPSIZE PARAMETER. * RI FO PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. * RI F CURRENT VALUE OF THE OBJECTIVE FUNCTION. * RI PO PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE. * II IPOM1 UPDATE SELECTION. * II IPOM2 METHOD SELECTION. * IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. * II IREST RESTART SPECIFICATION. IF IREST=0 DOES NOT HOLD THEN A * RESTART IS PERFORMED. * SUBROUTINE PUDRV1(R,FO,F,PO,IPOM1,IPOM2,NRED,IREST) IMPLICIT NONE INTEGER IPOM1,IPOM2,NRED,IREST DOUBLE PRECISION R,FO,F,PO DOUBLE PRECISION POM DOUBLE PRECISION CON1,CON2 PARAMETER(CON1=1.0D-1,CON2=1.0D-2) POM=(FO-F)/FO GO TO (1,2,3,4) IPOM2 1 IREST=1 IF (NRED.LE.0) THEN IPOM1=2 IREST=0 ELSE IPOM1=0 ENDIF GO TO 5 2 IREST=1 IF(POM.GE.CON2) THEN IPOM1=0 ELSE IF (F-FO.LE.R*PO) THEN IPOM1=0 ELSE IPOM1=1 IREST=0 ENDIF GO TO 5 3 IREST=1 IF (NRED.LE.0) THEN IF (IPOM1.NE.1) THEN IPOM1=2 IREST=0 ELSE IPOM1=0 ENDIF ELSE IF (POM.GE.CON2) THEN IPOM1=0 ELSE IF (IPOM1.NE.2) THEN IPOM1=1 IREST=0 ELSE IPOM1=0 ENDIF GO TO 5 4 IREST=1 IPOM1=0 5 CONTINUE RETURN 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, STAT) IMPLICIT NONE INTEGER NF,N,NC,IX(*),IC(*),ICA(*),KBF,KBC,INEW,IER,ITERM DOUBLE PRECISION X(*),XL(*),XU(*),CF(*),CFD(*),CL(*),CU(*), & CG(*),CR(*),CZ(*),H(*),S(*),R,EPS7,EPS9,GMAX,UMAX INTEGER I,J,K,L,IJ,IK,KC,KJ,KK,LL DOUBLE PRECISION DEN,TEMP c INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH INTEGER STAT(8) c COMMON /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH IF (KBC.GT.0) THEN IF (R.NE.0.0D 0) CALL PLDIRL(NC,CF,CFD,IC,R,KBC) IF (INEW.NE.0) THEN IF (KBF.GT.0) THEN DO 1 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, & STAT(4),IER) CALL MXVIND(IX,I,IER) IF (IER.LT.0) THEN ITERM=-15 RETURN ENDIF 1 CONTINUE ENDIF DO 2 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, & STAT(4),IER) CALL MXVIND(IC,KC,IER) IF (IER.LT.0) THEN ITERM=-15 RETURN ENDIF 2 CONTINUE ENDIF ELSE IF (KBF.GT.0) THEN K=0 DO 7 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.0D 0) THEN IJ=0 KJ=KK DO 4 J=1,N IF (J.LE.K) THEN KJ=KJ+1 ELSE KJ=KJ+J-1 ENDIF IF (J.NE.K) TEMP=H(KJ)/DEN IK=KK DO 3 I=1,J IF (I.LE.K) THEN IK=IK+1 ELSE IK=IK+I-1 ENDIF IJ=IJ+1 IF (I.NE.K.AND.J.NE.K) H(IJ)=H(IJ)+TEMP*H(IK) 3 CONTINUE 4 CONTINUE ENDIF LL=KK+K DO 6 I=K+1,N DO 5 J=1,I LL=LL+1 IF (J.NE.K) THEN KK=KK+1 H(KK)=H(LL) ENDIF 5 CONTINUE 6 CONTINUE N=N-1 ENDIF 7 CONTINUE ENDIF RETURN END * SUBROUTINE PYFUT1 ALL SYSTEMS 98/12/01 C PORTABILITY : ALL SYSTEMS C 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * TERMINATION CRITERIA AND TEST ON RESTART. * * PARAMETERS : * II N ACTUAL NUMBER OF VARIABLES. * RI F NEW VALUE OF THE OBJECTIVE FUNCTION. * RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. * RI UMAX MAXIMUN ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER. * RO GMAX NORM OF THE TRANSFORMED GRADIENT. * RI DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES. * RI TOLX LOWER BOUND FOR STEPLENGTH. * RI TOLF LOWER BOUND FOR FUNCTION DECREASE. * RI TOLB LOWER BOUND FOR FUNCTION VALUE. * RI TOLG LOWER BOUND FOR GRADIENT. * II KD DEGREE OF REQUIRED DERIVATIVES. * IU NIT ACTUAL NUMBER OF ITERATIONS. * II KIT NUMBER OF THE ITERATION AFTER RESTART. * II MIT MAXIMUM NUMBER OF ITERATIONS. * IU NFV ACTUAL NUMBER OF COMPUTED FUNCTION VALUES. * II MFV MAXIMUM NUMBER OF COMPUTED FUNCTION VALUES. * IU NFG ACTUAL NUMBER OF COMPUTED GRADIENT VALUES. * II MFG MAXIMUM NUMBER OF COMPUTED GRADIENT VALUES. * IU NTESX ACTUAL NUMBER OF TESTS ON STEPLENGTH. * II MTESX MAXIMUM NUMBER OF TESTS ON STEPLENGTH. * IU NTESF ACTUAL NUMBER OF TESTS ON FUNCTION DECREASE. * II MTESF MAXIMUM NUMBER OF TESTS ON FUNCTION DECREASE. * II ITES SYSTEM VARIBLE WHICH SPECIFIES TERMINATION. IF ITES=0 * THEN TERMINATION IS SUPPRESSED. * II IRES1 RESTART SPECIFICATION. RESTART IS PERFORMED AFTER * IRES1*N+IRES2 ITERATIONS. * II IRES2 RESTART SPECIFICATION. RESTART IS PERFORMED AFTER * IRES1*N+IRES2 ITERATIONS. * IU IREST RESTART INDICATOR. RESTART IS PERFORMED IF IREST>0. * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. * ITERS=0 FOR ZERO STEP. * IO ITERM TERMINATION INDICATOR. ITERM=1-TERMINATION AFTER MTESX * UNSUFFICIENT STEPLENGTHS. ITERM=2-TERMINATION AFTER MTESF * UNSUFFICIENT FUNCTION DECREASES. ITERM=3-TERMINATION ON LOWER * BOUND FOR FUNCTION VALUE. ITERM=4-TERMINATION ON LOWER BOUND * FOR GRADIENT. ITERM=11-TERMINATION AFTER MAXIMUM NUMBER OF * ITERATIONS. ITERM=12-TERMINATION AFTER MAXIMUM NUMBER OF * COMPUTED FUNCTION VALUES. * SUBROUTINE PYFUT1(N,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) IMPLICIT NONE INTEGER N,KD,NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF, & ITES,IRES1,IRES2,IREST,ITERS,ITERM DOUBLE PRECISION F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLG,TOLB DOUBLE PRECISION TEMP IF (ITERM.LT.0) RETURN IF (ITES .LE.0) GO TO 2 IF (ITERS.EQ.0) GO TO 1 IF (NIT.LE.0) FO=F+MIN(SQRT(ABS(F)),ABS(F)/1.0D 1) IF (F.LE.TOLB) THEN ITERM = 3 RETURN ENDIF IF (KD.GT.0) THEN IF (GMAX.LE.TOLG.AND.UMAX.LE.TOLG) THEN ITERM = 4 RETURN ENDIF ENDIF IF (NIT.LE.0) THEN NTESX = 0 NTESF = 0 ENDIF IF (DMAX.LE.TOLX) THEN ITERM = 1 NTESX = NTESX + 1 IF (NTESX .GE. MTESX) RETURN ELSE NTESX = 0 ENDIF TEMP=ABS(FO-F)/MAX(ABS(F),1.0D 0) IF (TEMP.LE.TOLF) THEN ITERM = 2 NTESF = NTESF+1 IF (NTESF.GE.MTESF) RETURN ELSE NTESF = 0 ENDIF 1 IF (NIT.GE.MIT) THEN ITERM = 11 RETURN ENDIF IF (NFV.GE.MFV) THEN ITERM = 12 RETURN ENDIF IF (NFG.GE.MFG) THEN ITERM = 13 RETURN ENDIF 2 ITERM = 0 IF (N.GT.0.AND.NIT-KIT.GE.IRES1*N+IRES2) THEN IREST=MAX(IREST,1) ENDIF NIT = NIT + 1 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, STAT) IMPLICIT NONE INTEGER NF,N,IX(*),IC(*),ICA(*),KBF,KBC,IOLD,KOLD,KREM,IER, $ ITERM DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),H(*),EPS8,UMAX, & GMAX INTEGER I,J,K,KC,L c INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH INTEGER STAT(8) c COMMON /STAT/ NRES,NDEC,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,STAT(3),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.0D 0,H(K+1)) H(K+N)=1.0D 0 KC=ICA(NF-N+1) IF (KC.GT.0) THEN IC(KC)=-IC(KC) ELSE K=-KC IX(K)=-IX(K) ENDIF ENDIF ELSE IOLD=0 ENDIF ELSE IF (KBF.GT.0) THEN IF (UMAX.GT.EPS8*GMAX) THEN IX(IOLD)=MIN(ABS(IX(IOLD)),3) DO 1 I=N,KOLD,-1 GN(I+1)=GN(I) 1 CONTINUE GN(KOLD)=G(IOLD) N=N+1 K=N*(N-1)/2 L=K+N DO 3 I=N,KOLD,-1 DO 2 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.0D 0 L=L-1 ELSE H(L)=0.0D 0 L=L-1 ENDIF 2 CONTINUE 3 CONTINUE ELSE IOLD=0 KOLD=0 ENDIF ENDIF 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) IMPLICIT NONE INTEGER NF,N,IX(*),ITERS,KBF,KBC DOUBLE PRECISION X(*),XO(*),G(*),GO(*),CZ(*),SN(*),R,F,FO,P,PO, & DMAX 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) ENDIF DMAX = 0.0D 0 IF (KBC.GT.0) THEN DO 1 I=1,NF DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0)) 1 CONTINUE IF (N.GT.0) THEN CALL MXVSCL(N,R,SN,XO) CALL MXVCOP(NF,GO,SN) CALL MXDRMM(NF,N,CZ,SN,GO) ENDIF ELSE IF (KBF.GT.0) THEN K=0 DO 2 I=1,NF IF (IX(I).LT.0) GO TO 2 K=K+1 DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0)) XO(K)=XO(I) GO(K)=GO(I) 2 CONTINUE ELSE DO 3 I=1,NF DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0)) 3 CONTINUE ENDIF 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) IMPLICIT NONE INTEGER NF,N,IX(*),IC(*),ICA(*),KBF,KBC,IOLD,KOLD DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),UMAX,GMAX DOUBLE PRECISION TEMP,MXVMAX,MXVDOT INTEGER NCA,NCZ,I,J,K,KC IOLD=0 KOLD=0 UMAX=0.0D 0 GMAX=0.0D 0 IF (KBC.GT.0) THEN IF (NF.GT.N) THEN NCA=NF-N NCZ=N*NF CALL MXVCOP(NF,G,GN) DO 1 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) ENDIF 1 CONTINUE CALL MXDPRB(NCA,CR,CZ(NCZ+1),0) DO 2 J=1,NCA TEMP=CZ(NCZ+J) KC=ICA(J) IF (KC.GT.0) THEN K=IC(KC) ELSE I=-KC K=IX(I) ENDIF IF (K.LE.-5) THEN ELSE IF ((K.EQ.-1.OR.K.EQ.-3).AND.UMAX+TEMP.GE.0.0D 0) THEN ELSE IF ((K.EQ.-2.OR.K.EQ.-4).AND.UMAX-TEMP.GE.0.0D 0) THEN ELSE IOLD=J UMAX=ABS(TEMP) ENDIF 2 CONTINUE ENDIF IF (N.GT.0) THEN CALL MXDRMM(NF,N,CZ,G,GN) GMAX=MXVMAX(N,GN) ENDIF ELSE IF (KBF.GT.0) THEN J=0 IOLD=0 KOLD=0 DO 3 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.0D 0) THEN ELSE IF ((K.EQ.-2.OR.K.EQ.-4).AND.UMAX-TEMP.GE.0.0D 0) THEN ELSE IOLD=I KOLD=J+1 UMAX=ABS(TEMP) ENDIF 3 CONTINUE N=J ELSE DO 4 I=1,NF TEMP=G(I) GMAX=MAX(GMAX,ABS(TEMP)) 4 CONTINUE N=NF ENDIF RETURN END * SUBROUTINE PYTRBH ALL SYSTEMS 98/12/01 C PORTABILITY : ALL SYSTEMS C 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * HESSIAN MATRIX OF THE OBJECTIVE FUNCTION OR ITS APPROXIMATION IS * SCALED AND REDUCED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * RI CR(NF*(NF+1)/2) TRIANGULAR DECOMPOSITION OF KERNEL OF THE * ORTHOGONAL PROJECTION. * RI CZ(NF*NF) MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE * CURRENT REDUCED SUBSPACE. * RI H(NF*(NF+1)/2) HESSIAN MATRIX OR ITS APPROXIMATION. * RA S(NF) AUXILIARY VECTOR. * II LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. * * SUBPROGRAMS USED : * S MXDSMM MATRIX VECTOR PRODUCT. * S MXVCOP COPYING OF A VECTOR. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * SUBROUTINE PYTRBH(NF,N,CR,CZ,H,S,LD,ITERS) IMPLICIT NONE INTEGER NF,N,LD,ITERS DOUBLE PRECISION CR(*),CZ(*),H(*),S(*) DOUBLE PRECISION MXVDOT INTEGER NCA,NCR,ICZ,JCZ,I,J,K IF (LD.NE.2.OR.ITERS.EQ.0) RETURN IF (N.LE.0) RETURN NCA=NF-N NCR=NCA*(NCA+1)/2 K=NCR JCZ=1 DO 4 J=1,N CALL MXDSMM(NF,H,CZ(JCZ),S) ICZ=1 DO 3 I=1,J K=K+1 CR(K)=MXVDOT(NF,CZ(ICZ),S) ICZ=ICZ+NF 3 CONTINUE JCZ=JCZ+NF 4 CONTINUE CALL MXVCOP(N*(N+1)/2,CR(NCR+1),H) 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) IMPLICIT NONE INTEGER NF,N,NC,IX(*),IC(*),KBF,KBC,KREM,INEW DOUBLE PRECISION X(*),XO(*),XL(*),XU(*),G(*),GO(*),CF(*),CFD(*), & CL(*),CU(*),CG(*),CZ(*),SN(*),S(*),RO,FP,FO,F,PO,P,RMAX INTEGER I,K FP = FO RO = 0.0D 0 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.0D 0,S) ENDIF ELSE IF (KBF.GT.0) THEN K=N+1 DO 1 I=NF,1,-1 IF (IX(I).LT.0) THEN S(I)=0.0D 0 ELSE K=K-1 S(I)=SN(K) ENDIF 1 CONTINUE INEW=0 CALL PLMAXS(NF,X,IX,XL,XU,S,RMAX,KBF,KREM,INEW) ENDIF 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) IMPLICIT NONE 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 * SUBROUTINE PYTRND ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD FOR MINIMAX * APPROXIMATION. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II N ACTUAL NUMBER OF VARIABLES. * II NC NUMBER OF CONSTRAINTS. * RI X(NF) VECTOR OF VARIABLES. * RI XN(NF) VECTOR OF SCALING FACTORS. * RO XO(NF) SAVED VECTOR OF VARIABLES. * II ICA(NF) VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS. * RI CG(NF*NC) MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR * CONSTRAINTS. * RO CZ(NF) VECTOR OF LAGRANGE MULTIPLIERS. * RO CZS(NF) SAVED VECTOR OF LAGRANGE MULTIPLIERS. * RI G(NF) GRADIENT OF THE LAGRANGIAN FUNCTION. * RI GO(NF) SAVED GRADIENT OF THE LAGRANGIAN FUNCTION. * RO R VALUE OF THE STEPSIZE PARAMETER. * RO F NEW VALUE OF THE OBJECTIVE FUNCTION. * RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. * RO P NEW VALUE OF THE DIRECTIONAL DERIVATIVE. * RI PO OLD VALUE OF THE DIRECTIONAL DERIVATIVE. * RI CMAX VALUE OF THE CONSTRAINT VIOLATION. * RO CMAXO SAVED VALUE OF THE CONSTRAINT VIOLATION. * RO DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES. * * COMMON DATA : * II NORMF SCALING SPECIFICATION. NORMF=0-NO SCALING PERFORMED. * NORMF=1-SCALING FACTORS ARE DETERMINED AUTOMATICALLY. * NORMF=2-SCALING FACTORS ARE SUPPLIED BY USER. * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. * ITERS=0 FOR ZERO STEP. * * SUBPROGRAMS USED : * S MXVCOP COPYING OF A VECTOR. * 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 PYTRND(NF,N,X,XO,ICA,CG,CZ,G,GO,R,F,FO,P,PO,CMAX,CMAXO, & DMAX,KD,LD,ITERS) IMPLICIT NONE INTEGER NF,N,KD,LD,ITERS INTEGER ICA(*) DOUBLE PRECISION X(*),XO(*),CG(*),CZ(*),G(*),GO(*),R,F,FO,P,PO, & CMAX,CMAXO,DMAX INTEGER I,J,L DO 2 J=1,NF-N L=ICA(J) IF (L.GT.0) THEN CALL MXVDIR(NF,-CZ(J),CG((L-1)*NF+1),G,G) ELSE L=-L G(L)=G(L)-CZ(J) ENDIF 2 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 F = FO P = PO CMAX = CMAXO CALL MXVSAV(NF,X,XO) CALL MXVSAV(NF,G,GO) LD = KD ENDIF DMAX = 0.0D0 DO 1 I=1,NF DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D0)) 1 CONTINUE N=NF RETURN END * SUBROUTINE PYTRUD 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 SCALED. TEST VALUE DMAX IS DETERMINED. * * PARAMETERS : * II NF DECLARED 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. * RO R VALUE OF THE STEPSIZE PARAMETER. * RO F NEW VALUE OF THE OBJECTIVE FUNCTION. * RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. * RO P NEW VALUE OF THE DIRECTIONAL DERIVATIVE. * RI PO OLD VALUE OF THE DIRECTIONAL DERIVATIVE. * RO DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES. * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. * ITERS=0 FOR ZERO STEP. * * SUBPROGRAMS USED : * S PYSET1 DEGREE DEFINITION OF THE COMPUTED DERIVATIVES. * S MXVDIF DIFFERENCE OF TWO VECTORS. * S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE * SUBSTRACTED ONE. * SUBROUTINE PYTRUD(NF,X,XO,G,GO,R,F,FO,P,PO,DMAX,KD,LD,ITERS) IMPLICIT NONE INTEGER NF,KD,LD,ITERS DOUBLE PRECISION X(*),XO(*),G(*),GO(*),R,F,FO,P,PO,DMAX INTEGER I 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) LD=KD ENDIF DMAX = 0.0D 0 DO 1 I=1,NF DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0)) 1 CONTINUE RETURN END * SUBROUTINE PYTRUF ALL SYSTEMS 98/12/01 * PORTABILITY : ALL SYSTEMS * 98/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VECTORS OF VARIABLES DIFFERENCE AND RIGHT HAND SIDES DIFFERENCE ARE * COMPUTED AND SCALED. TEST VALUE DMAX IS DETERMINED. * * PARAMETERS : * II NF DECLARED NUMBER OF VARIABLES. * II NA NUMBER OF APPROXIMATED FUNCTIONS. * RI X(NF) VECTOR OF VARIABLES. * RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. * RI AF(NA) VECTOR OF RIGHT HAND SIDES. * RI AFO(NA) VECTOR OF RIGHT HAND SIDES DIFFERENCE. * RO R VALUE OF THE STEPSIZE PARAMETER. * RO F NEW VALUE OF THE OBJECTIVE FUNCTION. * RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. * RO P NEW VALUE OF THE DIRECTIONAL DERIVATIVE. * RI PO OLD VALUE OF THE DIRECTIONAL DERIVATIVE. * RO DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES. * II KD DEGREE OF REQUIRED DERVATIVES. * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. * ITERS=0 FOR ZERO STEP. * * SUBPROGRAMS USED : * S PYSET1 DEGREE DEFINITION OF THE COMPUTED DERIVATIVES. * S MXVDIF DIFFERENCE OF TWO VECTORS. * S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE * SUBSTRACTED ONE. * SUBROUTINE PYTRUF(NF,NA,X,XO,AF,AFO,R,F,FO,P,PO,DMAX,KD,LD, & ITERS) IMPLICIT NONE INTEGER NF,NA,KD,LD,ITERS DOUBLE PRECISION X(*),XO(*),AF(*),AFO(*),R,F,FO,P,PO,DMAX INTEGER I IF (ITERS.GT.0) THEN CALL MXVDIF(NF,X,XO,XO) CALL MXVDIF(NA,AF,AFO,AFO) PO=R*PO P=R*P ELSE R = 0.0D 0 F = FO P = PO CALL MXVSAV(NF,X,XO) CALL MXVSAV(NA,AF,AFO) LD=KD ENDIF DMAX = 0.0D 0 DO 1 I=1,NF DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0)) 1 CONTINUE RETURN END * SUBROUTINE MXDCMM ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * MULTIPLICATION OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A * BY A VECTOR X. * * PARAMETERS : * II N NUMBER OF ROWS OF THE MATRIX A. * II M NUMBER OF COLUMNS OF THE MATRIX A. * RI A(N*M) RECTANGULAR MATRIX STORED COLUMNWISE IN THE * ONE-DIMENSIONAL ARRAY. * RI X(M) INPUT VECTOR. * RO Y(N) OUTPUT VECTOR EQUAL TO A*X. * * SUBPROGRAMS USED : * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE MXDCMM(N,M,A,X,Y) IMPLICIT NONE INTEGER N,M DOUBLE PRECISION A(*),X(*),Y(*) INTEGER J,K CALL MXVSET(N,0.0D 0,Y) K=0 DO 1 J=1,M CALL MXVDIR(N,X(J),A(K+1),Y,Y) K=K+N 1 CONTINUE RETURN END * SUBROUTINE MXDPGB ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A DENSE SYMMETRIC * POSITIVE DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L) * OBTAINED BY THE SUBROUTINE MXDPGF. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE MXDPGF. * RU X(N) ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR * EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR * EQUATIONS. * II JOB OPTION. IF JOB=0 THEN X:=(A+E)**(-1)*X. IF JOB>0 THEN * X:=L**(-1)*X. IF JOB<0 THEN X:=TRANS(L)**(-1)*X. * * METHOD : * BACK SUBSTITUTION * SUBROUTINE MXDPGB(N,A,X,JOB) IMPLICIT NONE INTEGER JOB,N DOUBLE PRECISION A(*),X(*) INTEGER I,II,IJ,J IF (JOB.GE.0) THEN * * PHASE 1 : X:=L**(-1)*X * IJ = 0 DO 20 I = 1,N DO 10 J = 1,I - 1 IJ = IJ + 1 X(I) = X(I) - A(IJ)*X(J) 10 CONTINUE IJ = IJ + 1 20 CONTINUE END IF IF (JOB.EQ.0) THEN * * PHASE 2 : X:=D**(-1)*X * II = 0 DO 30 I = 1,N II = II + I X(I) = X(I)/A(II) 30 CONTINUE END IF IF (JOB.LE.0) THEN * * PHASE 3 : X:=TRANS(L)**(-1)*X * II = N* (N-1)/2 DO 50 I = N - 1,1,-1 IJ = II DO 40 J = I + 1,N IJ = IJ + J - 1 X(I) = X(I) - A(IJ)*X(J) 40 CONTINUE II = II - I 50 CONTINUE END IF RETURN END * SUBROUTINE MXDPGF ALL SYSTEMS 89/12/01 * PORTABILITY : ALL SYSTEMS * 89/12/01 LU : ORIGINAL VERSION * * PURPOSE : * FACTORIZATION A+E=L*D*TRANS(L) OF A DENSE SYMMETRIC POSITIVE DEFINITE * MATRIX A+E WHERE D AND E ARE DIAGONAL POSITIVE DEFINITE MATRICES AND * L IS A LOWER TRIANGULAR MATRIX. IF A IS SUFFICIENTLY POSITIVE * DEFINITE THEN E=0. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) ON INPUT A GIVEN DENSE SYMMETRIC (USUALLY POSITIVE * DEFINITE) MATRIX A STORED IN THE PACKED FORM. ON OUTPUT THE * COMPUTED FACTORIZATION A+E=L*D*TRANS(L). * IO INF AN INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. IF * INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE AND E=0. IF * INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE AND E>0. IF * INF>0 THEN A IS INDEFINITE AND INF IS AN INDEX OF THE * MOST NEGATIVE DIAGONAL ELEMENT USED IN THE FACTORIZATION * PROCESS. * RU ALF ON INPUT A DESIRED TOLERANCE FOR POSITIVE DEFINITENESS. ON * OUTPUT THE MOST NEGATIVE DIAGONAL ELEMENT USED IN THE * FACTORIZATION PROCESS (IF INF>0). * RO TAU MAXIMUM DIAGONAL ELEMENT OF THE MATRIX E. * * METHOD : * P.E.GILL, W.MURRAY : NEWTON TYPE METHODS FOR UNCONSTRAINED AND * LINEARLY CONSTRAINED OPTIMIZATION, MATH. PROGRAMMING 28 (1974) * PP. 311-350. * SUBROUTINE MXDPGF(N,A,INF,ALF,TAU) IMPLICIT NONE DOUBLE PRECISION ALF,TAU INTEGER INF,N DOUBLE PRECISION A(*) DOUBLE PRECISION BET,DEL,GAM,RHO,SIG,TOL INTEGER I,IJ,IK,J,K,KJ,KK,L L = 0 INF = 0 TOL = ALF * * ESTIMATION OF THE MATRIX NORM * ALF = 0.0D0 BET = 0.0D0 GAM = 0.0D0 TAU = 0.0D0 KK = 0 DO 20 K = 1,N KK = KK + K BET = MAX(BET,ABS(A(KK))) KJ = KK DO 10 J = K + 1,N KJ = KJ + J - 1 GAM = MAX(GAM,ABS(A(KJ))) 10 CONTINUE 20 CONTINUE BET = MAX(TOL,BET,GAM/N) * DEL = TOL*BET DEL = TOL*MAX(BET,1.0D0) KK = 0 DO 60 K = 1,N KK = KK + K * * DETERMINATION OF A DIAGONAL CORRECTION * SIG = A(KK) IF (ALF.GT.SIG) THEN ALF = SIG L = K END IF GAM = 0.0D0 KJ = KK DO 30 J = K + 1,N KJ = KJ + J - 1 GAM = MAX(GAM,ABS(A(KJ))) 30 CONTINUE GAM = GAM*GAM RHO = MAX(ABS(SIG),GAM/BET,DEL) IF (TAU.LT.RHO-SIG) THEN TAU = RHO - SIG INF = -1 END IF * * GAUSSIAN ELIMINATION * A(KK) = RHO KJ = KK DO 50 J = K + 1,N KJ = KJ + J - 1 GAM = A(KJ) A(KJ) = GAM/RHO IK = KK IJ = KJ DO 40 I = K + 1,J IK = IK + I - 1 IJ = IJ + 1 A(IJ) = A(IJ) - A(IK)*GAM 40 CONTINUE 50 CONTINUE 60 CONTINUE IF (L.GT.0 .AND. ABS(ALF).GT.DEL) INF = L RETURN END * FUNCTION MXDPGP ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF THE NUMBER MXDPGP=TRANS(X)*D**(-1)*Y WHERE D IS A * DIAGONAL MATRIX IN THE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE MXDPGF. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE MXDPGF. * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RR MXDPGP COMPUTED NUMBER MXDPGP=TRANS(X)*D**(-1)*Y. * FUNCTION MXDPGP(N,A,X,Y) IMPLICIT NONE INTEGER N DOUBLE PRECISION A(*),X(*),Y(*),MXDPGP DOUBLE PRECISION TEMP INTEGER I,J J = 0 TEMP = 0.0D0 DO 10 I = 1,N J = J + I TEMP = TEMP + X(I)*Y(I)/A(J) 10 CONTINUE MXDPGP = TEMP RETURN END * SUBROUTINE MXDPGS ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SCALING OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX A+E USING THE * FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXDPGF. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE MXDPGF. * RI ALF SCALING FACTOR. * SUBROUTINE MXDPGS(N,A,ALF) IMPLICIT NONE DOUBLE PRECISION ALF INTEGER N DOUBLE PRECISION A(*) INTEGER I,J J = 0 DO 10 I = 1,N J = J + I A(J) = A(J)*ALF 10 CONTINUE RETURN END * SUBROUTINE MXDPGU ALL SYSTEMS 89/12/01 * PORTABILITY : ALL SYSTEMS * 89/12/01 LU : ORIGINAL VERSION * * PURPOSE : * CORRECTION OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX A+E IN THE * FACTORED FORM A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXDPGF. * THE CORRECTION IS DEFINED AS A+E:=A+E+ALF*X*TRANS(X) WHERE ALF IS A * GIVEN SCALING FACTOR AND X IS A GIVEN VECTOR. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE MXDPGF. * RI ALF SCALING FACTOR IN THE CORRECTION TERM. * RI X(N) VECTOR IN THE CORRECTION TERM. * RA Y(N) AUXILIARY VECTOR. * * METHOD : * P.E.GILL, W.MURRAY, M.SAUNDERS: METHODS FOR COMPUTING AND MODIFYING * THE LDV FACTORS OF A MATRIX, MATH. OF COMP. 29 (1974) PP. 1051-1077. * SUBROUTINE MXDPGU(N,A,ALF,X,Y) IMPLICIT NONE DOUBLE PRECISION ZERO,ONE,FOUR,CON PARAMETER (ZERO=0.0D0,ONE=1.0D0,FOUR=4.0D0,CON=1.0D-8) DOUBLE PRECISION ALF,ALFR INTEGER N DOUBLE PRECISION A(*),X(*),Y(*) DOUBLE PRECISION B,D,P,R,T,TO INTEGER I,II,IJ,J IF (ALF.GE.ZERO) THEN * * FORWARD CORRECTION IN CASE WHEN THE SCALING FACTOR IS NONNEGATIVE * ALFR = SQRT(ALF) CALL MXVSCL(N,ALFR,X,Y) TO = ONE II = 0 DO 30 I = 1,N II = II + I D = A(II) P = Y(I) T = TO + P*P/D R = TO/T A(II) = D/R B = P/ (D*T) IF (A(II).LE.FOUR*D) THEN * * AN EASY FORMULA FOR LIMITED DIAGONAL ELEMENT * IJ = II DO 10 J = I + 1,N IJ = IJ + J - 1 D = A(IJ) Y(J) = Y(J) - P*D A(IJ) = D + B*Y(J) 10 CONTINUE ELSE * * A MORE COMPLICATE BUT NUMERICALLY STABLE FORMULA FOR UNLIMITED * DIAGONAL ELEMENT * IJ = II DO 20 J = I + 1,N IJ = IJ + J - 1 D = A(IJ) A(IJ) = R*D + B*Y(J) Y(J) = Y(J) - P*D 20 CONTINUE END IF TO = T 30 CONTINUE ELSE * * BACKWARD CORRECTION IN CASE WHEN THE SCALING FACTOR IS NEGATIVE * ALFR = SQRT(-ALF) CALL MXVSCL(N,ALFR,X,Y) TO = ONE IJ = 0 DO 50 I = 1,N D = Y(I) DO 40 J = 1,I - 1 IJ = IJ + 1 D = D - A(IJ)*Y(J) 40 CONTINUE Y(I) = D IJ = IJ + 1 TO = TO - D*D/A(IJ) 50 CONTINUE IF (TO.LE.ZERO) TO = CON II = N* (N+1)/2 DO 70 I = N,1,-1 D = A(II) P = Y(I) T = TO + P*P/D A(II) = D*TO/T B = -P/ (D*TO) TO = T IJ = II DO 60 J = I + 1,N IJ = IJ + J - 1 D = A(IJ) A(IJ) = D + B*Y(J) Y(J) = Y(J) + P*D 60 CONTINUE II = II - I 70 CONTINUE END IF RETURN END * SUBROUTINE MXDPRB ALL SYSTEMS 89/12/01 * PORTABILITY : ALL SYSTEMS * 89/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A DENSE SYMMETRIC * POSITIVE DEFINITE MATRIX A USING THE FACTORIZATION A=TRANS(R)*R. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) FACTORIZATION A=TRANS(R)*R. * RU X(N) ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR * EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR * EQUATIONS. * II JOB OPTION. IF JOB=0 THEN X:=A**(-1)*X. IF JOB>0 THEN * X:=TRANS(R)**(-1)*X. IF JOB<0 THEN X:=R**(-1)*X. * * METHOD : * BACK SUBSTITUTION * SUBROUTINE MXDPRB(N,A,X,JOB) IMPLICIT NONE INTEGER JOB,N DOUBLE PRECISION A(*),X(*) INTEGER I,II,IJ,J IF (JOB.GE.0) THEN * * PHASE 1 : X:=TRANS(R)**(-1)*X * IJ = 0 DO 20 I = 1,N DO 10 J = 1,I - 1 IJ = IJ + 1 X(I) = X(I) - A(IJ)*X(J) 10 CONTINUE IJ = IJ + 1 X(I) = X(I)/A(IJ) 20 CONTINUE END IF IF (JOB.LE.0) THEN * * PHASE 2 : X:=R**(-1)*X * II = N* (N+1)/2 DO 40 I = N,1,-1 IJ = II DO 30 J = I + 1,N IJ = IJ + J - 1 X(I) = X(I) - A(IJ)*X(J) 30 CONTINUE X(I) = X(I)/A(II) II = II - I 40 CONTINUE END IF RETURN END * SUBROUTINE MXDPRC ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * CORRECTION OF A SINGULAR DENSE SYMMETRIC POSITIVE SEMIDEFINITE MATRIX * A DECOMPOSED AS A=TRANS(R)*R. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * IO INF AN INFORMATION OBTAINED IN THE CORRECTION PROCESS. IF * INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE. IF * INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE. * PROCESS. * RI TOL DESIRED TOLERANCE FOR POSITIVE DEFINITENESS. * SUBROUTINE MXDPRC(N,A,INF,TOL) IMPLICIT NONE INTEGER N,INF DOUBLE PRECISION A(N*(N+1)/2),TOL DOUBLE PRECISION TOL1,TEMP INTEGER L,I DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D 0) INF=0 TOL1=SQRT(TOL) TEMP=TOL1 DO 3 I=1,N*(N+1)/2 TEMP=MAX(TEMP,ABS(A(I))) 3 CONTINUE TEMP=TEMP*TOL1 L=0 DO 1 I=1,N L=L+I IF (ABS(A(L)).LE.TEMP) THEN A(L)=SIGN(TEMP,A(L)) INF=-1 ENDIF 1 CONTINUE RETURN END * SUBROUTINE MXDPRM ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * MULTIPLICATION OF A GIVEN VECTOR X BY A DENSE SYMMETRIC POSITIVE * DEFINITE MATRIX A USING THE FACTORIZATION A=TRANS(R)*R. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) FACTORIZATION A=TRANS(R)*R. * RU X(N) ON INPUT THE GIVEN VECTOR. ON OUTPUT THE RESULT OF * MULTIPLICATION. * II JOB OPTION. IF JOB=0 THEN X:=A*X. IF JOB>0 THEN X:=R*X. * IF JOB<0 THEN X:=TRANS(R)*X. * SUBROUTINE MXDPRM( N, A, X, JOB) IMPLICIT NONE INTEGER N, JOB DOUBLE PRECISION A(N*(N+1)/2), X(N) INTEGER I, J, II, IJ IF (JOB .GE. 0) THEN C C PHASE 1 : X:=R*X C II = 0 DO 3 I = 1, N II = II + I X(I) = A(II) * X(I) IJ = II DO 2 J = I+1, N IJ = IJ + J - 1 X(I) = X(I) + A(IJ) * X(J) 2 CONTINUE 3 CONTINUE ENDIF IF (JOB .LE. 0) THEN C C PHASE 2 : X:=TRANS(R)*X C IJ = N*(N+1)/2 DO 6 I = N, 1, -1 X(I) = A(IJ) * X(I) DO 5 J = I-1, 1, -1 IJ = IJ - 1 X(I) = X(I) + A(IJ) * X(J) 5 CONTINUE IJ = IJ - 1 6 CONTINUE ENDIF RETURN END * SUBROUTINE MXDRGR ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * PLANE ROTATION IS APPLIED TO A ROWWISE STORED DENSE RECTANGULAR * MATRIX A. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX A. * RU A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY. * II K FIRST INDEX OF THE PLANE ROTATION. * II L SECOND INDEX OF THE PLANE ROTATION. * RI CK DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * RI CL OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * II IER TYPE OF THE PLANE ROTATION. IER=0-GENERAL PLANE ROTATION. * IER=1-PERMUTATION. IER=2-TRANSFORMATION SUPPRESSED. * * SUBPROGRAMS USED : * S MXVROT PLANE ROTATION APPLIED TO TWO ELEMENTS. * SUBROUTINE MXDRGR(N,A,K,L,CK,CL,IER) IMPLICIT NONE INTEGER N,K,L,IER DOUBLE PRECISION A(*),CK,CL INTEGER I,IK,IL IF (IER.NE.0.AND.IER.NE.1) RETURN IK=(K-1)*N IL=(L-1)*N DO 1 I=1,N IK=IK+1 IL=IL+1 CALL MXVROT(A(IK),A(IL),CK,CL,IER) 1 CONTINUE RETURN END * SUBROUTINE MXDRMD ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A BY * A VECTOR X AND ADDITION OF A SCALED VECTOR ALF*Y. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX A. * II M NUMBER OF ROWS OF THE MATRIX A. * RI A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY. * RI X(N) INPUT VECTOR. * RI ALF SCALING FACTOR. * RI Y(M) INPUT VECTOR. * RO Z(M) OUTPUT VECTOR EQUAL TO A*X+ALF*Y. * SUBROUTINE MXDRMD(N,M,A,X,ALF,Y,Z) IMPLICIT NONE INTEGER N,M DOUBLE PRECISION A(M*N),X(N),ALF,Y(M),Z(M) DOUBLE PRECISION TEMP INTEGER I,J,K K=0 DO 2 J=1,M TEMP=ALF*Y(J) DO 1 I=1,N TEMP=TEMP+A(K+I)*X(I) 1 CONTINUE Z(J)=TEMP K=K+N 2 CONTINUE RETURN END * SUBROUTINE MXDRMI ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * ROWWISE STORED DENSE RECTANGULAR MATRIX A IS SET TO BE A PART OF THE * UNIT MATRIX. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX A. * II M NUMBER OF ROWS OF THE MATRIX A. * RO A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE ONE-DIMENSIONAL * ARRAY. THIS MATRIX IS SET TO TRANS([I,0]). * SUBROUTINE MXDRMI(N,M,A) IMPLICIT NONE INTEGER N,M DOUBLE PRECISION A(M*N) INTEGER I,J,K DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D 0,ONE=1.0D 0) K=0 DO 2 J=1,M DO 1 I=1,N A(I+K)=ZERO IF (I.EQ.J) A(I+K)=ONE 1 CONTINUE K=K+N 2 CONTINUE RETURN END * SUBROUTINE MXDRMM ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A BY * A VECTOR X. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX A. * II M NUMBER OF ROWS OF THE MATRIX A. * RI A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY. * RI X(N) INPUT VECTOR. * RO Y(M) OUTPUT VECTOR EQUAL TO A*X. * SUBROUTINE MXDRMM(N,M,A,X,Y) IMPLICIT NONE INTEGER N,M DOUBLE PRECISION A(*),X(*),Y(*) DOUBLE PRECISION TEMP INTEGER I,J,K K=0 DO 2 J=1,M TEMP=0.0D 0 DO 1 I=1,N TEMP=TEMP+A(K+I)*X(I) 1 CONTINUE Y(J)=TEMP K=K+N 2 CONTINUE RETURN END * FUNCTION MXDRMN ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * EUCLIDEAN NORM OF A PART OF THE I-TH COLUMN OF A ROWWISE STORED DENSE * RECTANGULAR MATRIX A IS COMPUTED. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX A. * II M NUMBER OF ROWS OF THE MATRIX A. * RI A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY. * II I INDEX OF THE COLUMN WHOSE NORM IS COMPUTED. * II J INDEX OF THE FIRST ELEMENT FROM WHICH THE NORM IS COMPUTED. * FUNCTION MXDRMN(N,M,A,I,J) IMPLICIT NONE INTEGER N,M,I,J DOUBLE PRECISION A(M*N),MXDRMN DOUBLE PRECISION POM,DEN INTEGER K,L DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D 0) DEN=ZERO L=(J-1)*N DO 1 K=J,M DEN=MAX(DEN,ABS(A(L+I))) L=L+N 1 CONTINUE POM=ZERO IF (DEN.GT.ZERO) THEN L=(J-1)*N DO 2 K=J,M POM=POM+(A(L+I)/DEN)**2 L=L+N 2 CONTINUE ENDIF MXDRMN=DEN*SQRT(POM) RETURN END * SUBROUTINE MXDRMV ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * K-TH COLUMN OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A IS COPIED * TO THE VECTOR X. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX A. * II M NUMBER OF ROWS OF THE MATRIX A. * RI A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY. * RO X(M) OUTPUT VECTOR SUCH THAT X(J)=A(J,K) FOR ALL J. * II K INDEX OF THE ROW BEING COPIED TO THE OUTPUT VECTOR. * SUBROUTINE MXDRMV(N,M,A,X,K) IMPLICIT NONE INTEGER N,M,K DOUBLE PRECISION A(*),X(*) INTEGER I,J IF (K.LT.1.OR.K.GT.N) RETURN I=K DO 1 J=1,M X(J)=A(I) I=I+N 1 CONTINUE RETURN END * SUBROUTINE MXDRQF ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * QR DECOMPOSITION OF ROWWISE STORED DENSE RECTANGULAR MATRIX Q USING * HOUSEHOLDER TRANSFORMATIONS WITHOUT PIVOTING. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX Q. * II M NUMBER OF ROWS OF THE MATRIX Q. * RU Q(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY. * RO R(N*(N+1)/2) UPPER TRIANGULAR MATRIX STORED IN THE PACKED FORM. * * SUBPROGRAMS USED : * S MXDRMN EUCLIDEAN NORM OF A PART OF THE ROWWISE STORED * RECTANGULAR MATRIX COLUMN. * * METHOD : * P.A.BUSSINGER, G.H.GOLUB : LINEAR LEAST SQUARES SOLUTION BY * HOUSEHOLDER TRANSFORMATION. NUMER. MATH. 7 (1965) 269-276. * SUBROUTINE MXDRQF(N,M,Q,R) IMPLICIT NONE INTEGER N,M DOUBLE PRECISION Q(M*N),R(N*(N+1)/2) DOUBLE PRECISION ALF,POM,MXDRMN INTEGER I,J,K,L,JP,KP,NM DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D 0,ONE=1.0D 0) NM=MIN(N,M) C C QR DECOMPOSITION C L=0 KP=0 DO 6 K=1,NM POM=MXDRMN(N,M,Q,K,K) IF (Q(KP+K).NE.ZERO) POM=SIGN(POM,Q(KP+K)) R(L+K)=-POM JP=0 DO 1 J=1,K-1 R(L+J)=Q(JP+K) Q(JP+K)=ZERO JP=JP+N 1 CONTINUE IF (POM.NE.ZERO) THEN C C HOUSEHOLDER TRANSFORMATION C DO 2 J=K,M Q(JP+K)=Q(JP+K)/POM JP=JP+N 2 CONTINUE Q(KP+K)=Q(KP+K)+ONE DO 5 I=K+1,N ALF=ZERO JP=KP DO 3 J=K,M ALF=ALF+Q(JP+K)*Q(JP+I) JP=JP+N 3 CONTINUE ALF=ALF/Q(KP+K) JP=KP DO 4 J=K,M Q(JP+I)=Q(JP+I)-ALF*Q(JP+K) JP=JP+N 4 CONTINUE 5 CONTINUE ENDIF L=L+K KP=KP+N 6 CONTINUE C C EXPLICIT FORMULATION OF THE ORTHOGONAL MATRIX C KP=N*N DO 11 K=N,1,-1 KP=KP-N IF (Q(KP+K).NE.ZERO) THEN DO 9 I=K+1,N ALF=ZERO JP=KP DO 7 J=K,M ALF=ALF+Q(JP+K)*Q(JP+I) JP=JP+N 7 CONTINUE ALF=ALF/Q(KP+K) JP=KP DO 8 J=K,M Q(JP+I)=Q(JP+I)-ALF*Q(JP+K) JP=JP+N 8 CONTINUE 9 CONTINUE JP=KP DO 10 J=K,M Q(JP+K)=-Q(JP+K) JP=JP+N 10 CONTINUE ENDIF Q(KP+K)=Q(KP+K)+ONE 11 CONTINUE RETURN END * SUBROUTINE MXDRQU ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * UPDATE OF A QR DECOMPOSITION. THIS QR DECOMPOSITION IS UPDATED * BY THE RULE Q*R:=Q*R+ALF*X*TRANS(Y). * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX Q. * II M NUMBER OF ROWS OF THE MATRIX Q. * RU Q(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY (PART OF THE ORTHOGONAL MATRIX). * RU R(N*(N+1)/2) UPPER TRIANGULAR MATRIX STORED IN A PACKED FORM. * RI ALF SCALAR PARAMETER. * RI X(M) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RA Z(N) AUXILIARY VECTOR. * IO INF INFORMATION. IF INF=0 THEN X LIES IN THE COLUMN SPACE OF Q. * IF INF=1 THEN X DOES NOT LIE IN THE COLUMN SPACE OF Q. * * SUBPROGRAMS USED : * RF MXVNOR EUCLIDEAN NORM OF A VECTOR. * S MXVORT COMPUTATION OF THE PLANE ROTATION MATRIX. * S MXVROT PLANE ROTATION IS APPLIED TO TWO NUMBERS. * * METHOD : * J.W.DANIEL, W.B.GRAGG, L.KAUFMAN, G.W.STEWARD : REORTHOGONALIZATION * AND STABLE ALGORITHMS FOR UPDATING THE GRAM-SCHMIDT QR FACTORIZATION. * MATHEMATICS OF COMPUTATION 30 (1976) 772-795. SUBROUTINE MXDRQU(N,M,Q,R,ALF,X,Y,Z,INF) IMPLICIT NONE INTEGER N,M,INF DOUBLE PRECISION Q(M*N),R(N*(N+1)/2),ALF,X(M),Y(N),Z(N) DOUBLE PRECISION CK,CL,ZK,ZL,MXVNOR INTEGER J,K,L,KJ,KK,IER DOUBLE PRECISION ONE,CON PARAMETER (ONE=1.0D 0,CON=1.0D-10) INF=0 KK=N*(N+1)/2 C C COMPUTATION OF THE VECTOR TRANS(Q)*X C CALL MXDCMM(N,M,Q,X,Z) IF (M.GT.N) THEN C C IF X DOES NOT LIE IN THE COLUMN SPACE OF Q WE HAVE TO USE C A SUBPROBLEM WHOSE DIMENSION IS BY ONE GREATER (INF=1). C ZK=MXVNOR(M,X) CALL MXDRMD(N,M,Q,Z,-ONE,X,X) ZL=MXVNOR(M,X) IF (ZL.GT.CON*ZK) THEN INF=1 CALL MXVSCL(M,-ONE/ZL,X,X) CALL MXVORT(Z(N),ZL,CK,CL,IER) IF (IER.EQ.0.OR.IER.EQ.1) THEN CALL MXVROT(R(KK),ZL,CK,CL,IER) KJ=N DO 1 J=1,M CALL MXVROT(Q(KJ),X(J),CK,CL,IER) KJ=KJ+N 1 CONTINUE ENDIF ENDIF ENDIF C C APPLICATION OF PLANE ROTATIONS TO THE VECTOR Z SO THAT C TRANS(Q1)*Z=E1 WHERE Q1 IS AN ORTHOGONAL MATRIX (ACCUMULATION OF C THE PLANE ROTATIONS) AND E1 IS THE FIRST COLUMN OF THE UNIT C MATRIX. AT THE SAME TIME BOTH THE UPPER HESSENBERG MATRIX C TRANS(Q1)*R AND THE ORTHOGONAL MATRIX Q*Q1 ARE CONSTRUCTED SO THAT C Q*Q1*R1=Q*Q1*(TRANS(Q1)*R+ALF*E1*TRANS(Y)) WHERE R1 IS AN UPPER C HESSENBERG MATRIX. C DO 4 L=N,2,-1 K=L-1 KK=KK-L CALL MXVORT(Z(K),Z(L),CK,CL,IER) IF (IER.EQ.0.OR.IER.EQ.1) THEN CALL MXVROT(R(KK),Z(L),CK,CL,IER) KJ=KK DO 2 J=L,N KJ=KJ+J-1 CALL MXVROT(R(KJ),R(KJ+1),CK,CL,IER) 2 CONTINUE KJ=K DO 3 J=1,M CALL MXVROT(Q(KJ),Q(KJ+1),CK,CL,IER) KJ=KJ+N 3 CONTINUE ENDIF 4 CONTINUE Z(1)=ALF*Z(1) KJ=1 DO 5 J=1,N R(KJ)=R(KJ)+Z(1)*Y(J) KJ=KJ+J 5 CONTINUE C C APPLICATION OF PLANE ROTATIONS TO THE UPPER HESSENBERG MATRIX R1 C GIVEN ABOVE SO THAT R2=TRANS(Q2)*R1 WHERE Q2 IS AN ORTHOGONAL C MATRIX (ACCUMULATION OF THE PLANE ROTATIONS) AND R2 IS AN UPPER C TRIANGULAR MATRIX. WE OBTAIN THE NEW QR DECOMPOSITION Q*Q1*Q2*R2. C KK=1 DO 8 L=2,N K=L-1 CALL MXVORT(R(KK),Z(L),CK,CL,IER) IF (IER.EQ.0.OR.IER.EQ.1) THEN KJ=KK DO 6 J=L,N KJ=KJ+J-1 CALL MXVROT(R(KJ),R(KJ+1),CK,CL,IER) 6 CONTINUE KJ=K DO 7 J=1,M CALL MXVROT(Q(KJ),Q(KJ+1),CK,CL,IER) KJ=KJ+N 7 CONTINUE ENDIF KK=KK+L 8 CONTINUE C C BACK TRANSFORMATION OF THE GREATER SUBPROBLEM IF INF=1. C IF (INF.EQ.1) THEN CALL MXVORT(R(KK),ZL,CK,CL,IER) IF (IER.EQ.0.OR.IER.EQ.1) THEN KJ=N DO 9 J=1,M CALL MXVROT(Q(KJ),X(J),CK,CL,IER) KJ=KJ+N 9 CONTINUE ENDIF ENDIF RETURN END * SUBROUTINE MXDSMI ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DENSE SYMMETRIC MATRIX A IS SET TO THE UNIT MATRIX WITH THE SAME * ORDER. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RO A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM * WHICH IS SET TO THE UNIT MATRIX (I.E. A:=I). * SUBROUTINE MXDSMI(N,A) IMPLICIT NONE INTEGER N DOUBLE PRECISION A(*) INTEGER I,M M = N* (N+1)/2 DO 10 I = 1,M A(I) = 0.0D0 10 CONTINUE M = 0 DO 20 I = 1,N M = M + I A(M) = 1.0D0 20 CONTINUE RETURN END * SUBROUTINE MXDSMM ALL SYSTEMS 89/12/01 * PORTABILITY : ALL SYSTEMS * 89/12/01 LU : ORIGINAL VERSION * * PURPOSE : * MULTIPLICATION OF A DENSE SYMMETRIC MATRIX A BY A VECTOR X. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * RI X(N) INPUT VECTOR. * RO Y(N) OUTPUT VECTOR EQUAL TO A*X. * SUBROUTINE MXDSMM(N,A,X,Y) IMPLICIT NONE INTEGER N DOUBLE PRECISION A(*),X(*),Y(*) DOUBLE PRECISION TEMP INTEGER I,J,K,L K = 0 DO 30 I = 1,N TEMP = 0.0D0 L = K DO 10 J = 1,I L = L + 1 TEMP = TEMP + A(L)*X(J) 10 CONTINUE DO 20 J = I + 1,N L = L + J - 1 TEMP = TEMP + A(L)*X(J) 20 CONTINUE Y(I) = TEMP K = K + I 30 CONTINUE RETURN END * SUBROUTINE MXDSMR ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * PLANE ROTATION IS APPLIED TO A DENSE SYMMETRIC MATRIX A. THE CASE * K=L+1 IS REQUIRED. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * II K FIRST INDEX OF PLANE ROTATION. * II L SECOND INDEX OF PLANE ROTATION. * RO CK DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * RO CL OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * IO IER INFORMATION ON THE TRANSFORMATION. IER<0-K OR L OUT OF * RANGE. IER=0-PLANE ROTATION. IER=1-PERMUTATION. * IER=2-TRANSFORMATION SUPPRESSED. * * SUBPROGRAMS USED : * S MXVROT PLANE ROTATION IS APPLIED TO TWO NUMBERS. * SUBROUTINE MXDSMR(N,A,K,L,CK,CL,IER) IMPLICIT NONE INTEGER N,K,L,IER DOUBLE PRECISION A(*),CK,CL DOUBLE PRECISION AKK,AKL,ALL,CKK,CKL,CLL INTEGER J,KJ,LJ,KK,KL,LL IF (IER.NE.0.AND.IER.NE.1) RETURN IF (K.NE.L+1) THEN IER=-1 RETURN ENDIF LJ=L*(L-1)/2 DO 1 J=1,N IF (J.LE.L) THEN LJ=LJ+1 KJ=LJ+L ELSE LJ=LJ+J-1 KJ=LJ+1 ENDIF IF (J.NE.K.AND.J.NE.L) THEN CALL MXVROT(A(KJ),A(LJ),CK,CL,IER) ENDIF 1 CONTINUE IF (IER.EQ.0) THEN CKK=CK**2 CKL=CK*CL CLL=CL**2 LL=L*(L+1)/2 KL=LL+L KK=LL+K AKL=(CKL+CKL)*A(KL) AKK=CKK*A(KK)+CLL*A(LL)+AKL ALL=CLL*A(KK)+CKK*A(LL)-AKL A(KL)=(CLL-CKK)*A(KL)+CKL*(A(KK)-A(LL)) A(KK)=AKK A(LL)=ALL ELSE LL=L*(L+1)/2 KK=LL+K AKK=A(KK) A(KK)=A(LL) A(LL)=AKK ENDIF RETURN END * SUBROUTINE MXDSMS ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SCALING OF A DENSE SYMMETRIC MATRIX. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM * WHICH IS SCALED BY THE VALUE ALF (I.E. A:=ALF*A). * RI ALF SCALING FACTOR. * SUBROUTINE MXDSMS( N, A, ALF) IMPLICIT NONE INTEGER N DOUBLE PRECISION A(N*(N+1)/2), ALF INTEGER I,M M = N * (N+1) / 2 DO 1 I = 1, M A(I) = A(I) * ALF 1 CONTINUE RETURN END * SUBROUTINE MXDSMU ALL SYSTEMS 89/12/01 C PORTABILITY : ALL SYSTEMS C 89/12/01 LU : ORIGINAL VERSION * * PURPOSE : * UPDATE OF A DENSE SYMMETRIC MATRIX A. THIS UPDATE IS DEFINED AS * A:=A+ALF*X*TRANS(X) WHERE ALF IS A GIVEN SCALING FACTOR AND X IS * A GIVEN VECTOR. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * RI ALF SCALING FACTOR IN THE CORRECTION TERM. * RI X(N) VECTOR IN THE CORRECTION TERM. * SUBROUTINE MXDSMU( N, A, ALF, X) IMPLICIT NONE INTEGER N DOUBLE PRECISION A(N*(N+1)/2), X(N), ALF DOUBLE PRECISION TEMP INTEGER I, J, K K = 0 DO 2 I = 1, N TEMP = ALF * X(I) DO 1 J = 1, I K = K + 1 A(K) = A(K) + TEMP * X(J) 1 CONTINUE 2 CONTINUE RETURN END * SUBROUTINE MXDSMV ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * K-TH ROW OF A DENSE SYMMETRIC MATRIX A IS COPIED TO THE VECTOR X. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * RO X(N) OUTPUT VECTOR. * II K INDEX OF COPIED ROW. * SUBROUTINE MXDSMV(N,A,X,K) IMPLICIT NONE INTEGER K,N DOUBLE PRECISION A(*),X(*) INTEGER I,L L = K* (K-1)/2 DO 10 I = 1,N IF (I.LE.K) THEN L = L + 1 ELSE L = L + I - 1 END IF X(I) = A(L) 10 CONTINUE RETURN END * SUBROUTINE MXVCOP ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COPYING OF A VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RO Y(N) OUTPUT VECTOR WHERE Y:= X. * SUBROUTINE MXVCOP(N,X,Y) IMPLICIT NONE INTEGER N DOUBLE PRECISION X(*),Y(*) INTEGER I DO 10 I = 1,N Y(I) = X(I) 10 CONTINUE RETURN END * SUBROUTINE MXVDIF ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VECTOR DIFFERENCE. * * PARAMETERS : * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RO Z(N) OUTPUT VECTOR WHERE Z:= X - Y. * SUBROUTINE MXVDIF(N,X,Y,Z) IMPLICIT NONE INTEGER N DOUBLE PRECISION X(*),Y(*),Z(*) INTEGER I DO 10 I = 1,N Z(I) = X(I) - Y(I) 10 CONTINUE RETURN END * SUBROUTINE MXVDIR ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VECTOR AUGMENTED BY THE SCALED VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI A SCALING FACTOR. * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RO Z(N) OUTPUT VECTOR WHERE Z:= Y + A*X. * SUBROUTINE MXVDIR(N,A,X,Y,Z) IMPLICIT NONE DOUBLE PRECISION A INTEGER N DOUBLE PRECISION X(*),Y(*),Z(*) INTEGER I DO 10 I = 1,N Z(I) = Y(I) + A*X(I) 10 CONTINUE RETURN END * FUNCTION MXVDOT ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DOT PRODUCT OF TWO VECTORS. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RR MXVDOT VALUE OF DOT PRODUCT MXVDOT=TRANS(X)*Y. * FUNCTION MXVDOT(N,X,Y) IMPLICIT NONE INTEGER N DOUBLE PRECISION X(*),Y(*),MXVDOT DOUBLE PRECISION TEMP INTEGER I TEMP = 0.0D0 DO 10 I = 1,N TEMP = TEMP + X(I)*Y(I) 10 CONTINUE MXVDOT = TEMP RETURN END * SUBROUTINE MXVINA ALL SYSTEMS 90/12/01 * PORTABILITY : ALL SYSTEMS * 90/12/01 LU : ORIGINAL VERSION * * PURPOSE : * ELEMENTS OF THE INTEGER VECTOR ARE REPLACED BY THEIR ABSOLUTE VALUES. * * PARAMETERS : * II N DIMENSION OF THE INTEGER VECTOR. * IU IX(N) INTEGER VECTOR WHICH IS UPDATED SO THAT IX(I):=ABS(IX(I)) * FOR ALL I. * SUBROUTINE MXVINA(N,IX) IMPLICIT NONE INTEGER N INTEGER IX(*) INTEGER I DO 10 I = 1,N IX(I) = ABS(IX(I)) IF (IX(I).GT.10) IX(I) = IX(I) - 10 10 CONTINUE RETURN END * SUBROUTINE MXVIND ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * CHANGE OF THE INTEGER VECTOR ELEMENT FOR THE CONSTRAINT ADDITION. * * PARAMETERS : * IU IX(N) INTEGER VECTOR. * II I INDEX OF THE CHANGED ELEMENT. * II JOB CHANGE SPECIFICATION. IS JOB.EQ.0 THEN IX(I)=10-IX(I). * SUBROUTINE MXVIND(IX,I,JOB) IMPLICIT NONE INTEGER IX(*),I,JOB IF (JOB.EQ.0) IX(I)=10-IX(I) RETURN END * SUBROUTINE MXVINV ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * CHANGE OF THE INTEGER VECTOR ELEMENT FOR THE CONSTRAINT ADDITION. * * PARAMETERS : * II N VECTOR DIMENSION. * IU IX(N) INTEGER VECTOR. * II I INDEX OF THE CHANGED ELEMENT. * II JOB CHANGE SPECIFICATION * SUBROUTINE MXVINV(IX,I,JOB) IMPLICIT NONE INTEGER I,JOB INTEGER IX(*) IF ((IX(I).EQ.3.OR.IX(I).EQ.5) .AND. JOB.LT.0) IX(I) = IX(I) + 1 IF ((IX(I).EQ.4.OR.IX(I).EQ.6) .AND. JOB.GT.0) IX(I) = IX(I) - 1 IX(I) = -IX(I) RETURN END * FUNCTION MXVMAX ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * L-INFINITY NORM OF A VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RR MXVMAX L-INFINITY NORM OF THE VECTOR X. * FUNCTION MXVMAX(N,X) IMPLICIT NONE INTEGER N DOUBLE PRECISION X(*),MXVMAX INTEGER I MXVMAX = 0.0D0 DO 10 I = 1,N MXVMAX = MAX(MXVMAX,ABS(X(I))) 10 CONTINUE RETURN END * SUBROUTINE MXVNEG ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * CHANGE THE SIGNS OF VECTOR ELEMENTS. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RO Y(N) OUTPUT VECTOR WHERE Y:= - X. * SUBROUTINE MXVNEG(N,X,Y) IMPLICIT NONE INTEGER N DOUBLE PRECISION X(*),Y(*) INTEGER I DO 10 I = 1,N Y(I) = -X(I) 10 CONTINUE RETURN END * FUNCTION MXVNOR ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * EUCLIDEAN NORM OF A VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RR MXVNOR EUCLIDEAN NORM OF X. * FUNCTION MXVNOR(N,X) IMPLICIT NONE INTEGER N DOUBLE PRECISION X(*),MXVNOR DOUBLE PRECISION POM,DEN INTEGER I DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D 0) DEN=ZERO DO 1 I=1,N DEN=MAX(DEN,ABS(X(I))) 1 CONTINUE POM=ZERO IF (DEN.GT.ZERO) THEN DO 2 I=1,N POM=POM+(X(I)/DEN)**2 2 CONTINUE ENDIF MXVNOR=DEN*SQRT(POM) RETURN END * SUBROUTINE MXVORT ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR PLANE ROTATION. * * PARAMETERS : * RU XK FIRST VALUE FOR PLANE ROTATION (XK IS TRANSFORMED TO * SQRT(XK**2+XL**2)) * RU XL SECOND VALUE FOR PLANE ROTATION (XL IS TRANSFORMED TO * ZERO) * RO CK DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * RO CL OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * IO IER INFORMATION ON THE TRANSFORMATION. IER=0-GENERAL PLANE * ROTATION. IER=1-PERMUTATION. IER=2-TRANSFORMATION SUPPRESSED. * SUBROUTINE MXVORT(XK,XL,CK,CL,IER) IMPLICIT NONE DOUBLE PRECISION CK,CL,XK,XL INTEGER IER DOUBLE PRECISION DEN,POM IF (XL.EQ.0.0D0) THEN IER = 2 ELSE IF (XK.EQ.0.0D0) THEN XK = XL XL = 0.0D0 IER = 1 ELSE IF (ABS(XK).GE.ABS(XL)) THEN POM = XL/XK DEN = SQRT(1.0D0+POM*POM) CK = 1.0D0/DEN CL = POM/DEN XK = XK*DEN ELSE POM = XK/XL DEN = SQRT(1.0D0+POM*POM) CL = 1.0D0/DEN CK = POM/DEN XK = XL*DEN END IF XL = 0.0D0 IER = 0 END IF RETURN END * SUBROUTINE MXVROT ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * PLANE ROTATION IS APPLIED TO TWO VALUES. * * PARAMETERS : * RU XK FIRST VALUE FOR PLANE ROTATION. * RU XL SECOND VALUE FOR PLANE ROTATION. * RI CK DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * RI CL OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * II IER INFORMATION ON THE TRANSFORMATION. IER=0-GENERAL PLANE * ROTATION. IER=1-PERMUTATION. IER=2-TRANSFORMATION SUPPRESSED. * SUBROUTINE MXVROT(XK,XL,CK,CL,IER) IMPLICIT NONE DOUBLE PRECISION CK,CL,XK,XL INTEGER IER DOUBLE PRECISION YK,YL IF (IER.EQ.0) THEN YK = XK YL = XL XK = CK*YK + CL*YL XL = CL*YK - CK*YL ELSE IF (IER.EQ.1) THEN YK = XK XK = XL XL = YK END IF RETURN END * SUBROUTINE MXVSAV ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DIFFERENCE OF TWO VECTORS RETURNED IN THE SUBSTRACTED ONE. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RU Y(N) UPDATE VECTOR WHERE Y:= X - Y. * SUBROUTINE MXVSAV(N,X,Y) IMPLICIT NONE INTEGER N DOUBLE PRECISION X(*),Y(*) DOUBLE PRECISION TEMP INTEGER I DO 10 I = 1,N TEMP = Y(I) Y(I) = X(I) - Y(I) X(I) = TEMP 10 CONTINUE RETURN END * SUBROUTINE MXVSCL ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SCALING OF A VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RI A SCALING FACTOR. * RO Y(N) OUTPUT VECTOR WHERE Y:= A*X. * SUBROUTINE MXVSCL(N,A,X,Y) IMPLICIT NONE DOUBLE PRECISION A INTEGER N DOUBLE PRECISION X(*),Y(*) INTEGER I DO 10 I = 1,N Y(I) = A*X(I) 10 CONTINUE RETURN END * SUBROUTINE MXVSET ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * A SCALAR IS SET TO ALL THE ELEMENTS OF A VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI A INITIAL VALUE. * RO X(N) OUTPUT VECTOR SUCH THAT X(I)=A FOR ALL I. * SUBROUTINE MXVSET(N,A,X) IMPLICIT NONE DOUBLE PRECISION A INTEGER N DOUBLE PRECISION X(*) INTEGER I DO 10 I = 1,N X(I) = A 10 CONTINUE RETURN END * SUBROUTINE MXVSUM ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SUM OF TWO VECTORS. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RO Z(N) OUTPUT VECTOR WHERE Z:= X + Y. * SUBROUTINE MXVSUM(N,X,Y,Z) IMPLICIT NONE INTEGER N DOUBLE PRECISION X(*),Y(*),Z(*) INTEGER I DO 10 I = 1,N Z(I) = X(I) + Y(I) 10 CONTINUE RETURN END