APBS  3.0.0
buildBd.c
1 
55 #include "buildBd.h"
56 
57 VPUBLIC 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 
119 VPUBLIC 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 
183 VPUBLIC 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