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 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; 16 17 PetscFunctionBegin; 18 k = PetscAbsInt(k); 19 ierr = PetscDTBinomialInt(N, k, &Nk);CHKERRQ(ierr); 20 ierr = PetscDTBinomialInt(M, k, &Mk);CHKERRQ(ierr); 21 if (negative) { 22 ierr = PetscMalloc1(Mk, &walloc);CHKERRQ(ierr); 23 ierr = PetscDTAltVStar(M, M - k, 1, w, walloc);CHKERRQ(ierr); 24 ww = walloc; 25 } else { 26 ww = w; 27 } 28 ierr = PetscMalloc2(Nk, &Lstarw, (M*k), &Lx);CHKERRQ(ierr); 29 ierr = PetscMalloc2(Nk * Mk, &Lstar, Nk, &Lstarwcheck);CHKERRQ(ierr); 30 ierr = PetscDTAltVPullback(N, M, L, negative ? -k : k, w, Lstarw);CHKERRQ(ierr); 31 ierr = PetscDTAltVPullbackMatrix(N, M, L, negative ? -k : k, Lstar);CHKERRQ(ierr); 32 if (negative) { 33 PetscReal *sLsw; 34 35 ierr = PetscMalloc1(Nk, &sLsw);CHKERRQ(ierr); 36 ierr = PetscDTAltVStar(N, N - k, 1, Lstarw, sLsw);CHKERRQ(ierr); 37 ierr = PetscDTAltVApply(N, k, sLsw, x, &Lstarwx);CHKERRQ(ierr); 38 ierr = PetscFree(sLsw);CHKERRQ(ierr); 39 } else { 40 ierr = PetscDTAltVApply(N, k, Lstarw, x, &Lstarwx);CHKERRQ(ierr); 41 } 42 for (l = 0; l < k; l++) { 43 for (i = 0; i < M; i++) { 44 PetscReal sum = 0.; 45 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 ierr = PetscViewerASCIIPrintf(viewer, "L:\n");CHKERRQ(ierr); 65 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 66 if (M*N > 0) {ierr = PetscRealView(M*N, L, viewer);CHKERRQ(ierr);} 67 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 68 69 ierr = PetscViewerASCIIPrintf(viewer, "L*:\n");CHKERRQ(ierr); 70 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 71 if (Nk * Mk > 0) {ierr = PetscRealView(Nk * Mk, Lstar, viewer);CHKERRQ(ierr);} 72 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 73 74 ierr = PetscViewerASCIIPrintf(viewer, "L*w:\n");CHKERRQ(ierr); 75 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 76 if (Nk > 0) {ierr = PetscRealView(Nk, Lstarw, viewer);CHKERRQ(ierr);} 77 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 78 } 79 ierr = PetscDTAltVApply(M, k, ww, Lx, &wLx);CHKERRQ(ierr); 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 application: 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 ierr = PetscFree2(Lstar, Lstarwcheck);CHKERRQ(ierr); 84 ierr = PetscFree2(Lstarw, Lx);CHKERRQ(ierr); 85 ierr = PetscFree(walloc);CHKERRQ(ierr); 86 PetscFunctionReturn(0); 87 } 88 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; 95 PetscErrorCode ierr; 96 97 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 98 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Options for exterior algebra tests","none");CHKERRQ(ierr); 99 ierr = PetscOptionsIntArray("-N", "Up to 5 vector space dimensions to test","ex7.c",n,&numTests,NULL);CHKERRQ(ierr); 100 ierr = PetscOptionsBool("-verbose", "Verbose test output","ex7.c",verbose,&verbose,NULL);CHKERRQ(ierr); 101 ierr = PetscOptionsEnd();CHKERRQ(ierr); 102 ierr = PetscRandomCreate(PETSC_COMM_SELF, &rand);CHKERRQ(ierr); 103 ierr = PetscRandomSetInterval(rand, -1., 1.);CHKERRQ(ierr); 104 ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr); 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]; 109 110 if (verbose) {ierr = PetscViewerASCIIPrintf(viewer, "N = %D:\n", N);CHKERRQ(ierr);} 111 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 112 113 if (verbose) { 114 PetscInt *perm; 115 PetscInt fac = 1; 116 117 ierr = PetscMalloc1(N, &perm);CHKERRQ(ierr); 118 119 for (k = 1; k <= N; k++) fac *= k; 120 ierr = PetscViewerASCIIPrintf(viewer, "Permutations of %D:\n", N);CHKERRQ(ierr); 121 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 122 for (k = 0; k < fac; k++) { 123 PetscBool isOdd, isOddCheck; 124 PetscInt j, kCheck; 125 126 ierr = PetscDTEnumPerm(N, k, perm, &isOdd);CHKERRQ(ierr); 127 ierr = PetscViewerASCIIPrintf(viewer, "%D:", k);CHKERRQ(ierr); 128 for (j = 0; j < N; j++) { 129 ierr = PetscPrintf(PETSC_COMM_WORLD, " %D", perm[j]);CHKERRQ(ierr); 130 } 131 ierr = PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even");CHKERRQ(ierr); 132 ierr = PetscDTPermIndex(N, perm, &kCheck, &isOddCheck);CHKERRQ(ierr); 133 if (kCheck != k || isOddCheck != isOdd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PetscDTEnumPerm / PetscDTPermIndex mismatch for (%D, %D)\n", N, k); 134 } 135 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 136 ierr = PetscFree(perm);CHKERRQ(ierr); 137 } 138 for (k = 0; k <= N; k++) { 139 PetscInt j, Nk, M; 140 PetscReal *w, *v, wv; 141 PetscInt *subset; 142 143 ierr = PetscDTBinomialInt(N, k, &Nk);CHKERRQ(ierr); 144 if (verbose) {ierr = PetscViewerASCIIPrintf(viewer, "k = %D:\n", k);CHKERRQ(ierr);} 145 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 146 if (verbose) {ierr = PetscViewerASCIIPrintf(viewer, "(%D choose %D): %D\n", N, k, Nk);CHKERRQ(ierr);} 147 148 /* Test subset and complement enumeration */ 149 ierr = PetscMalloc1(N, &subset);CHKERRQ(ierr); 150 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 151 for (j = 0; j < Nk; j++) { 152 PetscBool isOdd, isOddCheck; 153 PetscInt jCheck, kCheck; 154 155 ierr = PetscDTEnumSplit(N, k, j, subset, &isOdd);CHKERRQ(ierr); 156 ierr = PetscDTPermIndex(N, subset, &kCheck, &isOddCheck);CHKERRQ(ierr); 157 if (isOddCheck != isOdd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PetscDTEnumSplit sign does not mmatch PetscDTPermIndex sign"); 158 if (verbose) { 159 PetscInt l; 160 161 ierr = PetscViewerASCIIPrintf(viewer, "subset %D:", j);CHKERRQ(ierr); 162 for (l = 0; l < k; l++) { 163 ierr = PetscPrintf(PETSC_COMM_WORLD, " %D", subset[l]);CHKERRQ(ierr); 164 } 165 ierr = PetscPrintf(PETSC_COMM_WORLD, " |");CHKERRQ(ierr); 166 for (l = k; l < N; l++) { 167 ierr = PetscPrintf(PETSC_COMM_WORLD, " %D", subset[l]);CHKERRQ(ierr); 168 } 169 ierr = PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even");CHKERRQ(ierr); 170 } 171 ierr = PetscDTSubsetIndex(N, k, subset, &jCheck);CHKERRQ(ierr); 172 if (jCheck != j) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "jCheck (%D) != j (%D)", jCheck, j); 173 } 174 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 175 ierr = PetscFree(subset);CHKERRQ(ierr); 176 177 /* Make a random k form */ 178 ierr = PetscMalloc1(Nk, &w);CHKERRQ(ierr); 179 for (j = 0; j < Nk; j++) {ierr = PetscRandomGetValueReal(rand, &w[j]);CHKERRQ(ierr);} 180 /* Make a set of random vectors */ 181 ierr = PetscMalloc1(N*k, &v);CHKERRQ(ierr); 182 for (j = 0; j < N*k; j++) {ierr = PetscRandomGetValueReal(rand, &v[j]);CHKERRQ(ierr);} 183 184 ierr = PetscDTAltVApply(N, k, w, v, &wv);CHKERRQ(ierr); 185 186 if (verbose) { 187 ierr = PetscViewerASCIIPrintf(viewer, "w:\n");CHKERRQ(ierr); 188 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 189 if (Nk) {ierr = PetscRealView(Nk, w, viewer);CHKERRQ(ierr);} 190 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 191 ierr = PetscViewerASCIIPrintf(viewer, "v:\n");CHKERRQ(ierr); 192 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 193 if (N*k > 0) {ierr = PetscRealView(N*k, v, viewer);CHKERRQ(ierr);} 194 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 195 ierr = PetscViewerASCIIPrintf(viewer, "w(v): %g\n", (double) wv);CHKERRQ(ierr); 196 } 197 198 /* sanity checks */ 199 if (k == 1) { /* 1-forms are functionals (dot products) */ 200 PetscInt l; 201 PetscReal wvcheck = 0.; 202 PetscReal diff; 203 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; 210 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; 253 254 ierr = PetscMalloc2(N * k, &x, N * k, &axv);CHKERRQ(ierr); 255 ierr = PetscRandomGetValueReal(rand, &alpha);CHKERRQ(ierr); 256 ierr = PetscRandomSetInterval(rand, 0, k);CHKERRQ(ierr); 257 ierr = PetscRandomGetValueReal(rand, &rj);CHKERRQ(ierr); 258 j = (PetscInt) rj; 259 ierr = PetscRandomSetInterval(rand, -1., 1.);CHKERRQ(ierr); 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; 264 265 ierr = PetscRandomGetValueReal(rand, &val);CHKERRQ(ierr); 266 x[j * N + l] = val; 267 axv[j * N + l] += alpha * val; 268 } 269 ierr = PetscDTAltVApply(N, k, w, x, &wx);CHKERRQ(ierr); 270 ierr = PetscDTAltVApply(N, k, w, axv, &waxv);CHKERRQ(ierr); 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 ierr = PetscFree2(x,axv);CHKERRQ(ierr); 275 } 276 if (k > 1) { /* k-forms are antisymmetric */ 277 PetscReal rj, rl, *swapv, wswapv, diff; 278 PetscInt l, m; 279 280 ierr = PetscRandomSetInterval(rand, 0, k);CHKERRQ(ierr); 281 ierr = PetscRandomGetValueReal(rand, &rj);CHKERRQ(ierr); 282 j = (PetscInt) rj; 283 l = j; 284 while (l == j) { 285 ierr = PetscRandomGetValueReal(rand, &rl);CHKERRQ(ierr); 286 l = (PetscInt) rl; 287 } 288 ierr = PetscRandomSetInterval(rand, -1., 1.);CHKERRQ(ierr); 289 ierr = PetscMalloc1(N * k, &swapv);CHKERRQ(ierr); 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 ierr = PetscDTAltVApply(N, k, w, swapv, &wswapv);CHKERRQ(ierr); 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 ierr = PetscFree(swapv);CHKERRQ(ierr); 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; 304 305 if (verbose) {ierr = PetscViewerASCIIPrintf(viewer, "wedge j = %D:\n", j);CHKERRQ(ierr);} 306 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 307 ierr = PetscDTBinomialInt(N, j, &Nj);CHKERRQ(ierr); 308 ierr = PetscDTBinomialInt(N, j+k, &Njk);CHKERRQ(ierr); 309 ierr = PetscMalloc4(Nj, &u, Njk, &uWw, N*(j+k), &x, N*(j+k), &xsplit);CHKERRQ(ierr); 310 ierr = PetscMalloc1(j+k,&split);CHKERRQ(ierr); 311 for (l = 0; l < Nj; l++) {ierr = PetscRandomGetValueReal(rand, &u[l]);CHKERRQ(ierr);} 312 for (l = 0; l < N*(j+k); l++) {ierr = PetscRandomGetValueReal(rand, &x[l]);CHKERRQ(ierr);} 313 ierr = PetscDTAltVWedge(N, j, k, u, w, uWw);CHKERRQ(ierr); 314 ierr = PetscDTAltVApply(N, j+k, uWw, x, &uWwx);CHKERRQ(ierr); 315 if (verbose) { 316 ierr = PetscViewerASCIIPrintf(viewer, "u:\n");CHKERRQ(ierr); 317 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 318 if (Nj) {ierr = PetscRealView(Nj, u, viewer);CHKERRQ(ierr);} 319 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 320 ierr = PetscViewerASCIIPrintf(viewer, "u wedge w:\n");CHKERRQ(ierr); 321 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 322 if (Njk) {ierr = PetscRealView(Njk, uWw, viewer);CHKERRQ(ierr);} 323 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 324 ierr = PetscViewerASCIIPrintf(viewer, "x:\n");CHKERRQ(ierr); 325 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 326 if (N*(j+k) > 0) {ierr = PetscRealView(N*(j+k), x, viewer);CHKERRQ(ierr);} 327 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 328 ierr = PetscViewerASCIIPrintf(viewer, "u wedge w(x): %g\n", (double) uWwx);CHKERRQ(ierr); 329 } 330 /* verify wedge formula */ 331 uWwxcheck = 0.; 332 ierr = PetscDTBinomialInt(j+k, j, &JKj);CHKERRQ(ierr); 333 for (l = 0; l < JKj; l++) { 334 PetscBool isOdd; 335 PetscReal ux, wx; 336 PetscInt m, p; 337 338 ierr = PetscDTEnumSplit(j+k, j, l, split, &isOdd);CHKERRQ(ierr); 339 for (m = 0; m < j+k; m++) {for (p = 0; p < N; p++) {xsplit[m * N + p] = x[split[m] * N + p];}} 340 ierr = PetscDTAltVApply(N, j, u, xsplit, &ux);CHKERRQ(ierr); 341 ierr = PetscDTAltVApply(N, k, w, &xsplit[j*N], &wx);CHKERRQ(ierr); 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 ierr = PetscFree(split);CHKERRQ(ierr); 347 ierr = PetscMalloc2(Nk * Njk, &uWwmat, Njk, &uWwcheck);CHKERRQ(ierr); 348 ierr = PetscDTAltVWedgeMatrix(N, j, k, u, uWwmat);CHKERRQ(ierr); 349 if (verbose) { 350 ierr = PetscViewerASCIIPrintf(viewer, "(u wedge):\n");CHKERRQ(ierr); 351 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 352 if ((Nk * Njk) > 0) {ierr = PetscRealView(Nk * Njk, uWwmat, viewer);CHKERRQ(ierr);} 353 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 354 } 355 diff = 0.; 356 norm = 0.; 357 for (l = 0; l < Njk; l++) { 358 PetscInt m; 359 PetscReal sum = 0.; 360 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 application does not match wedge direct application"); 369 ierr = PetscFree2(uWwmat, uWwcheck);CHKERRQ(ierr); 370 ierr = PetscFree4(u, uWw, x, xsplit);CHKERRQ(ierr); 371 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 372 } 373 for (M = PetscMax(1,k); M <= N; M++) { /* pullback */ 374 PetscReal *L, *u, *x; 375 PetscInt Mk, l; 376 377 ierr = PetscDTBinomialInt(M, k, &Mk);CHKERRQ(ierr); 378 ierr = PetscMalloc3(M*N, &L, Mk, &u, M*k, &x);CHKERRQ(ierr); 379 for (l = 0; l < M*N; l++) {ierr = PetscRandomGetValueReal(rand, &L[l]);CHKERRQ(ierr);} 380 for (l = 0; l < Mk; l++) {ierr = PetscRandomGetValueReal(rand, &u[l]);CHKERRQ(ierr);} 381 for (l = 0; l < M*k; l++) {ierr = PetscRandomGetValueReal(rand, &x[l]);CHKERRQ(ierr);} 382 if (verbose) {ierr = PetscViewerASCIIPrintf(viewer, "pullback M = %D:\n", M);CHKERRQ(ierr);} 383 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 384 ierr = CheckPullback(M, N, L, k, w, x, verbose, viewer);CHKERRQ(ierr); 385 if (M != N) {ierr = CheckPullback(N, M, L, k, u, v, PETSC_FALSE, viewer);CHKERRQ(ierr);} 386 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 387 if ((k % N) && (N > 1)) { 388 if (verbose) {ierr = PetscViewerASCIIPrintf(viewer, "negative pullback M = %D:\n", M);CHKERRQ(ierr);} 389 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 390 ierr = CheckPullback(M, N, L, -k, w, x, verbose, viewer);CHKERRQ(ierr); 391 if (M != N) {ierr = CheckPullback(N, M, L, -k, u, v, PETSC_FALSE, viewer);CHKERRQ(ierr);} 392 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 393 } 394 ierr = PetscFree3(L, u, x);CHKERRQ(ierr); 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]; 401 402 ierr = PetscDTBinomialInt(N, k-1, &Nkm);CHKERRQ(ierr); 403 ierr = PetscMalloc5(Nkm, &wIntv0, Nkm, &wIntv0check, Nk * Nkm, &intv0mat, Nk * Nkm, &matcheck, Nk * k, &indices);CHKERRQ(ierr); 404 ierr = PetscDTAltVInterior(N, k, w, v, wIntv0);CHKERRQ(ierr); 405 ierr = PetscDTAltVInteriorMatrix(N, k, v, intv0mat);CHKERRQ(ierr); 406 ierr = PetscDTAltVInteriorPattern(N, k, indices);CHKERRQ(ierr); 407 if (verbose) { 408 ierr = PetscViewerASCIIPrintf(viewer, "interior product matrix pattern:\n");CHKERRQ(ierr); 409 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 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]; 414 415 ierr = PetscViewerASCIIPrintf(viewer,"intV[%D,%D] = %sV[%D]\n", row, col, x < 0 ? "-" : " ", x < 0 ? -(x + 1) : x);CHKERRQ(ierr); 416 } 417 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 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]; 424 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.; 444 445 for (m = 0; m < Nk; m++) sum += intv0mat[l * Nk + m] * w[m]; 446 wIntv0check[l] = sum; 447 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: application does not match matrix"); 454 if (verbose) { 455 ierr = PetscViewerASCIIPrintf(viewer, "(w int v_0):\n");CHKERRQ(ierr); 456 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 457 if (Nkm) {ierr = PetscRealView(Nkm, wIntv0, viewer);CHKERRQ(ierr);} 458 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 459 460 ierr = PetscViewerASCIIPrintf(viewer, "(int v_0):\n");CHKERRQ(ierr); 461 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 462 if (Nk * Nkm > 0) {ierr = PetscRealView(Nk * Nkm, intv0mat, viewer);CHKERRQ(ierr);} 463 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 464 } 465 ierr = PetscDTAltVApply(N, k - 1, wIntv0, &v[N], &wvcheck);CHKERRQ(ierr); 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 ierr = PetscFree5(wIntv0,wIntv0check,intv0mat,matcheck,indices);CHKERRQ(ierr); 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; 475 476 isOdd = (PetscBool) ((k * (N - k)) & 1); 477 ierr = PetscMalloc3(Nk, &u, Nk, &starw, Nk, &starstarw);CHKERRQ(ierr); 478 ierr = PetscDTAltVStar(N, k, 1, w, starw);CHKERRQ(ierr); 479 ierr = PetscDTAltVStar(N, N-k, 1, starw, starstarw);CHKERRQ(ierr); 480 if (verbose) { 481 ierr = PetscViewerASCIIPrintf(viewer, "star w:\n");CHKERRQ(ierr); 482 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 483 if (Nk) {ierr = PetscRealView(Nk, starw, viewer);CHKERRQ(ierr);} 484 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 485 486 ierr = PetscViewerASCIIPrintf(viewer, "star star w:\n");CHKERRQ(ierr); 487 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 488 if (Nk) {ierr = PetscRealView(Nk, starstarw, viewer);CHKERRQ(ierr);} 489 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 490 } 491 for (l = 0; l < Nk; l++) {ierr = PetscRandomGetValueReal(rand,&u[l]);CHKERRQ(ierr);} 492 ierr = PetscDTAltVWedge(N, k, N - k, w, u, &wu);CHKERRQ(ierr); 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); 497 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 ierr = PetscFree3(u, starw, starstarw);CHKERRQ(ierr); 508 } 509 ierr = PetscFree(v);CHKERRQ(ierr); 510 ierr = PetscFree(w);CHKERRQ(ierr); 511 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 512 } 513 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 514 } 515 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr); 516 ierr = PetscFinalize(); 517 return ierr; 518 } 519 520 /*TEST 521 test: 522 suffix: 1234 523 args: -verbose 524 test: 525 suffix: 56 526 args: -N 5,6 527 TEST*/ 528