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