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