53 int *ipc,
double *rpc,
54 double *ac,
double *cc,
double *fc,
55 double *x,
double *w1,
double *w2,
double *r,
56 int *
itmax,
int *iters,
57 double *
errtol,
double *omega,
58 int *iresid,
int *iadjoint) {
62 MAT2(ac, *
nx * *
ny * *
nz, 1);
65 numdia = VAT(ipc, 11);
69 RAT2(ac, 1,1), cc, fc,
70 RAT2(ac, 1,2), RAT2(ac, 1,3), RAT2(ac, 1,4),
73 }
else if (numdia == 27) {
76 RAT2(ac, 1, 1), cc, fc,
77 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
78 RAT2(ac, 1, 5), RAT2(ac, 1, 6),
79 RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
80 RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
84 Vnm_print(2,
"GSRB: invalid stencil type given...\n");
90 VPUBLIC
void Vgsrb7x(
int *
nx,
int *
ny,
int *
nz,
91 int *ipc,
double *rpc,
92 double *oC,
double *cc,
double *fc,
93 double *oE,
double *oN,
double *uC,
94 double *x,
double *w1,
double *w2,
double *r,
95 int *
itmax,
int *iters,
96 double *
errtol,
double *omega,
97 int *iresid,
int *iadjoint) {
113 for (*iters=1; *iters<=*
itmax; (*iters)++) {
116 #pragma omp parallel for private(i, j, k, ioff) 117 for (k=2; k<=*
nz-1; k++) {
118 for (j=2; j<=*
ny-1; j++) {
119 ioff = (1 - *iadjoint) * ( (j + k + 2) % 2)
120 + ( *iadjoint) * (1 - (j + k + 2) % 2);
121 for (i=2+ioff; i<=*
nx-1; i+=2) {
124 + VAT3(oN, i, j, k) * VAT3(x, i, j+1, k)
125 + VAT3(oN, i, j-1, k) * VAT3(x, i, j-1, k)
126 + VAT3(oE, i, j, k) * VAT3(x, i+1, j, k)
127 + VAT3(oE, i-1, j, k) * VAT3(x, i-1, j, k)
128 + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
129 + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
130 ) / (VAT3(oC, i, j, k) + VAT3(cc, i, j, k));
136 #pragma omp parallel for private(i, j, k, ioff) 137 for (k=2; k<=*
nz-1; k++) {
138 for (j=2; j<=*
ny-1; j++) {
139 ioff = ( *iadjoint) * ( (j + k + 2) % 2 )
140 + (1 - *iadjoint) * (1 - (j + k + 2) % 2 );
141 for (i=2+ioff;i<=*
nx-1; i+=2) {
144 + VAT3(oN, i, j, k) * VAT3(x, i,j+1, k)
145 + VAT3(oN, i, j-1, k) * VAT3(x, i,j-1, k)
146 + VAT3(oE, i, j, k) * VAT3(x, i+1, j, k)
147 + VAT3(oE, i-1, j, k) * VAT3(x, i-1, j, k)
148 + VAT3( uC, i, j, k-1) * VAT3(x, i, j,k-1)
149 + VAT3( uC, i, j, k) * VAT3(x, i, j,k+1)
150 ) / (VAT3(oC, i, j, k) + VAT3(cc, i, j, k));
157 Vmresid7_1s(
nx,
ny,
nz, ipc, rpc, oC, cc, fc, oE, oN, uC, x, r);
162 VPUBLIC
void Vgsrb27x(
int *
nx,
int *
ny,
int *
nz,
163 int *ipc,
double *rpc,
164 double *oC,
double *cc,
double *fc,
165 double *oE,
double *oN,
double *uC,
double *oNE,
double *oNW,
166 double *uE,
double *uW,
double *uN,
double *uS,
167 double *uNE,
double *uNW,
double *uSE,
double *uSW,
168 double *x,
double *w1,
double *w2,
double *r,
169 int *
itmax,
int *iters,
170 double *
errtol,
double *omega,
171 int *iresid,
int *iadjoint) {
179 double tmpO, tmpU, tmpD;
217 i1 = (1-*iadjoint) * 2 + *iadjoint * (*
nx-1);
218 i2 = *iadjoint * 2 + (1-*iadjoint) * (*
nx-1);
219 j1 = (1-*iadjoint) * 2 + *iadjoint * (*
ny-1);
220 j2 = *iadjoint * 2 + (1-*iadjoint) * (*
ny-1);
221 k1 = (1-*iadjoint) * 2 + *iadjoint * (*
nz-1);
222 k2 = *iadjoint * 2 + (1-*iadjoint) * (*
nz-1);
223 istep = *iadjoint*(-1) + (1-*iadjoint)*(1);
225 for (*iters=1; *iters<=*
itmax; (*iters)++) {
228 for (k=2; k<=*
nz-1; k++) {
230 for (j=2; j<=*
ny-1; j++) {
232 ioff = (1 - *iadjoint) * ( (j + k + 2) % 2)
233 + ( *iadjoint) * (1 - (j + k + 2) % 2);
235 for (i=2+ioff; i<=*
nx-1; i+=2) {
238 + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
239 + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
240 + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
241 + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
242 + VAT3( oNE, i, j, k) * VAT3(x, i+1, j+1, k)
243 + VAT3( oNW, i, j, k) * VAT3(x, i-1, j+1, k)
244 + VAT3( oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
245 + VAT3( oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
248 + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
249 + VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
250 + VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
251 + VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
252 + VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
253 + VAT3( uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
254 + VAT3( uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
255 + VAT3( uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
256 + VAT3( uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
259 + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
260 + VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
261 + VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
262 + VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
263 + VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
264 + VAT3( uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
265 + VAT3( uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
266 + VAT3( uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
267 + VAT3( uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
269 VAT3(x, i,j,k) = (VAT3(fc, i, j, k) + (tmpO + tmpU + tmpD))
270 / (VAT3(oC, i, j, k) + VAT3(cc, i, j, k));
277 for (k=2; k<=*
nz-1; k++) {
279 for (j=2; j<=*
ny-1; j++) {
281 ioff = ( *iadjoint) * ( (j + k + 2) % 2)
282 + (1 - *iadjoint) * (1 - (j + k + 2) % 2);
284 for (i=2+ioff; i<=*
nx-1; i+=2) {
287 + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
288 + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
289 + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
290 + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
291 + VAT3( oNE, i, j, k) * VAT3(x, i+1, j+1, k)
292 + VAT3( oNW, i, j, k) * VAT3(x, i-1, j+1, k)
293 + VAT3( oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
294 + VAT3( oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
297 + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
298 + VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
299 + VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
300 + VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
301 + VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
302 + VAT3( uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
303 + VAT3( uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
304 + VAT3( uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
305 + VAT3( uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
308 + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
309 + VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
310 + VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
311 + VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
312 + VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
313 + VAT3( uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
314 + VAT3( uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
315 + VAT3( uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
316 + VAT3( uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
318 VAT3(x, i,j,k) = (VAT3(fc, i, j, k) + (tmpO + tmpU + tmpD))
319 / (VAT3(oC, i, j, k) + VAT3(cc, i, j, k));
VPUBLIC void Vgsrb(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *x, double *w1, double *w2, double *r, int *itmax, int *iters, double *errtol, double *omega, int *iresid, int *iadjoint)
Guass-Seidel solver.