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