Actual source code: ex7.c
petsc-3.13.4 2020-08-01
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,
7: const PetscReal *w, PetscReal *x, PetscBool verbose, PetscViewer viewer)
8: {
9: PetscInt Nk, Mk, i, j, l;
10: PetscReal *Lstarw, *Lx, *Lstar, *Lstarwcheck, wLx, Lstarwx;
11: PetscReal diff, diffMat, normMat;
12: PetscReal *walloc = NULL;
13: const PetscReal *ww = NULL;
14: PetscBool negative = (PetscBool) (k < 0);
15: PetscErrorCode ierr;
18: k = PetscAbsInt(k);
19: PetscDTBinomialInt(N, k, &Nk);
20: PetscDTBinomialInt(M, k, &Mk);
21: if (negative) {
22: PetscMalloc1(Mk, &walloc);
23: PetscDTAltVStar(M, M - k, 1, w, walloc);
24: ww = walloc;
25: } else {
26: ww = w;
27: }
28: PetscMalloc2(Nk, &Lstarw, (M*k), &Lx);
29: PetscMalloc2(Nk * Mk, &Lstar, Nk, &Lstarwcheck);
30: PetscDTAltVPullback(N, M, L, negative ? -k : k, w, Lstarw);
31: PetscDTAltVPullbackMatrix(N, M, L, negative ? -k : k, Lstar);
32: if (negative) {
33: PetscReal *sLsw;
35: PetscMalloc1(Nk, &sLsw);
36: PetscDTAltVStar(N, N - k, 1, Lstarw, sLsw);
37: PetscDTAltVApply(N, k, sLsw, x, &Lstarwx);
38: PetscFree(sLsw);
39: } else {
40: PetscDTAltVApply(N, k, Lstarw, x, &Lstarwx);
41: }
42: for (l = 0; l < k; l++) {
43: for (i = 0; i < M; i++) {
44: PetscReal sum = 0.;
46: for (j = 0; j < N; j++) sum += L[i * N + j] * x[l * N + j];
47: Lx[l * M + i] = sum;
48: }
49: }
50: diffMat = 0.;
51: normMat = 0.;
52: for (i = 0; i < Nk; i++) {
53: PetscReal sum = 0.;
54: for (j = 0; j < Mk; j++) {
55: sum += Lstar[i * Mk + j] * w[j];
56: }
57: Lstarwcheck[i] = sum;
58: diffMat += PetscSqr(PetscAbsReal(Lstarwcheck[i] - Lstarw[i]));
59: normMat += PetscSqr(Lstarwcheck[i]) + PetscSqr(Lstarw[i]);
60: }
61: diffMat = PetscSqrtReal(diffMat);
62: normMat = PetscSqrtReal(normMat);
63: if (verbose) {
64: PetscViewerASCIIPrintf(viewer, "L:\n");
65: PetscViewerASCIIPushTab(viewer);
66: if (M*N > 0) {PetscRealView(M*N, L, viewer);}
67: PetscViewerASCIIPopTab(viewer);
69: PetscViewerASCIIPrintf(viewer, "L*:\n");
70: PetscViewerASCIIPushTab(viewer);
71: if (Nk * Mk > 0) {PetscRealView(Nk * Mk, Lstar, viewer);}
72: PetscViewerASCIIPopTab(viewer);
74: PetscViewerASCIIPrintf(viewer, "L*w:\n");
75: PetscViewerASCIIPushTab(viewer);
76: if (Nk > 0) {PetscRealView(Nk, Lstarw, viewer);}
77: PetscViewerASCIIPopTab(viewer);
78: }
79: PetscDTAltVApply(M, k, ww, Lx, &wLx);
80: diff = PetscAbsReal(wLx - Lstarwx);
81: if (diff > 10. * PETSC_SMALL * (PetscAbsReal(wLx) + PetscAbsReal(Lstarwx))) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "pullback check: pullback does not commute with Section 1.5 Writing Application Codes with PETSc: w(Lx)(%g) != (L* w)(x)(%g)", wLx, Lstarwx);
82: if (diffMat > PETSC_SMALL * normMat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "pullback check: pullback matrix does match matrix free result");
83: PetscFree2(Lstar, Lstarwcheck);
84: PetscFree2(Lstarw, Lx);
85: PetscFree(walloc);
86: return(0);
87: }
89: int main(int argc, char **argv)
90: {
91: PetscInt i, numTests = 5, n[5] = {0, 1, 2, 3, 4};
92: PetscBool verbose = PETSC_FALSE;
93: PetscRandom rand;
94: PetscViewer viewer;
97: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
98: PetscOptionsBegin(PETSC_COMM_WORLD,"","Options for exterior algebra tests","none");
99: PetscOptionsIntArray("-N", "Up to 5 vector space dimensions to test","ex6.c",n,&numTests,NULL);
100: PetscOptionsBool("-verbose", "Verbose test output","ex6.c",verbose,&verbose,NULL);
101: PetscOptionsEnd();
102: PetscRandomCreate(PETSC_COMM_SELF, &rand);
103: PetscRandomSetInterval(rand, -1., 1.);
104: PetscRandomSetFromOptions(rand);
105: if (!numTests) numTests = 5;
106: viewer = PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD);
107: for (i = 0; i < numTests; i++) {
108: PetscInt k, N = n[i];
110: if (verbose) {PetscViewerASCIIPrintf(viewer, "N = %D:\n", N);}
111: PetscViewerASCIIPushTab(viewer);
113: if (verbose) {
114: PetscInt *perm;
115: PetscInt fac = 1;
117: PetscMalloc1(N, &perm);
119: for (k = 1; k <= N; k++) fac *= k;
120: PetscViewerASCIIPrintf(viewer, "Permutations of %D:\n", N);
121: PetscViewerASCIIPushTab(viewer);
122: for (k = 0; k < fac; k++) {
123: PetscBool isOdd, isOddCheck;
124: PetscInt j, kCheck;
126: PetscDTEnumPerm(N, k, perm, &isOdd);
127: PetscViewerASCIIPrintf(viewer, "%D:", k);
128: for (j = 0; j < N; j++) {
129: PetscPrintf(PETSC_COMM_WORLD, " %D", perm[j]);
130: }
131: PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even");
132: PetscDTPermIndex(N, perm, &kCheck, &isOddCheck);
133: if (kCheck != k || isOddCheck != isOdd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PetscDTEnumPerm / PetscDTPermIndex mismatch for (%D, %D)\n", N, k);
134: }
135: PetscViewerASCIIPopTab(viewer);
136: PetscFree(perm);
137: }
138: for (k = 0; k <= N; k++) {
139: PetscInt j, Nk, M;
140: PetscReal *w, *v, wv;
141: PetscInt *subset;
143: PetscDTBinomialInt(N, k, &Nk);
144: if (verbose) {PetscViewerASCIIPrintf(viewer, "k = %D:\n", k);}
145: PetscViewerASCIIPushTab(viewer);
146: if (verbose) {PetscViewerASCIIPrintf(viewer, "(%D choose %D): %D\n", N, k, Nk);}
148: /* Test subset and complement enumeration */
149: PetscMalloc1(N, &subset);
150: PetscViewerASCIIPushTab(viewer);
151: for (j = 0; j < Nk; j++) {
152: PetscBool isOdd, isOddCheck;
153: PetscInt jCheck, kCheck;
155: PetscDTEnumSplit(N, k, j, subset, &isOdd);
156: PetscDTPermIndex(N, subset, &kCheck, &isOddCheck);
157: if (isOddCheck != isOdd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PetscDTEnumSplit sign does not mmatch PetscDTPermIndex sign");
158: if (verbose) {
159: PetscInt l;
161: PetscViewerASCIIPrintf(viewer, "subset %D:", j);
162: for (l = 0; l < k; l++) {
163: PetscPrintf(PETSC_COMM_WORLD, " %D", subset[l]);
164: }
165: PetscPrintf(PETSC_COMM_WORLD, " |");
166: for (l = k; l < N; l++) {
167: PetscPrintf(PETSC_COMM_WORLD, " %D", subset[l]);
168: }
169: PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even");
170: }
171: PetscDTSubsetIndex(N, k, subset, &jCheck);
172: if (jCheck != j) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "jCheck (%D) != j (%D)", jCheck, j);
173: }
174: PetscViewerASCIIPopTab(viewer);
175: PetscFree(subset);
177: /* Make a random k form */
178: PetscMalloc1(Nk, &w);
179: for (j = 0; j < Nk; j++) {PetscRandomGetValueReal(rand, &w[j]);}
180: /* Make a set of random vectors */
181: PetscMalloc1(N*k, &v);
182: for (j = 0; j < N*k; j++) {PetscRandomGetValueReal(rand, &v[j]);}
184: PetscDTAltVApply(N, k, w, v, &wv);
186: if (verbose) {
187: PetscViewerASCIIPrintf(viewer, "w:\n");
188: PetscViewerASCIIPushTab(viewer);
189: if (Nk) {PetscRealView(Nk, w, viewer);}
190: PetscViewerASCIIPopTab(viewer);
191: PetscViewerASCIIPrintf(viewer, "v:\n");
192: PetscViewerASCIIPushTab(viewer);
193: if (N*k > 0) {PetscRealView(N*k, v, viewer);}
194: PetscViewerASCIIPopTab(viewer);
195: PetscViewerASCIIPrintf(viewer, "w(v): %g\n", (double) wv);
196: }
198: /* sanity checks */
199: if (k == 1) { /* 1-forms are functionals (dot products) */
200: PetscInt l;
201: PetscReal wvcheck = 0.;
202: PetscReal diff;
204: for (l = 0; l < N; l++) wvcheck += w[l] * v[l];
205: diff = PetscSqrtReal(PetscSqr(wvcheck - wv));
206: if (diff >= PETSC_SMALL * (PetscAbsReal(wv) + PetscAbsReal(wvcheck))) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "1-form / dot product equivalence: wvcheck (%g) != wv (%g)", (double) wvcheck, (double) wv);
207: }
208: if (k == N && N < 5) { /* n-forms are scaled determinants */
209: PetscReal det, wvcheck, diff;
211: switch (k) {
212: case 0:
213: det = 1.;
214: break;
215: case 1:
216: det = v[0];
217: break;
218: case 2:
219: det = v[0] * v[3] - v[1] * v[2];
220: break;
221: case 3:
222: det = v[0] * (v[4] * v[8] - v[5] * v[7]) +
223: v[1] * (v[5] * v[6] - v[3] * v[8]) +
224: v[2] * (v[3] * v[7] - v[4] * v[6]);
225: break;
226: case 4:
227: det = v[0] * (v[5] * (v[10] * v[15] - v[11] * v[14]) +
228: v[6] * (v[11] * v[13] - v[ 9] * v[15]) +
229: v[7] * (v[ 9] * v[14] - v[10] * v[13])) -
230: v[1] * (v[4] * (v[10] * v[15] - v[11] * v[14]) +
231: v[6] * (v[11] * v[12] - v[ 8] * v[15]) +
232: v[7] * (v[ 8] * v[14] - v[10] * v[12])) +
233: v[2] * (v[4] * (v[ 9] * v[15] - v[11] * v[13]) +
234: v[5] * (v[11] * v[12] - v[ 8] * v[15]) +
235: v[7] * (v[ 8] * v[13] - v[ 9] * v[12])) -
236: v[3] * (v[4] * (v[ 9] * v[14] - v[10] * v[13]) +
237: v[5] * (v[10] * v[12] - v[ 8] * v[14]) +
238: v[6] * (v[ 8] * v[13] - v[ 9] * v[12]));
239: break;
240: default:
241: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "invalid k");
242: }
243: wvcheck = det * w[0];
244: diff = PetscSqrtReal(PetscSqr(wvcheck - wv));
245: if (diff >= PETSC_SMALL * (PetscAbsReal(wv) + PetscAbsReal(wvcheck))) SETERRQ3(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "n-form / determinant equivalence: wvcheck (%g) != wv (%g) %g", (double) wvcheck, (double) wv, (double) diff);
246: }
247: if (k > 0) { /* k-forms are linear in each component */
248: PetscReal alpha;
249: PetscReal *x, *axv, wx, waxv, waxvcheck;
250: PetscReal diff;
251: PetscReal rj;
252: PetscInt l;
254: PetscMalloc2(N * k, &x, N * k, &axv);
255: PetscRandomGetValueReal(rand, &alpha);
256: PetscRandomSetInterval(rand, 0, k);
257: PetscRandomGetValueReal(rand, &rj);
258: j = (PetscInt) rj;
259: PetscRandomSetInterval(rand, -1., 1.);
260: for (l = 0; l < N*k; l++) x[l] = v[l];
261: for (l = 0; l < N*k; l++) axv[l] = v[l];
262: for (l = 0; l < N; l++) {
263: PetscReal val;
265: PetscRandomGetValueReal(rand, &val);
266: x[j * N + l] = val;
267: axv[j * N + l] += alpha * val;
268: }
269: PetscDTAltVApply(N, k, w, x, &wx);
270: PetscDTAltVApply(N, k, w, axv, &waxv);
271: waxvcheck = alpha * wx + wv;
272: diff = waxv - waxvcheck;
273: if (PetscAbsReal(diff) > 10. * PETSC_SMALL * (PetscAbsReal(waxv) + PetscAbsReal(waxvcheck))) SETERRQ3(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "linearity check: component %D, waxvcheck (%g) != waxv (%g)", j, (double) waxvcheck, (double) waxv);
274: PetscFree2(x,axv);
275: }
276: if (k > 1) { /* k-forms are antisymmetric */
277: PetscReal rj, rl, *swapv, wswapv, diff;
278: PetscInt l, m;
280: PetscRandomSetInterval(rand, 0, k);
281: PetscRandomGetValueReal(rand, &rj);
282: j = (PetscInt) rj;
283: l = j;
284: while (l == j) {
285: PetscRandomGetValueReal(rand, &rl);
286: l = (PetscInt) rl;
287: }
288: PetscRandomSetInterval(rand, -1., 1.);
289: PetscMalloc1(N * k, &swapv);
290: for (m = 0; m < N * k; m++) swapv[m] = v[m];
291: for (m = 0; m < N; m++) {
292: swapv[j * N + m] = v[l * N + m];
293: swapv[l * N + m] = v[j * N + m];
294: }
295: PetscDTAltVApply(N, k, w, swapv, &wswapv);
296: diff = PetscAbsReal(wswapv + wv);
297: if (diff > PETSC_SMALL * (PetscAbsReal(wswapv) + PetscAbsReal(wv))) SETERRQ4(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "antisymmetry check: components %D & %D, wswapv (%g) != -wv (%g)", j, l, (double) wswapv, (double) wv);
298: PetscFree(swapv);
299: }
300: for (j = 0; j <= k && j + k <= N; j++) { /* wedge product */
301: PetscInt Nj, Njk, l, JKj;
302: PetscReal *u, *uWw, *uWwcheck, *uWwmat, *x, *xsplit, uWwx, uWwxcheck, diff, norm;
303: PetscInt *split;
305: if (verbose) {PetscViewerASCIIPrintf(viewer, "wedge j = %D:\n", j);}
306: PetscViewerASCIIPushTab(viewer);
307: PetscDTBinomialInt(N, j, &Nj);
308: PetscDTBinomialInt(N, j+k, &Njk);
309: PetscMalloc4(Nj, &u, Njk, &uWw, N*(j+k), &x, N*(j+k), &xsplit);
310: PetscMalloc1(j+k,&split);
311: for (l = 0; l < Nj; l++) {PetscRandomGetValueReal(rand, &u[l]);}
312: for (l = 0; l < N*(j+k); l++) {PetscRandomGetValueReal(rand, &x[l]);}
313: PetscDTAltVWedge(N, j, k, u, w, uWw);
314: PetscDTAltVApply(N, j+k, uWw, x, &uWwx);
315: if (verbose) {
316: PetscViewerASCIIPrintf(viewer, "u:\n");
317: PetscViewerASCIIPushTab(viewer);
318: if (Nj) {PetscRealView(Nj, u, viewer);}
319: PetscViewerASCIIPopTab(viewer);
320: PetscViewerASCIIPrintf(viewer, "u wedge w:\n");
321: PetscViewerASCIIPushTab(viewer);
322: if (Njk) {PetscRealView(Njk, uWw, viewer);}
323: PetscViewerASCIIPopTab(viewer);
324: PetscViewerASCIIPrintf(viewer, "x:\n");
325: PetscViewerASCIIPushTab(viewer);
326: if (N*(j+k) > 0) {PetscRealView(N*(j+k), x, viewer);}
327: PetscViewerASCIIPopTab(viewer);
328: PetscViewerASCIIPrintf(viewer, "u wedge w(x): %g\n", (double) uWwx);
329: }
330: /* verify wedge formula */
331: uWwxcheck = 0.;
332: PetscDTBinomialInt(j+k, j, &JKj);
333: for (l = 0; l < JKj; l++) {
334: PetscBool isOdd;
335: PetscReal ux, wx;
336: PetscInt m, p;
338: PetscDTEnumSplit(j+k, j, l, split, &isOdd);
339: for (m = 0; m < j+k; m++) {for (p = 0; p < N; p++) {xsplit[m * N + p] = x[split[m] * N + p];}}
340: PetscDTAltVApply(N, j, u, xsplit, &ux);
341: PetscDTAltVApply(N, k, w, &xsplit[j*N], &wx);
342: uWwxcheck += isOdd ? -(ux * wx) : (ux * wx);
343: }
344: diff = PetscAbsReal(uWwx - uWwxcheck);
345: if (diff > 10. * PETSC_SMALL * (PetscAbsReal(uWwx) + PetscAbsReal(uWwxcheck))) SETERRQ4(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "wedge check: forms %D & %D, uWwxcheck (%g) != uWwx (%g)", j, k, (double) uWwxcheck, (double) uWwx);
346: PetscFree(split);
347: PetscMalloc2(Nk * Njk, &uWwmat, Njk, &uWwcheck);
348: PetscDTAltVWedgeMatrix(N, j, k, u, uWwmat);
349: if (verbose) {
350: PetscViewerASCIIPrintf(viewer, "(u wedge):\n");
351: PetscViewerASCIIPushTab(viewer);
352: if ((Nk * Njk) > 0) {PetscRealView(Nk * Njk, uWwmat, viewer);}
353: PetscViewerASCIIPopTab(viewer);
354: }
355: diff = 0.;
356: norm = 0.;
357: for (l = 0; l < Njk; l++) {
358: PetscInt m;
359: PetscReal sum = 0.;
361: for (m = 0; m < Nk; m++) sum += uWwmat[l * Nk + m] * w[m];
362: uWwcheck[l] = sum;
363: diff += PetscSqr(uWwcheck[l] - uWw[l]);
364: norm += PetscSqr(uWwcheck[l]) + PetscSqr(uWw[l]);
365: }
366: diff = PetscSqrtReal(diff);
367: norm = PetscSqrtReal(norm);
368: if (diff > PETSC_SMALL * norm) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "wedge matrix check: wedge matrix Section 1.5 Writing Application Codes with PETSc does not match wedge direct Section 1.5 Writing Application Codes with PETSc");
369: PetscFree2(uWwmat, uWwcheck);
370: PetscFree4(u, uWw, x, xsplit);
371: PetscViewerASCIIPopTab(viewer);
372: }
373: for (M = PetscMax(1,k); M <= N; M++) { /* pullback */
374: PetscReal *L, *u, *x;
375: PetscInt Mk, l;
377: PetscDTBinomialInt(M, k, &Mk);
378: PetscMalloc3(M*N, &L, Mk, &u, M*k, &x);
379: for (l = 0; l < M*N; l++) {PetscRandomGetValueReal(rand, &L[l]);}
380: for (l = 0; l < Mk; l++) {PetscRandomGetValueReal(rand, &u[l]);}
381: for (l = 0; l < M*k; l++) {PetscRandomGetValueReal(rand, &x[l]);}
382: if (verbose) {PetscViewerASCIIPrintf(viewer, "pullback M = %D:\n", M);}
383: PetscViewerASCIIPushTab(viewer);
384: CheckPullback(M, N, L, k, w, x, verbose, viewer);
385: if (M != N) {CheckPullback(N, M, L, k, u, v, PETSC_FALSE, viewer);}
386: PetscViewerASCIIPopTab(viewer);
387: if ((k % N) && (N > 1)) {
388: if (verbose) {PetscViewerASCIIPrintf(viewer, "negative pullback M = %D:\n", M);}
389: PetscViewerASCIIPushTab(viewer);
390: CheckPullback(M, N, L, -k, w, x, verbose, viewer);
391: if (M != N) {CheckPullback(N, M, L, -k, u, v, PETSC_FALSE, viewer);}
392: PetscViewerASCIIPopTab(viewer);
393: }
394: PetscFree3(L, u, x);
395: }
396: if (k > 0) { /* Interior */
397: PetscInt Nkm, l, m;
398: PetscReal *wIntv0, *wIntv0check, wvcheck, diff, diffMat, normMat;
399: PetscReal *intv0mat, *matcheck;
400: PetscInt (*indices)[3];
402: PetscDTBinomialInt(N, k-1, &Nkm);
403: PetscMalloc5(Nkm, &wIntv0, Nkm, &wIntv0check, Nk * Nkm, &intv0mat, Nk * Nkm, &matcheck, Nk * k, &indices);
404: PetscDTAltVInterior(N, k, w, v, wIntv0);
405: PetscDTAltVInteriorMatrix(N, k, v, intv0mat);
406: PetscDTAltVInteriorPattern(N, k, indices);
407: if (verbose) {
408: PetscViewerASCIIPrintf(viewer, "interior product matrix pattern:\n");
409: PetscViewerASCIIPushTab(viewer);
410: for (l = 0; l < Nk * k; l++) {
411: PetscInt row = indices[l][0];
412: PetscInt col = indices[l][1];
413: PetscInt x = indices[l][2];
415: PetscViewerASCIIPrintf(viewer,"intV[%D,%D] = %sV[%D]\n", row, col, x < 0 ? "-" : " ", x < 0 ? -(x + 1) : x);
416: }
417: PetscViewerASCIIPopTab(viewer);
418: }
419: for (l = 0; l < Nkm * Nk; l++) matcheck[l] = 0.;
420: for (l = 0; l < Nk * k; l++) {
421: PetscInt row = indices[l][0];
422: PetscInt col = indices[l][1];
423: PetscInt x = indices[l][2];
425: if (x < 0) {
426: matcheck[row * Nk + col] = -v[-(x+1)];
427: } else {
428: matcheck[row * Nk + col] = v[x];
429: }
430: }
431: diffMat = 0.;
432: normMat = 0.;
433: for (l = 0; l < Nkm * Nk; l++) {
434: diffMat += PetscSqr(PetscAbsReal(matcheck[l] - intv0mat[l]));
435: normMat += PetscSqr(matcheck[l]) + PetscSqr(intv0mat[l]);
436: }
437: diffMat = PetscSqrtReal(diffMat);
438: normMat = PetscSqrtReal(normMat);
439: if (diffMat > PETSC_SMALL * normMat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior product check: matrix pattern does not match matrix");
440: diffMat = 0.;
441: normMat = 0.;
442: for (l = 0; l < Nkm; l++) {
443: PetscReal sum = 0.;
445: for (m = 0; m < Nk; m++) sum += intv0mat[l * Nk + m] * w[m];
446: wIntv0check[l] = sum;
448: diffMat += PetscSqr(PetscAbsReal(wIntv0check[l] - wIntv0[l]));
449: normMat += PetscSqr(wIntv0check[l]) + PetscSqr(wIntv0[l]);
450: }
451: diffMat = PetscSqrtReal(diffMat);
452: normMat = PetscSqrtReal(normMat);
453: if (diffMat > PETSC_SMALL * normMat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior product check: Section 1.5 Writing Application Codes with PETSc does not match matrix");
454: if (verbose) {
455: PetscViewerASCIIPrintf(viewer, "(w int v_0):\n");
456: PetscViewerASCIIPushTab(viewer);
457: if (Nkm) {PetscRealView(Nkm, wIntv0, viewer);}
458: PetscViewerASCIIPopTab(viewer);
460: PetscViewerASCIIPrintf(viewer, "(int v_0):\n");
461: PetscViewerASCIIPushTab(viewer);
462: if (Nk * Nkm > 0) {PetscRealView(Nk * Nkm, intv0mat, viewer);}
463: PetscViewerASCIIPopTab(viewer);
464: }
465: PetscDTAltVApply(N, k - 1, wIntv0, &v[N], &wvcheck);
466: diff = PetscSqrtReal(PetscSqr(wvcheck - wv));
467: if (diff >= PETSC_SMALL * (PetscAbsReal(wv) + PetscAbsReal(wvcheck))) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior product check: (w Int v0)(v_rem) (%g) != w(v) (%g)", (double) wvcheck, (double) wv);
468: PetscFree5(wIntv0,wIntv0check,intv0mat,matcheck,indices);
469: }
470: if (k >= N - k) { /* Hodge star */
471: PetscReal *u, *starw, *starstarw, wu, starwdotu;
472: PetscReal diff, norm;
473: PetscBool isOdd;
474: PetscInt l;
476: isOdd = (PetscBool) ((k * (N - k)) & 1);
477: PetscMalloc3(Nk, &u, Nk, &starw, Nk, &starstarw);
478: PetscDTAltVStar(N, k, 1, w, starw);
479: PetscDTAltVStar(N, N-k, 1, starw, starstarw);
480: if (verbose) {
481: PetscViewerASCIIPrintf(viewer, "star w:\n");
482: PetscViewerASCIIPushTab(viewer);
483: if (Nk) {PetscRealView(Nk, starw, viewer);}
484: PetscViewerASCIIPopTab(viewer);
486: PetscViewerASCIIPrintf(viewer, "star star w:\n");
487: PetscViewerASCIIPushTab(viewer);
488: if (Nk) {PetscRealView(Nk, starstarw, viewer);}
489: PetscViewerASCIIPopTab(viewer);
490: }
491: for (l = 0; l < Nk; l++) {PetscRandomGetValueReal(rand,&u[l]);}
492: PetscDTAltVWedge(N, k, N - k, w, u, &wu);
493: starwdotu = 0.;
494: for (l = 0; l < Nk; l++) starwdotu += starw[l] * u[l];
495: diff = PetscAbsReal(wu - starwdotu);
496: if (diff > PETSC_SMALL * (PetscAbsReal(wu) + PetscAbsReal(starwdotu))) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Hodge star check: (star w, u) (%g) != (w wedge u) (%g)", (double) starwdotu, (double) wu);
498: diff = 0.;
499: norm = 0.;
500: for (l = 0; l < Nk; l++) {
501: diff += PetscSqr(w[l] - (isOdd ? -starstarw[l] : starstarw[l]));
502: norm += PetscSqr(w[l]) + PetscSqr(starstarw[l]);
503: }
504: diff = PetscSqrtReal(diff);
505: norm = PetscSqrtReal(norm);
506: if (diff > PETSC_SMALL * norm) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Hodge star check: star(star(w)) != (-1)^(N*(N-k)) w");
507: PetscFree3(u, starw, starstarw);
508: }
509: PetscFree(v);
510: PetscFree(w);
511: PetscViewerASCIIPopTab(viewer);
512: }
513: PetscViewerASCIIPopTab(viewer);
514: }
515: PetscRandomDestroy(&rand);
516: PetscFinalize();
517: return ierr;
518: }
520: /*TEST
521: test:
522: suffix: 1234
523: args: -verbose
524: test:
525: suffix: 56
526: args: -N 5,6
527: TEST*/