Actual source code: ex7.c

petsc-3.13.4 2020-08-01
Report Typos and Errors
  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*/