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: PetscCall(PetscViewerASCIIPrintf(viewer, "subset %" PetscInt_FMT ":", j));
154: for (PetscInt l = 0; l < k; l++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, subset[l]));
155: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " |"));
156: for (PetscInt l = k; l < N; l++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, subset[l]));
157: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even"));
158: }
159: PetscCall(PetscDTSubsetIndex(N, k, subset, &jCheck));
160: PetscCheck(jCheck == j, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "jCheck (%" PetscInt_FMT ") != j (%" PetscInt_FMT ")", jCheck, j);
161: }
162: PetscCall(PetscViewerASCIIPopTab(viewer));
163: PetscCall(PetscFree(subset));
165: /* Make a random k form */
166: PetscCall(PetscMalloc1(Nk, &w));
167: for (j = 0; j < Nk; j++) PetscCall(PetscRandomGetValueReal(rand, &w[j]));
168: /* Make a set of random vectors */
169: PetscCall(PetscMalloc1(N * k, &v));
170: for (j = 0; j < N * k; j++) PetscCall(PetscRandomGetValueReal(rand, &v[j]));
172: PetscCall(PetscDTAltVApply(N, k, w, v, &wv));
174: if (verbose) {
175: PetscCall(PetscViewerASCIIPrintf(viewer, "w:\n"));
176: PetscCall(PetscViewerASCIIPushTab(viewer));
177: if (Nk) PetscCall(PetscRealView(Nk, w, viewer));
178: PetscCall(PetscViewerASCIIPopTab(viewer));
179: PetscCall(PetscViewerASCIIPrintf(viewer, "v:\n"));
180: PetscCall(PetscViewerASCIIPushTab(viewer));
181: if (N * k > 0) PetscCall(PetscRealView(N * k, v, viewer));
182: PetscCall(PetscViewerASCIIPopTab(viewer));
183: PetscCall(PetscViewerASCIIPrintf(viewer, "w(v): %g\n", (double)wv));
184: }
186: /* sanity checks */
187: if (k == 1) { /* 1-forms are functionals (dot products) */
188: PetscInt l;
189: PetscReal wvcheck = 0.;
190: PetscReal diff;
192: for (l = 0; l < N; l++) wvcheck += w[l] * v[l];
193: diff = PetscSqrtReal(PetscSqr(wvcheck - wv));
194: 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);
195: }
196: if (k == N && N < 5) { /* n-forms are scaled determinants */
197: PetscReal det, wvcheck, diff;
199: switch (k) {
200: case 0:
201: det = 1.;
202: break;
203: case 1:
204: det = v[0];
205: break;
206: case 2:
207: det = v[0] * v[3] - v[1] * v[2];
208: break;
209: case 3:
210: 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]);
211: break;
212: case 4:
213: 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]));
214: break;
215: default:
216: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "invalid k");
217: }
218: wvcheck = det * w[0];
219: diff = PetscSqrtReal(PetscSqr(wvcheck - wv));
220: 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);
221: }
222: if (k > 0) { /* k-forms are linear in each component */
223: PetscReal alpha;
224: PetscReal *x, *axv, wx, waxv, waxvcheck;
225: PetscReal diff;
226: PetscReal rj;
227: PetscInt l;
229: PetscCall(PetscMalloc2(N * k, &x, N * k, &axv));
230: PetscCall(PetscRandomGetValueReal(rand, &alpha));
231: PetscCall(PetscRandomSetInterval(rand, 0, k));
232: PetscCall(PetscRandomGetValueReal(rand, &rj));
233: j = (PetscInt)rj;
234: PetscCall(PetscRandomSetInterval(rand, -1., 1.));
235: for (l = 0; l < N * k; l++) x[l] = v[l];
236: for (l = 0; l < N * k; l++) axv[l] = v[l];
237: for (l = 0; l < N; l++) {
238: PetscReal val;
240: PetscCall(PetscRandomGetValueReal(rand, &val));
241: x[j * N + l] = val;
242: axv[j * N + l] += alpha * val;
243: }
244: PetscCall(PetscDTAltVApply(N, k, w, x, &wx));
245: PetscCall(PetscDTAltVApply(N, k, w, axv, &waxv));
246: waxvcheck = alpha * wx + wv;
247: diff = waxv - waxvcheck;
248: 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);
249: PetscCall(PetscFree2(x, axv));
250: }
251: if (k > 1) { /* k-forms are antisymmetric */
252: PetscReal rj, rl, *swapv, wswapv, diff;
253: PetscInt l, m;
255: PetscCall(PetscRandomSetInterval(rand, 0, k));
256: PetscCall(PetscRandomGetValueReal(rand, &rj));
257: j = (PetscInt)rj;
258: l = j;
259: while (l == j) {
260: PetscCall(PetscRandomGetValueReal(rand, &rl));
261: l = (PetscInt)rl;
262: }
263: PetscCall(PetscRandomSetInterval(rand, -1., 1.));
264: PetscCall(PetscMalloc1(N * k, &swapv));
265: for (m = 0; m < N * k; m++) swapv[m] = v[m];
266: for (m = 0; m < N; m++) {
267: swapv[j * N + m] = v[l * N + m];
268: swapv[l * N + m] = v[j * N + m];
269: }
270: PetscCall(PetscDTAltVApply(N, k, w, swapv, &wswapv));
271: diff = PetscAbsReal(wswapv + wv);
272: 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);
273: PetscCall(PetscFree(swapv));
274: }
275: for (j = 0; j <= k && j + k <= N; j++) { /* wedge product */
276: PetscInt Nj, Njk, l, JKj;
277: PetscReal *u, *uWw, *uWwcheck, *uWwmat, *x, *xsplit, uWwx, uWwxcheck, diff, norm;
278: PetscInt *split;
280: if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "wedge j = %" PetscInt_FMT ":\n", j));
281: PetscCall(PetscViewerASCIIPushTab(viewer));
282: PetscCall(PetscDTBinomialInt(N, j, &Nj));
283: PetscCall(PetscDTBinomialInt(N, j + k, &Njk));
284: PetscCall(PetscMalloc4(Nj, &u, Njk, &uWw, N * (j + k), &x, N * (j + k), &xsplit));
285: PetscCall(PetscMalloc1(j + k, &split));
286: for (l = 0; l < Nj; l++) PetscCall(PetscRandomGetValueReal(rand, &u[l]));
287: for (l = 0; l < N * (j + k); l++) PetscCall(PetscRandomGetValueReal(rand, &x[l]));
288: PetscCall(PetscDTAltVWedge(N, j, k, u, w, uWw));
289: PetscCall(PetscDTAltVApply(N, j + k, uWw, x, &uWwx));
290: if (verbose) {
291: PetscCall(PetscViewerASCIIPrintf(viewer, "u:\n"));
292: PetscCall(PetscViewerASCIIPushTab(viewer));
293: if (Nj) PetscCall(PetscRealView(Nj, u, viewer));
294: PetscCall(PetscViewerASCIIPopTab(viewer));
295: PetscCall(PetscViewerASCIIPrintf(viewer, "u wedge w:\n"));
296: PetscCall(PetscViewerASCIIPushTab(viewer));
297: if (Njk) PetscCall(PetscRealView(Njk, uWw, viewer));
298: PetscCall(PetscViewerASCIIPopTab(viewer));
299: PetscCall(PetscViewerASCIIPrintf(viewer, "x:\n"));
300: PetscCall(PetscViewerASCIIPushTab(viewer));
301: if (N * (j + k) > 0) PetscCall(PetscRealView(N * (j + k), x, viewer));
302: PetscCall(PetscViewerASCIIPopTab(viewer));
303: PetscCall(PetscViewerASCIIPrintf(viewer, "u wedge w(x): %g\n", (double)uWwx));
304: }
305: /* verify wedge formula */
306: uWwxcheck = 0.;
307: PetscCall(PetscDTBinomialInt(j + k, j, &JKj));
308: for (l = 0; l < JKj; l++) {
309: PetscBool isOdd;
310: PetscReal ux, wx;
311: PetscInt m;
313: PetscCall(PetscDTEnumSplit(j + k, j, l, split, &isOdd));
314: for (m = 0; m < j + k; m++) {
315: for (PetscInt p = 0; p < N; p++) xsplit[m * N + p] = x[split[m] * N + p];
316: }
317: PetscCall(PetscDTAltVApply(N, j, u, xsplit, &ux));
318: PetscCall(PetscDTAltVApply(N, k, w, PetscSafePointerPlusOffset(xsplit, j * N), &wx));
319: uWwxcheck += isOdd ? -(ux * wx) : (ux * wx);
320: }
321: diff = PetscAbsReal(uWwx - uWwxcheck);
322: 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);
323: PetscCall(PetscFree(split));
324: PetscCall(PetscMalloc2(Nk * Njk, &uWwmat, Njk, &uWwcheck));
325: PetscCall(PetscDTAltVWedgeMatrix(N, j, k, u, uWwmat));
326: if (verbose) {
327: PetscCall(PetscViewerASCIIPrintf(viewer, "(u wedge):\n"));
328: PetscCall(PetscViewerASCIIPushTab(viewer));
329: if ((Nk * Njk) > 0) PetscCall(PetscRealView(Nk * Njk, uWwmat, viewer));
330: PetscCall(PetscViewerASCIIPopTab(viewer));
331: }
332: diff = 0.;
333: norm = 0.;
334: for (l = 0; l < Njk; l++) {
335: PetscReal sum = 0.;
337: for (PetscInt m = 0; m < Nk; m++) sum += uWwmat[l * Nk + m] * w[m];
338: uWwcheck[l] = sum;
339: diff += PetscSqr(uWwcheck[l] - uWw[l]);
340: norm += PetscSqr(uWwcheck[l]) + PetscSqr(uWw[l]);
341: }
342: diff = PetscSqrtReal(diff);
343: norm = PetscSqrtReal(norm);
344: PetscCheck(diff <= PETSC_SMALL * norm, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "wedge matrix check: wedge matrix application does not match wedge direct application");
345: PetscCall(PetscFree2(uWwmat, uWwcheck));
346: PetscCall(PetscFree4(u, uWw, x, xsplit));
347: PetscCall(PetscViewerASCIIPopTab(viewer));
348: }
349: for (M = PetscMax(1, k); M <= N; M++) { /* pullback */
350: PetscReal *L, *u, *x;
351: PetscInt Mk, l;
353: PetscCall(PetscDTBinomialInt(M, k, &Mk));
354: PetscCall(PetscMalloc3(M * N, &L, Mk, &u, M * k, &x));
355: for (l = 0; l < M * N; l++) PetscCall(PetscRandomGetValueReal(rand, &L[l]));
356: for (l = 0; l < Mk; l++) PetscCall(PetscRandomGetValueReal(rand, &u[l]));
357: for (l = 0; l < M * k; l++) PetscCall(PetscRandomGetValueReal(rand, &x[l]));
358: if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "pullback M = %" PetscInt_FMT ":\n", M));
359: PetscCall(PetscViewerASCIIPushTab(viewer));
360: PetscCall(CheckPullback(M, N, L, k, w, x, verbose, viewer));
361: if (M != N) PetscCall(CheckPullback(N, M, L, k, u, v, PETSC_FALSE, viewer));
362: PetscCall(PetscViewerASCIIPopTab(viewer));
363: if ((k % N) && (N > 1)) {
364: if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "negative pullback M = %" PetscInt_FMT ":\n", M));
365: PetscCall(PetscViewerASCIIPushTab(viewer));
366: PetscCall(CheckPullback(M, N, L, -k, w, x, verbose, viewer));
367: if (M != N) PetscCall(CheckPullback(N, M, L, -k, u, v, PETSC_FALSE, viewer));
368: PetscCall(PetscViewerASCIIPopTab(viewer));
369: }
370: PetscCall(PetscFree3(L, u, x));
371: }
372: if (k > 0) { /* Interior */
373: PetscInt Nkm, l, m;
374: PetscReal *wIntv0, *wIntv0check, wvcheck, diff, diffMat, normMat;
375: PetscReal *intv0mat, *matcheck;
376: PetscInt (*indices)[3];
378: PetscCall(PetscDTBinomialInt(N, k - 1, &Nkm));
379: PetscCall(PetscMalloc5(Nkm, &wIntv0, Nkm, &wIntv0check, Nk * Nkm, &intv0mat, Nk * Nkm, &matcheck, Nk * k, &indices));
380: PetscCall(PetscDTAltVInterior(N, k, w, v, wIntv0));
381: PetscCall(PetscDTAltVInteriorMatrix(N, k, v, intv0mat));
382: PetscCall(PetscDTAltVInteriorPattern(N, k, indices));
383: if (verbose) {
384: PetscCall(PetscViewerASCIIPrintf(viewer, "interior product matrix pattern:\n"));
385: PetscCall(PetscViewerASCIIPushTab(viewer));
386: for (l = 0; l < Nk * k; l++) {
387: PetscInt row = indices[l][0];
388: PetscInt col = indices[l][1];
389: PetscInt x = indices[l][2];
391: PetscCall(PetscViewerASCIIPrintf(viewer, "intV[%" PetscInt_FMT ",%" PetscInt_FMT "] = %sV[%" PetscInt_FMT "]\n", row, col, x < 0 ? "-" : " ", x < 0 ? -(x + 1) : x));
392: }
393: PetscCall(PetscViewerASCIIPopTab(viewer));
394: }
395: for (l = 0; l < Nkm * Nk; l++) matcheck[l] = 0.;
396: for (l = 0; l < Nk * k; l++) {
397: PetscInt row = indices[l][0];
398: PetscInt col = indices[l][1];
399: PetscInt x = indices[l][2];
401: if (x < 0) {
402: matcheck[row * Nk + col] = -v[-(x + 1)];
403: } else {
404: matcheck[row * Nk + col] = v[x];
405: }
406: }
407: diffMat = 0.;
408: normMat = 0.;
409: for (l = 0; l < Nkm * Nk; l++) {
410: diffMat += PetscSqr(PetscAbsReal(matcheck[l] - intv0mat[l]));
411: normMat += PetscSqr(matcheck[l]) + PetscSqr(intv0mat[l]);
412: }
413: diffMat = PetscSqrtReal(diffMat);
414: normMat = PetscSqrtReal(normMat);
415: PetscCheck(diffMat <= PETSC_SMALL * normMat, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior product check: matrix pattern does not match matrix");
416: diffMat = 0.;
417: normMat = 0.;
418: for (l = 0; l < Nkm; l++) {
419: PetscReal sum = 0.;
421: for (m = 0; m < Nk; m++) sum += intv0mat[l * Nk + m] * w[m];
422: wIntv0check[l] = sum;
424: diffMat += PetscSqr(PetscAbsReal(wIntv0check[l] - wIntv0[l]));
425: normMat += PetscSqr(wIntv0check[l]) + PetscSqr(wIntv0[l]);
426: }
427: diffMat = PetscSqrtReal(diffMat);
428: normMat = PetscSqrtReal(normMat);
429: PetscCheck(diffMat <= PETSC_SMALL * normMat, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior product check: application does not match matrix");
430: if (verbose) {
431: PetscCall(PetscViewerASCIIPrintf(viewer, "(w int v_0):\n"));
432: PetscCall(PetscViewerASCIIPushTab(viewer));
433: if (Nkm) PetscCall(PetscRealView(Nkm, wIntv0, viewer));
434: PetscCall(PetscViewerASCIIPopTab(viewer));
436: PetscCall(PetscViewerASCIIPrintf(viewer, "(int v_0):\n"));
437: PetscCall(PetscViewerASCIIPushTab(viewer));
438: if (Nk * Nkm > 0) PetscCall(PetscRealView(Nk * Nkm, intv0mat, viewer));
439: PetscCall(PetscViewerASCIIPopTab(viewer));
440: }
441: PetscCall(PetscDTAltVApply(N, k - 1, wIntv0, &v[N], &wvcheck));
442: diff = PetscSqrtReal(PetscSqr(wvcheck - wv));
443: 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);
444: PetscCall(PetscFree5(wIntv0, wIntv0check, intv0mat, matcheck, indices));
445: }
446: if (k >= N - k) { /* Hodge star */
447: PetscReal *u, *starw, *starstarw, wu, starwdotu;
448: PetscReal diff, norm;
449: PetscBool isOdd;
451: isOdd = (PetscBool)((k * (N - k)) & 1);
452: PetscCall(PetscMalloc3(Nk, &u, Nk, &starw, Nk, &starstarw));
453: PetscCall(PetscDTAltVStar(N, k, 1, w, starw));
454: PetscCall(PetscDTAltVStar(N, N - k, 1, starw, starstarw));
455: if (verbose) {
456: PetscCall(PetscViewerASCIIPrintf(viewer, "star w:\n"));
457: PetscCall(PetscViewerASCIIPushTab(viewer));
458: if (Nk) PetscCall(PetscRealView(Nk, starw, viewer));
459: PetscCall(PetscViewerASCIIPopTab(viewer));
461: PetscCall(PetscViewerASCIIPrintf(viewer, "star star w:\n"));
462: PetscCall(PetscViewerASCIIPushTab(viewer));
463: if (Nk) PetscCall(PetscRealView(Nk, starstarw, viewer));
464: PetscCall(PetscViewerASCIIPopTab(viewer));
465: }
466: for (PetscInt l = 0; l < Nk; l++) PetscCall(PetscRandomGetValueReal(rand, &u[l]));
467: PetscCall(PetscDTAltVWedge(N, k, N - k, w, u, &wu));
468: starwdotu = 0.;
469: for (PetscInt l = 0; l < Nk; l++) starwdotu += starw[l] * u[l];
470: diff = PetscAbsReal(wu - starwdotu);
471: 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);
473: diff = 0.;
474: norm = 0.;
475: for (PetscInt l = 0; l < Nk; l++) {
476: diff += PetscSqr(w[l] - (isOdd ? -starstarw[l] : starstarw[l]));
477: norm += PetscSqr(w[l]) + PetscSqr(starstarw[l]);
478: }
479: diff = PetscSqrtReal(diff);
480: norm = PetscSqrtReal(norm);
481: PetscCheck(diff <= PETSC_SMALL * norm, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Hodge star check: star(star(w)) != (-1)^(N*(N-k)) w");
482: PetscCall(PetscFree3(u, starw, starstarw));
483: }
484: PetscCall(PetscFree(v));
485: PetscCall(PetscFree(w));
486: PetscCall(PetscViewerASCIIPopTab(viewer));
487: }
488: PetscCall(PetscViewerASCIIPopTab(viewer));
489: }
490: PetscCall(PetscRandomDestroy(&rand));
491: PetscCall(PetscFinalize());
492: return 0;
493: }
495: /*TEST
496: test:
497: suffix: 1234
498: args: -verbose
499: test:
500: suffix: 56
501: args: -N 5,6
502: output_file: output/empty.out
503: TEST*/