APBS  3.0.0
mgsubd.c
1 
55 #include "mgsubd.h"
56 
57 VPUBLIC 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 
257 VPUBLIC 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 
341 VPUBLIC 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 
377 VPUBLIC 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 
400 VPUBLIC 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 
422 VPUBLIC 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 
448 VPUBLIC void Vfiner(int *nold, int *nnew) {
449  // Find the coarser grid size ***
450  *nnew = (*nold - 1) * 2 + 1;
451 }
452 
453 
454 
455 VPUBLIC 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 
471 VPUBLIC 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 
503 VPUBLIC 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 
555 VPUBLIC 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 
595 VEXTERNC 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 
718 VPUBLIC 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 
904 VPUBLIC 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