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:   k = PetscAbsInt(k);
 16:   PetscDTBinomialInt(N, k, &Nk);
 17:   PetscDTBinomialInt(M, k, &Mk);
 18:   if (negative) {
 19:     PetscMalloc1(Mk, &walloc);
 20:     PetscDTAltVStar(M, M - k, 1, w, walloc);
 21:     ww = walloc;
 22:   } else {
 23:     ww = w;
 24:   }
 25:   PetscMalloc2(Nk, &Lstarw, (M * k), &Lx);
 26:   PetscMalloc2(Nk * Mk, &Lstar, Nk, &Lstarwcheck);
 27:   PetscDTAltVPullback(N, M, L, negative ? -k : k, w, Lstarw);
 28:   PetscDTAltVPullbackMatrix(N, M, L, negative ? -k : k, Lstar);
 29:   if (negative) {
 30:     PetscReal *sLsw;

 32:     PetscMalloc1(Nk, &sLsw);
 33:     PetscDTAltVStar(N, N - k, 1, Lstarw, sLsw);
 34:     PetscDTAltVApply(N, k, sLsw, x, &Lstarwx);
 35:     PetscFree(sLsw);
 36:   } else {
 37:     PetscDTAltVApply(N, k, Lstarw, x, &Lstarwx);
 38:   }
 39:   for (l = 0; l < k; l++) {
 40:     for (i = 0; i < M; i++) {
 41:       PetscReal sum = 0.;

 43:       for (j = 0; j < N; j++) sum += L[i * N + j] * x[l * N + j];
 44:       Lx[l * M + i] = sum;
 45:     }
 46:   }
 47:   diffMat = 0.;
 48:   normMat = 0.;
 49:   for (i = 0; i < Nk; i++) {
 50:     PetscReal sum = 0.;
 51:     for (j = 0; j < Mk; j++) sum += Lstar[i * Mk + j] * w[j];
 52:     Lstarwcheck[i] = sum;
 53:     diffMat += PetscSqr(PetscAbsReal(Lstarwcheck[i] - Lstarw[i]));
 54:     normMat += PetscSqr(Lstarwcheck[i]) + PetscSqr(Lstarw[i]);
 55:   }
 56:   diffMat = PetscSqrtReal(diffMat);
 57:   normMat = PetscSqrtReal(normMat);
 58:   if (verbose) {
 59:     PetscViewerASCIIPrintf(viewer, "L:\n");
 60:     PetscViewerASCIIPushTab(viewer);
 61:     if (M * N > 0) PetscRealView(M * N, L, viewer);
 62:     PetscViewerASCIIPopTab(viewer);

 64:     PetscViewerASCIIPrintf(viewer, "L*:\n");
 65:     PetscViewerASCIIPushTab(viewer);
 66:     if (Nk * Mk > 0) PetscRealView(Nk * Mk, Lstar, viewer);
 67:     PetscViewerASCIIPopTab(viewer);

 69:     PetscViewerASCIIPrintf(viewer, "L*w:\n");
 70:     PetscViewerASCIIPushTab(viewer);
 71:     if (Nk > 0) PetscRealView(Nk, Lstarw, viewer);
 72:     PetscViewerASCIIPopTab(viewer);
 73:   }
 74:   PetscDTAltVApply(M, k, ww, Lx, &wLx);
 75:   diff = PetscAbsReal(wLx - Lstarwx);
 78:   PetscFree2(Lstar, Lstarwcheck);
 79:   PetscFree2(Lstarw, Lx);
 80:   PetscFree(walloc);
 81:   return 0;
 82: }

 84: int main(int argc, char **argv)
 85: {
 86:   PetscInt    i, numTests = 5, n[5] = {0, 1, 2, 3, 4};
 87:   PetscBool   verbose = PETSC_FALSE;
 88:   PetscRandom rand;
 89:   PetscViewer viewer;

 92:   PetscInitialize(&argc, &argv, NULL, help);
 93:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for exterior algebra tests", "none");
 94:   PetscOptionsIntArray("-N", "Up to 5 vector space dimensions to test", "ex7.c", n, &numTests, NULL);
 95:   PetscOptionsBool("-verbose", "Verbose test output", "ex7.c", verbose, &verbose, NULL);
 96:   PetscOptionsEnd();
 97:   PetscRandomCreate(PETSC_COMM_SELF, &rand);
 98:   PetscRandomSetInterval(rand, -1., 1.);
 99:   PetscRandomSetFromOptions(rand);
100:   if (!numTests) numTests = 5;
101:   viewer = PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD);
102:   for (i = 0; i < numTests; i++) {
103:     PetscInt k, N = n[i];

105:     if (verbose) PetscViewerASCIIPrintf(viewer, "N = %" PetscInt_FMT ":\n", N);
106:     PetscViewerASCIIPushTab(viewer);

108:     if (verbose) {
109:       PetscInt *perm;
110:       PetscInt  fac = 1;

112:       PetscMalloc1(N, &perm);

114:       for (k = 1; k <= N; k++) fac *= k;
115:       PetscViewerASCIIPrintf(viewer, "Permutations of %" PetscInt_FMT ":\n", N);
116:       PetscViewerASCIIPushTab(viewer);
117:       for (k = 0; k < fac; k++) {
118:         PetscBool isOdd, isOddCheck;
119:         PetscInt  j, kCheck;

121:         PetscDTEnumPerm(N, k, perm, &isOdd);
122:         PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ":", k);
123:         for (j = 0; j < N; j++) PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, perm[j]);
124:         PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even");
125:         PetscDTPermIndex(N, perm, &kCheck, &isOddCheck);
127:       }
128:       PetscViewerASCIIPopTab(viewer);
129:       PetscFree(perm);
130:     }
131:     for (k = 0; k <= N; k++) {
132:       PetscInt   j, Nk, M;
133:       PetscReal *w, *v, wv;
134:       PetscInt  *subset;

136:       PetscDTBinomialInt(N, k, &Nk);
137:       if (verbose) PetscViewerASCIIPrintf(viewer, "k = %" PetscInt_FMT ":\n", k);
138:       PetscViewerASCIIPushTab(viewer);
139:       if (verbose) PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT " choose %" PetscInt_FMT "): %" PetscInt_FMT "\n", N, k, Nk);

141:       /* Test subset and complement enumeration */
142:       PetscMalloc1(N, &subset);
143:       PetscViewerASCIIPushTab(viewer);
144:       for (j = 0; j < Nk; j++) {
145:         PetscBool isOdd, isOddCheck;
146:         PetscInt  jCheck, kCheck;

148:         PetscDTEnumSplit(N, k, j, subset, &isOdd);
149:         PetscDTPermIndex(N, subset, &kCheck, &isOddCheck);
151:         if (verbose) {
152:           PetscInt l;

154:           PetscViewerASCIIPrintf(viewer, "subset %" PetscInt_FMT ":", j);
155:           for (l = 0; l < k; l++) PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, subset[l]);
156:           PetscPrintf(PETSC_COMM_WORLD, " |");
157:           for (l = k; l < N; l++) PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, subset[l]);
158:           PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even");
159:         }
160:         PetscDTSubsetIndex(N, k, subset, &jCheck);
162:       }
163:       PetscViewerASCIIPopTab(viewer);
164:       PetscFree(subset);

166:       /* Make a random k form */
167:       PetscMalloc1(Nk, &w);
168:       for (j = 0; j < Nk; j++) PetscRandomGetValueReal(rand, &w[j]);
169:       /* Make a set of random vectors */
170:       PetscMalloc1(N * k, &v);
171:       for (j = 0; j < N * k; j++) PetscRandomGetValueReal(rand, &v[j]);

173:       PetscDTAltVApply(N, k, w, v, &wv);

175:       if (verbose) {
176:         PetscViewerASCIIPrintf(viewer, "w:\n");
177:         PetscViewerASCIIPushTab(viewer);
178:         if (Nk) PetscRealView(Nk, w, viewer);
179:         PetscViewerASCIIPopTab(viewer);
180:         PetscViewerASCIIPrintf(viewer, "v:\n");
181:         PetscViewerASCIIPushTab(viewer);
182:         if (N * k > 0) PetscRealView(N * k, v, viewer);
183:         PetscViewerASCIIPopTab(viewer);
184:         PetscViewerASCIIPrintf(viewer, "w(v): %g\n", (double)wv);
185:       }

187:       /* sanity checks */
188:       if (k == 1) { /* 1-forms are functionals (dot products) */
189:         PetscInt  l;
190:         PetscReal wvcheck = 0.;
191:         PetscReal diff;

193:         for (l = 0; l < N; l++) wvcheck += w[l] * v[l];
194:         diff = PetscSqrtReal(PetscSqr(wvcheck - wv));
196:       }
197:       if (k == N && N < 5) { /* n-forms are scaled determinants */
198:         PetscReal det, wvcheck, diff;

200:         switch (k) {
201:         case 0:
202:           det = 1.;
203:           break;
204:         case 1:
205:           det = v[0];
206:           break;
207:         case 2:
208:           det = v[0] * v[3] - v[1] * v[2];
209:           break;
210:         case 3:
211:           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]);
212:           break;
213:         case 4:
214:           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]));
215:           break;
216:         default:
217:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "invalid k");
218:         }
219:         wvcheck = det * w[0];
220:         diff    = PetscSqrtReal(PetscSqr(wvcheck - wv));
222:       }
223:       if (k > 0) { /* k-forms are linear in each component */
224:         PetscReal  alpha;
225:         PetscReal *x, *axv, wx, waxv, waxvcheck;
226:         PetscReal  diff;
227:         PetscReal  rj;
228:         PetscInt   l;

230:         PetscMalloc2(N * k, &x, N * k, &axv);
231:         PetscRandomGetValueReal(rand, &alpha);
232:         PetscRandomSetInterval(rand, 0, k);
233:         PetscRandomGetValueReal(rand, &rj);
234:         j = (PetscInt)rj;
235:         PetscRandomSetInterval(rand, -1., 1.);
236:         for (l = 0; l < N * k; l++) x[l] = v[l];
237:         for (l = 0; l < N * k; l++) axv[l] = v[l];
238:         for (l = 0; l < N; l++) {
239:           PetscReal val;

241:           PetscRandomGetValueReal(rand, &val);
242:           x[j * N + l] = val;
243:           axv[j * N + l] += alpha * val;
244:         }
245:         PetscDTAltVApply(N, k, w, x, &wx);
246:         PetscDTAltVApply(N, k, w, axv, &waxv);
247:         waxvcheck = alpha * wx + wv;
248:         diff      = waxv - waxvcheck;
250:         PetscFree2(x, axv);
251:       }
252:       if (k > 1) { /* k-forms are antisymmetric */
253:         PetscReal rj, rl, *swapv, wswapv, diff;
254:         PetscInt  l, m;

256:         PetscRandomSetInterval(rand, 0, k);
257:         PetscRandomGetValueReal(rand, &rj);
258:         j = (PetscInt)rj;
259:         l = j;
260:         while (l == j) {
261:           PetscRandomGetValueReal(rand, &rl);
262:           l = (PetscInt)rl;
263:         }
264:         PetscRandomSetInterval(rand, -1., 1.);
265:         PetscMalloc1(N * k, &swapv);
266:         for (m = 0; m < N * k; m++) swapv[m] = v[m];
267:         for (m = 0; m < N; m++) {
268:           swapv[j * N + m] = v[l * N + m];
269:           swapv[l * N + m] = v[j * N + m];
270:         }
271:         PetscDTAltVApply(N, k, w, swapv, &wswapv);
272:         diff = PetscAbsReal(wswapv + wv);
274:         PetscFree(swapv);
275:       }
276:       for (j = 0; j <= k && j + k <= N; j++) { /* wedge product */
277:         PetscInt   Nj, Njk, l, JKj;
278:         PetscReal *u, *uWw, *uWwcheck, *uWwmat, *x, *xsplit, uWwx, uWwxcheck, diff, norm;
279:         PetscInt  *split;

281:         if (verbose) PetscViewerASCIIPrintf(viewer, "wedge j = %" PetscInt_FMT ":\n", j);
282:         PetscViewerASCIIPushTab(viewer);
283:         PetscDTBinomialInt(N, j, &Nj);
284:         PetscDTBinomialInt(N, j + k, &Njk);
285:         PetscMalloc4(Nj, &u, Njk, &uWw, N * (j + k), &x, N * (j + k), &xsplit);
286:         PetscMalloc1(j + k, &split);
287:         for (l = 0; l < Nj; l++) PetscRandomGetValueReal(rand, &u[l]);
288:         for (l = 0; l < N * (j + k); l++) PetscRandomGetValueReal(rand, &x[l]);
289:         PetscDTAltVWedge(N, j, k, u, w, uWw);
290:         PetscDTAltVApply(N, j + k, uWw, x, &uWwx);
291:         if (verbose) {
292:           PetscViewerASCIIPrintf(viewer, "u:\n");
293:           PetscViewerASCIIPushTab(viewer);
294:           if (Nj) PetscRealView(Nj, u, viewer);
295:           PetscViewerASCIIPopTab(viewer);
296:           PetscViewerASCIIPrintf(viewer, "u wedge w:\n");
297:           PetscViewerASCIIPushTab(viewer);
298:           if (Njk) PetscRealView(Njk, uWw, viewer);
299:           PetscViewerASCIIPopTab(viewer);
300:           PetscViewerASCIIPrintf(viewer, "x:\n");
301:           PetscViewerASCIIPushTab(viewer);
302:           if (N * (j + k) > 0) PetscRealView(N * (j + k), x, viewer);
303:           PetscViewerASCIIPopTab(viewer);
304:           PetscViewerASCIIPrintf(viewer, "u wedge w(x): %g\n", (double)uWwx);
305:         }
306:         /* verify wedge formula */
307:         uWwxcheck = 0.;
308:         PetscDTBinomialInt(j + k, j, &JKj);
309:         for (l = 0; l < JKj; l++) {
310:           PetscBool isOdd;
311:           PetscReal ux, wx;
312:           PetscInt  m, p;

314:           PetscDTEnumSplit(j + k, j, l, split, &isOdd);
315:           for (m = 0; m < j + k; m++) {
316:             for (p = 0; p < N; p++) xsplit[m * N + p] = x[split[m] * N + p];
317:           }
318:           PetscDTAltVApply(N, j, u, xsplit, &ux);
319:           PetscDTAltVApply(N, k, w, &xsplit[j * N], &wx);
320:           uWwxcheck += isOdd ? -(ux * wx) : (ux * wx);
321:         }
322:         diff = PetscAbsReal(uWwx - uWwxcheck);
324:         PetscFree(split);
325:         PetscMalloc2(Nk * Njk, &uWwmat, Njk, &uWwcheck);
326:         PetscDTAltVWedgeMatrix(N, j, k, u, uWwmat);
327:         if (verbose) {
328:           PetscViewerASCIIPrintf(viewer, "(u wedge):\n");
329:           PetscViewerASCIIPushTab(viewer);
330:           if ((Nk * Njk) > 0) PetscRealView(Nk * Njk, uWwmat, viewer);
331:           PetscViewerASCIIPopTab(viewer);
332:         }
333:         diff = 0.;
334:         norm = 0.;
335:         for (l = 0; l < Njk; l++) {
336:           PetscInt  m;
337:           PetscReal sum = 0.;

339:           for (m = 0; m < Nk; m++) sum += uWwmat[l * Nk + m] * w[m];
340:           uWwcheck[l] = sum;
341:           diff += PetscSqr(uWwcheck[l] - uWw[l]);
342:           norm += PetscSqr(uWwcheck[l]) + PetscSqr(uWw[l]);
343:         }
344:         diff = PetscSqrtReal(diff);
345:         norm = PetscSqrtReal(norm);
347:         PetscFree2(uWwmat, uWwcheck);
348:         PetscFree4(u, uWw, x, xsplit);
349:         PetscViewerASCIIPopTab(viewer);
350:       }
351:       for (M = PetscMax(1, k); M <= N; M++) { /* pullback */
352:         PetscReal *L, *u, *x;
353:         PetscInt   Mk, l;

355:         PetscDTBinomialInt(M, k, &Mk);
356:         PetscMalloc3(M * N, &L, Mk, &u, M * k, &x);
357:         for (l = 0; l < M * N; l++) PetscRandomGetValueReal(rand, &L[l]);
358:         for (l = 0; l < Mk; l++) PetscRandomGetValueReal(rand, &u[l]);
359:         for (l = 0; l < M * k; l++) PetscRandomGetValueReal(rand, &x[l]);
360:         if (verbose) PetscViewerASCIIPrintf(viewer, "pullback M = %" PetscInt_FMT ":\n", M);
361:         PetscViewerASCIIPushTab(viewer);
362:         CheckPullback(M, N, L, k, w, x, verbose, viewer);
363:         if (M != N) CheckPullback(N, M, L, k, u, v, PETSC_FALSE, viewer);
364:         PetscViewerASCIIPopTab(viewer);
365:         if ((k % N) && (N > 1)) {
366:           if (verbose) PetscViewerASCIIPrintf(viewer, "negative pullback M = %" PetscInt_FMT ":\n", M);
367:           PetscViewerASCIIPushTab(viewer);
368:           CheckPullback(M, N, L, -k, w, x, verbose, viewer);
369:           if (M != N) CheckPullback(N, M, L, -k, u, v, PETSC_FALSE, viewer);
370:           PetscViewerASCIIPopTab(viewer);
371:         }
372:         PetscFree3(L, u, x);
373:       }
374:       if (k > 0) { /* Interior */
375:         PetscInt   Nkm, l, m;
376:         PetscReal *wIntv0, *wIntv0check, wvcheck, diff, diffMat, normMat;
377:         PetscReal *intv0mat, *matcheck;
378:         PetscInt(*indices)[3];

380:         PetscDTBinomialInt(N, k - 1, &Nkm);
381:         PetscMalloc5(Nkm, &wIntv0, Nkm, &wIntv0check, Nk * Nkm, &intv0mat, Nk * Nkm, &matcheck, Nk * k, &indices);
382:         PetscDTAltVInterior(N, k, w, v, wIntv0);
383:         PetscDTAltVInteriorMatrix(N, k, v, intv0mat);
384:         PetscDTAltVInteriorPattern(N, k, indices);
385:         if (verbose) {
386:           PetscViewerASCIIPrintf(viewer, "interior product matrix pattern:\n");
387:           PetscViewerASCIIPushTab(viewer);
388:           for (l = 0; l < Nk * k; l++) {
389:             PetscInt row = indices[l][0];
390:             PetscInt col = indices[l][1];
391:             PetscInt x   = indices[l][2];

393:             PetscViewerASCIIPrintf(viewer, "intV[%" PetscInt_FMT ",%" PetscInt_FMT "] = %sV[%" PetscInt_FMT "]\n", row, col, x < 0 ? "-" : " ", x < 0 ? -(x + 1) : x);
394:           }
395:           PetscViewerASCIIPopTab(viewer);
396:         }
397:         for (l = 0; l < Nkm * Nk; l++) matcheck[l] = 0.;
398:         for (l = 0; l < Nk * k; l++) {
399:           PetscInt row = indices[l][0];
400:           PetscInt col = indices[l][1];
401:           PetscInt x   = indices[l][2];

403:           if (x < 0) {
404:             matcheck[row * Nk + col] = -v[-(x + 1)];
405:           } else {
406:             matcheck[row * Nk + col] = v[x];
407:           }
408:         }
409:         diffMat = 0.;
410:         normMat = 0.;
411:         for (l = 0; l < Nkm * Nk; l++) {
412:           diffMat += PetscSqr(PetscAbsReal(matcheck[l] - intv0mat[l]));
413:           normMat += PetscSqr(matcheck[l]) + PetscSqr(intv0mat[l]);
414:         }
415:         diffMat = PetscSqrtReal(diffMat);
416:         normMat = PetscSqrtReal(normMat);
418:         diffMat = 0.;
419:         normMat = 0.;
420:         for (l = 0; l < Nkm; l++) {
421:           PetscReal sum = 0.;

423:           for (m = 0; m < Nk; m++) sum += intv0mat[l * Nk + m] * w[m];
424:           wIntv0check[l] = sum;

426:           diffMat += PetscSqr(PetscAbsReal(wIntv0check[l] - wIntv0[l]));
427:           normMat += PetscSqr(wIntv0check[l]) + PetscSqr(wIntv0[l]);
428:         }
429:         diffMat = PetscSqrtReal(diffMat);
430:         normMat = PetscSqrtReal(normMat);
432:         if (verbose) {
433:           PetscViewerASCIIPrintf(viewer, "(w int v_0):\n");
434:           PetscViewerASCIIPushTab(viewer);
435:           if (Nkm) PetscRealView(Nkm, wIntv0, viewer);
436:           PetscViewerASCIIPopTab(viewer);

438:           PetscViewerASCIIPrintf(viewer, "(int v_0):\n");
439:           PetscViewerASCIIPushTab(viewer);
440:           if (Nk * Nkm > 0) PetscRealView(Nk * Nkm, intv0mat, viewer);
441:           PetscViewerASCIIPopTab(viewer);
442:         }
443:         PetscDTAltVApply(N, k - 1, wIntv0, &v[N], &wvcheck);
444:         diff = PetscSqrtReal(PetscSqr(wvcheck - wv));
446:         PetscFree5(wIntv0, wIntv0check, intv0mat, matcheck, indices);
447:       }
448:       if (k >= N - k) { /* Hodge star */
449:         PetscReal *u, *starw, *starstarw, wu, starwdotu;
450:         PetscReal  diff, norm;
451:         PetscBool  isOdd;
452:         PetscInt   l;

454:         isOdd = (PetscBool)((k * (N - k)) & 1);
455:         PetscMalloc3(Nk, &u, Nk, &starw, Nk, &starstarw);
456:         PetscDTAltVStar(N, k, 1, w, starw);
457:         PetscDTAltVStar(N, N - k, 1, starw, starstarw);
458:         if (verbose) {
459:           PetscViewerASCIIPrintf(viewer, "star w:\n");
460:           PetscViewerASCIIPushTab(viewer);
461:           if (Nk) PetscRealView(Nk, starw, viewer);
462:           PetscViewerASCIIPopTab(viewer);

464:           PetscViewerASCIIPrintf(viewer, "star star w:\n");
465:           PetscViewerASCIIPushTab(viewer);
466:           if (Nk) PetscRealView(Nk, starstarw, viewer);
467:           PetscViewerASCIIPopTab(viewer);
468:         }
469:         for (l = 0; l < Nk; l++) PetscRandomGetValueReal(rand, &u[l]);
470:         PetscDTAltVWedge(N, k, N - k, w, u, &wu);
471:         starwdotu = 0.;
472:         for (l = 0; l < Nk; l++) starwdotu += starw[l] * u[l];
473:         diff = PetscAbsReal(wu - starwdotu);

476:         diff = 0.;
477:         norm = 0.;
478:         for (l = 0; l < Nk; l++) {
479:           diff += PetscSqr(w[l] - (isOdd ? -starstarw[l] : starstarw[l]));
480:           norm += PetscSqr(w[l]) + PetscSqr(starstarw[l]);
481:         }
482:         diff = PetscSqrtReal(diff);
483:         norm = PetscSqrtReal(norm);
485:         PetscFree3(u, starw, starstarw);
486:       }
487:       PetscFree(v);
488:       PetscFree(w);
489:       PetscViewerASCIIPopTab(viewer);
490:     }
491:     PetscViewerASCIIPopTab(viewer);
492:   }
493:   PetscRandomDestroy(&rand);
494:   PetscFinalize();
495:   return 0;
496: }

498: /*TEST
499:   test:
500:     suffix: 1234
501:     args: -verbose
502:   test:
503:     suffix: 56
504:     args: -N 5,6
505: TEST*/