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