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*/