APBS 3.0.0
Loading...
Searching...
No Matches
mikpckd.c
1
55#include "mikpckd.h"
56
57VPUBLIC 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
75VPUBLIC 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
91VPUBLIC 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
112VPUBLIC 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
131VPUBLIC 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
152VPUBLIC 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
173VPUBLIC 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
195VPUBLIC 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
209VPUBLIC 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
258VPUBLIC 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
291VPUBLIC 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
318VPUBLIC 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
332VPUBLIC 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
361VPUBLIC 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
393VEXTERNC 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
445VPUBLIC 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