APBS 3.0.0
Loading...
Searching...
No Matches
mgsubd.c
1
55#include "mgsubd.h"
56
57VPUBLIC void Vbuildops(
58 int *nx, int *ny, int *nz,
59 int *nlev, int *ipkey, int *iinfo,
60 int *ido, int *iz,
61 int *mgprol, int *mgcoar, int *mgsolv, int *mgdisc, int *ipc,
62 double *rpc, double *pc, double *ac, double *cc, double *fc,
63 double *xf, double *yf, double *zf,
64 double *gxcf, double *gycf, double *gzcf,
65 double *a1cf, double *a2cf, double *a3cf,
66 double *ccf, double *fcf, double *tcf
67 ) {
68
69 // @todo Document this function
70 int lev = 0;
71 int nxx = 0;
72 int nyy = 0;
73 int nzz = 0;
74 int nxold = 0;
75 int nyold = 0;
76 int nzold = 0;
77 int numdia = 0;
78 int key = 0;
79
80 // Utility variables
81 int i;
82
83 MAT2(iz, 50, *nlev);
84
85 // Setup
86 nxx = *nx;
87 nyy = *ny;
88 nzz = *nz;
89
90 // Build the operator a on the finest level
91 if (*ido == 0 || *ido == 2) {
92
93 lev = 1;
94
95 // Some i/o
96 if (*iinfo > 0)
97 VMESSAGE3("Fine: (%03d, %03d, %03d)", nxx, nyy, nzz);
98
99 // Finest level discretization
100 VbuildA(&nxx, &nyy, &nzz,
101 ipkey, mgdisc, &numdia,
102 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
103 RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
104 RAT(xf, VAT2(iz, 8,lev)), RAT(yf, VAT2(iz, 9,lev)), RAT(zf, VAT2(iz, 10,lev)),
105 RAT(gxcf, VAT2(iz, 2,lev)), RAT(gycf, VAT2(iz, 3,lev)), RAT(gzcf, VAT2(iz, 4,lev)),
106 RAT(a1cf, VAT2(iz, 1,lev)), RAT(a2cf, VAT2(iz, 1,lev)), RAT(a3cf, VAT2(iz, 1,lev)),
107 RAT(ccf, VAT2(iz, 1,lev)), RAT(fcf, VAT2(iz, 1,lev)));
108
109 VMESSAGE2("Operator stencil (lev, numdia) = (%d, %d)", lev, numdia);
110
111 // Now initialize the differential operator offset
112 VAT2(iz, 7, lev+1) = VAT2(iz, 7, lev) + numdia * nxx * nyy * nzz;
113
114 // Debug
115 if (*iinfo > 7) {
116 Vprtmatd(&nxx, &nyy, &nzz,
117 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)), RAT(ac, VAT2(iz, 7,lev)));
118 }
119 }
120
121 // Build the (nlev-1) level operators
122 if (*ido == 1 || *ido == 2 || *ido == 3) {
123
124 for (lev=2; lev<=*nlev; lev++) {
125 nxold = nxx;
126 nyold = nyy;
127 nzold = nzz;
128 i = 1;
129
130
131 Vmkcors(&i, &nxold, &nyold, &nzold, &nxx, &nyy, &nzz);
132 if (*ido != 3) {
133
134 // Build the interpolation operator on this level
135 VbuildP(&nxold, &nyold, &nzold,
136 &nxx, &nyy, &nzz,
137 mgprol,
138 RAT(ipc, VAT2(iz, 5,lev-1)), RAT(rpc, VAT2(iz, 6,lev-1)),
139 RAT(pc, VAT2(iz, 11,lev-1)), RAT(ac, VAT2(iz, 7,lev-1)),
140 RAT(xf, VAT2(iz, 8,lev-1)), RAT(yf, VAT2(iz, 9,lev-1)), RAT(zf, VAT2(iz, 10,lev-1)));
141
142 // Differential operator this level with standard disc.
143 if (*mgcoar == 0) {
144
145 // Some i/o
146 if (*iinfo > 0)
147 VMESSAGE3("Stand: (%03d, %03d, %03d)", nxx, nyy, nzz);
148
149
150
151 Vbuildcopy0(&nxx, &nyy, &nzz,
152 &nxold, &nyold, &nzold,
153 RAT(xf, VAT2(iz, 8,lev )), RAT(yf, VAT2(iz, 9,lev )), RAT(zf, VAT2(iz, 10,lev )),
154 RAT(gxcf, VAT2(iz, 2,lev )), RAT(gycf, VAT2(iz, 3,lev )), RAT(gzcf, VAT2(iz, 4,lev )),
155 RAT(a1cf, VAT2(iz, 1,lev )), RAT(a2cf, VAT2(iz, 1,lev )), RAT(a3cf, VAT2(iz, 1,lev )),
156 RAT(ccf, VAT2(iz, 1,lev )), RAT(fcf, VAT2(iz, 1,lev )), RAT(tcf, VAT2(iz, 1,lev )),
157 RAT(xf, VAT2(iz, 8,lev-1)), RAT(yf, VAT2(iz, 9,lev-1)), RAT(zf, VAT2(iz, 10,lev-1)),
158 RAT(gxcf, VAT2(iz, 2,lev-1)), RAT(gycf, VAT2(iz, 3,lev-1)), RAT(gzcf, VAT2(iz, 4,lev-1)),
159 RAT(a1cf, VAT2(iz, 1,lev-1)), RAT(a2cf, VAT2(iz, 1,lev-1)), RAT(a3cf, VAT2(iz, 1,lev-1)),
160 RAT(ccf, VAT2(iz, 1,lev-1)), RAT(fcf, VAT2(iz, 1,lev-1)), RAT(tcf, VAT2(iz, 1,lev-1)));
161
162 VbuildA(&nxx, &nyy, &nzz,
163 ipkey, mgdisc, &numdia,
164 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
165 RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
166 RAT(xf, VAT2(iz, 8,lev)), RAT(yf, VAT2(iz, 9,lev)), RAT(zf, VAT2(iz, 10,lev)),
167 RAT(gxcf, VAT2(iz, 2,lev)), RAT(gycf, VAT2(iz, 3,lev)), RAT(gzcf, VAT2(iz, 4,lev)),
168 RAT(a1cf, VAT2(iz, 1,lev)), RAT(a2cf, VAT2(iz, 1,lev)), RAT(a3cf, VAT2(iz, 1,lev)),
169 RAT(ccf, VAT2(iz, 1,lev)), RAT(fcf, VAT2(iz, 1,lev)));
170 }
171
172 // Differential operator this level with harmonic disc.
173 else if (*mgcoar == 1) {
174
175 // Some i/o
176 if (*iinfo > 0)
177 VMESSAGE3("Harm0: (%03d, %03d, %03d)", nxx, nyy, nzz);
178
179 Vbuildharm0(&nxx, &nyy, &nzz, &nxold, &nyold, &nzold,
180 RAT(xf, VAT2(iz, 8, lev )), RAT(yf, VAT2(iz, 9, lev )), RAT(zf, VAT2(iz, 10, lev )),
181 RAT(gxcf, VAT2(iz, 2, lev )), RAT(gycf, VAT2(iz, 3, lev )), RAT(gzcf, VAT2(iz, 4, lev )),
182 RAT(a1cf, VAT2(iz, 1, lev )), RAT(a2cf, VAT2(iz, 1, lev )), RAT(a3cf, VAT2(iz, 1, lev )),
183 RAT(ccf, VAT2(iz, 1, lev )), RAT(fcf, VAT2(iz, 1, lev )), RAT(tcf, VAT2(iz, 1, lev )),
184 RAT(xf, VAT2(iz, 8, lev-1)), RAT(yf, VAT2(iz, 9, lev-1)), RAT(zf, VAT2(iz, 10, lev-1)),
185 RAT(gxcf, VAT2(iz, 2, lev-1)), RAT(gycf, VAT2(iz, 3, lev-1)), RAT(gzcf, VAT2(iz, 4, lev-1)),
186 RAT(a1cf, VAT2(iz, 1, lev-1)), RAT(a2cf, VAT2(iz, 1, lev-1)), RAT(a3cf, VAT2(iz, 1, lev-1)),
187 RAT(ccf, VAT2(iz, 1, lev-1)), RAT(fcf, VAT2(iz, 1, lev-1)), RAT(tcf, VAT2(iz, 1, lev-1)));
188
189 VbuildA(&nxx, &nyy, &nzz,
190 ipkey, mgdisc, &numdia,
191 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
192 RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
193 RAT(xf, VAT2(iz, 8,lev)), RAT(yf, VAT2(iz, 9,lev)), RAT(zf, VAT2(iz, 10,lev)),
194 RAT(gxcf, VAT2(iz, 2,lev)), RAT(gycf, VAT2(iz, 3,lev)), RAT(gzcf, VAT2(iz, 4,lev)),
195 RAT(a1cf, VAT2(iz, 1,lev)), RAT(a2cf, VAT2(iz, 1,lev)), RAT(a3cf, VAT2(iz, 1,lev)),
196 RAT(ccf, VAT2(iz, 1,lev)), RAT(fcf, VAT2(iz, 1,lev)));
197 }
198
199 // Differential operator with galerkin formulation ***
200 else if (*mgcoar == 2) {
201
202 // Some i/o
203 if (*iinfo > 0)
204 VMESSAGE3("Galer: (%03d, %03d, %03d)", nxx, nyy, nzz);
205
206 Vbuildgaler0(&nxold, &nyold, &nzold,
207 &nxx, &nyy, &nzz,
208 ipkey, &numdia,
209 RAT(pc, VAT2(iz, 11,lev-1)),
210 RAT(ipc, VAT2(iz, 5,lev-1)), RAT(rpc, VAT2(iz, 6,lev-1)),
211 RAT(ac, VAT2(iz, 7,lev-1)), RAT(cc, VAT2(iz, 1,lev-1)), RAT(fc, VAT2(iz, 1,lev-1)),
212 RAT(ipc, VAT2(iz, 5,lev )), RAT(rpc, VAT2(iz, 6,lev )),
213 RAT(ac, VAT2(iz, 7,lev )), RAT(cc, VAT2(iz, 1,lev )), RAT(fc, VAT2(iz, 1,lev )));
214
215
216
217 Vextrac(&nxold, &nyold, &nzold,
218 &nxx, &nyy, &nzz,
219 RAT(tcf, VAT2(iz, 1,lev-1)), RAT(tcf, VAT2(iz, 1,lev)));
220 }
221 else {
222 VABORT_MSG1("Bad mgcoar value given: %d", *mgcoar);
223 }
224
225 // Now initialize the differential operator offset
226 // Vnm_print(0, "BUILDOPS: operator stencil (lev,numdia) = (%d, %d)\n",
227 // lev,numdia);
228 VAT2(iz, 7, lev+1) = VAT2(iz, 7,lev) + numdia * nxx * nyy * nzz;
229
230 // Debug
231 if (*iinfo > 8) {
232
233 Vprtmatd(&nxx, &nyy, &nzz,
234 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)), RAT(ac, VAT2(iz, 7,lev)));
235 }
236 }
237 }
238
239 // Build a sparse format coarse grid operator
240 if (*mgsolv == 1) {
241 lev = *nlev;
242
243 Vbuildband(&key, &nxx, &nyy, &nzz,
244 RAT(ipc, VAT2(iz, 5,lev )), RAT(rpc, VAT2(iz, 6,lev )), RAT(ac, VAT2(iz, 7,lev )),
245 RAT(ipc, VAT2(iz, 5,lev+1)), RAT(rpc, VAT2(iz, 6,lev+1)), RAT(ac, VAT2(iz, 7,lev+1)));
246
247 if (key == 1) {
248 VERRMSG0("Changing your mgsolv to iterative");
249 *mgsolv = 0;
250 }
251 }
252 }
253}
254
255
256
257VPUBLIC void Vbuildstr(int *nx, int *ny, int *nz, int *nlev, int *iz) {
258
259 int nxold, nyold, nzold;
260 int nxnew, nynew, nznew;
261 int n, lev;
262
263 // Utility variable
264 int numlev;
265
266 MAT2(iz, 50, *nlev);
267
268 // Setup
269 nxnew = *nx;
270 nynew = *ny;
271 nznew = *nz;
272 n = nxnew * nynew * nznew;
273
274 // Start with level 1
275 lev = 1;
276
277 // Mark beginning of everything at level 1 ***
278 VAT2(iz, 1, lev) = 1;
279 VAT2(iz, 2, lev) = 1;
280 VAT2(iz, 3, lev) = 1;
281 VAT2(iz, 4, lev) = 1;
282 VAT2(iz, 5, lev) = 1;
283 VAT2(iz, 6, lev) = 1;
284 VAT2(iz, 7, lev) = 1;
285 VAT2(iz, 8, lev) = 1;
286 VAT2(iz, 9, lev) = 1;
287 VAT2(iz, 10, lev) = 1;
288 VAT2(iz, 11, lev) = 1;
289
290 // Mark beginning of everything at level 2
291 VAT2(iz, 1, lev + 1) = VAT2(iz, 1, lev) + n;
292 VAT2(iz, 2, lev + 1) = VAT2(iz, 2, lev) + 4 * nynew * nznew;
293 VAT2(iz, 3, lev + 1) = VAT2(iz, 3, lev) + 4 * nxnew * nznew;
294 VAT2(iz, 4, lev + 1) = VAT2(iz, 4, lev) + 4 * nxnew * nynew;
295 VAT2(iz, 5, lev + 1) = VAT2(iz, 5, lev) + 100;
296 VAT2(iz, 6, lev + 1) = VAT2(iz, 6, lev) + 100;
297 VAT2(iz, 8, lev + 1) = VAT2(iz, 8, lev) + nxnew;
298 VAT2(iz, 9, lev + 1) = VAT2(iz, 9, lev) + nynew;
299 VAT2(iz, 10, lev + 1) = VAT2(iz, 10, lev) + nznew;
300
301 /* ******************************************************************
302 * ***NOTE: we mark operator offsets as we build the operators ***
303 * ***VAT2(iz, 7,lev+1) = VAT2(iz, 7, lev) + 4*n ***
304 * ******************************************************************
305 * ***NOTE: we mark prolongation operator offsets lagging a level ***
306 * ***VAT2(iz, 11, lev) = VAT2(iz, 11,lev-1) + 27*nsmall ***
307 * ******************************************************************/
308
309 // Mark the beginning of everything at (nlev-1) more
310 for (lev=2; lev<=*nlev; lev++) {
311 nxold = nxnew;
312 nyold = nynew;
313 nzold = nznew;
314 numlev = 1;
315 Vmkcors(&numlev, &nxold, &nyold, &nzold, &nxnew, &nynew, &nznew);
316 n = nxnew * nynew * nznew;
317
318 // Mark the beginning of everything at level (lev+1)
319 VAT2(iz, 1, lev + 1) = VAT2(iz, 1, lev) + n;
320 VAT2(iz, 2, lev + 1) = VAT2(iz, 2, lev) + 4 * nynew * nznew;
321 VAT2(iz, 3, lev + 1) = VAT2(iz, 3, lev) + 4 * nxnew * nznew;
322 VAT2(iz, 4, lev + 1) = VAT2(iz, 4, lev) + 4 * nxnew * nynew;
323 VAT2(iz, 5, lev + 1) = VAT2(iz, 5, lev) + 100;
324 VAT2(iz, 6, lev + 1) = VAT2(iz, 6, lev) + 100;
325 VAT2(iz, 7, lev + 1) = VAT2(iz, 7, lev) + 4 * n;
326 VAT2(iz, 8, lev + 1) = VAT2(iz, 8, lev) + nxnew;
327 VAT2(iz, 9, lev + 1) = VAT2(iz, 9, lev) + nynew;
328 VAT2(iz, 10, lev + 1) = VAT2(iz, 10, lev) + nznew;
329
330 // Mark prolongation operator storage for previous level
331 VAT2(iz, 11, lev) = VAT2(iz, 11, lev - 1) + 27 * n;
332
333 /* ****************************************************************
334 * *** NOTE: we mark operator offsets as we build the operators ***
335 * *** VAT2(iz, 7, lev + 1) = VAT2(iz, 7, lev) + 4 * n ***
336 ******************************************************************/
337 }
338}
339
340
341VPUBLIC void Vbuildgaler0(int *nxf, int *nyf, int *nzf,
342 int *nxc, int *nyc, int *nzc,
343 int *ipkey, int *numdia,
344 double *pcFF, int *ipcFF, double *rpcFF,
345 double *acFF, double *ccFF, double *fcFF,
346 int *ipc, double *rpc,
347 double *ac, double *cc, double *fc) {
348
349 int numdia_loc;
350
351 // Call the algebraic galerkin routine
352 numdia_loc = VAT(ipcFF, 11);
353 VbuildG(nxf, nyf, nzf,
354 nxc, nyc, nzc,
355 &numdia_loc,
356 pcFF, acFF, ac);
357
358 // Note how many nonzeros in this new discretization stencil
359 VAT(ipc, 11) = 27;
360 *numdia = 14;
361
362 // Save the problem key with this new operator
363 VAT(ipc, 10) = *ipkey;
364
365 // Restrict the helmholtz term and source function
366 Vrestrc(nxf, nyf, nzf,
367 nxc, nyc, nzc,
368 ccFF, cc, pcFF);
369
370 Vrestrc(nxf, nyf, nzf,
371 nxc, nyc, nzc,
372 fcFF, fc, pcFF);
373}
374
375
376
377VPUBLIC void Vmkcors(int *numlev,
378 int *nxold, int *nyold, int *nzold,
379 int *nxnew, int *nynew, int *nznew) {
380 int nxtmp, nytmp, nztmp; // Temporary variables to hold current x,y,z values
381 int i; // Index used in for loops
382
383 // Determine the coarser grid
384 *nxnew = *nxold;
385 *nynew = *nyold;
386 *nznew = *nzold;
387
388 for (i=1; i<=*numlev; i++) {
389 nxtmp = *nxnew;
390 nytmp = *nynew;
391 nztmp = *nznew;
392 Vcorsr(&nxtmp,nxnew);
393 Vcorsr(&nytmp,nynew);
394 Vcorsr(&nztmp,nznew);
395 }
396}
397
398
399
400VPUBLIC void Vcorsr(int *nold, int *nnew) {
401
402 // Find the coarser grid size ***
403 *nnew = (*nold - 1) / 2 + 1;
404
405 // Check a few things
406 if ((*nnew - 1) * 2 != *nold - 1) {
407 Vnm_print(2, "Vcorsr: WARNING! The grid dimensions are not\n");
408 Vnm_print(2, "Vcorsr: consistent with nlev! Your\n");
409 Vnm_print(2, "Vcorsr: calculation will only work if you\n");
410 Vnm_print(2, "Vcorsr: are performing a mg-dummy run.\n");
411 }
412 if (*nnew < 1) {
414 Vnm_print(2, "Vcorsr: ERROR! The grid dimensions are not\n");
415 Vnm_print(2, "Vcorsr: consistent with nlev!\n");
416 Vnm_print(2, "Vcorsr: Grid coarsened below zero.\n");
417 }
418}
419
420
421
422VPUBLIC void Vmkfine(int *numlev,
423 int *nxold, int *nyold, int *nzold,
424 int *nxnew, int *nynew, int *nznew) {
425
426 int nxtmp, nytmp, nztmp, i;
427
428 // Determine the finer grid
429 *nxnew = *nxold;
430 *nynew = *nyold;
431 *nznew = *nzold;
432
433 for (i=1; i<=*numlev; i++) {
434
435 nxtmp = *nxnew;
436 nytmp = *nynew;
437 nztmp = *nznew;
438
439 Vfiner(&nxtmp, nxnew);
440 Vfiner(&nytmp, nynew);
441 Vfiner(&nztmp, nznew);
442 }
443
444}
445
446
447
448VPUBLIC void Vfiner(int *nold, int *nnew) {
449 // Find the coarser grid size ***
450 *nnew = (*nold - 1) * 2 + 1;
451}
452
453
454
455VPUBLIC int Vivariv(int *nu, int *level) {
456
458 int ivariv;
459
460 // Variable V-cycle
461 // ivariv = *nu * VPOW(2, level - 1);
462
463 // Standard V-cycle
464 ivariv = *nu;
465
466 return ivariv;
467}
468
469
470
471VPUBLIC int Vmaxlev(int n1, int n2, int n3) {
472
473 int n1c;
474 int n2c;
475 int n3c;
476 int lev;
477 int iden;
478 int idone;
479
480 // Fine the common level
481 idone = 0;
482 lev = 0;
483 do {
484 lev += 1;
485 iden = (int)VPOW(2, lev - 1);
486 n1c = (n1 - 1) / iden + 1;
487 n2c = (n2 - 1) / iden + 1;
488 n3c = (n3 - 1) / iden + 1;
489 if (((n1c - 1) * iden != (n1 - 1)) || (n1c <= 2) )
490 idone = 1;
491 if (((n2c - 1) * iden != (n2 - 1)) || (n2c <= 2) )
492 idone = 1;
493 if (((n3c - 1) * iden != (n3 - 1)) || (n3c <= 2) )
494 idone = 1;
495
496 } while (idone != 1);
497
498 return lev - 1;
499}
500
501
502
503VPUBLIC void Vprtstp(int iok, int iters,
504 double rsnrm, double rsden, double orsnrm) {
505
506 double relres = 0.0;
507 double contrac = 0.0;
508
509 // Initializing timer
510 if (iters == -99) {
511 // Vnm_tstart(40, "MG iteration");
512 cputme = 0.0;
513 return;
514 }
515
516 // Setup for the iteration
517 else if (iters == -1) {
518 Vnm_tstop(40, "MG iteration");
519 return;
520 }
521
522 // During the iteration
523 else {
524
525 // Stop the timer
526 // Vnm_tstop(40, "MG iteration");
527
528 // Relative residual
529 if (rsden == 0.0) {
530 relres = 1.0e6;
531 VERRMSG0("Vprtstp: avoided division by zero\n");
532 } else {
533 relres = rsnrm / rsden;
534 }
535
536 // Contraction number
537 if (orsnrm == 0.0) {
538 contrac = 1.0e6;
539 VERRMSG0("avoided division by zero\n");
540 } else {
541 contrac = rsnrm / orsnrm;
542 }
543
544 // The i/o
545 if (iok == 1 || iok == 2) {
546 VMESSAGE1("iteration = %d", iters);
547 VMESSAGE1("relative residual = %e", relres);
548 VMESSAGE1("contraction number = %e", contrac);
549 }
550 }
551}
552
553
554
555VPUBLIC void Vpackmg(int *iparm, double *rparm, size_t *nrwk, int *niwk,
556 int *nx, int *ny, int *nz, int *nlev, int *nu1, int *nu2, int *mgkey,
557 int *itmax, int *istop, int *ipcon, int *nonlin, int *mgsmoo, int *mgprol,
558 int *mgcoar, int *mgsolv, int *mgdisc, int *iinfo, double *errtol,
559 int *ipkey, double *omegal, double *omegan, int *irite, int *iperf) {
560
562
563 // Encode iparm parameters ***
564 VAT(iparm, 1) = *nrwk;
565 VAT(iparm, 2) = *niwk;
566 VAT(iparm, 3) = *nx;
567 VAT(iparm, 4) = *ny;
568 VAT(iparm, 5) = *nz;
569 VAT(iparm, 6) = *nlev;
570 VAT(iparm, 7) = *nu1;
571 VAT(iparm, 8) = *nu2;
572 VAT(iparm, 9) = *mgkey;
573 VAT(iparm, 10) = *itmax;
574 VAT(iparm, 11) = *istop;
575 VAT(iparm, 12) = *iinfo;
576 VAT(iparm, 13) = *irite;
577 VAT(iparm, 14) = *ipkey;
578 VAT(iparm, 15) = *ipcon;
579 VAT(iparm, 16) = *nonlin;
580 VAT(iparm, 17) = *mgprol;
581 VAT(iparm, 18) = *mgcoar;
582 VAT(iparm, 19) = *mgdisc;
583 VAT(iparm, 20) = *mgsmoo;
584 VAT(iparm, 21) = *mgsolv;
585 VAT(iparm, 22) = *iperf;
586
587 // Encode rparm parameters
588 VAT(rparm, 1) = *errtol;
589 VAT(rparm, 9) = *omegal;
590 VAT(rparm, 10) = *omegan;
591}
592
593
594
595VEXTERNC void Vbuildcopy0(int *nx, int *ny, int *nz,
596 int *nxf, int *nyf, int *nzf,
597 double *xc, double *yc, double *zc,
598 double *gxc, double *gyc, double *gzc,
599 double *a1c, double *a2c, double *a3c,
600 double *cc, double *fc, double *tc,
601 double *xf, double *yf, double *zf,
602 double *gxcf, double *gycf, double *gzcf,
603 double *a1cf, double *a2cf, double *a3cf,
604 double *ccf, double *fcf, double *tcf) {
605
606
607 int i, j, k;
608 int ii, jj, kk;
609 int iadd, jadd, kadd;
610
611 MAT3( gxc, *ny, *nz, 4);
612 MAT3( gyc, *nx, *nz, 4);
613 MAT3( gzc, *nx, *ny, 4);
614 MAT3( a1c, *nx, *ny, *nz);
615 MAT3( a2c, *nx, *ny, *nz);
616 MAT3( a3c, *nx, *ny, *nz);
617 MAT3( cc, *nx, *ny, *nz);
618 MAT3( fc, *nx, *ny, *nz);
619 MAT3( tc, *nx, *ny, *nz);
620 MAT3(gxcf, *nyf, *nzf, 4);
621 MAT3(gycf, *nxf, *nzf, 4);
622 MAT3(gzcf, *nxf, *nyf, 4);
623 MAT3(a1cf, *nxf, *nyf, *nzf);
624 MAT3(a2cf, *nxf, *nyf, *nzf);
625 MAT3(a3cf, *nxf, *nyf, *nzf);
626 MAT3( tcf, *nxf, *nyf, *nzf);
627 MAT3( ccf, *nxf, *nyf, *nzf);
628 MAT3( fcf, *nxf, *nyf, *nzf);
629
630 WARN_UNTESTED;
631
632 // How far to step into the coefficient arrays
633 iadd = (*nxf - 1) / (*nx - 1);
634 jadd = (*nyf - 1) / (*ny - 1);
635 kadd = (*nzf - 1) / (*nz - 1);
636
637 if (iadd != 2 || jadd != 2 || kadd != 2) {
638 Vnm_print(2, "Vbuildcopy0: Problem with grid dimensions...\n");
639 }
640
641 // Compute the coarse grid pde coefficients
642 for (k=1; k<=*nz; k++) {
643 kk = 2 * k - 1;
644 VAT(zc, k) = VAT(zf, kk);
645
646 for (j=1; j<=*ny; j++) {
647 jj = 2 * j - 1;
648 VAT(yc, j) = VAT(yf, jj);
649
650 for (i=1; i<=*nx; i++){
651 ii = 2 * i - 1;
652 VAT(xc, i) = VAT(xf, ii);
653
654 // True solution
655 VAT3( tc, i, j, k) = VAT3( tcf, ii, jj, kk);
656
657 // Helmholtz coefficient
658 VAT3( cc, i, j, k) = VAT3( ccf, ii, jj, kk);
659
660 // Source function
661 VAT3( fc, i, j, k) = VAT3( fcf, ii, jj, kk);
662
663 // East/West neighbor
664 VAT3(a1c, i, j, k) = VAT3(a1cf, ii, jj, kk);
665
666 // North/South neighbor
667 VAT3(a2c, i, j, k) = VAT3(a2cf, ii, jj, kk);
668
669 // Up/Down neighbor
670 VAT3(a3c, i, j, k) = VAT3(a3cf, ii, jj, kk);
671 }
672 }
673 }
674
675 // The (i=1) and (i=nx) boundaries
676 for (k=1; k<=*nz; k++) {
677 kk = 2 * k - 1;
678
679 for (j=1; j<=*ny; j++) {
680 jj = 2 * j - 1;
681
682 VAT3(gxc, j, k, 1) = VAT3(gxcf, jj, kk, 1);
683 VAT3(gxc, j, k, 2) = VAT3(gxcf, jj, kk, 2);
684 VAT3(gxc, j, k, 3) = VAT3(gxcf, jj, kk, 3);
685 VAT3(gxc, j, k, 4) = VAT3(gxcf, jj, kk, 4);
686 }
687 }
688
689 // The (j=1) and (j=ny) boundaries
690 for (k=1; k<=*nz; k++) {
691 kk = 2 * k - 1;
692
693 for (i=1; i<=*nx; i++) {
694 ii = 2 * i - 1;
695
696 VAT3(gyc, i, k, 1) = VAT3(gycf, ii, kk, 1);
697 VAT3(gyc, i, k, 2) = VAT3(gycf, ii, kk, 2);
698 VAT3(gyc, i, k, 3) = VAT3(gycf, ii, kk, 3);
699 VAT3(gyc, i, k, 4) = VAT3(gycf, ii, kk, 4);
700 }
701 }
702
703 // The (k=1) and (k=nz) boundaries
704 for (j=1; j<=*ny; j++) {
705 jj = 2 * j - 1;
706
707 for (i=1; i<=*nx; i++) {
708 ii = 2 * i - 1;
709
710 VAT3(gzc, i, j, 1) = VAT3(gzcf, ii, jj, 1);
711 VAT3(gzc, i, j, 2) = VAT3(gzcf, ii, jj, 2);
712 VAT3(gzc, i, j, 3) = VAT3(gzcf, ii, jj, 3);
713 VAT3(gzc, i, j, 4) = VAT3(gzcf, ii, jj, 4);
714 }
715 }
716}
717
718VPUBLIC void Vbuildharm0(int *nx, int *ny, int *nz,
719 int *nxf, int *nyf, int *nzf,
720 double *xc, double *yc, double *zc,
721 double *gxc, double *gyc, double *gzc,
722 double *a1c, double *a2c, double *a3c,
723 double *cc, double *fc, double *tc,
724 double *xf, double *yf, double *zf,
725 double *gxcf, double *gycf, double *gzcf,
726 double *a1cf, double *a2cf, double *a3cf,
727 double *ccf, double *fcf, double *tcf) {
728#if 1
729 Vnm_print(2, "WARNING: FUNCTION IS NOT FULLY IMPLEMENTED YET!!!");
730#else
731 int i, j, k;
732 int ii, jj, kk;
733 int iadd, jadd, kadd;
734
735 MAT3( gxc, *ny, *nz, 4);
736 MAT3( gyc, *nx, *nz, 4);
737 MAT3( gzc, *nx, *ny, 4);
738
739 MAT3( a1c, *nx, *ny, *nz);
740 MAT3( a2c, *nx, *ny, *nz);
741 MAT3( a3c, *nx, *ny, *nz);
742
743 MAT3( cc, *nx, *ny, *nz);
744 MAT3( fc, *nx, *ny, *nz);
745 MAT3( tc, *nx, *ny, *nz);
746
747 MAT3(gxcf, *nyf, *nzf, 4);
748 MAT3(gycf, *nxf, *nzf, 4);
749 MAT3(gzcf, *nxf, *nyf, 4);
750 MAT3(a1cf, *nxf, *nyf, *nzf);
751 MAT3(a2cf, *nxf, *nyf, *nzf);
752 MAT3(a3cf, *nxf, *nyf, *nzf);
753 MAT3( tcf, *nxf, *nyf, *nzf);
754 MAT3( ccf, *nxf, *nyf, *nzf);
755 MAT3( fcf, *nxf, *nyf, *nzf);
756
757 // Statement functions
759 double a, b, c, d, e, f, g, h;
760
761 // How far to step into the coefficient arrays
762 iadd = (*nxf - 1) / (*nx - 1);
763 jadd = (*nyf - 1) / (*ny - 1);
764 kadd = (*nzf - 1) / (*nz - 1);
765 if (iadd !=2 || jadd != 2 || kadd != 2) {
766 Vnm_print(2, "BUILDHARM0: problem with grid dimensions...\n");
767 }
768
769 // Compute the coarse grid pde coefficients
770 for (k=1; k<=*nz; k++) {
771 kk = 2 * k - 1;
772 VAT(zc, k) = VAT(zf, kk);
773
774 for (j=1; j<=*ny; j++) {
775 jj = 2 * j - 1;
776 VAT(yc, j) = VAT(yf,jj);
777
778 for (i=1; i<=*nx; i++) {
779 ii = 2 * i - 1;
780 VAT(xc, i) = VAT(xf, ii);
781
782 // True solution
783 VAT3(tc, i, j, k) = VAT3(tcf, ii, jj, kk);
784
785 // Helmholtz coefficient
786 VAT3(cc, i, j, k) = VAT3(ccf, ii, jj, kk);
787
788 /* Commented out in original fortran code
789 cc(i,j,k) = (
790 +0.5e0 * ccf(ii,jj,kk)
791 +0.5e0 * ARITH6( ccf(max0(1,ii-1),jj,kk),
792 ccf(min0(nxf,ii+1),jj,kk),
793 ccf(ii,max0(1,jj-1),kk),
794 ccf(ii,min0(nyf,jj+1),kk),
795 ccf(ii,jj,max0(nzf,kk-1)),
796 ccf(ii,jj,min0(nzf,kk+1)) )
797 )
798 */
799
800 // Source function
801 VAT3(fc, i, j, k) = VAT3(fcf, ii, jj, kk);
802 /*
803 fc(i,j,k) = (
804 +0.5e0 * fcf(ii,jj,kk)
805 +0.5e0 * ARITH6( fcf(max0(1,ii-1),jj,kk),
806 fcf(min0(nxf,ii+1),jj,kk),
807 fcf(ii,max0(1,jj-1),kk),
808 fcf(ii,min0(nyf,jj+1),kk),
809 fcf(ii,jj,max0(nzf,kk-1)),
810 fcf(ii,jj,min0(nzf,kk+1)) )
811 )
812 */
813
814 // East/West neighbor
815 VAT3(a1c, i, j, k) = (
816 +0.500 * HARMO2(VAT3(a1cf, ii, jj, kk),
817 VAT3(a1cf, VMIN2(*nxf, ii+1), jj, kk))
818 +0.125 * HARMO2(VAT3(a1cf, ii, jj, VMAX2(1, kk-1)),
819 VAT3(a1cf, VMIN2(*nxf, ii+1), jj, VMAX2(1, kk-1)))
820 +0.125 * HARMO2(VAT3(a1cf, ii, jj, VMIN2(*nzf, kk+1)),
821 VAT3(a1cf, VMIN2(*nxf, ii+1), jj, VMIN2(*nzf, kk+1)))
822 +0.125 * HARMO2(VAT3(a1cf, ii, VMAX2(1, jj-1), kk),
823 VAT3(a1cf, VMIN2(*nxf, ii+1), VMAX2(1, jj-1), kk))
824 +0.125 * HARMO2(VAT3(a1cf, ii, VMIN2(*nyf, jj+1), kk),
825 VAT3(a1cf, VMIN2(*nxf, ii+1), VMIN2(*nyf, jj+1), kk))
826 );
827
828 // North/South neighbor
829 VAT3(a2c, i, j, k) = (
830 +0.500 * HARMO2(VAT3(a2cf, ii, jj, kk),
831 VAT3(a2cf, ii, VMIN2(*nyf, jj+1), kk))
832 +0.125 * HARMO2(VAT3(a2cf, ii, jj, VMAX2(1, kk-1)),
833 VAT3(a2cf, ii, VMIN2(*nyf, jj+1), VMAX2(1, kk-1)))
834 +0.125 * HARMO2(VAT3(a2cf, ii, jj, VMIN2(*nzf, kk+1)),
835 VAT3(a2cf, ii, VMIN2(*nyf, jj+1), VMIN2(*nzf, kk+1)))
836 +0.125 * HARMO2(VAT3(a2cf, VMAX2(1, ii-1), jj, kk),
837 VAT3(a2cf, VMAX2(1, ii-1), VMIN2(nyf, jj+1), kk))
838 +0.125 * HARMO2(VAT3(a2cf, VMIN2(*nxf, ii+1), jj, kk),
839 VAT3(a2cf, VMIN2(*nxf, ii+1), VMIN2(*nyf, jj+1), kk))
840 );
841
842 // Up/Down neighbor
843 VAT3(a3c, i, j, k) = (
844 +0.500 * HARMO2(VAT3(a3cf, ii, jj, kk),
845 VAT3(a3cf, ii, jj, VMIN2(*nzf, kk+1)))
846 +0.125 * HARMO2(VAT3(a3cf, ii, VMAX2(1, jj-1), kk),
847 VAT3(a3cf, ii, VMAX2(1, jj-1), VMIN2(*nzf, kk+1)))
848 +0.125 * HARMO2(VAT3(a3cf, ii, VMIN2(*nyf, jj+1), kk),
849 VAT3(a3cf, ii, VMIN2(*nyf, jj+1), VMIN2(*nzf, kk+1)))
850 +0.125 * HARMO2(VAT3(a3cf, VMAX2(1, ii-1), jj, kk),
851 VAT3(a3cf, VMAX2(1, ii-1), jj, VMIN2(*nzf, kk+1)))
852 +0.125 * HARMO2(VAT3(a3cf, VMIN2(*nxf, ii+1), jj, kk),
853 VAT3(a3cf, VMIN2(*nxf, ii+1), jj, VMIN2(*nzf, kk+1)))
854 );
855 }
856 }
857 }
858
859 // The (i=1) and (i=nx) boundaries ***
860 for (k=1; k<=*nz; k++) {
861 kk = 2 * k - 1;
862
863 for (j=1; j<=*ny; j++) {
864 jj = 2 * j - 1;
865
866 VAT3(gxc, j, k, 1) = VAT3(gxcf, jj, kk, 1);
867 VAT3(gxc, j, k, 2) = VAT3(gxcf, jj, kk, 2);
868 VAT3(gxc, j, k, 3) = VAT3(gxcf, jj, kk, 3);
869 VAT3(gxc, j, k, 4) = VAT3(gxcf, jj, kk, 4);
870 }
871 }
872
873 // The (j=1) and (j=ny) boundaries
874 for (k=1; k<=*nz; k++) {
875 kk = 2 * k - 1;
876
877 for (i=1; i<=*nx; i++) {
878 ii = 2 * i - 1;
879 VAT3(gyc, i, k, 1) = VAT3(gycf, ii, kk, 1);
880 VAT3(gyc, i, k, 2) = VAT3(gycf, ii, kk, 2);
881 VAT3(gyc, i, k, 3) = VAT3(gycf, ii, kk, 3);
882 VAT3(gyc, i, k, 4) = VAT3(gycf, ii, kk, 4);
883 }
884 }
885
886 // The (k=1) and (k=nz) boundaries
887 for (j=1; j<=*ny; j++) {
888 jj = 2 * j - 1;
889
890 for (i=1; i<=*nx; i++) {
891 ii = 2 * i - 1;
892
893 VAT3(gzc, i, j, 1) = VAT3(gzcf, ii, jj, 1);
894 VAT3(gzc, i, j, 2) = VAT3(gzcf, ii, jj, 2);
895 VAT3(gzc, i, j, 3) = VAT3(gzcf, ii, jj, 3);
896 VAT3(gzc, i, j, 4) = VAT3(gzcf, ii, jj, 4);
897 }
898 }
899#endif
900}
901
902
903
904VPUBLIC void Vbuildalg(int *nx, int *ny, int *nz,
905 int *mode, int *nlev, int *iz,
906 int *ipc, double *rpc,
907 double *ac, double *cc, double *fc,
908 double *x, double *y, double *tmp) {
909
910 int nxx, nyy, nzz;
911 int nxold, nyold, nzold;
912 int lev, numlev;
913
914 MAT2(iz, 50, *nlev);
915
916 // Setup
917 nxx = *nx;
918 nyy = *ny;
919 nzz = *nz;
920
921 // Build the rhs the finest level
922 lev = 1;
923 if ((*mode == 1) || (*mode == 2)) {
924 Vnmatvec(&nxx, &nyy, &nzz,
925 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
926 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
927 RAT( x, VAT2(iz, 1, lev)), RAT( y, VAT2(iz, 1, lev)),
928 tmp);
929 } else {
930 Vmatvec(&nxx, &nyy, &nzz,
931 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
932 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
933 RAT( x, VAT2(iz, 1, lev)), RAT( y, VAT2(iz, 1, lev)));
934 }
935
936 // Build the (nlev-1) level rhs function
937 for (lev=2; lev <= *nlev; lev++) {
938 nxold = nxx;
939 nyold = nyy;
940 nzold = nzz;
941
942 numlev = 1;
943 Vmkcors(&numlev, &nxold, &nyold, &nzold, &nxx, &nyy, &nzz);
944
945 // Build the rhs on this level
946 if ((*mode == 1) || (*mode == 2)) {
947 Vnmatvec(&nxx, &nyy, &nzz,
948 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
949 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
950 RAT( x, VAT2(iz, 1, lev)), RAT( y, VAT2(iz, 1, lev)),
951 tmp);
952 } else {
953 Vmatvec(&nxx, &nyy, &nzz,
954 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
955 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
956 RAT( x, VAT2(iz, 1, lev)), RAT( y, VAT2(iz, 1, lev)));
957 }
958 }
959}
VPUBLIC void Vbuildops(int *nx, int *ny, int *nz, int *nlev, int *ipkey, int *iinfo, int *ido, int *iz, int *mgprol, int *mgcoar, int *mgsolv, int *mgdisc, int *ipc, double *rpc, double *pc, double *ac, double *cc, double *fc, double *xf, double *yf, double *zf, double *gxcf, double *gycf, double *gzcf, double *a1cf, double *a2cf, double *a3cf, double *ccf, double *fcf, double *tcf)
Build operators, boundary arrays, modify affine vectors ido==0: do only fine level ido==1: do only co...
Definition mgsubd.c:57
VPUBLIC void Vbuildband(int *key, int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, int *ipcB, double *rpcB, double *acB)
Banded matrix builder.
Definition buildBd.c:57
#define HARMO2(a, b)
Multigrid subroutines.
Definition mgsubd.h:65
VPUBLIC void Vmatvec(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *x, double *y)
Matrix-vector multiplication routines.
Definition matvecd.c:57
VPUBLIC void VbuildA(int *nx, int *ny, int *nz, int *ipkey, int *mgdisc, int *numdia, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *xf, double *yf, double *zf, double *gxcf, double *gycf, double *gzcf, double *a1cf, double *a2cf, double *a3cf, double *ccf, double *fcf)
Build the Laplacian.
Definition buildAd.c:57
VPUBLIC void Vprtmatd(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac)
Definition mikpckd.c:332
VPUBLIC void VbuildP(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, int *mgprol, int *ipc, double *rpc, double *pc, double *ac, double *xf, double *yf, double *zf)
Builds prolongation matrix.
Definition buildPd.c:57
VPUBLIC void Vbuildstr(int *nx, int *ny, int *nz, int *nlev, int *iz)
Build the nexted operator framework in the array iz.
Definition mgsubd.c:257
VPUBLIC void Vbuildgaler0(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, int *ipkey, int *numdia, double *pcFF, int *ipcFF, double *rpcFF, double *acFF, double *ccFF, double *fcFF, int *ipc, double *rpc, double *ac, double *cc, double *fc)
Form the Galerkin coarse grid system.
Definition mgsubd.c:341
VPUBLIC void VbuildG(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, int *numdia, double *pcFF, double *acFF, double *ac)
Build Galerkin matrix structures.
Definition buildGd.c:57
VPUBLIC void Vrestrc(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, double *xin, double *xout, double *pc)
Apply the restriction operator.
Definition matvecd.c:782
VPUBLIC void Vextrac(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, double *xin, double *xout)
Simple injection of a fine grid function into coarse grid.
Definition matvecd.c:1078
VEXTERNC void Vnmatvec(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *x, double *y, double *w1)
Break the matrix data-structure into diagonals and then call the matrix-vector routine.
Definition matvecd.c:232
VPUBLIC void Vpackmg(int *iparm, double *rparm, size_t *nrwk, int *niwk, int *nx, int *ny, int *nz, int *nlev, int *nu1, int *nu2, int *mgkey, int *itmax, int *istop, int *ipcon, int *nonlin, int *mgsmoo, int *mgprol, int *mgcoar, int *mgsolv, int *mgdisc, int *iinfo, double *errtol, int *ipkey, double *omegal, double *omegan, int *irite, int *iperf)
Print out a column-compressed sparse matrix in Harwell-Boeing format.
Definition mgsubd.c:555