1 static char help[] = "Test the PetscDTAltV interface for k-forms (alternating k-linear maps).\n\n";
2
3 #include <petscviewer.h>
4 #include <petscdt.h>
5
CheckPullback(PetscInt N,PetscInt M,const PetscReal * L,PetscInt k,const PetscReal * w,PetscReal * x,PetscBool verbose,PetscViewer viewer)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);
14
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;
32
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.;
43
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));
64
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));
69
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 }
84
main(int argc,char ** argv)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;
91
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];
105
106 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "N = %" PetscInt_FMT ":\n", N));
107 PetscCall(PetscViewerASCIIPushTab(viewer));
108
109 if (verbose) {
110 PetscInt *perm;
111 PetscInt fac = 1;
112
113 PetscCall(PetscMalloc1(N, &perm));
114
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;
121
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;
136
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));
141
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;
148
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;
154
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));
166
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]));
173
174 PetscCall(PetscDTAltVApply(N, k, w, v, &wv));
175
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 }
187
188 /* sanity checks */
189 if (k == 1) { /* 1-forms are functionals (dot products) */
190 PetscInt l;
191 PetscReal wvcheck = 0.;
192 PetscReal diff;
193
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;
200
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;
230
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;
241
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;
256
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;
281
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;
314
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, PetscSafePointerPlusOffset(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.;
339
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;
355
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];
380
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];
393
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];
403
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.;
423
424 for (m = 0; m < Nk; m++) sum += intv0mat[l * Nk + m] * w[m];
425 wIntv0check[l] = sum;
426
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));
438
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;
454
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));
464
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);
476
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 }
498
499 /*TEST
500 test:
501 suffix: 1234
502 args: -verbose
503 test:
504 suffix: 56
505 args: -N 5,6
506 output_file: output/empty.out
507 TEST*/
508