APBS 3.0.0
Loading...
Searching...
No Matches
buildBd.c
1
55#include "buildBd.h"
56
57VPUBLIC void Vbuildband(int *key, int *nx, int *ny, int *nz,
58 int *ipc, double *rpc, double *ac,
59 int *ipcB, double *rpcB, double *acB) {
60
61 int numdia;
62 int n, m;
63 int lda, info;
64
65 MAT2(ac, *nx * *ny * *nz, 1);
66
67 // Do in one step
68 numdia = VAT(ipc, 11);
69 if (numdia == 7) {
70
71 n = (*nx - 2) * (*ny - 2) * (*nz - 2);
72 m = (*nx - 2) * (*ny - 2);
73 lda = m + 1;
74
76 (nx, ny, nz,
77 ipc, rpc,
78 RAT2(ac, 1, 1), RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
79 ipcB, rpcB, acB,
80 &n, &m, &lda);
81
82 } else if (numdia == 27) {
83
84 n = (*nx - 2) * (*ny - 2) * (*nz - 2);
85 m = (*nx - 2) * (*ny - 2) + (*nx - 2) + 1;
86 lda = m + 1;
87
89 (nx, ny, nz,
90 ipc, rpc,
91 RAT2(ac, 1, 1), RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
92 RAT2(ac, 1, 5), RAT2(ac, 1, 6),
93 RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1, 10),
94 RAT2(ac, 1, 11), RAT2(ac, 1, 12), RAT2(ac, 1, 13), RAT2(ac, 1, 14),
95 ipcB, rpcB, acB,
96 &n, &m, &lda);
97 } else {
98 Vnm_print(2, "Vbuildband: invalid stencil type given...");
99 }
100
101 // Factor the system
102 *key = 0;
103 info = 0;
104
105 Vdpbfa(acB, &lda, &n, &m, &info);
106 VAT(ipcB, 4) = 1;
107
108 if (info != 0) {
109
110 Vnm_print(2, "Vbuildband: dpbfa problem: %d\n", info);
111 Vnm_print(2, "Vbuildband: leading principle minor not PD...\n");
112
113 *key = 1;
114 }
115}
116
117
118
119VPUBLIC void Vbuildband1_7(int *nx, int *ny, int *nz,
120 int *ipc, double *rpc,
121 double *oC, double *oE, double *oN, double *uC,
122 int *ipcB, double *rpcB, double *acB,
123 int *n, int *m, int *lda) {
124
125 int i, j, k;
126 int ii, jj, kk;
127
128 MAT2(acB, *lda, *ny-1);
129
130 MAT3(oC, *nx, *ny, *nz);
131 MAT3(oE, *nx, *ny, *nz);
132 MAT3(oN, *nx, *ny, *nz);
133 MAT3(uC, *nx, *ny, *nz);
134
135 WARN_UNTESTED;
136
137 // Do it
138 VAT(ipcB, 1) = *n;
139 VAT(ipcB, 2) = *m;
140 VAT(ipcB, 3) = *lda;
141 VAT(ipcB, 4) = 0;
142
143 jj = 0;
144
145 //fprintf(data, "%s\n", PRINT_FUNC);
146
147 for (k=2; k<=*nz-1; k++) {
148
149 for (j=2; j<=*ny-1; j++) {
150
151 for (i=2; i<=*nx-1; i++) {
152 jj++;
153
154 // Diagonal term
155 ii = jj;
156 kk = ii - jj + *m + 1;
157
158 VAT2(acB, kk, jj) = VAT3(oC, i, j, k);
159
160 // East neighbor
161 ii = jj - 1;
162 kk = ii - jj + *m + 1;
163 VAT2(acB, kk, jj) = - VAT3(oE, i-1, j, k);
164
165 // North neighbor
166 ii = jj - (*nx - 2);
167 kk = ii - jj + *m + 1;
168 VAT2(acB, kk, jj) = - VAT3(oN, i, j-1, k);
169
170 // Up neighbor ***
171 ii = jj - (*nx - 2) * (*ny - 2);
172 kk = ii - jj + *m + 1;
173 VAT2(acB, kk, jj) = - VAT3(uC, i, j, k-1);
174
175 //fprintf(data, "%19.12E\n", VAT2(acB, kk, jj));
176 }
177 }
178 }
179}
180
181
182
183VPUBLIC void Vbuildband1_27(int *nx, int *ny, int *nz,
184 int *ipc, double *rpc,
185 double *oC, double *oE, double *oN, double *uC,
186 double *oNE, double *oNW,
187 double *uE, double *uW, double *uN, double *uS,
188 double *uNE, double *uNW, double *uSE, double *uSW,
189 int *ipcB, double *rpcB, double *acB,
190 int *n, int *m, int *lda) {
191
192 int i, j, k;
193 int ii, jj, kk;
194
195 MAT2(acB, *lda, *ny-1);
196
197 MAT3( oC, *nx, *ny, *nz);
198 MAT3( oE, *nx, *ny, *nz);
199 MAT3( oN, *nx, *ny, *nz);
200 MAT3( uC, *nx, *ny, *nz);
201
202 MAT3(oNE, *nx, *ny, *nz);
203 MAT3(oNW, *nx, *ny, *nz);
204
205 MAT3( uE, *nx, *ny, *nz);
206 MAT3( uW, *nx, *ny, *nz);
207 MAT3( uN, *nx, *ny, *nz);
208 MAT3( uS, *nx, *ny, *nz);
209
210 MAT3(uNE, *nx, *ny, *nz);
211 MAT3(uNW, *nx, *ny, *nz);
212 MAT3(uSE, *nx, *ny, *nz);
213 MAT3(uSW, *nx, *ny, *nz);
214
215 // Do it
216 VAT(ipcB, 1) = *n;
217 VAT(ipcB, 2) = *m;
218 VAT(ipcB, 3) = *lda;
219 VAT(ipcB, 4) = 0;
220
221 jj = 0;
222
223 //fprintf(data, "%s\n", PRINT_FUNC);
224
225 for (k=2; k<=*nz-1; k++) {
226
227 for (j=2; j<=*ny-1; j++) {
228
229 for (i=2; i<=*nx-1; i++) {
230 jj++;
231
232 // Diagonal term
233 ii = jj;
234 kk = ii - jj + *m + 1;
235 VAT2(acB, kk, jj) = VAT3(oC, i, j, k);
236
237 // East neighbor
238 ii = jj - 1;
239 kk = ii - jj + *m + 1;
240 VAT2(acB, kk, jj) = - VAT3(oE, i-1, j, k);
241
242 // North neighbor
243 ii = jj - (*nx - 2);
244 kk = ii - jj + *m + 1;
245 VAT2(acB, kk, jj) = - VAT3(oN, i, j-1, k);
246
247 // North-east neighbor
248 ii = jj - (*nx - 2) + 1;
249 kk = ii - jj + *m + 1;
250 VAT2(acB, kk, jj) = - VAT3(oNE, i, j-1, k);
251
252 // North-west neighbor
253 ii = jj - (*nx - 2) - 1;
254 kk = ii - jj + *m + 1;
255 VAT2(acB, kk, jj) = - VAT3(oNW, i, j-1, k);
256
257 // Up neighbor
258 ii = jj - (*nx - 2) * (*ny - 2);
259 kk = ii - jj + *m + 1;
260 VAT2(acB, kk, jj) = - VAT3(uC, i, j, k-1);
261
262 // Up-east neighbor
263 ii = jj - (*nx - 2) * (*ny - 2) +1;
264 kk = ii - jj + *m + 1;
265 VAT2(acB, kk, jj) = - VAT3(uE, i, j, k-1);
266
267 // Up-west neighbor
268 ii = jj - (*nx - 2) * (*ny - 2) - 1;
269 kk = ii - jj + *m + 1;
270 VAT2(acB, kk, jj) = - VAT3(uW, i, j, k-1);
271
272 // Up-north neighbor
273 ii = jj - (*nx - 2) * (*ny - 2) + (*nx - 2);
274 kk = ii - jj + *m + 1;
275 VAT2(acB, kk, jj) = - VAT3(uN, i, j, k-1);
276
277 // Up-south neighbor
278 ii = jj - (*nx - 2) * (*ny - 2) - (*nx - 2);
279 kk = ii - jj + *m + 1;
280 VAT2(acB, kk, jj) = - VAT3(uS, i, j, k-1);
281
282 // Up-north-east neighbor
283 ii = jj - (*nx - 2) * (*ny - 2) + (*nx - 2) + 1;
284 kk = ii - jj + *m + 1;
285 VAT2(acB, kk, jj) = - VAT3(uNE, i, j, k-1);
286
287 // Up-north-west neighbor
288 ii = jj - (*nx - 2) * (*ny - 2) + (*nx - 2) - 1;
289 kk = ii - jj + *m + 1;
290 VAT2(acB, kk, jj) = - VAT3(uNW, i, j, k-1);
291
292 // Up-south-east neighbor
293 ii = jj - (*nx - 2) * (*ny - 2) - (*nx - 2) + 1;
294 kk = ii - jj + *m + 1;
295 VAT2(acB, kk, jj) = - VAT3(uSE, i, j, k-1);
296
297 // Up-south-west neighbor
298 ii = jj - (*nx - 2) * (*ny - 2) - (*nx - 2) - 1;
299 kk = ii - jj + *m + 1;
300 VAT2(acB, kk, jj) = - VAT3(uSW, i, j, k-1);
301
302 //fprintf(data, "%19.12E\n", VAT2(acB, kk, jj));
303 }
304 }
305 }
306}
VPUBLIC void Vbuildband1_7(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *oC, double *oE, double *oN, double *uC, int *ipcB, double *rpcB, double *acB, int *n, int *m, int *lda)
Build the operator in banded form given the 7-diagonal form.
Definition buildBd.c:119
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
VPUBLIC void Vbuildband1_27(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, int *ipcB, double *rpcB, double *acB, int *n, int *m, int *lda)
Build the operator in banded form given the 27-diagonal form.
Definition buildBd.c:183