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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
|
*DECK MCGACA
SUBROUTINE MCGACA(LFORW,PACA,N,NG,NFIRST,NGEFF,M,LC,NGIND,NGINDV,
1 NCONV,KPSYS,JPMACR,NZON,IPERM,IM,MCU,JU,XIN,
2 LC0,IM0,MCU0,XOUT)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the product of the left-hand side ACA matrix in its multigroup
* form with a vector and apply group per group left preconditioner.
*
*Copyright:
* Copyright (C) 2002 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): R. Le Tellier
*
*Parameters: input
* LFORW flag set to .false. to transpose the coefficient matrix.
* PACA type of preconditioner to solve the ACA corrective system.
* N number of unknowns per group.
* NG total number of groups.
* NFIRST first group to proceed.
* NGEFF number of unconverged groups.
* M number of material mixtures.
* LC dimension of profiled matrices MCU and CQ.
* NGIND index of the groups to process.
* NGINDV index to pass from "NGEFF format" to "NG format".
* NCONV logical array of convergence status for each group (.TRUE.
* not converged).
* KPSYS pointer array for each group properties.
* JPMACR pointer to the macrolib LCM object ('GROUP' directory).
* NZON index-number of the mixture type assigned to each volume.
* IPERM permutation array for ACA.
* IM connection matrix.
* MCU connection matrix.
* JU used for ilu0 preconditioner.
* XIN undefined.
* LC0 used in ILU0-ACA acceleration.
* IM0 used in ILU0-ACA acceleration.
* MCU0 used in ILU0-ACA acceleration.
*
*Parameters: output
* XOUT product.
*
*-----------------------------------------------------------------------
*
USE GANLIB
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
TYPE(C_PTR) KPSYS(NGEFF),JPMACR
INTEGER PACA,N,NFIRST,NGEFF,NG,M,LC,NGIND(NGEFF),NGINDV(NG),
1 NZON(N),IPERM(N),IM(N+1),MCU(LC),JU(N),LC0,IM0(*),MCU0(*)
DOUBLE PRECISION XIN(N,NGEFF),XOUT(N,NGEFF)
LOGICAL LFORW,NCONV(NGEFF)
*----
* LOCAL VARIABLES
*----
TYPE(C_PTR) JPSYS,KPMACR
INTEGER I,J,II,IG,JG,JJ,JND,IBM
REAL, TARGET, SAVE, DIMENSION(1) :: DUMMY
*----
* ALLOCATABLE ARRAYS
*----
INTEGER, ALLOCATABLE, DIMENSION(:) :: NJJ,IJJ,IPOS
REAL, ALLOCATABLE, DIMENSION(:) :: XSCAT
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: TEMP
*
TYPE(C_PTR) DIAGF_PTR,CF_PTR,LUDF_PTR,LUCF_PTR,DIAGQ_PTR,CQ_PTR
REAL, POINTER, DIMENSION(:) :: DIAGF,CF,LUDF,LUCF,DIAGQ,CQ
*----
* INITIALIZE POINTERS
*----
LUDF=>DUMMY
LUCF=>DUMMY
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(NJJ(0:M),IJJ(0:M),IPOS(0:M),XSCAT(0:M*NG),TEMP(N))
*
DO II=NFIRST,NGEFF
IF(NCONV(II)) THEN
* compute temp=sum_{g' ne g} Sigma_s^{g<-g'} XIN^{g'}
IG=NGIND(II)
JPSYS=KPSYS(II)
CALL LCMGPD(JPSYS,'DIAGF$MCCG',DIAGF_PTR)
CALL LCMGPD(JPSYS,'CF$MCCG',CF_PTR)
CALL C_F_POINTER(DIAGF_PTR,DIAGF,(/ N /))
CALL C_F_POINTER(CF_PTR,CF,(/ LC /))
IF(PACA.GE.2) THEN
CALL LCMGPD(JPSYS,'ILUDF$MCCG',LUDF_PTR)
CALL C_F_POINTER(LUDF_PTR,LUDF,(/ N /))
IF(PACA.LT.4) THEN
CALL LCMGPD(JPSYS,'ILUCF$MCCG',LUCF_PTR)
CALL C_F_POINTER(LUCF_PTR,LUCF,(/ LC /))
ENDIF
ENDIF
KPMACR=LCMGIL(JPMACR,IG)
CALL LCMGPD(JPSYS,'DIAGQ$MCCG',DIAGQ_PTR)
CALL LCMGPD(JPSYS,'CQ$MCCG',CQ_PTR)
CALL C_F_POINTER(DIAGQ_PTR,DIAGQ,(/ N /))
CALL C_F_POINTER(CQ_PTR,CQ,(/ LC /))
CALL LCMGET(KPMACR,'NJJS00',NJJ(1))
CALL LCMGET(KPMACR,'IJJS00',IJJ(1))
CALL LCMGET(KPMACR,'IPOS00',IPOS(1))
CALL LCMGET(KPMACR,'SCAT00',XSCAT(1))
DO I=1,N
TEMP(I)=0.0D0
J=IPERM(I)
IBM=NZON(J)
IF(IBM.GT.0) THEN
JG=IJJ(IBM)
DO 10 JND=1,NJJ(IBM)
IF(JG.NE.IG) THEN
JJ=NGINDV(JG)
IF(JJ.GE.NFIRST) THEN
TEMP(I)=TEMP(I)+XSCAT(IPOS(IBM)+JND-1)*XIN(I,JJ)
ENDIF
ENDIF
JG=JG-1
10 CONTINUE
ENDIF
ENDDO
* compute E^{g}*temp
CALL MCGPRA(LFORW,1,PACA,.FALSE.,N,LC,IM,MCU,JU,DIAGQ,CQ,
1 LUDF,LUCF,DIAGF,TEMP(1),XOUT(1,II),LC0,IM0,MCU0,CF)
* compute D^{g}*XIN^{g}
CALL MCGPRA(LFORW,1,PACA,.FALSE.,N,LC,IM,MCU,JU,DIAGF,CF,
1 LUDF,LUCF,DIAGF,XIN(1,II),TEMP(1),LC0,IM0,MCU0,CF)
* temp=D^{g}*XIN^{g}-E^{g}*sum_{g' ne g} Sigma_s^{g<-g'} XIN^{g'}
DO I=1,N
TEMP(I)=TEMP(I)-XOUT(I,II)
ENDDO
* apply single-group preconditioner XOUT^{g}=P^{g}*temp
CALL MCGPRA(LFORW,2,PACA,.TRUE.,N,LC,IM,MCU,JU,DIAGF,CF,
1 LUDF,LUCF,DIAGF,XOUT(1,II),TEMP(1),LC0,IM0,MCU0,CF)
ENDIF
ENDDO
*----
* SCRATCH STORAGE DEALLOCATION
*----
DEALLOCATE(TEMP,XSCAT,IPOS,IJJ,NJJ)
RETURN
END
|
