1 /* ./src_f77/gsd.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__9 = 9;
11 static integer c__1 = 1;
12
13 /* * /////////////////////////////////////////////////////////////////////////// */
14 /* * @file gsd.f */
15 /* * @author Michael Holst */
16 /* * @brief Gauss-Seidel iteration. */
17 /* * @version $Id: gsd.c,v 1.3 2010/08/12 05:45:21 fetk Exp $ */
18 /* * @attention */
19 /* * @verbatim */
20 /* * */
21 /* * PMG -- Parallel algebraic MultiGrid */
22 /* * Copyright (C) 1994-- Michael Holst. */
23 /* * */
24 /* * Michael Holst <mholst@math.ucsd.edu> */
25 /* * University of California, San Diego */
26 /* * Department of Mathematics, 5739 AP&M */
27 /* * 9500 Gilman Drive, Dept. 0112 */
28 /* * La Jolla, CA 92093-0112 USA */
29 /* * http://math.ucsd.edu/~mholst */
30 /* * */
31 /* * This file is part of PMG. */
32 /* * */
33 /* * PMG is free software; you can redistribute it and/or modify */
34 /* * it under the terms of the GNU General Public License as published by */
35 /* * the Free Software Foundation; either version 2 of the License, or */
36 /* * (at your option) any later version. */
37 /* * */
38 /* * PMG is distributed in the hope that it will be useful, */
39 /* * but WITHOUT ANY WARRANTY; without even the implied warranty of */
40 /* * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
41 /* * GNU General Public License for more details. */
42 /* * */
43 /* * You should have received a copy of the GNU General Public License */
44 /* * along with PMG; if not, write to the Free Software */
45 /* * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
46 /* * */
47 /* * Linking PMG statically or dynamically with other modules is making a */
48 /* * combined work based on PMG. Thus, the terms and conditions of the GNU */
49 /* * General Public License cover the whole combination. */
50 /* * */
51 /* * SPECIAL GPL EXCEPTION */
52 /* * In addition, as a special exception, the copyright holders of PMG */
53 /* * give you permission to combine the PMG program with free software */
54 /* * programs and libraries that are released under the GNU LGPL or with */
55 /* * code included in releases of ISIM, PMV, PyMOL, SMOL, VMD, and Vision. */
56 /* * Such combined software may be linked with PMG and redistributed together */
57 /* * in original or modified form as mere aggregation without requirement that */
58 /* * the entire work be under the scope of the GNU General Public License. */
59 /* * This special exception permission is also extended to any software listed */
60 /* * in the SPECIAL GPL EXCEPTION clauses by the FEtk and APBS libraries. */
61 /* * */
62 /* * Note that people who make modified versions of PMG are not obligated */
63 /* * to grant this special exception for their modified versions; it is */
64 /* * their choice whether to do so. The GNU General Public License gives */
65 /* * permission to release a modified version without this exception; this */
66 /* * exception also makes it possible to release a modified version which */
67 /* * carries forward this exception. */
68 /* * */
69 /* * @endverbatim */
70 /* * /////////////////////////////////////////////////////////////////////////// */
gsrb_(integer * nx,integer * ny,integer * nz,integer * ipc,doublereal * rpc,doublereal * ac,doublereal * cc,doublereal * fc,doublereal * x,doublereal * w1,doublereal * w2,doublereal * r__,integer * itmax,integer * iters,doublereal * errtol,doublereal * omega,integer * iresid,integer * iadjoint)71 /* Subroutine */ int gsrb_(integer *nx, integer *ny, integer *nz, integer *
72 ipc, doublereal *rpc, doublereal *ac, doublereal *cc, doublereal *fc,
73 doublereal *x, doublereal *w1, doublereal *w2, doublereal *r__,
74 integer *itmax, integer *iters, doublereal *errtol, doublereal *omega,
75 integer *iresid, integer *iadjoint)
76 {
77 /* System generated locals */
78 integer ac_dim1, ac_offset, cc_dim1, cc_dim2, cc_offset, fc_dim1, fc_dim2,
79 fc_offset, x_dim1, x_dim2, x_offset, w1_dim1, w1_dim2, w1_offset,
80 w2_dim1, w2_dim2, w2_offset, r_dim1, r_dim2, r_offset;
81
82 /* Builtin functions */
83 integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
84 e_wsle(void);
85
86 /* Local variables */
87 extern /* Subroutine */ int gsrb7_(integer *, integer *, integer *,
88 integer *, doublereal *, doublereal *, doublereal *, doublereal *,
89 doublereal *, doublereal *, doublereal *, doublereal *,
90 doublereal *, doublereal *, doublereal *, integer *, integer *,
91 doublereal *, doublereal *, integer *, integer *), gsrb27_(
92 integer *, integer *, integer *, integer *, doublereal *,
93 doublereal *, doublereal *, doublereal *, doublereal *,
94 doublereal *, doublereal *, doublereal *, doublereal *,
95 doublereal *, doublereal *, doublereal *, doublereal *,
96 doublereal *, doublereal *, doublereal *, doublereal *,
97 doublereal *, doublereal *, doublereal *, doublereal *, integer *,
98 integer *, doublereal *, doublereal *, integer *, integer *);
99 static integer numdia;
100
101 /* Fortran I/O blocks */
102 static cilist io___2 = { 0, 6, 0, 0, 0 };
103
104
105 /* * ********************************************************************* */
106 /* * purpose: */
107 /* * */
108 /* * call the fast diagonal iterative method. */
109 /* * */
110 /* * author: michael holst */
111 /* * ********************************************************************* */
112 /* * */
113 /* mdir 0 0 */
114 /* * */
115 /* * *** do in one step *** */
116 /* Parameter adjustments */
117 r_dim1 = *nx;
118 r_dim2 = *ny;
119 r_offset = 1 + r_dim1 * (1 + r_dim2);
120 r__ -= r_offset;
121 w2_dim1 = *nx;
122 w2_dim2 = *ny;
123 w2_offset = 1 + w2_dim1 * (1 + w2_dim2);
124 w2 -= w2_offset;
125 w1_dim1 = *nx;
126 w1_dim2 = *ny;
127 w1_offset = 1 + w1_dim1 * (1 + w1_dim2);
128 w1 -= w1_offset;
129 x_dim1 = *nx;
130 x_dim2 = *ny;
131 x_offset = 1 + x_dim1 * (1 + x_dim2);
132 x -= x_offset;
133 fc_dim1 = *nx;
134 fc_dim2 = *ny;
135 fc_offset = 1 + fc_dim1 * (1 + fc_dim2);
136 fc -= fc_offset;
137 cc_dim1 = *nx;
138 cc_dim2 = *ny;
139 cc_offset = 1 + cc_dim1 * (1 + cc_dim2);
140 cc -= cc_offset;
141 ac_dim1 = *nx * *ny * *nz;
142 ac_offset = 1 + ac_dim1;
143 ac -= ac_offset;
144 --ipc;
145 --rpc;
146
147 /* Function Body */
148 numdia = ipc[11];
149 if (numdia == 7) {
150 gsrb7_(nx, ny, nz, &ipc[1], &rpc[1], &ac[ac_dim1 + 1], &cc[cc_offset],
151 &fc[fc_offset], &ac[(ac_dim1 << 1) + 1], &ac[ac_dim1 * 3 + 1]
152 , &ac[(ac_dim1 << 2) + 1], &x[x_offset], &w1[w1_offset], &w2[
153 w2_offset], &r__[r_offset], itmax, iters, errtol, omega,
154 iresid, iadjoint);
155 } else if (numdia == 27) {
156 gsrb27_(nx, ny, nz, &ipc[1], &rpc[1], &ac[ac_dim1 + 1], &cc[cc_offset]
157 , &fc[fc_offset], &ac[(ac_dim1 << 1) + 1], &ac[ac_dim1 * 3 +
158 1], &ac[(ac_dim1 << 2) + 1], &ac[ac_dim1 * 5 + 1], &ac[
159 ac_dim1 * 6 + 1], &ac[ac_dim1 * 7 + 1], &ac[(ac_dim1 << 3) +
160 1], &ac[ac_dim1 * 9 + 1], &ac[ac_dim1 * 10 + 1], &ac[ac_dim1 *
161 11 + 1], &ac[ac_dim1 * 12 + 1], &ac[ac_dim1 * 13 + 1], &ac[
162 ac_dim1 * 14 + 1], &x[x_offset], &w1[w1_offset], &w2[
163 w2_offset], &r__[r_offset], itmax, iters, errtol, omega,
164 iresid, iadjoint);
165 } else {
166 s_wsle(&io___2);
167 do_lio(&c__9, &c__1, "% GSRB: invalid stencil type given...", (ftnlen)
168 37);
169 e_wsle();
170 }
171 /* * */
172 /* * *** return and end *** */
173 return 0;
174 } /* gsrb_ */
175
gsrb7_(integer * nx,integer * ny,integer * nz,integer * ipc,doublereal * rpc,doublereal * oc,doublereal * cc,doublereal * fc,doublereal * oe,doublereal * on,doublereal * uc,doublereal * x,doublereal * w1,doublereal * w2,doublereal * r__,integer * itmax,integer * iters,doublereal * errtol,doublereal * omega,integer * iresid,integer * iadjoint)176 /* Subroutine */ int gsrb7_(integer *nx, integer *ny, integer *nz, integer *
177 ipc, doublereal *rpc, doublereal *oc, doublereal *cc, doublereal *fc,
178 doublereal *oe, doublereal *on, doublereal *uc, doublereal *x,
179 doublereal *w1, doublereal *w2, doublereal *r__, integer *itmax,
180 integer *iters, doublereal *errtol, doublereal *omega, integer *
181 iresid, integer *iadjoint)
182 {
183 /* System generated locals */
184 integer oe_dim1, oe_dim2, oe_offset, on_dim1, on_dim2, on_offset, uc_dim1,
185 uc_dim2, uc_offset, fc_dim1, fc_dim2, fc_offset, oc_dim1,
186 oc_dim2, oc_offset, cc_dim1, cc_dim2, cc_offset, x_dim1, x_dim2,
187 x_offset, w1_dim1, w1_dim2, w1_offset, w2_dim1, w2_dim2,
188 w2_offset, r_dim1, r_dim2, r_offset, i__1, i__2, i__3, i__4, i__5,
189 i__6, i__7;
190
191 /* Local variables */
192 static integer i__, j, k, i1, i2, j1, j2, k1, k2;
193 extern /* Subroutine */ int mresid7_1s__(integer *, integer *, integer *,
194 integer *, doublereal *, doublereal *, doublereal *, doublereal *,
195 doublereal *, doublereal *, doublereal *, doublereal *,
196 doublereal *);
197 static integer istep;
198
199 /* * ********************************************************************* */
200 /* * purpose: */
201 /* * */
202 /* * 7 diagonal gauss-seidel routine. */
203 /* * */
204 /* * this routine applies the gauss-seidel operator or its */
205 /* * adjoint depending on the flag iadjoint. */
206 /* * */
207 /* * author: michael holst */
208 /* * ********************************************************************* */
209 /* * */
210 /* mdir 0 0 */
211 /* * */
212 /* * *** do the gauss-seidel iteration itmax times *** */
213 /* Parameter adjustments */
214 r_dim1 = *nx;
215 r_dim2 = *ny;
216 r_offset = 1 + r_dim1 * (1 + r_dim2);
217 r__ -= r_offset;
218 w2_dim1 = *nx;
219 w2_dim2 = *ny;
220 w2_offset = 1 + w2_dim1 * (1 + w2_dim2);
221 w2 -= w2_offset;
222 w1_dim1 = *nx;
223 w1_dim2 = *ny;
224 w1_offset = 1 + w1_dim1 * (1 + w1_dim2);
225 w1 -= w1_offset;
226 x_dim1 = *nx;
227 x_dim2 = *ny;
228 x_offset = 1 + x_dim1 * (1 + x_dim2);
229 x -= x_offset;
230 uc_dim1 = *nx;
231 uc_dim2 = *ny;
232 uc_offset = 1 + uc_dim1 * (1 + uc_dim2);
233 uc -= uc_offset;
234 on_dim1 = *nx;
235 on_dim2 = *ny;
236 on_offset = 1 + on_dim1 * (1 + on_dim2);
237 on -= on_offset;
238 oe_dim1 = *nx;
239 oe_dim2 = *ny;
240 oe_offset = 1 + oe_dim1 * (1 + oe_dim2);
241 oe -= oe_offset;
242 fc_dim1 = *nx;
243 fc_dim2 = *ny;
244 fc_offset = 1 + fc_dim1 * (1 + fc_dim2);
245 fc -= fc_offset;
246 cc_dim1 = *nx;
247 cc_dim2 = *ny;
248 cc_offset = 1 + cc_dim1 * (1 + cc_dim2);
249 cc -= cc_offset;
250 oc_dim1 = *nx;
251 oc_dim2 = *ny;
252 oc_offset = 1 + oc_dim1 * (1 + oc_dim2);
253 oc -= oc_offset;
254 --ipc;
255 --rpc;
256
257 /* Function Body */
258 i1 = (1 - *iadjoint << 1) + *iadjoint * (*nx - 1);
259 i2 = (*iadjoint << 1) + (1 - *iadjoint) * (*nx - 1);
260 j1 = (1 - *iadjoint << 1) + *iadjoint * (*ny - 1);
261 j2 = (*iadjoint << 1) + (1 - *iadjoint) * (*ny - 1);
262 k1 = (1 - *iadjoint << 1) + *iadjoint * (*nz - 1);
263 k2 = (*iadjoint << 1) + (1 - *iadjoint) * (*nz - 1);
264 istep = -(*iadjoint) + (1 - *iadjoint);
265 i__1 = *itmax;
266 for (*iters = 1; *iters <= i__1; ++(*iters)) {
267 /* * */
268 /* * *** do the red points *** */
269 /* mdir 3 1 */
270 i__2 = k2;
271 i__3 = istep;
272 for (k = k1; i__3 < 0 ? k >= i__2 : k <= i__2; k += i__3) {
273 /* mdir 3 2 */
274 i__4 = j2;
275 i__5 = istep;
276 for (j = j1; i__5 < 0 ? j >= i__4 : j <= i__4; j += i__5) {
277 /* mdir 3 3 */
278 i__6 = i2;
279 i__7 = istep;
280 for (i__ = i1; i__7 < 0 ? i__ >= i__6 : i__ <= i__6; i__ +=
281 i__7) {
282 x[i__ + (j + k * x_dim2) * x_dim1] = (fc[i__ + (j + k *
283 fc_dim2) * fc_dim1] + (on[i__ + (j + k * on_dim2)
284 * on_dim1] * x[i__ + (j + 1 + k * x_dim2) *
285 x_dim1] + on[i__ + (j - 1 + k * on_dim2) *
286 on_dim1] * x[i__ + (j - 1 + k * x_dim2) * x_dim1]
287 + oe[i__ + (j + k * oe_dim2) * oe_dim1] * x[i__ +
288 1 + (j + k * x_dim2) * x_dim1] + oe[i__ - 1 + (j
289 + k * oe_dim2) * oe_dim1] * x[i__ - 1 + (j + k *
290 x_dim2) * x_dim1] + uc[i__ + (j + (k - 1) *
291 uc_dim2) * uc_dim1] * x[i__ + (j + (k - 1) *
292 x_dim2) * x_dim1] + uc[i__ + (j + k * uc_dim2) *
293 uc_dim1] * x[i__ + (j + (k + 1) * x_dim2) *
294 x_dim1])) / (oc[i__ + (j + k * oc_dim2) * oc_dim1]
295 + cc[i__ + (j + k * cc_dim2) * cc_dim1]);
296 /* L12: */
297 }
298 /* L11: */
299 }
300 /* L10: */
301 }
302 /* * */
303 /* * *** main loop *** */
304 /* L30: */
305 }
306 /* * */
307 /* * *** if specified, return the new residual as well *** */
308 if (*iresid == 1) {
309 mresid7_1s__(nx, ny, nz, &ipc[1], &rpc[1], &oc[oc_offset], &cc[
310 cc_offset], &fc[fc_offset], &oe[oe_offset], &on[on_offset], &
311 uc[uc_offset], &x[x_offset], &r__[r_offset]);
312 }
313 /* * */
314 /* * *** return and end *** */
315 return 0;
316 } /* gsrb7_ */
317
gsrb27_(integer * nx,integer * ny,integer * nz,integer * ipc,doublereal * rpc,doublereal * oc,doublereal * cc,doublereal * fc,doublereal * oe,doublereal * on,doublereal * uc,doublereal * one,doublereal * onw,doublereal * ue,doublereal * uw,doublereal * un,doublereal * us,doublereal * une,doublereal * unw,doublereal * use,doublereal * usw,doublereal * x,doublereal * w1,doublereal * w2,doublereal * r__,integer * itmax,integer * iters,doublereal * errtol,doublereal * omega,integer * iresid,integer * iadjoint)318 /* Subroutine */ int gsrb27_(integer *nx, integer *ny, integer *nz, integer *
319 ipc, doublereal *rpc, doublereal *oc, doublereal *cc, doublereal *fc,
320 doublereal *oe, doublereal *on, doublereal *uc, doublereal *one,
321 doublereal *onw, doublereal *ue, doublereal *uw, doublereal *un,
322 doublereal *us, doublereal *une, doublereal *unw, doublereal *use,
323 doublereal *usw, doublereal *x, doublereal *w1, doublereal *w2,
324 doublereal *r__, integer *itmax, integer *iters, doublereal *errtol,
325 doublereal *omega, integer *iresid, integer *iadjoint)
326 {
327 /* System generated locals */
328 integer oe_dim1, oe_dim2, oe_offset, on_dim1, on_dim2, on_offset, uc_dim1,
329 uc_dim2, uc_offset, one_dim1, one_dim2, one_offset, onw_dim1,
330 onw_dim2, onw_offset, ue_dim1, ue_dim2, ue_offset, uw_dim1,
331 uw_dim2, uw_offset, un_dim1, un_dim2, un_offset, us_dim1, us_dim2,
332 us_offset, une_dim1, une_dim2, une_offset, unw_dim1, unw_dim2,
333 unw_offset, use_dim1, use_dim2, use_offset, usw_dim1, usw_dim2,
334 usw_offset, fc_dim1, fc_dim2, fc_offset, oc_dim1, oc_dim2,
335 oc_offset, cc_dim1, cc_dim2, cc_offset, x_dim1, x_dim2, x_offset,
336 w1_dim1, w1_dim2, w1_offset, w2_dim1, w2_dim2, w2_offset, r_dim1,
337 r_dim2, r_offset, i__1, i__2, i__3, i__4, i__5, i__6, i__7;
338
339 /* Local variables */
340 static integer i__, j, k, i1, i2, j1, j2, k1, k2;
341 extern /* Subroutine */ int mresid27_1s__(integer *, integer *, integer *,
342 integer *, doublereal *, doublereal *, doublereal *, doublereal *
343 , doublereal *, doublereal *, doublereal *, doublereal *,
344 doublereal *, doublereal *, doublereal *, doublereal *,
345 doublereal *, doublereal *, doublereal *, doublereal *,
346 doublereal *, doublereal *, doublereal *);
347 static doublereal tmpd, tmpo, tmpu;
348 static integer istep;
349
350 /* * ********************************************************************* */
351 /* * purpose: */
352 /* * */
353 /* * 27 diagonal gauss-seidel routine. */
354 /* * */
355 /* * this routine applies the gauss-seidel operator or its */
356 /* * adjoint depending on the flag iadjoint. */
357 /* * */
358 /* * author: michael holst */
359 /* * ********************************************************************* */
360 /* * */
361 /* mdir 0 0 */
362 /* * */
363 /* * *** do the gauss-seidel iteration itmax times *** */
364 /* Parameter adjustments */
365 r_dim1 = *nx;
366 r_dim2 = *ny;
367 r_offset = 1 + r_dim1 * (1 + r_dim2);
368 r__ -= r_offset;
369 w2_dim1 = *nx;
370 w2_dim2 = *ny;
371 w2_offset = 1 + w2_dim1 * (1 + w2_dim2);
372 w2 -= w2_offset;
373 w1_dim1 = *nx;
374 w1_dim2 = *ny;
375 w1_offset = 1 + w1_dim1 * (1 + w1_dim2);
376 w1 -= w1_offset;
377 x_dim1 = *nx;
378 x_dim2 = *ny;
379 x_offset = 1 + x_dim1 * (1 + x_dim2);
380 x -= x_offset;
381 usw_dim1 = *nx;
382 usw_dim2 = *ny;
383 usw_offset = 1 + usw_dim1 * (1 + usw_dim2);
384 usw -= usw_offset;
385 use_dim1 = *nx;
386 use_dim2 = *ny;
387 use_offset = 1 + use_dim1 * (1 + use_dim2);
388 use -= use_offset;
389 unw_dim1 = *nx;
390 unw_dim2 = *ny;
391 unw_offset = 1 + unw_dim1 * (1 + unw_dim2);
392 unw -= unw_offset;
393 une_dim1 = *nx;
394 une_dim2 = *ny;
395 une_offset = 1 + une_dim1 * (1 + une_dim2);
396 une -= une_offset;
397 us_dim1 = *nx;
398 us_dim2 = *ny;
399 us_offset = 1 + us_dim1 * (1 + us_dim2);
400 us -= us_offset;
401 un_dim1 = *nx;
402 un_dim2 = *ny;
403 un_offset = 1 + un_dim1 * (1 + un_dim2);
404 un -= un_offset;
405 uw_dim1 = *nx;
406 uw_dim2 = *ny;
407 uw_offset = 1 + uw_dim1 * (1 + uw_dim2);
408 uw -= uw_offset;
409 ue_dim1 = *nx;
410 ue_dim2 = *ny;
411 ue_offset = 1 + ue_dim1 * (1 + ue_dim2);
412 ue -= ue_offset;
413 onw_dim1 = *nx;
414 onw_dim2 = *ny;
415 onw_offset = 1 + onw_dim1 * (1 + onw_dim2);
416 onw -= onw_offset;
417 one_dim1 = *nx;
418 one_dim2 = *ny;
419 one_offset = 1 + one_dim1 * (1 + one_dim2);
420 one -= one_offset;
421 uc_dim1 = *nx;
422 uc_dim2 = *ny;
423 uc_offset = 1 + uc_dim1 * (1 + uc_dim2);
424 uc -= uc_offset;
425 on_dim1 = *nx;
426 on_dim2 = *ny;
427 on_offset = 1 + on_dim1 * (1 + on_dim2);
428 on -= on_offset;
429 oe_dim1 = *nx;
430 oe_dim2 = *ny;
431 oe_offset = 1 + oe_dim1 * (1 + oe_dim2);
432 oe -= oe_offset;
433 fc_dim1 = *nx;
434 fc_dim2 = *ny;
435 fc_offset = 1 + fc_dim1 * (1 + fc_dim2);
436 fc -= fc_offset;
437 cc_dim1 = *nx;
438 cc_dim2 = *ny;
439 cc_offset = 1 + cc_dim1 * (1 + cc_dim2);
440 cc -= cc_offset;
441 oc_dim1 = *nx;
442 oc_dim2 = *ny;
443 oc_offset = 1 + oc_dim1 * (1 + oc_dim2);
444 oc -= oc_offset;
445 --ipc;
446 --rpc;
447
448 /* Function Body */
449 i1 = (1 - *iadjoint << 1) + *iadjoint * (*nx - 1);
450 i2 = (*iadjoint << 1) + (1 - *iadjoint) * (*nx - 1);
451 j1 = (1 - *iadjoint << 1) + *iadjoint * (*ny - 1);
452 j2 = (*iadjoint << 1) + (1 - *iadjoint) * (*ny - 1);
453 k1 = (1 - *iadjoint << 1) + *iadjoint * (*nz - 1);
454 k2 = (*iadjoint << 1) + (1 - *iadjoint) * (*nz - 1);
455 istep = -(*iadjoint) + (1 - *iadjoint);
456 i__1 = *itmax;
457 for (*iters = 1; *iters <= i__1; ++(*iters)) {
458 /* * */
459 /* * *** do all of the points *** */
460 /* mdir 3 1 */
461 i__2 = k2;
462 i__3 = istep;
463 for (k = k1; i__3 < 0 ? k >= i__2 : k <= i__2; k += i__3) {
464 /* mdir 3 2 */
465 i__4 = j2;
466 i__5 = istep;
467 for (j = j1; i__5 < 0 ? j >= i__4 : j <= i__4; j += i__5) {
468 /* mdir 3 3 */
469 i__6 = i2;
470 i__7 = istep;
471 for (i__ = i1; i__7 < 0 ? i__ >= i__6 : i__ <= i__6; i__ +=
472 i__7) {
473 tmpo = on[i__ + (j + k * on_dim2) * on_dim1] * x[i__ + (j
474 + 1 + k * x_dim2) * x_dim1] + on[i__ + (j - 1 + k
475 * on_dim2) * on_dim1] * x[i__ + (j - 1 + k *
476 x_dim2) * x_dim1] + oe[i__ + (j + k * oe_dim2) *
477 oe_dim1] * x[i__ + 1 + (j + k * x_dim2) * x_dim1]
478 + oe[i__ - 1 + (j + k * oe_dim2) * oe_dim1] * x[
479 i__ - 1 + (j + k * x_dim2) * x_dim1] + one[i__ + (
480 j + k * one_dim2) * one_dim1] * x[i__ + 1 + (j +
481 1 + k * x_dim2) * x_dim1] + onw[i__ + (j + k *
482 onw_dim2) * onw_dim1] * x[i__ - 1 + (j + 1 + k *
483 x_dim2) * x_dim1] + onw[i__ + 1 + (j - 1 + k *
484 onw_dim2) * onw_dim1] * x[i__ + 1 + (j - 1 + k *
485 x_dim2) * x_dim1] + one[i__ - 1 + (j - 1 + k *
486 one_dim2) * one_dim1] * x[i__ - 1 + (j - 1 + k *
487 x_dim2) * x_dim1];
488 tmpu = uc[i__ + (j + k * uc_dim2) * uc_dim1] * x[i__ + (j
489 + (k + 1) * x_dim2) * x_dim1] + un[i__ + (j + k *
490 un_dim2) * un_dim1] * x[i__ + (j + 1 + (k + 1) *
491 x_dim2) * x_dim1] + us[i__ + (j + k * us_dim2) *
492 us_dim1] * x[i__ + (j - 1 + (k + 1) * x_dim2) *
493 x_dim1] + ue[i__ + (j + k * ue_dim2) * ue_dim1] *
494 x[i__ + 1 + (j + (k + 1) * x_dim2) * x_dim1] + uw[
495 i__ + (j + k * uw_dim2) * uw_dim1] * x[i__ - 1 + (
496 j + (k + 1) * x_dim2) * x_dim1] + une[i__ + (j +
497 k * une_dim2) * une_dim1] * x[i__ + 1 + (j + 1 + (
498 k + 1) * x_dim2) * x_dim1] + unw[i__ + (j + k *
499 unw_dim2) * unw_dim1] * x[i__ - 1 + (j + 1 + (k +
500 1) * x_dim2) * x_dim1] + use[i__ + (j + k *
501 use_dim2) * use_dim1] * x[i__ + 1 + (j - 1 + (k +
502 1) * x_dim2) * x_dim1] + usw[i__ + (j + k *
503 usw_dim2) * usw_dim1] * x[i__ - 1 + (j - 1 + (k +
504 1) * x_dim2) * x_dim1];
505 tmpd = uc[i__ + (j + (k - 1) * uc_dim2) * uc_dim1] * x[
506 i__ + (j + (k - 1) * x_dim2) * x_dim1] + us[i__ +
507 (j + 1 + (k - 1) * us_dim2) * us_dim1] * x[i__ + (
508 j + 1 + (k - 1) * x_dim2) * x_dim1] + un[i__ + (j
509 - 1 + (k - 1) * un_dim2) * un_dim1] * x[i__ + (j
510 - 1 + (k - 1) * x_dim2) * x_dim1] + uw[i__ + 1 + (
511 j + (k - 1) * uw_dim2) * uw_dim1] * x[i__ + 1 + (
512 j + (k - 1) * x_dim2) * x_dim1] + ue[i__ - 1 + (j
513 + (k - 1) * ue_dim2) * ue_dim1] * x[i__ - 1 + (j
514 + (k - 1) * x_dim2) * x_dim1] + usw[i__ + 1 + (j
515 + 1 + (k - 1) * usw_dim2) * usw_dim1] * x[i__ + 1
516 + (j + 1 + (k - 1) * x_dim2) * x_dim1] + use[i__
517 - 1 + (j + 1 + (k - 1) * use_dim2) * use_dim1] *
518 x[i__ - 1 + (j + 1 + (k - 1) * x_dim2) * x_dim1]
519 + unw[i__ + 1 + (j - 1 + (k - 1) * unw_dim2) *
520 unw_dim1] * x[i__ + 1 + (j - 1 + (k - 1) * x_dim2)
521 * x_dim1] + une[i__ - 1 + (j - 1 + (k - 1) *
522 une_dim2) * une_dim1] * x[i__ - 1 + (j - 1 + (k -
523 1) * x_dim2) * x_dim1];
524 x[i__ + (j + k * x_dim2) * x_dim1] = (fc[i__ + (j + k *
525 fc_dim2) * fc_dim1] + (tmpo + tmpu + tmpd)) / (oc[
526 i__ + (j + k * oc_dim2) * oc_dim1] + cc[i__ + (j
527 + k * cc_dim2) * cc_dim1]);
528 /* L12: */
529 }
530 /* L11: */
531 }
532 /* L10: */
533 }
534 /* * */
535 /* * *** main loop *** */
536 /* L30: */
537 }
538 /* * */
539 /* * *** if specified, return the new residual as well *** */
540 if (*iresid == 1) {
541 mresid27_1s__(nx, ny, nz, &ipc[1], &rpc[1], &oc[oc_offset], &cc[
542 cc_offset], &fc[fc_offset], &oe[oe_offset], &on[on_offset], &
543 uc[uc_offset], &one[one_offset], &onw[onw_offset], &ue[
544 ue_offset], &uw[uw_offset], &un[un_offset], &us[us_offset], &
545 une[une_offset], &unw[unw_offset], &use[use_offset], &usw[
546 usw_offset], &x[x_offset], &r__[r_offset]);
547 }
548 /* * */
549 /* * *** return and end *** */
550 return 0;
551 } /* gsrb27_ */
552
gsrb7x_(integer * nx,integer * ny,integer * nz,integer * ipc,doublereal * rpc,doublereal * oc,doublereal * cc,doublereal * fc,doublereal * oe,doublereal * on,doublereal * uc,doublereal * x,doublereal * w1,doublereal * w2,doublereal * r__,integer * itmax,integer * iters,doublereal * errtol,doublereal * omega,integer * iresid,integer * iadjoint)553 /* Subroutine */ int gsrb7x_(integer *nx, integer *ny, integer *nz, integer *
554 ipc, doublereal *rpc, doublereal *oc, doublereal *cc, doublereal *fc,
555 doublereal *oe, doublereal *on, doublereal *uc, doublereal *x,
556 doublereal *w1, doublereal *w2, doublereal *r__, integer *itmax,
557 integer *iters, doublereal *errtol, doublereal *omega, integer *
558 iresid, integer *iadjoint)
559 {
560 /* System generated locals */
561 integer oe_dim1, oe_dim2, oe_offset, on_dim1, on_dim2, on_offset, uc_dim1,
562 uc_dim2, uc_offset, fc_dim1, fc_dim2, fc_offset, oc_dim1,
563 oc_dim2, oc_offset, cc_dim1, cc_dim2, cc_offset, x_dim1, x_dim2,
564 x_offset, w1_dim1, w1_dim2, w1_offset, w2_dim1, w2_dim2,
565 w2_offset, r_dim1, r_dim2, r_offset, i__1, i__2, i__3, i__4;
566
567 /* Local variables */
568 static integer i__, j, k;
569 extern /* Subroutine */ int mresid7_1s__(integer *, integer *, integer *,
570 integer *, doublereal *, doublereal *, doublereal *, doublereal *,
571 doublereal *, doublereal *, doublereal *, doublereal *,
572 doublereal *);
573 static integer ioff;
574
575 /* * ********************************************************************* */
576 /* * purpose: */
577 /* * */
578 /* * fast 7 diagonal red/black gauss-seidel routine. */
579 /* * */
580 /* * this routine applies the red/black gauss-seidel operator or its */
581 /* * adjoint depending on the flag iadjoint. note that the adjoint */
582 /* * in this case is simply doing the red or black points first. */
583 /* * */
584 /* * note that the interior grid points from 2,...,nx-1, etc. */
585 /* * we then begin coloring with a red point, and do a point red/black */
586 /* * coloring. */
587 /* * */
588 /* * the red points are: */
589 /* * */
590 /* * if ((j even) and (k even)) then */
591 /* * begin row at first point, i=2 (i even), or ioff = 0 */
592 /* * else if ((j odd) and (k even)) then */
593 /* * begin row at second point, i=3 (i odd), or ioff = 1 */
594 /* * else if ((j even) and (k odd)) then */
595 /* * begin row at second point, i=3 (i odd), or ioff = 1 */
596 /* * else if ((j odd) and (k odd)) then */
597 /* * begin row at first point, i=2 (i even), or ioff = 0 */
598 /* * endif */
599 /* * then: begin row at: i=2+ioff */
600 /* * */
601 /* * the appropriate ioff function for the red points is then: */
602 /* * ioff = dabs(mod(j,2)-mod(k,2)) */
603 /* * */
604 /* * the appropriate ioff function for the black points is then: */
605 /* * ioff = 1 - dabs(mod(j,2)-mod(k,2)) */
606 /* * */
607 /* * */
608 /* * alternatively, the red points are: */
609 /* * */
610 /* * those whose indices add up to an even number. */
611 /* * to see this, consider that all surrounding points are */
612 /* * only one index different, hence the sum will differ by one. */
613 /* * thus, if a given point has an even sum, then the surrounding */
614 /* * points will have an odd sum. */
615 /* * */
616 /* * thus, the black points are: */
617 /* * */
618 /* * therefore those whose indices add up to an odd number. */
619 /* * */
620 /* * author: michael holst */
621 /* * ********************************************************************* */
622 /* * */
623 /* mdir 0 0 */
624 /* * */
625 /* * *** do the gauss-seidel iteration itmax times *** */
626 /* Parameter adjustments */
627 r_dim1 = *nx;
628 r_dim2 = *ny;
629 r_offset = 1 + r_dim1 * (1 + r_dim2);
630 r__ -= r_offset;
631 w2_dim1 = *nx;
632 w2_dim2 = *ny;
633 w2_offset = 1 + w2_dim1 * (1 + w2_dim2);
634 w2 -= w2_offset;
635 w1_dim1 = *nx;
636 w1_dim2 = *ny;
637 w1_offset = 1 + w1_dim1 * (1 + w1_dim2);
638 w1 -= w1_offset;
639 x_dim1 = *nx;
640 x_dim2 = *ny;
641 x_offset = 1 + x_dim1 * (1 + x_dim2);
642 x -= x_offset;
643 uc_dim1 = *nx;
644 uc_dim2 = *ny;
645 uc_offset = 1 + uc_dim1 * (1 + uc_dim2);
646 uc -= uc_offset;
647 on_dim1 = *nx;
648 on_dim2 = *ny;
649 on_offset = 1 + on_dim1 * (1 + on_dim2);
650 on -= on_offset;
651 oe_dim1 = *nx;
652 oe_dim2 = *ny;
653 oe_offset = 1 + oe_dim1 * (1 + oe_dim2);
654 oe -= oe_offset;
655 fc_dim1 = *nx;
656 fc_dim2 = *ny;
657 fc_offset = 1 + fc_dim1 * (1 + fc_dim2);
658 fc -= fc_offset;
659 cc_dim1 = *nx;
660 cc_dim2 = *ny;
661 cc_offset = 1 + cc_dim1 * (1 + cc_dim2);
662 cc -= cc_offset;
663 oc_dim1 = *nx;
664 oc_dim2 = *ny;
665 oc_offset = 1 + oc_dim1 * (1 + oc_dim2);
666 oc -= oc_offset;
667 --ipc;
668 --rpc;
669
670 /* Function Body */
671 i__1 = *itmax;
672 for (*iters = 1; *iters <= i__1; ++(*iters)) {
673 /* * */
674 /* * *** do the red points *** */
675 /* mdir 3 1 */
676 i__2 = *nz - 1;
677 for (k = 2; k <= i__2; ++k) {
678 /* mdir 3 2 */
679 i__3 = *ny - 1;
680 for (j = 2; j <= i__3; ++j) {
681 /* ZZZ ioff = mod((j+k+2),2) */
682 ioff = (1 - *iadjoint) * ((j + k + 2) % 2) + *iadjoint * (1 -
683 (j + k + 2) % 2);
684 /* mdir 3 3 */
685 i__4 = *nx - 1;
686 for (i__ = ioff + 2; i__ <= i__4; i__ += 2) {
687 x[i__ + (j + k * x_dim2) * x_dim1] = (fc[i__ + (j + k *
688 fc_dim2) * fc_dim1] + (on[i__ + (j + k * on_dim2)
689 * on_dim1] * x[i__ + (j + 1 + k * x_dim2) *
690 x_dim1] + on[i__ + (j - 1 + k * on_dim2) *
691 on_dim1] * x[i__ + (j - 1 + k * x_dim2) * x_dim1]
692 + oe[i__ + (j + k * oe_dim2) * oe_dim1] * x[i__ +
693 1 + (j + k * x_dim2) * x_dim1] + oe[i__ - 1 + (j
694 + k * oe_dim2) * oe_dim1] * x[i__ - 1 + (j + k *
695 x_dim2) * x_dim1] + uc[i__ + (j + (k - 1) *
696 uc_dim2) * uc_dim1] * x[i__ + (j + (k - 1) *
697 x_dim2) * x_dim1] + uc[i__ + (j + k * uc_dim2) *
698 uc_dim1] * x[i__ + (j + (k + 1) * x_dim2) *
699 x_dim1])) / (oc[i__ + (j + k * oc_dim2) * oc_dim1]
700 + cc[i__ + (j + k * cc_dim2) * cc_dim1]);
701 /* L12: */
702 }
703 /* L11: */
704 }
705 /* L10: */
706 }
707 /* * */
708 /* * *** do the black points *** */
709 /* mdir 3 1 */
710 i__2 = *nz - 1;
711 for (k = 2; k <= i__2; ++k) {
712 /* mdir 3 2 */
713 i__3 = *ny - 1;
714 for (j = 2; j <= i__3; ++j) {
715 /* ZZZ ioff = 1-mod((j+k+2),2) */
716 ioff = *iadjoint * ((j + k + 2) % 2) + (1 - *iadjoint) * (1 -
717 (j + k + 2) % 2);
718 /* mdir 3 3 */
719 i__4 = *nx - 1;
720 for (i__ = ioff + 2; i__ <= i__4; i__ += 2) {
721 x[i__ + (j + k * x_dim2) * x_dim1] = (fc[i__ + (j + k *
722 fc_dim2) * fc_dim1] + (on[i__ + (j + k * on_dim2)
723 * on_dim1] * x[i__ + (j + 1 + k * x_dim2) *
724 x_dim1] + on[i__ + (j - 1 + k * on_dim2) *
725 on_dim1] * x[i__ + (j - 1 + k * x_dim2) * x_dim1]
726 + oe[i__ + (j + k * oe_dim2) * oe_dim1] * x[i__ +
727 1 + (j + k * x_dim2) * x_dim1] + oe[i__ - 1 + (j
728 + k * oe_dim2) * oe_dim1] * x[i__ - 1 + (j + k *
729 x_dim2) * x_dim1] + uc[i__ + (j + (k - 1) *
730 uc_dim2) * uc_dim1] * x[i__ + (j + (k - 1) *
731 x_dim2) * x_dim1] + uc[i__ + (j + k * uc_dim2) *
732 uc_dim1] * x[i__ + (j + (k + 1) * x_dim2) *
733 x_dim1])) / (oc[i__ + (j + k * oc_dim2) * oc_dim1]
734 + cc[i__ + (j + k * cc_dim2) * cc_dim1]);
735 /* L22: */
736 }
737 /* L21: */
738 }
739 /* L20: */
740 }
741 /* * */
742 /* * *** main loop *** */
743 /* L30: */
744 }
745 /* * */
746 /* * *** if specified, return the new residual as well *** */
747 if (*iresid == 1) {
748 mresid7_1s__(nx, ny, nz, &ipc[1], &rpc[1], &oc[oc_offset], &cc[
749 cc_offset], &fc[fc_offset], &oe[oe_offset], &on[on_offset], &
750 uc[uc_offset], &x[x_offset], &r__[r_offset]);
751 }
752 /* * */
753 /* * *** return and end *** */
754 return 0;
755 } /* gsrb7x_ */
756
757