1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
|
*DECK GRA001
SUBROUTINE GRA001(IPFLX,IPGPT,NVAR,NCST,DERIV)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute new gradients of system characteristics (part 2).
*
*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
* IPFLX pointer of a LCM object containing a set solutions of
* fixed-source eigenvalue problems.
* IPGPT pointer of a LCM object containing a set of fixed sources.
* NVAR number of control variables.
* NCST number of constraints with indirect effects (can be zero).
*
*Parameters: output
* DERIV gradient matrix.
*
*-----------------------------------------------------------------------
*
USE GANLIB
*----
* SUBROUTINE ARGUMENTS
*----
TYPE(C_PTR) IPFLX,IPGPT
INTEGER NVAR,NCST
DOUBLE PRECISION DERIV(NVAR,NCST+1)
*----
* LOCAL VARIABLES
*----
TYPE(C_PTR) JPFLX,KPFLX,JPGPT,KPGPT
PARAMETER (NSTATE=40)
INTEGER ISTATE(NSTATE)
DOUBLE PRECISION SUM
REAL, ALLOCATABLE, DIMENSION(:) :: DFLUX,SOUR
*
CALL LCMGET(IPFLX,'STATE-VECTOR',ISTATE)
NGRP=ISTATE(1)
NUN=ISTATE(2)
ITYPE=ISTATE(3)
NGPT=ISTATE(5)
CALL LCMGET(IPGPT,'STATE-VECTOR',ISTATE)
IF(ISTATE(1).NE.NGRP) CALL XABORT('GRA001: INVALID NGRP')
IF(ISTATE(2).NE.NUN) CALL XABORT('GRA001: INVALID NUN')
ND=ISTATE(3)
NA=ISTATE(4)
*
ALLOCATE(SOUR(NUN),DFLUX(NUN))
DERIV(:NVAR,:NCST+1)=0.0D0
IF(ITYPE.EQ.100) THEN
* EXPLICIT APPROACH
IF(NVAR.NE.NGPT) CALL XABORT('GRA001: INVALID NGPT(1)')
IF(NCST+1.NE.NA) CALL XABORT('GRA001: INVALID NA(1)')
JPFLX=LCMGID(IPFLX,'DFLUX')
JPGPT=LCMGID(IPGPT,'ASOUR')
DO 25 IVAR=1,NVAR
KPFLX=LCMGIL(JPFLX,IVAR)
DO 20 ICST=1,NCST+1
KPGPT=LCMGIL(JPGPT,ICST)
SUM=0.0D0
DO 15 IGR=1,NGRP
CALL LCMGDL(KPGPT,IGR,SOUR)
CALL LCMGDL(KPFLX,IGR,DFLUX)
DO 10 IUN=1,NUN
SUM=SUM+SOUR(IUN)*DFLUX(IUN)
10 CONTINUE
15 CONTINUE
DERIV(IVAR,ICST)=SUM
20 CONTINUE
25 CONTINUE
ELSE IF(ITYPE.EQ.1000) THEN
* IMPLICIT APPROACH
IF(NVAR.NE.ND) CALL XABORT('GRA001: INVALID ND(2)')
IF(NCST+1.NE.NGPT) CALL XABORT('GRA001: INVALID NGPT(2)')
JPFLX=LCMGID(IPFLX,'ADFLUX')
JPGPT=LCMGID(IPGPT,'DSOUR')
DO 45 ICST=1,NCST+1
KPFLX=LCMGIL(JPFLX,ICST)
DO 40 IVAR=1,NVAR
KPGPT=LCMGIL(JPGPT,IVAR)
SUM=0.0D0
DO 35 IGR=1,NGRP
CALL LCMGDL(KPGPT,IGR,SOUR)
CALL LCMGDL(KPFLX,IGR,DFLUX)
DO 30 IUN=1,NUN
SUM=SUM+SOUR(IUN)*DFLUX(IUN)
30 CONTINUE
35 CONTINUE
DERIV(IVAR,ICST)=SUM
40 CONTINUE
45 CONTINUE
ELSE
CALL XABORT('GRA001: INVALID FLUX OBJECT')
ENDIF
DEALLOCATE(DFLUX,SOUR)
RETURN
END
|
