Actual source code: ex7.c

  1: static char help[] = "Test the PetscDTAltV interface for k-forms (alternating k-linear maps).\n\n";

  3: #include <petscviewer.h>
  4: #include <petscdt.h>

  6: static PetscErrorCode CheckPullback(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, const PetscReal *w, PetscReal *x, PetscBool verbose, PetscViewer viewer)
  7: {
  8:   PetscInt         Nk, Mk, i, j, l;
  9:   PetscReal       *Lstarw, *Lx, *Lstar, *Lstarwcheck, wLx, Lstarwx;
 10:   PetscReal        diff, diffMat, normMat;
 11:   PetscReal       *walloc   = NULL;
 12:   const PetscReal *ww       = NULL;
 13:   PetscBool        negative = (PetscBool)(k < 0);

 15:   PetscFunctionBegin;
 16:   k = PetscAbsInt(k);
 17:   PetscCall(PetscDTBinomialInt(N, k, &Nk));
 18:   PetscCall(PetscDTBinomialInt(M, k, &Mk));
 19:   if (negative) {
 20:     PetscCall(PetscMalloc1(Mk, &walloc));
 21:     PetscCall(PetscDTAltVStar(M, M - k, 1, w, walloc));
 22:     ww = walloc;
 23:   } else {
 24:     ww = w;
 25:   }
 26:   PetscCall(PetscMalloc2(Nk, &Lstarw, (M * k), &Lx));
 27:   PetscCall(PetscMalloc2(Nk * Mk, &Lstar, Nk, &Lstarwcheck));
 28:   PetscCall(PetscDTAltVPullback(N, M, L, negative ? -k : k, w, Lstarw));
 29:   PetscCall(PetscDTAltVPullbackMatrix(N, M, L, negative ? -k : k, Lstar));
 30:   if (negative) {
 31:     PetscReal *sLsw;

 33:     PetscCall(PetscMalloc1(Nk, &sLsw));
 34:     PetscCall(PetscDTAltVStar(N, N - k, 1, Lstarw, sLsw));
 35:     PetscCall(PetscDTAltVApply(N, k, sLsw, x, &Lstarwx));
 36:     PetscCall(PetscFree(sLsw));
 37:   } else {
 38:     PetscCall(PetscDTAltVApply(N, k, Lstarw, x, &Lstarwx));
 39:   }
 40:   for (l = 0; l < k; l++) {
 41:     for (i = 0; i < M; i++) {
 42:       PetscReal sum = 0.;

 44:       for (j = 0; j < N; j++) sum += L[i * N + j] * x[l * N + j];
 45:       Lx[l * M + i] = sum;
 46:     }
 47:   }
 48:   diffMat = 0.;
 49:   normMat = 0.;
 50:   for (i = 0; i < Nk; i++) {
 51:     PetscReal sum = 0.;
 52:     for (j = 0; j < Mk; j++) sum += Lstar[i * Mk + j] * w[j];
 53:     Lstarwcheck[i] = sum;
 54:     diffMat += PetscSqr(PetscAbsReal(Lstarwcheck[i] - Lstarw[i]));
 55:     normMat += PetscSqr(Lstarwcheck[i]) + PetscSqr(Lstarw[i]);
 56:   }
 57:   diffMat = PetscSqrtReal(diffMat);
 58:   normMat = PetscSqrtReal(normMat);
 59:   if (verbose) {
 60:     PetscCall(PetscViewerASCIIPrintf(viewer, "L:\n"));
 61:     PetscCall(PetscViewerASCIIPushTab(viewer));
 62:     if (M * N > 0) PetscCall(PetscRealView(M * N, L, viewer));
 63:     PetscCall(PetscViewerASCIIPopTab(viewer));

 65:     PetscCall(PetscViewerASCIIPrintf(viewer, "L*:\n"));
 66:     PetscCall(PetscViewerASCIIPushTab(viewer));
 67:     if (Nk * Mk > 0) PetscCall(PetscRealView(Nk * Mk, Lstar, viewer));
 68:     PetscCall(PetscViewerASCIIPopTab(viewer));

 70:     PetscCall(PetscViewerASCIIPrintf(viewer, "L*w:\n"));
 71:     PetscCall(PetscViewerASCIIPushTab(viewer));
 72:     if (Nk > 0) PetscCall(PetscRealView(Nk, Lstarw, viewer));
 73:     PetscCall(PetscViewerASCIIPopTab(viewer));
 74:   }
 75:   PetscCall(PetscDTAltVApply(M, k, ww, Lx, &wLx));
 76:   diff = PetscAbsReal(wLx - Lstarwx);
 77:   PetscCheck(diff <= 10. * PETSC_SMALL * (PetscAbsReal(wLx) + PetscAbsReal(Lstarwx)), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "pullback check: pullback does not commute with application: w(Lx)(%g) != (L* w)(x)(%g)", (double)wLx, (double)Lstarwx);
 78:   PetscCheck(diffMat <= PETSC_SMALL * normMat, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "pullback check: pullback matrix does match matrix-free result");
 79:   PetscCall(PetscFree2(Lstar, Lstarwcheck));
 80:   PetscCall(PetscFree2(Lstarw, Lx));
 81:   PetscCall(PetscFree(walloc));
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

 85: int main(int argc, char **argv)
 86: {
 87:   PetscInt    i, numTests = 5, n[5] = {0, 1, 2, 3, 4};
 88:   PetscBool   verbose = PETSC_FALSE;
 89:   PetscRandom rand;
 90:   PetscViewer viewer;

 92:   PetscFunctionBeginUser;
 93:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 94:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for exterior algebra tests", "none");
 95:   PetscCall(PetscOptionsIntArray("-N", "Up to 5 vector space dimensions to test", "ex7.c", n, &numTests, NULL));
 96:   PetscCall(PetscOptionsBool("-verbose", "Verbose test output", "ex7.c", verbose, &verbose, NULL));
 97:   PetscOptionsEnd();
 98:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand));
 99:   PetscCall(PetscRandomSetInterval(rand, -1., 1.));
100:   PetscCall(PetscRandomSetFromOptions(rand));
101:   if (!numTests) numTests = 5;
102:   viewer = PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD);
103:   for (i = 0; i < numTests; i++) {
104:     PetscInt k, N = n[i];

106:     if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "N = %" PetscInt_FMT ":\n", N));
107:     PetscCall(PetscViewerASCIIPushTab(viewer));

109:     if (verbose) {
110:       PetscInt *perm;
111:       PetscInt  fac = 1;

113:       PetscCall(PetscMalloc1(N, &perm));

115:       for (k = 1; k <= N; k++) fac *= k;
116:       PetscCall(PetscViewerASCIIPrintf(viewer, "Permutations of %" PetscInt_FMT ":\n", N));
117:       PetscCall(PetscViewerASCIIPushTab(viewer));
118:       for (k = 0; k < fac; k++) {
119:         PetscBool isOdd, isOddCheck;
120:         PetscInt  j, kCheck;

122:         PetscCall(PetscDTEnumPerm(N, k, perm, &isOdd));
123:         PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ":", k));
124:         for (j = 0; j < N; j++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, perm[j]));
125:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even"));
126:         PetscCall(PetscDTPermIndex(N, perm, &kCheck, &isOddCheck));
127:         PetscCheck(kCheck == k && isOddCheck == isOdd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PetscDTEnumPerm / PetscDTPermIndex mismatch for (%" PetscInt_FMT ", %" PetscInt_FMT ")", N, k);
128:       }
129:       PetscCall(PetscViewerASCIIPopTab(viewer));
130:       PetscCall(PetscFree(perm));
131:     }
132:     for (k = 0; k <= N; k++) {
133:       PetscInt   j, Nk, M;
134:       PetscReal *w, *v, wv;
135:       PetscInt  *subset;

137:       PetscCall(PetscDTBinomialInt(N, k, &Nk));
138:       if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "k = %" PetscInt_FMT ":\n", k));
139:       PetscCall(PetscViewerASCIIPushTab(viewer));
140:       if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT " choose %" PetscInt_FMT "): %" PetscInt_FMT "\n", N, k, Nk));

142:       /* Test subset and complement enumeration */
143:       PetscCall(PetscMalloc1(N, &subset));
144:       PetscCall(PetscViewerASCIIPushTab(viewer));
145:       for (j = 0; j < Nk; j++) {
146:         PetscBool isOdd, isOddCheck;
147:         PetscInt  jCheck, kCheck;

149:         PetscCall(PetscDTEnumSplit(N, k, j, subset, &isOdd));
150:         PetscCall(PetscDTPermIndex(N, subset, &kCheck, &isOddCheck));
151:         PetscCheck(isOddCheck == isOdd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PetscDTEnumSplit sign does not mmatch PetscDTPermIndex sign");
152:         if (verbose) {
153:           PetscInt l;

155:           PetscCall(PetscViewerASCIIPrintf(viewer, "subset %" PetscInt_FMT ":", j));
156:           for (l = 0; l < k; l++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, subset[l]));
157:           PetscCall(PetscPrintf(PETSC_COMM_WORLD, " |"));
158:           for (l = k; l < N; l++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, subset[l]));
159:           PetscCall(PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even"));
160:         }
161:         PetscCall(PetscDTSubsetIndex(N, k, subset, &jCheck));
162:         PetscCheck(jCheck == j, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "jCheck (%" PetscInt_FMT ") != j (%" PetscInt_FMT ")", jCheck, j);
163:       }
164:       PetscCall(PetscViewerASCIIPopTab(viewer));
165:       PetscCall(PetscFree(subset));

167:       /* Make a random k form */
168:       PetscCall(PetscMalloc1(Nk, &w));
169:       for (j = 0; j < Nk; j++) PetscCall(PetscRandomGetValueReal(rand, &w[j]));
170:       /* Make a set of random vectors */
171:       PetscCall(PetscMalloc1(N * k, &v));
172:       for (j = 0; j < N * k; j++) PetscCall(PetscRandomGetValueReal(rand, &v[j]));

174:       PetscCall(PetscDTAltVApply(N, k, w, v, &wv));

176:       if (verbose) {
177:         PetscCall(PetscViewerASCIIPrintf(viewer, "w:\n"));
178:         PetscCall(PetscViewerASCIIPushTab(viewer));
179:         if (Nk) PetscCall(PetscRealView(Nk, w, viewer));
180:         PetscCall(PetscViewerASCIIPopTab(viewer));
181:         PetscCall(PetscViewerASCIIPrintf(viewer, "v:\n"));
182:         PetscCall(PetscViewerASCIIPushTab(viewer));
183:         if (N * k > 0) PetscCall(PetscRealView(N * k, v, viewer));
184:         PetscCall(PetscViewerASCIIPopTab(viewer));
185:         PetscCall(PetscViewerASCIIPrintf(viewer, "w(v): %g\n", (double)wv));
186:       }

188:       /* sanity checks */
189:       if (k == 1) { /* 1-forms are functionals (dot products) */
190:         PetscInt  l;
191:         PetscReal wvcheck = 0.;
192:         PetscReal diff;

194:         for (l = 0; l < N; l++) wvcheck += w[l] * v[l];
195:         diff = PetscSqrtReal(PetscSqr(wvcheck - wv));
196:         PetscCheck(diff < PETSC_SMALL * (PetscAbsReal(wv) + PetscAbsReal(wvcheck)), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "1-form / dot product equivalence: wvcheck (%g) != wv (%g)", (double)wvcheck, (double)wv);
197:       }
198:       if (k == N && N < 5) { /* n-forms are scaled determinants */
199:         PetscReal det, wvcheck, diff;

201:         switch (k) {
202:         case 0:
203:           det = 1.;
204:           break;
205:         case 1:
206:           det = v[0];
207:           break;
208:         case 2:
209:           det = v[0] * v[3] - v[1] * v[2];
210:           break;
211:         case 3:
212:           det = v[0] * (v[4] * v[8] - v[5] * v[7]) + v[1] * (v[5] * v[6] - v[3] * v[8]) + v[2] * (v[3] * v[7] - v[4] * v[6]);
213:           break;
214:         case 4:
215:           det = v[0] * (v[5] * (v[10] * v[15] - v[11] * v[14]) + v[6] * (v[11] * v[13] - v[9] * v[15]) + v[7] * (v[9] * v[14] - v[10] * v[13])) - v[1] * (v[4] * (v[10] * v[15] - v[11] * v[14]) + v[6] * (v[11] * v[12] - v[8] * v[15]) + v[7] * (v[8] * v[14] - v[10] * v[12])) + v[2] * (v[4] * (v[9] * v[15] - v[11] * v[13]) + v[5] * (v[11] * v[12] - v[8] * v[15]) + v[7] * (v[8] * v[13] - v[9] * v[12])) - v[3] * (v[4] * (v[9] * v[14] - v[10] * v[13]) + v[5] * (v[10] * v[12] - v[8] * v[14]) + v[6] * (v[8] * v[13] - v[9] * v[12]));
216:           break;
217:         default:
218:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "invalid k");
219:         }
220:         wvcheck = det * w[0];
221:         diff    = PetscSqrtReal(PetscSqr(wvcheck - wv));
222:         PetscCheck(diff < PETSC_SMALL * (PetscAbsReal(wv) + PetscAbsReal(wvcheck)), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "n-form / determinant equivalence: wvcheck (%g) != wv (%g) %g", (double)wvcheck, (double)wv, (double)diff);
223:       }
224:       if (k > 0) { /* k-forms are linear in each component */
225:         PetscReal  alpha;
226:         PetscReal *x, *axv, wx, waxv, waxvcheck;
227:         PetscReal  diff;
228:         PetscReal  rj;
229:         PetscInt   l;

231:         PetscCall(PetscMalloc2(N * k, &x, N * k, &axv));
232:         PetscCall(PetscRandomGetValueReal(rand, &alpha));
233:         PetscCall(PetscRandomSetInterval(rand, 0, k));
234:         PetscCall(PetscRandomGetValueReal(rand, &rj));
235:         j = (PetscInt)rj;
236:         PetscCall(PetscRandomSetInterval(rand, -1., 1.));
237:         for (l = 0; l < N * k; l++) x[l] = v[l];
238:         for (l = 0; l < N * k; l++) axv[l] = v[l];
239:         for (l = 0; l < N; l++) {
240:           PetscReal val;

242:           PetscCall(PetscRandomGetValueReal(rand, &val));
243:           x[j * N + l] = val;
244:           axv[j * N + l] += alpha * val;
245:         }
246:         PetscCall(PetscDTAltVApply(N, k, w, x, &wx));
247:         PetscCall(PetscDTAltVApply(N, k, w, axv, &waxv));
248:         waxvcheck = alpha * wx + wv;
249:         diff      = waxv - waxvcheck;
250:         PetscCheck(PetscAbsReal(diff) <= 10. * PETSC_SMALL * (PetscAbsReal(waxv) + PetscAbsReal(waxvcheck)), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "linearity check: component %" PetscInt_FMT ", waxvcheck (%g) != waxv (%g)", j, (double)waxvcheck, (double)waxv);
251:         PetscCall(PetscFree2(x, axv));
252:       }
253:       if (k > 1) { /* k-forms are antisymmetric */
254:         PetscReal rj, rl, *swapv, wswapv, diff;
255:         PetscInt  l, m;

257:         PetscCall(PetscRandomSetInterval(rand, 0, k));
258:         PetscCall(PetscRandomGetValueReal(rand, &rj));
259:         j = (PetscInt)rj;
260:         l = j;
261:         while (l == j) {
262:           PetscCall(PetscRandomGetValueReal(rand, &rl));
263:           l = (PetscInt)rl;
264:         }
265:         PetscCall(PetscRandomSetInterval(rand, -1., 1.));
266:         PetscCall(PetscMalloc1(N * k, &swapv));
267:         for (m = 0; m < N * k; m++) swapv[m] = v[m];
268:         for (m = 0; m < N; m++) {
269:           swapv[j * N + m] = v[l * N + m];
270:           swapv[l * N + m] = v[j * N + m];
271:         }
272:         PetscCall(PetscDTAltVApply(N, k, w, swapv, &wswapv));
273:         diff = PetscAbsReal(wswapv + wv);
274:         PetscCheck(diff <= PETSC_SMALL * (PetscAbsReal(wswapv) + PetscAbsReal(wv)), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "antisymmetry check: components %" PetscInt_FMT " & %" PetscInt_FMT ", wswapv (%g) != -wv (%g)", j, l, (double)wswapv, (double)wv);
275:         PetscCall(PetscFree(swapv));
276:       }
277:       for (j = 0; j <= k && j + k <= N; j++) { /* wedge product */
278:         PetscInt   Nj, Njk, l, JKj;
279:         PetscReal *u, *uWw, *uWwcheck, *uWwmat, *x, *xsplit, uWwx, uWwxcheck, diff, norm;
280:         PetscInt  *split;

282:         if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "wedge j = %" PetscInt_FMT ":\n", j));
283:         PetscCall(PetscViewerASCIIPushTab(viewer));
284:         PetscCall(PetscDTBinomialInt(N, j, &Nj));
285:         PetscCall(PetscDTBinomialInt(N, j + k, &Njk));
286:         PetscCall(PetscMalloc4(Nj, &u, Njk, &uWw, N * (j + k), &x, N * (j + k), &xsplit));
287:         PetscCall(PetscMalloc1(j + k, &split));
288:         for (l = 0; l < Nj; l++) PetscCall(PetscRandomGetValueReal(rand, &u[l]));
289:         for (l = 0; l < N * (j + k); l++) PetscCall(PetscRandomGetValueReal(rand, &x[l]));
290:         PetscCall(PetscDTAltVWedge(N, j, k, u, w, uWw));
291:         PetscCall(PetscDTAltVApply(N, j + k, uWw, x, &uWwx));
292:         if (verbose) {
293:           PetscCall(PetscViewerASCIIPrintf(viewer, "u:\n"));
294:           PetscCall(PetscViewerASCIIPushTab(viewer));
295:           if (Nj) PetscCall(PetscRealView(Nj, u, viewer));
296:           PetscCall(PetscViewerASCIIPopTab(viewer));
297:           PetscCall(PetscViewerASCIIPrintf(viewer, "u wedge w:\n"));
298:           PetscCall(PetscViewerASCIIPushTab(viewer));
299:           if (Njk) PetscCall(PetscRealView(Njk, uWw, viewer));
300:           PetscCall(PetscViewerASCIIPopTab(viewer));
301:           PetscCall(PetscViewerASCIIPrintf(viewer, "x:\n"));
302:           PetscCall(PetscViewerASCIIPushTab(viewer));
303:           if (N * (j + k) > 0) PetscCall(PetscRealView(N * (j + k), x, viewer));
304:           PetscCall(PetscViewerASCIIPopTab(viewer));
305:           PetscCall(PetscViewerASCIIPrintf(viewer, "u wedge w(x): %g\n", (double)uWwx));
306:         }
307:         /* verify wedge formula */
308:         uWwxcheck = 0.;
309:         PetscCall(PetscDTBinomialInt(j + k, j, &JKj));
310:         for (l = 0; l < JKj; l++) {
311:           PetscBool isOdd;
312:           PetscReal ux, wx;
313:           PetscInt  m, p;

315:           PetscCall(PetscDTEnumSplit(j + k, j, l, split, &isOdd));
316:           for (m = 0; m < j + k; m++) {
317:             for (p = 0; p < N; p++) xsplit[m * N + p] = x[split[m] * N + p];
318:           }
319:           PetscCall(PetscDTAltVApply(N, j, u, xsplit, &ux));
320:           PetscCall(PetscDTAltVApply(N, k, w, &xsplit[j * N], &wx));
321:           uWwxcheck += isOdd ? -(ux * wx) : (ux * wx);
322:         }
323:         diff = PetscAbsReal(uWwx - uWwxcheck);
324:         PetscCheck(diff <= 10. * PETSC_SMALL * (PetscAbsReal(uWwx) + PetscAbsReal(uWwxcheck)), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "wedge check: forms %" PetscInt_FMT " & %" PetscInt_FMT ", uWwxcheck (%g) != uWwx (%g)", j, k, (double)uWwxcheck, (double)uWwx);
325:         PetscCall(PetscFree(split));
326:         PetscCall(PetscMalloc2(Nk * Njk, &uWwmat, Njk, &uWwcheck));
327:         PetscCall(PetscDTAltVWedgeMatrix(N, j, k, u, uWwmat));
328:         if (verbose) {
329:           PetscCall(PetscViewerASCIIPrintf(viewer, "(u wedge):\n"));
330:           PetscCall(PetscViewerASCIIPushTab(viewer));
331:           if ((Nk * Njk) > 0) PetscCall(PetscRealView(Nk * Njk, uWwmat, viewer));
332:           PetscCall(PetscViewerASCIIPopTab(viewer));
333:         }
334:         diff = 0.;
335:         norm = 0.;
336:         for (l = 0; l < Njk; l++) {
337:           PetscInt  m;
338:           PetscReal sum = 0.;

340:           for (m = 0; m < Nk; m++) sum += uWwmat[l * Nk + m] * w[m];
341:           uWwcheck[l] = sum;
342:           diff += PetscSqr(uWwcheck[l] - uWw[l]);
343:           norm += PetscSqr(uWwcheck[l]) + PetscSqr(uWw[l]);
344:         }
345:         diff = PetscSqrtReal(diff);
346:         norm = PetscSqrtReal(norm);
347:         PetscCheck(diff <= PETSC_SMALL * norm, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "wedge matrix check: wedge matrix application does not match wedge direct application");
348:         PetscCall(PetscFree2(uWwmat, uWwcheck));
349:         PetscCall(PetscFree4(u, uWw, x, xsplit));
350:         PetscCall(PetscViewerASCIIPopTab(viewer));
351:       }
352:       for (M = PetscMax(1, k); M <= N; M++) { /* pullback */
353:         PetscReal *L, *u, *x;
354:         PetscInt   Mk, l;

356:         PetscCall(PetscDTBinomialInt(M, k, &Mk));
357:         PetscCall(PetscMalloc3(M * N, &L, Mk, &u, M * k, &x));
358:         for (l = 0; l < M * N; l++) PetscCall(PetscRandomGetValueReal(rand, &L[l]));
359:         for (l = 0; l < Mk; l++) PetscCall(PetscRandomGetValueReal(rand, &u[l]));
360:         for (l = 0; l < M * k; l++) PetscCall(PetscRandomGetValueReal(rand, &x[l]));
361:         if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "pullback M = %" PetscInt_FMT ":\n", M));
362:         PetscCall(PetscViewerASCIIPushTab(viewer));
363:         PetscCall(CheckPullback(M, N, L, k, w, x, verbose, viewer));
364:         if (M != N) PetscCall(CheckPullback(N, M, L, k, u, v, PETSC_FALSE, viewer));
365:         PetscCall(PetscViewerASCIIPopTab(viewer));
366:         if ((k % N) && (N > 1)) {
367:           if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "negative pullback M = %" PetscInt_FMT ":\n", M));
368:           PetscCall(PetscViewerASCIIPushTab(viewer));
369:           PetscCall(CheckPullback(M, N, L, -k, w, x, verbose, viewer));
370:           if (M != N) PetscCall(CheckPullback(N, M, L, -k, u, v, PETSC_FALSE, viewer));
371:           PetscCall(PetscViewerASCIIPopTab(viewer));
372:         }
373:         PetscCall(PetscFree3(L, u, x));
374:       }
375:       if (k > 0) { /* Interior */
376:         PetscInt   Nkm, l, m;
377:         PetscReal *wIntv0, *wIntv0check, wvcheck, diff, diffMat, normMat;
378:         PetscReal *intv0mat, *matcheck;
379:         PetscInt(*indices)[3];

381:         PetscCall(PetscDTBinomialInt(N, k - 1, &Nkm));
382:         PetscCall(PetscMalloc5(Nkm, &wIntv0, Nkm, &wIntv0check, Nk * Nkm, &intv0mat, Nk * Nkm, &matcheck, Nk * k, &indices));
383:         PetscCall(PetscDTAltVInterior(N, k, w, v, wIntv0));
384:         PetscCall(PetscDTAltVInteriorMatrix(N, k, v, intv0mat));
385:         PetscCall(PetscDTAltVInteriorPattern(N, k, indices));
386:         if (verbose) {
387:           PetscCall(PetscViewerASCIIPrintf(viewer, "interior product matrix pattern:\n"));
388:           PetscCall(PetscViewerASCIIPushTab(viewer));
389:           for (l = 0; l < Nk * k; l++) {
390:             PetscInt row = indices[l][0];
391:             PetscInt col = indices[l][1];
392:             PetscInt x   = indices[l][2];

394:             PetscCall(PetscViewerASCIIPrintf(viewer, "intV[%" PetscInt_FMT ",%" PetscInt_FMT "] = %sV[%" PetscInt_FMT "]\n", row, col, x < 0 ? "-" : " ", x < 0 ? -(x + 1) : x));
395:           }
396:           PetscCall(PetscViewerASCIIPopTab(viewer));
397:         }
398:         for (l = 0; l < Nkm * Nk; l++) matcheck[l] = 0.;
399:         for (l = 0; l < Nk * k; l++) {
400:           PetscInt row = indices[l][0];
401:           PetscInt col = indices[l][1];
402:           PetscInt x   = indices[l][2];

404:           if (x < 0) {
405:             matcheck[row * Nk + col] = -v[-(x + 1)];
406:           } else {
407:             matcheck[row * Nk + col] = v[x];
408:           }
409:         }
410:         diffMat = 0.;
411:         normMat = 0.;
412:         for (l = 0; l < Nkm * Nk; l++) {
413:           diffMat += PetscSqr(PetscAbsReal(matcheck[l] - intv0mat[l]));
414:           normMat += PetscSqr(matcheck[l]) + PetscSqr(intv0mat[l]);
415:         }
416:         diffMat = PetscSqrtReal(diffMat);
417:         normMat = PetscSqrtReal(normMat);
418:         PetscCheck(diffMat <= PETSC_SMALL * normMat, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior product check: matrix pattern does not match matrix");
419:         diffMat = 0.;
420:         normMat = 0.;
421:         for (l = 0; l < Nkm; l++) {
422:           PetscReal sum = 0.;

424:           for (m = 0; m < Nk; m++) sum += intv0mat[l * Nk + m] * w[m];
425:           wIntv0check[l] = sum;

427:           diffMat += PetscSqr(PetscAbsReal(wIntv0check[l] - wIntv0[l]));
428:           normMat += PetscSqr(wIntv0check[l]) + PetscSqr(wIntv0[l]);
429:         }
430:         diffMat = PetscSqrtReal(diffMat);
431:         normMat = PetscSqrtReal(normMat);
432:         PetscCheck(diffMat <= PETSC_SMALL * normMat, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior product check: application does not match matrix");
433:         if (verbose) {
434:           PetscCall(PetscViewerASCIIPrintf(viewer, "(w int v_0):\n"));
435:           PetscCall(PetscViewerASCIIPushTab(viewer));
436:           if (Nkm) PetscCall(PetscRealView(Nkm, wIntv0, viewer));
437:           PetscCall(PetscViewerASCIIPopTab(viewer));

439:           PetscCall(PetscViewerASCIIPrintf(viewer, "(int v_0):\n"));
440:           PetscCall(PetscViewerASCIIPushTab(viewer));
441:           if (Nk * Nkm > 0) PetscCall(PetscRealView(Nk * Nkm, intv0mat, viewer));
442:           PetscCall(PetscViewerASCIIPopTab(viewer));
443:         }
444:         PetscCall(PetscDTAltVApply(N, k - 1, wIntv0, &v[N], &wvcheck));
445:         diff = PetscSqrtReal(PetscSqr(wvcheck - wv));
446:         PetscCheck(diff < PETSC_SMALL * (PetscAbsReal(wv) + PetscAbsReal(wvcheck)), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior product check: (w Int v0)(v_rem) (%g) != w(v) (%g)", (double)wvcheck, (double)wv);
447:         PetscCall(PetscFree5(wIntv0, wIntv0check, intv0mat, matcheck, indices));
448:       }
449:       if (k >= N - k) { /* Hodge star */
450:         PetscReal *u, *starw, *starstarw, wu, starwdotu;
451:         PetscReal  diff, norm;
452:         PetscBool  isOdd;
453:         PetscInt   l;

455:         isOdd = (PetscBool)((k * (N - k)) & 1);
456:         PetscCall(PetscMalloc3(Nk, &u, Nk, &starw, Nk, &starstarw));
457:         PetscCall(PetscDTAltVStar(N, k, 1, w, starw));
458:         PetscCall(PetscDTAltVStar(N, N - k, 1, starw, starstarw));
459:         if (verbose) {
460:           PetscCall(PetscViewerASCIIPrintf(viewer, "star w:\n"));
461:           PetscCall(PetscViewerASCIIPushTab(viewer));
462:           if (Nk) PetscCall(PetscRealView(Nk, starw, viewer));
463:           PetscCall(PetscViewerASCIIPopTab(viewer));

465:           PetscCall(PetscViewerASCIIPrintf(viewer, "star star w:\n"));
466:           PetscCall(PetscViewerASCIIPushTab(viewer));
467:           if (Nk) PetscCall(PetscRealView(Nk, starstarw, viewer));
468:           PetscCall(PetscViewerASCIIPopTab(viewer));
469:         }
470:         for (l = 0; l < Nk; l++) PetscCall(PetscRandomGetValueReal(rand, &u[l]));
471:         PetscCall(PetscDTAltVWedge(N, k, N - k, w, u, &wu));
472:         starwdotu = 0.;
473:         for (l = 0; l < Nk; l++) starwdotu += starw[l] * u[l];
474:         diff = PetscAbsReal(wu - starwdotu);
475:         PetscCheck(diff <= PETSC_SMALL * (PetscAbsReal(wu) + PetscAbsReal(starwdotu)), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Hodge star check: (star w, u) (%g) != (w wedge u) (%g)", (double)starwdotu, (double)wu);

477:         diff = 0.;
478:         norm = 0.;
479:         for (l = 0; l < Nk; l++) {
480:           diff += PetscSqr(w[l] - (isOdd ? -starstarw[l] : starstarw[l]));
481:           norm += PetscSqr(w[l]) + PetscSqr(starstarw[l]);
482:         }
483:         diff = PetscSqrtReal(diff);
484:         norm = PetscSqrtReal(norm);
485:         PetscCheck(diff <= PETSC_SMALL * norm, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Hodge star check: star(star(w)) != (-1)^(N*(N-k)) w");
486:         PetscCall(PetscFree3(u, starw, starstarw));
487:       }
488:       PetscCall(PetscFree(v));
489:       PetscCall(PetscFree(w));
490:       PetscCall(PetscViewerASCIIPopTab(viewer));
491:     }
492:     PetscCall(PetscViewerASCIIPopTab(viewer));
493:   }
494:   PetscCall(PetscRandomDestroy(&rand));
495:   PetscCall(PetscFinalize());
496:   return 0;
497: }

499: /*TEST
500:   test:
501:     suffix: 1234
502:     args: -verbose
503:   test:
504:     suffix: 56
505:     args: -N 5,6
506: TEST*/