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) {
65 MAT2(ac, *nx * *ny * *nz, 1);
68 numdia = VAT(ipc, 11);
71 n = (*nx - 2) * (*ny - 2) * (*nz - 2);
72 m = (*nx - 2) * (*ny - 2);
78 RAT2(ac, 1, 1), RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
82 }
else if (numdia == 27) {
84 n = (*nx - 2) * (*ny - 2) * (*nz - 2);
85 m = (*nx - 2) * (*ny - 2) + (*nx - 2) + 1;
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),
98 Vnm_print(2,
"Vbuildband: invalid stencil type given...");
105 Vdpbfa(acB, &lda, &n, &m, &info);
110 Vnm_print(2,
"Vbuildband: dpbfa problem: %d\n", info);
111 Vnm_print(2,
"Vbuildband: leading principle minor not PD...\n");
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) {
195 MAT2(acB, *lda, *ny-1);
197 MAT3( oC, *nx, *ny, *nz);
198 MAT3( oE, *nx, *ny, *nz);
199 MAT3( oN, *nx, *ny, *nz);
200 MAT3( uC, *nx, *ny, *nz);
202 MAT3(oNE, *nx, *ny, *nz);
203 MAT3(oNW, *nx, *ny, *nz);
205 MAT3( uE, *nx, *ny, *nz);
206 MAT3( uW, *nx, *ny, *nz);
207 MAT3( uN, *nx, *ny, *nz);
208 MAT3( uS, *nx, *ny, *nz);
210 MAT3(uNE, *nx, *ny, *nz);
211 MAT3(uNW, *nx, *ny, *nz);
212 MAT3(uSE, *nx, *ny, *nz);
213 MAT3(uSW, *nx, *ny, *nz);
225 for (k=2; k<=*nz-1; k++) {
227 for (j=2; j<=*ny-1; j++) {
229 for (i=2; i<=*nx-1; i++) {
234 kk = ii - jj + *m + 1;
235 VAT2(acB, kk, jj) = VAT3(oC, i, j, k);
239 kk = ii - jj + *m + 1;
240 VAT2(acB, kk, jj) = - VAT3(oE, i-1, j, k);
244 kk = ii - jj + *m + 1;
245 VAT2(acB, kk, jj) = - VAT3(oN, i, j-1, k);
248 ii = jj - (*nx - 2) + 1;
249 kk = ii - jj + *m + 1;
250 VAT2(acB, kk, jj) = - VAT3(oNE, i, j-1, k);
253 ii = jj - (*nx - 2) - 1;
254 kk = ii - jj + *m + 1;
255 VAT2(acB, kk, jj) = - VAT3(oNW, i, j-1, k);
258 ii = jj - (*nx - 2) * (*ny - 2);
259 kk = ii - jj + *m + 1;
260 VAT2(acB, kk, jj) = - VAT3(uC, i, j, k-1);
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);
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);
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);
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);
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);
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);
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);
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);
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.
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.