*DECK GRAD SUBROUTINE GRAD(NENTRY,HENTRY,IENTRY,JENTRY,KENTRY) * *----------------------------------------------------------------------- * *Purpose: * Compute gradients of system characteristics. * *Copyright: * Copyright (C) 2012 Ecole Polytechnique de Montreal * This library is 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 * *Author(s): * A. Hebert * *Parameters: input * NENTRY number of data structures transfered to this module. * HENTRY name of the data structures. * IENTRY data structure type where: * IENTRY=1 for LCM memory object; * IENTRY=2 for XSM file; * IENTRY=3 for sequential binary file; * IENTRY=4 for sequential ASCII file. * JENTRY access permission for the data structure where: * JENTRY=0 for a data structure in creation mode; * JENTRY=1 for a data structure in modifications mode; * JENTRY=2 for a data structure in read-only mode. * KENTRY data structure pointer. * *Comments: * The calling specifications are: * OPTIM := GRAD: [ OPTIM ] DFLUX GPT :: (grad\_data) ; * where * OPTIM : name of the \emph{optimize} object (L\_OPTIMIZE signature) * containing the optimization informations. Object OPTIM must appear on the * RHS to be able to updated the previous values. * DFLUX : name of the \emph{flux} object (L\_FLUX signature) containing a set * of solutions of fixed-source eigenvalue problems. * GPT : name of the \emph{gpt} object (L\_GPT signature) containing a set * of direct or adjoint sources. * (grad\_data) : structure containing the data to the module GRAD:. * *----------------------------------------------------------------------- * USE GANLIB IMPLICIT DOUBLE PRECISION (A-H,O-Z) *---- * SUBROUTINE ARGUMENTS *---- INTEGER NENTRY,IENTRY(NENTRY),JENTRY(NENTRY) TYPE(C_PTR) KENTRY(NENTRY) CHARACTER HENTRY(NENTRY)*12 *---- * LOCAL VARIABLES *---- PARAMETER (NSTATE=40) TYPE(C_PTR) IPFLX,IPGPT,IPGRAD CHARACTER HSIGN*12,TEXT4*4,TEXT12*12,TEXT16*16 INTEGER ISTATE(NSTATE) REAL FLOTT DOUBLE PRECISION DFLOTT,SR,EPS1,EPS2,EPS3,EPS4 DOUBLE PRECISION OPTPRR(NSTATE) *---- * ALLOCATABLE ARRAYS *---- INTEGER, ALLOCATABLE, DIMENSION(:) :: IREL DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: VARV,CSTV,RHS, 1 DERIV,DERIV0 *---- * PARAMETER VALIDATION. *---- IF((IENTRY(1).NE.1).AND.(IENTRY(1).NE.2)) CALL XABORT('GRAD: LCM' 1 //' OBJECT EXPECTED AT LHS.') IF(JENTRY(1).EQ.2) CALL XABORT('GRAD: OPTIMIZE ENTRY IN CREATE O' 1 //'R MODIFICATION MODE EXPECTED.') DO I=2,NENTRY TEXT12=HENTRY(I) IF((JENTRY(I).NE.2).OR.((IENTRY(I).NE.1).AND.(IENTRY(I).NE.2))) 1 CALL XABORT('GRAD: LCM OBJECT IN READ-ONLY MODE EXPECTED AT R' 2 //'HS ('//TEXT12//').') ENDDO IPGRAD=KENTRY(1) IPFLX=C_NULL_PTR IPGPT=C_NULL_PTR *---- * RECOVER THE ACTUAL FLUX SOLUTION AND CORRESPONDING TRACKING. *---- NVAR0=0 NCST0=0 ITYPE=0 IF(NENTRY.EQ.3) THEN CALL LCMGTC(KENTRY(2),'SIGNATURE',12,HSIGN) IF(HSIGN.NE.'L_FLUX') THEN TEXT12=HENTRY(2) CALL XABORT('GRAD: SIGNATURE OF '//TEXT12//' IS '//HSIGN// 1 '. L_FLUX EXPECTED.') ENDIF IPFLX=KENTRY(2) CALL LCMGET(IPFLX,'STATE-VECTOR',ISTATE) ITYPE=ISTATE(3) NGPT=ISTATE(5) IF(NGPT.EQ.0) CALL XABORT('GRAD: MISSING FIXED-SOURCE EIGENVA' 1 //'LUE SOLUTION') IPGPT=KENTRY(3) CALL LCMGTC(IPGPT,'SIGNATURE',12,HSIGN) IF(HSIGN.NE.'L_SOURCE') THEN TEXT12=HENTRY(3) CALL XABORT('GRAD: SIGNATURE OF '//TEXT12//' IS '//HSIGN// 1 '. L_SOURCE EXPECTED.') ENDIF CALL LCMGET(IPGPT,'STATE-VECTOR',ISTATE) ND=ISTATE(3) NA=ISTATE(4) *---- * COMPUTE THE NUMBER OF CONSTRAINTS AND OF CONTROL VARIABLES *---- IF(ITYPE.EQ.100) THEN NVAR0=NGPT NCST0=NA-1 ELSE IF(ITYPE.EQ.1000) THEN NVAR0=ND NCST0=NGPT-1 ELSE CALL XABORT('GRAD: INVALID FLUX OBJECT') ENDIF ENDIF *---- * READ INPUT PARAMETERS *---- IPRINT=1 IOPT=1 ICONV=0 IEXT=0 IEDSTP=2 IHESS=0 ISEARC=0 IMETH=2 ISTEP=0 JCONV=0 SR=1.0D0 EPS1=0.1D0 EPS2=1.0D-4 EPS3=1.0D-4 EPS4=1.0D-4 IF(JENTRY(1).EQ.0) THEN HSIGN='L_OPTIMIZE' CALL LCMPTC(IPGRAD,'SIGNATURE',12,HSIGN) ELSE IF (JENTRY(1).EQ.1) THEN CALL LCMGTC(IPGRAD,'SIGNATURE',12,HSIGN) IF(HSIGN.NE.'L_OPTIMIZE') THEN TEXT12=HENTRY(3) CALL XABORT('GRAD: SIGNATURE OF '//TEXT12//' IS '//HSIGN// 1 '. L_OPTIMIZE EXPECTED.') ENDIF CALL LCMGET(IPGRAD,'STATE-VECTOR',ISTATE) NVAR=ISTATE(1) NCST=ISTATE(2) IOPT=ISTATE(3) ICONV=ISTATE(4) IEXT=ISTATE(5) IEDSTP=ISTATE(6) IHESS=ISTATE(7) ISEARC=ISTATE(8) IMETH=ISTATE(9) MAXEXT=ISTATE(12) NSTART=ISTATE(13) CALL LCMGET(IPGRAD,'OPT-PARAM-R',OPTPRR) SR=OPTPRR(1) EPS1=OPTPRR(2) EPS2=OPTPRR(3) EPS3=OPTPRR(4) EPS4=OPTPRR(5) ENDIF 10 CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT) IF(INDIC.EQ.10) GO TO 20 IF(INDIC.NE.3) CALL XABORT('GRAD: CHARACTER DATA EXPECTED.') IF(TEXT12(:4).EQ.'EDIT') THEN CALL REDGET(INDIC,IPRINT,FLOTT,TEXT12,DFLOTT) IF(INDIC.NE.1) CALL XABORT('GRAD: INTEGER DATA EXPECTED FOR IP' 1 //'RINT.') ELSE IF(TEXT12(:8).EQ.'MINIMIZE') THEN IOPT=1 ELSE IF(TEXT12(:8).EQ.'MAXIMIZE') THEN IOPT=-1 ELSE IF(TEXT12.EQ.'OUT-STEP-LIM') THEN CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT) IF(INDIC.NE.2) CALL XABORT('GRAD: REAL DATA EXPECTED(1)') SR=FLOTT ELSE IF((TEXT12(:9).EQ.'VAR-VALUE').OR. 1 (TEXT12(:10).EQ.'VAR-WEIGHT')) THEN ALLOCATE(VARV(NVAR)) DO IVAR=1,NVAR CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT) IF(INDIC.NE.2) CALL XABORT('GRAD: REAL DATA EXPECTED(2)') VARV(IVAR)=FLOTT ENDDO CALL LCMPUT(IPGRAD,TEXT12,NVAR,4,VARV) DEALLOCATE(VARV) ELSE IF((TEXT12(:11).EQ.'VAR-VAL-MIN').OR. 1 (TEXT12(:11).EQ.'VAR-VAL-MAX')) THEN ALLOCATE(VARV(NVAR)) CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT) IF(INDIC.EQ.2) THEN VARV=FLOTT DO IVAR=2,NVAR CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT) IF(INDIC.NE.2) CALL XABORT('GRAD: REAL DATA EXPECTED(3)') VARV(IVAR)=FLOTT ENDDO ELSE IF((INDIC.EQ.3).AND.(TEXT4.EQ.'ALL')) THEN CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT) IF(INDIC.NE.2) CALL XABORT('GRAD: REAL DATA EXPECTED(4)') DO IVAR=1,NVAR VARV(IVAR)=FLOTT ENDDO ELSE CALL XABORT('GRAD: REAL DATA OR ALL KEYWORD EXPECTED') ENDIF CALL LCMPUT(IPGRAD,TEXT12,NVAR,4,VARV) DEALLOCATE(VARV) ELSE IF(TEXT12.EQ.'FOBJ-CST-VAL') THEN ALLOCATE(CSTV(NCST+1)) DO ICST=1,NCST+1 CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT) IF(INDIC.NE.2) CALL XABORT('GRAD: REAL DATA EXPECTED(5)') CSTV(ICST)=FLOTT ENDDO CALL LCMPUT(IPGRAD,'FOBJ-CST-VAL',NCST+1,4,CSTV) OBJNEW=CSTV(1) DEALLOCATE(CSTV) ELSE IF(TEXT12(:8).EQ.'CST-TYPE') THEN IF(NCST.EQ.0) CALL XABORT('GRAD: CST-TYPE KEYWORD FORBIDDEN') ALLOCATE(IREL(NCST)) DO ICST=1,NCST CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT) IF(INDIC.NE.1) THEN CALL XABORT('GRAD: INTEGER DATA EXPECTED') ELSE IF((NITMA.LT.-1).OR.(NITMA.GT.1)) THEN CALL XABORT('GRAD: -1, 0 or 1 EXPECTED') ENDIF IREL(ICST)=NITMA ENDDO CALL LCMPUT(IPGRAD,'CST-TYPE',NCST,1,IREL) DEALLOCATE(IREL) ELSE IF(TEXT12(:7).EQ.'CST-OBJ') THEN IF(NCST.EQ.0) CALL XABORT('GRAD: CST-OBJ KEYWORD FORBIDDEN') CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT) IF(INDIC.EQ.2) THEN ALLOCATE(RHS(NCST)) RHS(1)=FLOTT DO ICST=2,NCST CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT) IF(INDIC.NE.2) CALL XABORT('GRAD: REAL DATA EXPECTED(6)') RHS(ICST)=FLOTT ENDDO CALL LCMPUT(IPGRAD,'CST-OBJ',NCST,4,RHS) DEALLOCATE(RHS) ELSE CALL XABORT('GRAD: REAL DATA OR KEEP KEYWORD EXPECTED') ENDIF ELSE IF(TEXT12(:10).EQ.'CST-WEIGHT') THEN IF(NCST.EQ.0) CALL XABORT('GRAD: CST-WEIGHT KEYWORD FORBIDDEN') ALLOCATE(RHS(NCST)) DO ICST=1,NCST CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DFLOTT) IF(INDIC.NE.2) CALL XABORT('GRAD: REAL DATA EXPECTED(7)') RHS(ICST)=FLOTT ENDDO CALL LCMPUT(IPGRAD,'CST-WEIGHT',NCST,4,RHS) DEALLOCATE(RHS) ELSE IF(TEXT12(:1).EQ.';') THEN GO TO 20 ELSE CALL XABORT('GRAD: '//TEXT12//' IS AN INVALID KEYWORD') ENDIF GO TO 10 *---- * CALCULATION OF THE NEW GRADIENT *---- 20 IF(IPRINT.GT.0) THEN IF(ITYPE.EQ.100) THEN WRITE(6,'(/25H GRAD: EXPLICIT APPROACH.)') ELSE IF(ITYPE.EQ.1000) THEN WRITE(6,'(/25H GRAD: IMPLICIT APPROACH.)') ENDIF ENDIF ALLOCATE(DERIV(NVAR*(NCST+1))) DERIV(:NVAR*(NCST+1))=0.0D0 IF(C_ASSOCIATED(IPFLX).AND.C_ASSOCIATED(IPGPT)) THEN IF(NVAR0.NE.NVAR) CALL XABORT('GRAD: INCONSISTENT NVAR.') IF(NCST0.GT.NCST) CALL XABORT('GRAD: INCONSISTENT NCST.') * ------------------------------------------ CALL GRA001(IPFLX,IPGPT,NVAR0,NCST0,DERIV) * ------------------------------------------ ENDIF CALL LCMLEN(IPGRAD,'GRADIENT-DIR',LENGTH,ITYLCM) IF(LENGTH.EQ.NVAR*(NCST+1)) THEN ALLOCATE(DERIV0(NVAR*(NCST+1))) CALL LCMGET(IPGRAD,'GRADIENT-DIR',DERIV0) DO I=1,NVAR*(NCST+1) DERIV(I)=DERIV(I)+DERIV0(I) ENDDO DEALLOCATE(DERIV0) ENDIF CALL LCMPUT(IPGRAD,'GRADIENT',NVAR*(NCST+1),4,DERIV) DEALLOCATE(DERIV) *---- * PRINT INFORMATION *---- IF(IPRINT.GT.0) THEN WRITE(6,'(/31H GRAD: INFORMATION AT ITERATION,I5)') IEXT+1 CALL LCMLEN(IPGRAD,'VAR-VALUE',ILONG,ITYLCM) IF(ILONG.GT.0) THEN ALLOCATE(VARV(NVAR)) CALL LCMGET(IPGRAD,'VAR-VALUE',VARV) WRITE(6,100) 'CONTROL VARIABLES:',(VARV(IVAR),IVAR=1,NVAR) DEALLOCATE(VARV) ENDIF IF(IPRINT.GT.1) THEN ALLOCATE(DERIV(NVAR*(NCST+1))) CALL LCMGET(IPGRAD,'GRADIENT',DERIV) WRITE(6,'(/29H GRADIENTS-------------------)') WRITE(6,100) 'OBJECTIVE FUNCTION:',(DERIV(IVAR),IVAR=1,NVAR) IF(IPRINT.GT.2) THEN DO 60 ICST=1,NCST WRITE(TEXT16,'(10HCONSTRAINT,I4,1H:)') ICST WRITE(6,100) TEXT16,(DERIV(ICST*NVAR+IVAR),IVAR=1,NVAR) 60 CONTINUE ENDIF DEALLOCATE(DERIV) ENDIF ENDIF *---- * SAVE THE STATE VECTORS *---- ISTATE(1)=NVAR ISTATE(2)=NCST ISTATE(3)=IOPT ISTATE(4)=ICONV ISTATE(5)=IEXT ISTATE(6)=IEDSTP ISTATE(7)=IHESS ISTATE(8)=ISEARC ISTATE(9)=IMETH ISTATE(10)=ISTEP ISTATE(11)=JCONV ISTATE(14)=0 IF(IPRINT.GT.0) WRITE(6,110) (ISTATE(I),I=1,9) CALL LCMPUT(IPGRAD,'STATE-VECTOR',NSTATE,1,ISTATE) OPTPRR(:NSTATE)=0.0D0 OPTPRR(1)=SR OPTPRR(2)=EPS1 OPTPRR(3)=EPS2 OPTPRR(4)=EPS3 OPTPRR(5)=EPS4 IF(IPRINT.GT.0) WRITE(6,120) (OPTPRR(I),I=1,5) CALL LCMPUT(IPGRAD,'OPT-PARAM-R',NSTATE,4,OPTPRR) IF(IPRINT.GT.2) CALL LCMLIB(IPGRAD) RETURN * 100 FORMAT(1X,A28,1P,8E12.4/(29X,8E12.4)) 110 FORMAT(/8H OPTIONS/8H -------/ 1 7H NVAR ,I8,32H (NUMBER OF CONTROL VARIABLES)/ 2 7H NCST ,I8,26H (NUMBER OF CONSTRAINTS)/ 3 7H IOPT ,I8,37H (=1/-1: MINIMIZATION/MAXIMIZATION)/ 4 7H ICONV ,I8,43H (=0/1: EXTERNAL NOT CONVERGED/CONVERGED)/ 5 7H IEXT ,I8,32H (INDEX OF EXTERNAL ITERATION)/ 6 7H IEDSTP,I8,43H (=1/2: HALF REDUCTION/PARABOLIC FORMULA)/ 7 7H IHESS ,I8,29H (=0/1/2: STEEPEST/CG/BFGS)/ 8 7H ISEARC,I8,35H (=0/1/2: NO SEARCH/OPTEX/NEWTON)/ 9 7H IMETH ,I8,42H (=1/2/3: SIMPLEX-LEMKE/LEMKE-LEMKE/MAP)) 120 FORMAT(/ 1 12H REAL PARAM:,1P/12H -----------/ 2 7H SR ,D12.4,39H (RADIUS OF THE QUADRATIC CONSTRAINT)/ 3 7H EPS1 ,D12.4,13H (NOT USED)/ 4 7H EPS2 ,D12.4,31H (EXTERNAL CONVERGENCE LIMIT)/ 5 7H EPS3 ,D12.4,31H (INTERNAL CONVERGENCE LIMIT)/ 6 7H EPS4 ,D12.4,43H (QUADRATIC CONSTRAINT CONVERGENCE LIMIT)) END