Initial commit from Polytechnique Montreal

1 files changed, 382 insertions, 0 deletions

diff --git a/Donjon/src/GRAD.f b/Donjon/src/GRAD.f
new file mode 100644
index 0000000..26ca070
--- /dev/null
+++ b/Donjon/src/GRAD.f
@@ -0,0 +1,382 @@
+*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