APBS  3.0.0
mikpckd.c
1 
55 #include "mikpckd.h"
56 
57 VPUBLIC void Vxcopy(int *nx, int *ny, int *nz, double *x, double *y) {
58 
59  MAT3(x, *nx, *ny, *nz);
60  MAT3(y, *nx, *ny, *nz);
61 
62  // The indices used to traverse the matrices
63  int i, j, k;
64 
66  #pragma omp parallel for private(i,j,k)
67  for(k=2; k<=*nz-1; k++)
68  for(j=2; j<=*ny-1; j++)
69  for(i=2; i<=*nx-1; i++)
70  VAT3(y, i, j, k) = VAT3(x, i, j, k);
71 }
72 
73 
74 
75 VPUBLIC void Vxcopy_small(int *nx, int *ny, int *nz, double *x, double *y) {
76 
77  MAT3(x, *nx, *ny, *nz);
78  MAT3(y, *nx - 2, *ny - 2, *nz - 2);
79 
80  // The indices used to traverse the matrices
81  int i, j, k;
82 
83  for(k=2; k<=*nz-1; k++)
84  for(j=2; j<=*ny-1; j++)
85  for(i=2; i<=*nx-1; i++)
86  VAT3(y, i - 1, j - 1, k - 1) = VAT3(x, i, j, k);
87 }
88 
89 
90 
91 VPUBLIC void Vxcopy_large(int *nx, int *ny, int *nz, double *x, double *y) {
92 
98  MAT3(x, *nx - 2, *ny - 2, *nz - 2);
99  MAT3(y, *nx, *ny, *nz);
100 
101  // The indices used to traverse the matrices
102  int i, j, k;
103 
104  for(k=2; k<=*nz-1; k++)
105  for(j=2; j<=*ny-1; j++)
106  for(i=2; i<=*nx-1; i++)
107  VAT3(y, i, j, k) = VAT3(x, i - 1, j - 1, k - 1);
108 }
109 
110 
111 
112 VPUBLIC void Vxaxpy(int *nx, int *ny, int *nz,
113  double *alpha, double *x, double *y) {
114 
115  // Create the wrappers
116  MAT3(x, *nx, *ny, *nz);
117  MAT3(y, *nx, *ny, *nz);
118 
119  // The indices used to traverse the matrices
120  int i, j, k;
121 
123  for(k=2; k<=*nz-1; k++)
124  for(j=2; j<=*ny-1; j++)
125  for(i=2; i<=*nx-1; i++)
126  VAT3(y, i, j, k) += *alpha * VAT3(x, i, j, k);
127 }
128 
129 
130 
131 VPUBLIC double Vxnrm1(int *nx, int *ny, int *nz,
132  double *x) {
133 
134  double xnrm1 = 0.0;
135 
136  MAT3(x, *nx, *ny, *nz);
137 
138  // The indices used to traverse the matrices
139  int i, j, k;
140 
142  for(k=2; k<=*nz-1; k++)
143  for(j=2; j<=*ny-1; j++)
144  for(i=2; i<=*nx-1; i++)
145  xnrm1 += VABS(VAT3(x, i, j, k));
146 
147  return xnrm1;
148 }
149 
150 
151 
152 VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz,
153  double *x) {
154 
155  double xnrm2 = 0.0;
156 
157  MAT3(x, *nx, *ny, *nz);
158 
159  // The indices used to traverse the matrices
160  int i, j, k;
161 
163  for(k=2; k<=*nz-1; k++)
164  for(j=2; j<=*ny-1; j++)
165  for(i=2; i<=*nx-1; i++)
166  xnrm2 += VAT3(x, i, j, k) * VAT3(x, i, j, k);
167 
168  return VSQRT(xnrm2);
169 }
170 
171 
172 
173 VPUBLIC double Vxdot(int *nx, int *ny, int *nz,
174  double *x, double *y) {
175 
176  int i, j, k;
177 
178  // Initialize
179  double xdot = 0.0;
180 
181  MAT3(x, *nx, *ny, *nz);
182  MAT3(y, *nx, *ny, *nz);
183 
184  // Do it
185  for(k=2; k<=*nz-1; k++)
186  for(j=2; j<=*ny-1; j++)
187  for(i=2; i<=*nx-1; i++)
188  xdot += VAT3(x, i, j, k) * VAT3(y, i, j, k);
189 
190  return xdot;
191 }
192 
193 
194 
195 VPUBLIC void Vazeros(int *nx, int *ny, int *nz, double *x) {
196 
197  int i, n;
198  int nproc = 1;
199 
200  n = *nx * *ny * *nz;
201 
202  #pragma omp parallel for private(i)
203  for (i=1; i<=n; i++)
204  VAT(x, i) = 0.0;
205 }
206 
207 
208 
209 VPUBLIC void VfboundPMG(int *ibound, int *nx, int *ny, int *nz,
210  double *x, double *gxc, double *gyc, double *gzc) {
211 
212  int i, j, k;
213 
214  // Create and bind the wrappers for the source data
215  MAT3( x, *nx, *ny, *nz);
216  MAT3(gxc, *ny, *nz, 2);
217  MAT3(gyc, *nx, *nz, 2);
218  MAT3(gzc, *nx, *ny, 2);
219 
220  // Dirichlet test
221  if (ibound == 0) {
222 
223  // Dero dirichlet
224  VfboundPMG00(nx, ny, nz, x);
225 
226  } else {
227 
228  // Nonzero dirichlet
229 
230  // The (i=1) and (i=nx) boundaries
231  for (k=1; k<=*nz; k++) {
232  for (j=1; j<=*ny; j++) {
233  VAT3(x, 1, j, k) = VAT3(gxc, j, k, 1);
234  VAT3(x, *nx, j, k) = VAT3(gxc, j, k, 2);
235  }
236  }
237 
238  // The (j=1) and (j=ny) boundaries
239  for (k=1; k<=*nz; k++) {
240  for (i=1; i<=*nx; i++) {
241  VAT3(x, i, 1 ,k) = VAT3(gyc, i, k, 1);
242  VAT3(x, i, *ny, k) = VAT3(gyc, i, k, 2);
243  }
244  }
245 
246  // The (k=1) and (k=nz) boundaries
247  for (j=1; j<=*ny; j++) {
248  for (i=1; i<=*nx; i++) {
249  VAT3(x, i, j, 1) = VAT3(gzc, i, j, 1);
250  VAT3(x, i, j, *nz) = VAT3(gzc, i, j, 2);
251  }
252  }
253  }
254 }
255 
256 
257 
258 VPUBLIC void VfboundPMG00(int *nx, int *ny, int *nz, double *x) {
259 
260  int i, j, k;
261 
262  MAT3( x, *nx, *ny, *nz);
263 
264  // The (i=1) and (i=nx) boundaries
265  for (k=1; k<=*nz; k++) {
266  for (j=1; j<=*ny; j++) {
267  VAT3(x, 1, j, k) = 0.0;
268  VAT3(x, *nx, j, k) = 0.0;
269  }
270  }
271 
272  // The (j=1) and (j=ny) boundaries
273  for (k=1; k<=*nz; k++) {
274  for(i=1; i<=*nx; i++) {
275  VAT3(x, i, 1, k) = 0.0;
276  VAT3(x, i, *ny, k) = 0.0;
277  }
278  }
279 
280  // The (k=1) and (k=nz) boundaries
281  for (j=1; j<=*ny; j++) {
282  for (i=1; i<=*nx; i++) {
283  VAT3(x, i, j, 1) = 0.0;
284  VAT3(x, i, j, *nz) = 0.0;
285  }
286  }
287 }
288 
289 
290 
291 VPUBLIC void Vaxrand(int *nx, int *ny, int *nz, double *x) {
292 
293  int n, i, ii, ipara, ivect, iflag;
294  int nproc = 1;
295  double xdum;
296 
297  WARN_UNTESTED;
298 
299  // Find parallel loops (ipara), remainder (ivect)
300  n = *nx * *ny * *nz;
301  ipara = n / nproc;
302  ivect = n % nproc;
303  iflag = 1;
304  xdum = (double)(VRAND);
305 
306  // Do parallel loops
307  for (ii=1; ii<=nproc; ii++)
308  for (i=1+(ipara*(ii-1)); i<=ipara*ii; i++)
309  VAT(x, i) = (double)(VRAND);
310 
311  // Do vector loops
312  for (i=ipara*nproc+1; i<=n; i++)
313  VAT(x, i) = (double)(VRAND);
314 }
315 
316 
317 
318 VPUBLIC void Vxscal(int *nx, int *ny, int *nz, double *fac, double *x) {
319 
320  int i, j, k;
321 
322  MAT3(x, *nx, *ny, *nz);
323 
324  for (k=2; k<=*nz-1; k++)
325  for (j=2; j<=*ny-1; j++)
326  for (i=2; i<=*nx-1; i++)
327  VAT3(x, i, j, k) *= *fac;
328 }
329 
330 
331 
332 VPUBLIC void Vprtmatd(int *nx, int *ny, int *nz,
333  int *ipc, double *rpc, double *ac) {
334 
335  int numdia;
336 
337  MAT2(ac, *nx * *ny * *nz, 1);
338 
339  WARN_UNTESTED;
340 
341  // Do the printing
342  numdia = VAT(ipc, 11);
343  if (numdia == 7) {
344  Vprtmatd7(nx, ny, nz,
345  ipc, rpc,
346  RAT2(ac, 1, 1), RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4));
347  } else if (numdia == 27) {
348  Vprtmatd27(nx, ny, nz,
349  ipc, rpc,
350  RAT2(ac, 1, 1), RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
351  RAT2(ac, 1, 5), RAT2(ac, 1, 6),
352  RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1, 10),
353  RAT2(ac, 1, 11), RAT2(ac, 1, 12), RAT2(ac, 1, 13), RAT2(ac, 1, 14));
354  } else {
355  Vnm_print(2, "Vprtmatd: invalid stencil type given: %d\n", numdia);
356  }
357 }
358 
359 
360 
361 VPUBLIC void Vprtmatd7(int *nx, int *ny, int *nz,
362  int *ipc, double *rpc,
363  double *oC, double *oE, double *oN, double *uC) {
364 
365  int n, i, j, k;
366 
367  MAT3(oC, *nx, *ny, *nz);
368  MAT3(oE, *nx, *ny, *nz);
369  MAT3(oN, *nx, *ny, *nz);
370  MAT3(uC, *nx, *ny, *nz);
371 
372  WARN_UNTESTED;
373 
374  // Recover matrix dimension
375  n = (*nx - 2) * (*ny - 2) * (*nz - 2);
376 
377  Vnm_print(2, "Vprtmatd7: Dimension of matrix = %d\n", n);
378  Vnm_print(2, "Vprtmatd7: Begin diagonal matrix\n");
379 
380  // Handle first block
381  for (k=2; k<=*nz-1; k++)
382  for (j=2; j<=*ny-1; j++)
383  for (i=2; i<=*nx-1; i++)
384  Vnm_print(2, "Vprtmatd7: (%d,%d,%d) - oC=%g, oE=%g, oN=%g, uC=%g\n",
385  VAT3(oC,i,j,k), VAT3(oE,i,j,k), VAT3(oN,i,j,k), VAT3(uC,i,j,k));
386 
387  // Finish up
388  Vnm_print(2, "Vprtmatd7: End diagonal matrix\n");
389 }
390 
391 
392 
393 VEXTERNC void Vprtmatd27(int *nx, int *ny, int *nz,
394  int *ipc, double *rpc,
395  double *oC, double *oE, double *oN, double *uC,
396  double *oNE, double *oNW,
397  double *uE, double *uW, double *uN, double *uS,
398  double *uNE, double *uNW, double *uSE, double *uSW) {
399 
400  int n, i, j, k;
401 
402  MAT3( oC, *nx, *ny, *nz);
403  MAT3( oE, *nx, *ny, *nz);
404  MAT3( oN, *nx, *ny, *nz);
405  MAT3( uC, *nx, *ny, *nz);
406  MAT3(oNE, *nx, *ny, *nz);
407  MAT3(oNW, *nx, *ny, *nz);
408  MAT3( uE, *nx, *ny, *nz);
409  MAT3( uW, *nx, *ny, *nz);
410  MAT3( uN, *nx, *ny, *nz);
411  MAT3( uS, *nx, *ny, *nz);
412  MAT3(uNE, *nx, *ny, *nz);
413  MAT3(uNW, *nx, *ny, *nz);
414  MAT3(uSE, *nx, *ny, *nz);
415  MAT3(uSW, *nx, *ny, *nz);
416 
417  WARN_UNTESTED;
418 
419  // Recover matrix dimension
420  n = (*nx - 2) * (*ny - 2) * (*nz - 2);
421 
422  Vnm_print(2, "Vprtmatd27: Dimension of matrix = %d\n", n);
423  Vnm_print(2, "Vprtmatd27: Begin diagonal matrix\n");
424 
425  // Handle first block
426  for (k=2; k<=*nz-1; k++)
427  for (j=2; j<=*ny-1; j++)
428  for (i=2; i<=*nx-1; i++)
429  Vnm_print(2, "Vprtmatd7: (%d,%d,%d) - oC = %g, oE = %g, "
430  "oNW = %g, oN = %g, oNE = %g, uSW = %g, uS = %g, "
431  "uSE = %g, uW = %g, uC = %g, uE = %g, uNW = %g, "
432  "uN = %g, uNE = %g\n",
433  VAT3( oC, i, j, k), VAT3( oE, i, j, k),
434  VAT3(oNW, i, j, k), VAT3( oN, i, j, k),
435  VAT3(oNE, i, j, k), VAT3(uSW, i, j, k),
436  VAT3( uS, i, j, k), VAT3(uSE, i, j, k),
437  VAT3( uW, i, j, k), VAT3( uC, i, j, k),
438  VAT3( uE, i, j, k), VAT3(uNW, i, j, k),
439  VAT3( uN, i, j, k), VAT3(uNE, i, j, k) );
440 
441  // Finish up
442  Vnm_print(2, "Vprtmatd27: End diagonal matrix\n");
443 }
444 
445 VPUBLIC void Vlinesearch(int *nx, int *ny, int *nz,
446  double *alpha,
447  int *ipc, double *rpc,
448  double *ac, double *cc, double *fc,
449  double *p, double *x, double *r,
450  double *ap, double *zk, double *zkp1) {
451  VABORT_MSG0("Not translated yet");
452 }
VPUBLIC void VfboundPMG(int *ibound, int *nx, int *ny, int *nz, double *x, double *gxc, double *gyc, double *gzc)
Initialize a grid function to have a certain boundary value,.
Definition: mikpckd.c:209
VPUBLIC void Vaxrand(int *nx, int *ny, int *nz, double *x)
Fill grid function with random values, including boundary values.
Definition: mikpckd.c:291
VPUBLIC double Vxdot(int *nx, int *ny, int *nz, double *x, double *y)
Inner product operation for a grid function with boundary values.
Definition: mikpckd.c:173
VPUBLIC void Vprtmatd(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac)
Definition: mikpckd.c:332
VPUBLIC void Vxcopy_small(int *nx, int *ny, int *nz, double *x, double *y)
Copy operation for a grid function with boundary values. Quite simply copies one 3d matrix to another...
Definition: mikpckd.c:75
VPUBLIC void Vxscal(int *nx, int *ny, int *nz, double *fac, double *x)
Scale operation for a grid function with boundary values.
Definition: mikpckd.c:318
VPUBLIC void Vxcopy(int *nx, int *ny, int *nz, double *x, double *y)
A collection of useful low-level routines (timing, etc).
Definition: mikpckd.c:57
VPUBLIC void VfboundPMG00(int *nx, int *ny, int *nz, double *x)
Initialize a grid function to have a zero boundary value.
Definition: mikpckd.c:258
VPUBLIC void Vxaxpy(int *nx, int *ny, int *nz, double *alpha, double *x, double *y)
saxpy operation for a grid function with boundary values.
Definition: mikpckd.c:112
VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
Definition: mikpckd.c:152
VPUBLIC void Vazeros(int *nx, int *ny, int *nz, double *x)
Zero out operation for a grid function, including boundary values.
Definition: mikpckd.c:195
VPUBLIC double Vxnrm1(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
Definition: mikpckd.c:131
VPUBLIC void Vxcopy_large(int *nx, int *ny, int *nz, double *x, double *y)
Copy operation for a grid function with boundary values. Quite simply copies one 3d matrix to another...
Definition: mikpckd.c:91