53 int *ipc,
double *rpc,
double *ac,
54 int *ipcB,
double *rpcB,
double *acB) {
60 MAT2(ac, *nx * *ny * *nz, 1);
63 numdia = VAT(ipc, 11);
66 n = (*nx - 2) * (*ny - 2) * (*nz - 2);
67 m = (*nx - 2) * (*ny - 2);
73 RAT2(ac, 1, 1), RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
77 }
else if (numdia == 27) {
79 n = (*nx - 2) * (*ny - 2) * (*nz - 2);
80 m = (*nx - 2) * (*ny - 2) + (*nx - 2) + 1;
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),
93 Vnm_print(2,
"Vbuildband: invalid stencil type given...");
100 Vdpbfa(acB, &lda, &n, &m, &info);
105 Vnm_print(2,
"Vbuildband: dpbfa problem: %d\n", info);
106 Vnm_print(2,
"Vbuildband: leading principle minor not PD...\n");
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) {
123 MAT2(acB, *lda, *ny-1);
125 MAT3(oC, *nx, *ny, *nz);
126 MAT3(oE, *nx, *ny, *nz);
127 MAT3(oN, *nx, *ny, *nz);
128 MAT3(uC, *nx, *ny, *nz);
142 for (k=2; k<=*nz-1; k++) {
144 for (j=2; j<=*ny-1; j++) {
146 for (i=2; i<=*nx-1; i++) {
151 kk = ii - jj + *m + 1;
153 VAT2(acB, kk, jj) = VAT3(oC, i, j, k);
157 kk = ii - jj + *m + 1;
158 VAT2(acB, kk, jj) = - VAT3(oE, i-1, j, k);
162 kk = ii - jj + *m + 1;
163 VAT2(acB, kk, jj) = - VAT3(oN, i, j-1, k);
166 ii = jj - (*nx - 2) * (*ny - 2);
167 kk = ii - jj + *m + 1;
168 VAT2(acB, kk, jj) = - VAT3(uC, i, j, k-1);
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) {
190 MAT2(acB, *lda, *ny-1);
192 MAT3( oC, *nx, *ny, *nz);
193 MAT3( oE, *nx, *ny, *nz);
194 MAT3( oN, *nx, *ny, *nz);
195 MAT3( uC, *nx, *ny, *nz);
197 MAT3(oNE, *nx, *ny, *nz);
198 MAT3(oNW, *nx, *ny, *nz);
200 MAT3( uE, *nx, *ny, *nz);
201 MAT3( uW, *nx, *ny, *nz);
202 MAT3( uN, *nx, *ny, *nz);
203 MAT3( uS, *nx, *ny, *nz);
205 MAT3(uNE, *nx, *ny, *nz);
206 MAT3(uNW, *nx, *ny, *nz);
207 MAT3(uSE, *nx, *ny, *nz);
208 MAT3(uSW, *nx, *ny, *nz);
220 for (k=2; k<=*nz-1; k++) {
222 for (j=2; j<=*ny-1; j++) {
224 for (i=2; i<=*nx-1; i++) {
229 kk = ii - jj + *m + 1;
230 VAT2(acB, kk, jj) = VAT3(oC, i, j, k);
234 kk = ii - jj + *m + 1;
235 VAT2(acB, kk, jj) = - VAT3(oE, i-1, j, k);
239 kk = ii - jj + *m + 1;
240 VAT2(acB, kk, jj) = - VAT3(oN, i, j-1, k);
243 ii = jj - (*nx - 2) + 1;
244 kk = ii - jj + *m + 1;
245 VAT2(acB, kk, jj) = - VAT3(oNE, i, j-1, k);
248 ii = jj - (*nx - 2) - 1;
249 kk = ii - jj + *m + 1;
250 VAT2(acB, kk, jj) = - VAT3(oNW, i, j-1, k);
253 ii = jj - (*nx - 2) * (*ny - 2);
254 kk = ii - jj + *m + 1;
255 VAT2(acB, kk, jj) = - VAT3(uC, i, j, k-1);
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);
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);
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);
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);
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);
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);
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);
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);
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.
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.
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.