1 /**
2 * @ingroup PMGC
3 * @author Tucker Beck [fortran ->c translation], Michael Holst [original]
4 * @brief
5 * @version $Id:
6 *
7 * @attention
8 * @verbatim
9 *
10 * APBS -- Adaptive Poisson-Boltzmann Solver
11 *
12 * Nathan A. Baker (nathan.baker@pnl.gov)
13 * Pacific Northwest National Laboratory
14 *
15 * Additional contributing authors listed in the code documentation.
16 *
17 * Copyright (c) 2010-2014 Battelle Memorial Institute. Developed at the Pacific Northwest National Laboratory, operated by Battelle Memorial Institute, Pacific Northwest Division for the U.S. Department Energy. Portions Copyright (c) 2002-2010, Washington University in St. Louis. Portions Copyright (c) 2002-2010, Nathan A. Baker. Portions Copyright (c) 1999-2002, The Regents of the University of California. Portions Copyright (c) 1995, Michael Holst.
18 * All rights reserved.
19 *
20 *
21 * Redistribution and use in source and binary forms, with or without
22 * modification, are permitted provided that the following conditions are met:
23 *
24 * - Redistributions of source code must retain the above copyright notice, this
25 * list of conditions and the following disclaimer.
26 *
27 * - Redistributions in binary form must reproduce the above copyright notice,
28 * this list of conditions and the following disclaimer in the documentation
29 * and/or other materials provided with the distribution.
30 *
31 * - Neither the name of Washington University in St. Louis nor the names of its
32 * contributors may be used to endorse or promote products derived from this
33 * software without specific prior written permission.
34 *
35 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
36 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
37 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
38 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
39 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
40 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
41 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
42 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
43 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
44 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
45 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
46 *
47 * @endverbatim
48 */
49
50 #include "mikpckd.h"
51
Vxcopy(int * nx,int * ny,int * nz,double * x,double * y)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
60 /// @todo: Once the refactor begins, this will need to be unrolled
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
Vxcopy_small(int * nx,int * ny,int * nz,double * x,double * y)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
Vxcopy_large(int * nx,int * ny,int * nz,double * x,double * y)86 VPUBLIC void Vxcopy_large(int *nx, int *ny, int *nz, double *x, double *y) {
87
88 /** @note This function is exactly equivalent to calling xcopy_small with
89 * the matrix arguments reversed.
90 * @todo Replace this. it's worthless
91 */
92
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
Vxaxpy(int * nx,int * ny,int * nz,double * alpha,double * x,double * y)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
117 /// @todo parallel optimization
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
Vxnrm1(int * nx,int * ny,int * nz,double * x)126 VPUBLIC double Vxnrm1(int *nx, int *ny, int *nz,
127 double *x) {
128
129 double xnrm1 = 0.0; ///< Accumulates the calculated normal value
130
131 MAT3(x, *nx, *ny, *nz);
132
133 // The indices used to traverse the matrices
134 int i, j, k;
135
136 /// @todo parallel optimization
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
Vxnrm2(int * nx,int * ny,int * nz,double * x)147 VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz,
148 double *x) {
149
150 double xnrm2 = 0.0; ///< Accumulates the calculated normal value
151
152 MAT3(x, *nx, *ny, *nz);
153
154 // The indices used to traverse the matrices
155 int i, j, k;
156
157 /// @todo parallel optimization
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
Vxdot(int * nx,int * ny,int * nz,double * x,double * y)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
Vazeros(int * nx,int * ny,int * nz,double * x)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
VfboundPMG(int * ibound,int * nx,int * ny,int * nz,double * x,double * gxc,double * gyc,double * gzc)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
VfboundPMG00(int * nx,int * ny,int * nz,double * x)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
Vaxrand(int * nx,int * ny,int * nz,double * x)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
Vxscal(int * nx,int * ny,int * nz,double * fac,double * x)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
Vprtmatd(int * nx,int * ny,int * nz,int * ipc,double * rpc,double * ac)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
Vprtmatd7(int * nx,int * ny,int * nz,int * ipc,double * rpc,double * oC,double * oE,double * oN,double * uC)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
Vprtmatd27(int * nx,int * ny,int * nz,int * ipc,double * rpc,double * oC,double * oE,double * oN,double * uC,double * oNE,double * oNW,double * uE,double * uW,double * uN,double * uS,double * uNE,double * uNW,double * uSE,double * uSW)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
Vlinesearch(int * nx,int * ny,int * nz,double * alpha,int * ipc,double * rpc,double * ac,double * cc,double * fc,double * p,double * x,double * r,double * ap,double * zk,double * zkp1)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 }
448