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 16 PetscFunctionBegin; 17 k = PetscAbsInt(k); 18 PetscCall(PetscDTBinomialInt(N, k, &Nk)); 19 PetscCall(PetscDTBinomialInt(M, k, &Mk)); 20 if (negative) { 21 PetscCall(PetscMalloc1(Mk, &walloc)); 22 PetscCall(PetscDTAltVStar(M, M - k, 1, w, walloc)); 23 ww = walloc; 24 } else { 25 ww = w; 26 } 27 PetscCall(PetscMalloc2(Nk, &Lstarw, (M*k), &Lx)); 28 PetscCall(PetscMalloc2(Nk * Mk, &Lstar, Nk, &Lstarwcheck)); 29 PetscCall(PetscDTAltVPullback(N, M, L, negative ? -k : k, w, Lstarw)); 30 PetscCall(PetscDTAltVPullbackMatrix(N, M, L, negative ? -k : k, Lstar)); 31 if (negative) { 32 PetscReal *sLsw; 33 34 PetscCall(PetscMalloc1(Nk, &sLsw)); 35 PetscCall(PetscDTAltVStar(N, N - k, 1, Lstarw, sLsw)); 36 PetscCall(PetscDTAltVApply(N, k, sLsw, x, &Lstarwx)); 37 PetscCall(PetscFree(sLsw)); 38 } else { 39 PetscCall(PetscDTAltVApply(N, k, Lstarw, x, &Lstarwx)); 40 } 41 for (l = 0; l < k; l++) { 42 for (i = 0; i < M; i++) { 43 PetscReal sum = 0.; 44 45 for (j = 0; j < N; j++) sum += L[i * N + j] * x[l * N + j]; 46 Lx[l * M + i] = sum; 47 } 48 } 49 diffMat = 0.; 50 normMat = 0.; 51 for (i = 0; i < Nk; i++) { 52 PetscReal sum = 0.; 53 for (j = 0; j < Mk; j++) { 54 sum += Lstar[i * Mk + j] * w[j]; 55 } 56 Lstarwcheck[i] = sum; 57 diffMat += PetscSqr(PetscAbsReal(Lstarwcheck[i] - Lstarw[i])); 58 normMat += PetscSqr(Lstarwcheck[i]) + PetscSqr(Lstarw[i]); 59 } 60 diffMat = PetscSqrtReal(diffMat); 61 normMat = PetscSqrtReal(normMat); 62 if (verbose) { 63 PetscCall(PetscViewerASCIIPrintf(viewer, "L:\n")); 64 PetscCall(PetscViewerASCIIPushTab(viewer)); 65 if (M*N > 0) PetscCall(PetscRealView(M*N, L, viewer)); 66 PetscCall(PetscViewerASCIIPopTab(viewer)); 67 68 PetscCall(PetscViewerASCIIPrintf(viewer, "L*:\n")); 69 PetscCall(PetscViewerASCIIPushTab(viewer)); 70 if (Nk * Mk > 0) PetscCall(PetscRealView(Nk * Mk, Lstar, viewer)); 71 PetscCall(PetscViewerASCIIPopTab(viewer)); 72 73 PetscCall(PetscViewerASCIIPrintf(viewer, "L*w:\n")); 74 PetscCall(PetscViewerASCIIPushTab(viewer)); 75 if (Nk > 0) PetscCall(PetscRealView(Nk, Lstarw, viewer)); 76 PetscCall(PetscViewerASCIIPopTab(viewer)); 77 } 78 PetscCall(PetscDTAltVApply(M, k, ww, Lx, &wLx)); 79 diff = PetscAbsReal(wLx - Lstarwx); 80 PetscCheckFalse(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)", wLx, Lstarwx); 81 PetscCheckFalse(diffMat > PETSC_SMALL * normMat,PETSC_COMM_WORLD, PETSC_ERR_PLIB, "pullback check: pullback matrix does match matrix free result"); 82 PetscCall(PetscFree2(Lstar, Lstarwcheck)); 83 PetscCall(PetscFree2(Lstarw, Lx)); 84 PetscCall(PetscFree(walloc)); 85 PetscFunctionReturn(0); 86 } 87 88 int main(int argc, char **argv) 89 { 90 PetscInt i, numTests = 5, n[5] = {0, 1, 2, 3, 4}; 91 PetscBool verbose = PETSC_FALSE; 92 PetscRandom rand; 93 PetscViewer viewer; 94 95 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 96 PetscOptionsBegin(PETSC_COMM_WORLD,"","Options for exterior algebra tests","none"); 97 PetscCall(PetscOptionsIntArray("-N", "Up to 5 vector space dimensions to test","ex7.c",n,&numTests,NULL)); 98 PetscCall(PetscOptionsBool("-verbose", "Verbose test output","ex7.c",verbose,&verbose,NULL)); 99 PetscOptionsEnd(); 100 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand)); 101 PetscCall(PetscRandomSetInterval(rand, -1., 1.)); 102 PetscCall(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]; 107 108 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "N = %D:\n", N)); 109 PetscCall(PetscViewerASCIIPushTab(viewer)); 110 111 if (verbose) { 112 PetscInt *perm; 113 PetscInt fac = 1; 114 115 PetscCall(PetscMalloc1(N, &perm)); 116 117 for (k = 1; k <= N; k++) fac *= k; 118 PetscCall(PetscViewerASCIIPrintf(viewer, "Permutations of %D:\n", N)); 119 PetscCall(PetscViewerASCIIPushTab(viewer)); 120 for (k = 0; k < fac; k++) { 121 PetscBool isOdd, isOddCheck; 122 PetscInt j, kCheck; 123 124 PetscCall(PetscDTEnumPerm(N, k, perm, &isOdd)); 125 PetscCall(PetscViewerASCIIPrintf(viewer, "%D:", k)); 126 for (j = 0; j < N; j++) { 127 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %D", perm[j])); 128 } 129 PetscCall(PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even")); 130 PetscCall(PetscDTPermIndex(N, perm, &kCheck, &isOddCheck)); 131 PetscCheckFalse(kCheck != k || isOddCheck != isOdd,PETSC_COMM_SELF, PETSC_ERR_PLIB, "PetscDTEnumPerm / PetscDTPermIndex mismatch for (%D, %D)", N, k); 132 } 133 PetscCall(PetscViewerASCIIPopTab(viewer)); 134 PetscCall(PetscFree(perm)); 135 } 136 for (k = 0; k <= N; k++) { 137 PetscInt j, Nk, M; 138 PetscReal *w, *v, wv; 139 PetscInt *subset; 140 141 PetscCall(PetscDTBinomialInt(N, k, &Nk)); 142 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "k = %D:\n", k)); 143 PetscCall(PetscViewerASCIIPushTab(viewer)); 144 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "(%D choose %D): %D\n", N, k, Nk)); 145 146 /* Test subset and complement enumeration */ 147 PetscCall(PetscMalloc1(N, &subset)); 148 PetscCall(PetscViewerASCIIPushTab(viewer)); 149 for (j = 0; j < Nk; j++) { 150 PetscBool isOdd, isOddCheck; 151 PetscInt jCheck, kCheck; 152 153 PetscCall(PetscDTEnumSplit(N, k, j, subset, &isOdd)); 154 PetscCall(PetscDTPermIndex(N, subset, &kCheck, &isOddCheck)); 155 PetscCheck(isOddCheck == isOdd,PETSC_COMM_SELF, PETSC_ERR_PLIB, "PetscDTEnumSplit sign does not mmatch PetscDTPermIndex sign"); 156 if (verbose) { 157 PetscInt l; 158 159 PetscCall(PetscViewerASCIIPrintf(viewer, "subset %D:", j)); 160 for (l = 0; l < k; l++) { 161 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %D", subset[l])); 162 } 163 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " |")); 164 for (l = k; l < N; l++) { 165 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %D", subset[l])); 166 } 167 PetscCall(PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even")); 168 } 169 PetscCall(PetscDTSubsetIndex(N, k, subset, &jCheck)); 170 PetscCheck(jCheck == j,PETSC_COMM_WORLD, PETSC_ERR_PLIB, "jCheck (%D) != j (%D)", jCheck, j); 171 } 172 PetscCall(PetscViewerASCIIPopTab(viewer)); 173 PetscCall(PetscFree(subset)); 174 175 /* Make a random k form */ 176 PetscCall(PetscMalloc1(Nk, &w)); 177 for (j = 0; j < Nk; j++) PetscCall(PetscRandomGetValueReal(rand, &w[j])); 178 /* Make a set of random vectors */ 179 PetscCall(PetscMalloc1(N*k, &v)); 180 for (j = 0; j < N*k; j++) PetscCall(PetscRandomGetValueReal(rand, &v[j])); 181 182 PetscCall(PetscDTAltVApply(N, k, w, v, &wv)); 183 184 if (verbose) { 185 PetscCall(PetscViewerASCIIPrintf(viewer, "w:\n")); 186 PetscCall(PetscViewerASCIIPushTab(viewer)); 187 if (Nk) PetscCall(PetscRealView(Nk, w, viewer)); 188 PetscCall(PetscViewerASCIIPopTab(viewer)); 189 PetscCall(PetscViewerASCIIPrintf(viewer, "v:\n")); 190 PetscCall(PetscViewerASCIIPushTab(viewer)); 191 if (N*k > 0) PetscCall(PetscRealView(N*k, v, viewer)); 192 PetscCall(PetscViewerASCIIPopTab(viewer)); 193 PetscCall(PetscViewerASCIIPrintf(viewer, "w(v): %g\n", (double) wv)); 194 } 195 196 /* sanity checks */ 197 if (k == 1) { /* 1-forms are functionals (dot products) */ 198 PetscInt l; 199 PetscReal wvcheck = 0.; 200 PetscReal diff; 201 202 for (l = 0; l < N; l++) wvcheck += w[l] * v[l]; 203 diff = PetscSqrtReal(PetscSqr(wvcheck - wv)); 204 PetscCheckFalse(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); 205 } 206 if (k == N && N < 5) { /* n-forms are scaled determinants */ 207 PetscReal det, wvcheck, diff; 208 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)); 243 PetscCheckFalse(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); 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; 251 252 PetscCall(PetscMalloc2(N * k, &x, N * k, &axv)); 253 PetscCall(PetscRandomGetValueReal(rand, &alpha)); 254 PetscCall(PetscRandomSetInterval(rand, 0, k)); 255 PetscCall(PetscRandomGetValueReal(rand, &rj)); 256 j = (PetscInt) rj; 257 PetscCall(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; 262 263 PetscCall(PetscRandomGetValueReal(rand, &val)); 264 x[j * N + l] = val; 265 axv[j * N + l] += alpha * val; 266 } 267 PetscCall(PetscDTAltVApply(N, k, w, x, &wx)); 268 PetscCall(PetscDTAltVApply(N, k, w, axv, &waxv)); 269 waxvcheck = alpha * wx + wv; 270 diff = waxv - waxvcheck; 271 PetscCheckFalse(PetscAbsReal(diff) > 10. * PETSC_SMALL * (PetscAbsReal(waxv) + PetscAbsReal(waxvcheck)),PETSC_COMM_WORLD, PETSC_ERR_PLIB, "linearity check: component %D, waxvcheck (%g) != waxv (%g)", j, (double) waxvcheck, (double) waxv); 272 PetscCall(PetscFree2(x,axv)); 273 } 274 if (k > 1) { /* k-forms are antisymmetric */ 275 PetscReal rj, rl, *swapv, wswapv, diff; 276 PetscInt l, m; 277 278 PetscCall(PetscRandomSetInterval(rand, 0, k)); 279 PetscCall(PetscRandomGetValueReal(rand, &rj)); 280 j = (PetscInt) rj; 281 l = j; 282 while (l == j) { 283 PetscCall(PetscRandomGetValueReal(rand, &rl)); 284 l = (PetscInt) rl; 285 } 286 PetscCall(PetscRandomSetInterval(rand, -1., 1.)); 287 PetscCall(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 PetscCall(PetscDTAltVApply(N, k, w, swapv, &wswapv)); 294 diff = PetscAbsReal(wswapv + wv); 295 PetscCheckFalse(diff > PETSC_SMALL * (PetscAbsReal(wswapv) + PetscAbsReal(wv)),PETSC_COMM_WORLD, PETSC_ERR_PLIB, "antisymmetry check: components %D & %D, wswapv (%g) != -wv (%g)", j, l, (double) wswapv, (double) wv); 296 PetscCall(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; 302 303 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "wedge j = %D:\n", j)); 304 PetscCall(PetscViewerASCIIPushTab(viewer)); 305 PetscCall(PetscDTBinomialInt(N, j, &Nj)); 306 PetscCall(PetscDTBinomialInt(N, j+k, &Njk)); 307 PetscCall(PetscMalloc4(Nj, &u, Njk, &uWw, N*(j+k), &x, N*(j+k), &xsplit)); 308 PetscCall(PetscMalloc1(j+k,&split)); 309 for (l = 0; l < Nj; l++) PetscCall(PetscRandomGetValueReal(rand, &u[l])); 310 for (l = 0; l < N*(j+k); l++) PetscCall(PetscRandomGetValueReal(rand, &x[l])); 311 PetscCall(PetscDTAltVWedge(N, j, k, u, w, uWw)); 312 PetscCall(PetscDTAltVApply(N, j+k, uWw, x, &uWwx)); 313 if (verbose) { 314 PetscCall(PetscViewerASCIIPrintf(viewer, "u:\n")); 315 PetscCall(PetscViewerASCIIPushTab(viewer)); 316 if (Nj) PetscCall(PetscRealView(Nj, u, viewer)); 317 PetscCall(PetscViewerASCIIPopTab(viewer)); 318 PetscCall(PetscViewerASCIIPrintf(viewer, "u wedge w:\n")); 319 PetscCall(PetscViewerASCIIPushTab(viewer)); 320 if (Njk) PetscCall(PetscRealView(Njk, uWw, viewer)); 321 PetscCall(PetscViewerASCIIPopTab(viewer)); 322 PetscCall(PetscViewerASCIIPrintf(viewer, "x:\n")); 323 PetscCall(PetscViewerASCIIPushTab(viewer)); 324 if (N*(j+k) > 0) PetscCall(PetscRealView(N*(j+k), x, viewer)); 325 PetscCall(PetscViewerASCIIPopTab(viewer)); 326 PetscCall(PetscViewerASCIIPrintf(viewer, "u wedge w(x): %g\n", (double) uWwx)); 327 } 328 /* verify wedge formula */ 329 uWwxcheck = 0.; 330 PetscCall(PetscDTBinomialInt(j+k, j, &JKj)); 331 for (l = 0; l < JKj; l++) { 332 PetscBool isOdd; 333 PetscReal ux, wx; 334 PetscInt m, p; 335 336 PetscCall(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 PetscCall(PetscDTAltVApply(N, j, u, xsplit, &ux)); 339 PetscCall(PetscDTAltVApply(N, k, w, &xsplit[j*N], &wx)); 340 uWwxcheck += isOdd ? -(ux * wx) : (ux * wx); 341 } 342 diff = PetscAbsReal(uWwx - uWwxcheck); 343 PetscCheckFalse(diff > 10. * PETSC_SMALL * (PetscAbsReal(uWwx) + PetscAbsReal(uWwxcheck)),PETSC_COMM_WORLD, PETSC_ERR_PLIB, "wedge check: forms %D & %D, uWwxcheck (%g) != uWwx (%g)", j, k, (double) uWwxcheck, (double) uWwx); 344 PetscCall(PetscFree(split)); 345 PetscCall(PetscMalloc2(Nk * Njk, &uWwmat, Njk, &uWwcheck)); 346 PetscCall(PetscDTAltVWedgeMatrix(N, j, k, u, uWwmat)); 347 if (verbose) { 348 PetscCall(PetscViewerASCIIPrintf(viewer, "(u wedge):\n")); 349 PetscCall(PetscViewerASCIIPushTab(viewer)); 350 if ((Nk * Njk) > 0) PetscCall(PetscRealView(Nk * Njk, uWwmat, viewer)); 351 PetscCall(PetscViewerASCIIPopTab(viewer)); 352 } 353 diff = 0.; 354 norm = 0.; 355 for (l = 0; l < Njk; l++) { 356 PetscInt m; 357 PetscReal sum = 0.; 358 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); 366 PetscCheckFalse(diff > PETSC_SMALL * norm,PETSC_COMM_WORLD, PETSC_ERR_PLIB, "wedge matrix check: wedge matrix application does not match wedge direct application"); 367 PetscCall(PetscFree2(uWwmat, uWwcheck)); 368 PetscCall(PetscFree4(u, uWw, x, xsplit)); 369 PetscCall(PetscViewerASCIIPopTab(viewer)); 370 } 371 for (M = PetscMax(1,k); M <= N; M++) { /* pullback */ 372 PetscReal *L, *u, *x; 373 PetscInt Mk, l; 374 375 PetscCall(PetscDTBinomialInt(M, k, &Mk)); 376 PetscCall(PetscMalloc3(M*N, &L, Mk, &u, M*k, &x)); 377 for (l = 0; l < M*N; l++) PetscCall(PetscRandomGetValueReal(rand, &L[l])); 378 for (l = 0; l < Mk; l++) PetscCall(PetscRandomGetValueReal(rand, &u[l])); 379 for (l = 0; l < M*k; l++) PetscCall(PetscRandomGetValueReal(rand, &x[l])); 380 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "pullback M = %D:\n", M)); 381 PetscCall(PetscViewerASCIIPushTab(viewer)); 382 PetscCall(CheckPullback(M, N, L, k, w, x, verbose, viewer)); 383 if (M != N) PetscCall(CheckPullback(N, M, L, k, u, v, PETSC_FALSE, viewer)); 384 PetscCall(PetscViewerASCIIPopTab(viewer)); 385 if ((k % N) && (N > 1)) { 386 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "negative pullback M = %D:\n", M)); 387 PetscCall(PetscViewerASCIIPushTab(viewer)); 388 PetscCall(CheckPullback(M, N, L, -k, w, x, verbose, viewer)); 389 if (M != N) PetscCall(CheckPullback(N, M, L, -k, u, v, PETSC_FALSE, viewer)); 390 PetscCall(PetscViewerASCIIPopTab(viewer)); 391 } 392 PetscCall(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]; 399 400 PetscCall(PetscDTBinomialInt(N, k-1, &Nkm)); 401 PetscCall(PetscMalloc5(Nkm, &wIntv0, Nkm, &wIntv0check, Nk * Nkm, &intv0mat, Nk * Nkm, &matcheck, Nk * k, &indices)); 402 PetscCall(PetscDTAltVInterior(N, k, w, v, wIntv0)); 403 PetscCall(PetscDTAltVInteriorMatrix(N, k, v, intv0mat)); 404 PetscCall(PetscDTAltVInteriorPattern(N, k, indices)); 405 if (verbose) { 406 PetscCall(PetscViewerASCIIPrintf(viewer, "interior product matrix pattern:\n")); 407 PetscCall(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]; 412 413 PetscCall(PetscViewerASCIIPrintf(viewer,"intV[%D,%D] = %sV[%D]\n", row, col, x < 0 ? "-" : " ", x < 0 ? -(x + 1) : x)); 414 } 415 PetscCall(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]; 422 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); 437 PetscCheckFalse(diffMat > PETSC_SMALL * normMat,PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior product check: matrix pattern does not match matrix"); 438 diffMat = 0.; 439 normMat = 0.; 440 for (l = 0; l < Nkm; l++) { 441 PetscReal sum = 0.; 442 443 for (m = 0; m < Nk; m++) sum += intv0mat[l * Nk + m] * w[m]; 444 wIntv0check[l] = sum; 445 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); 451 PetscCheckFalse(diffMat > PETSC_SMALL * normMat,PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior product check: application does not match matrix"); 452 if (verbose) { 453 PetscCall(PetscViewerASCIIPrintf(viewer, "(w int v_0):\n")); 454 PetscCall(PetscViewerASCIIPushTab(viewer)); 455 if (Nkm) PetscCall(PetscRealView(Nkm, wIntv0, viewer)); 456 PetscCall(PetscViewerASCIIPopTab(viewer)); 457 458 PetscCall(PetscViewerASCIIPrintf(viewer, "(int v_0):\n")); 459 PetscCall(PetscViewerASCIIPushTab(viewer)); 460 if (Nk * Nkm > 0) PetscCall(PetscRealView(Nk * Nkm, intv0mat, viewer)); 461 PetscCall(PetscViewerASCIIPopTab(viewer)); 462 } 463 PetscCall(PetscDTAltVApply(N, k - 1, wIntv0, &v[N], &wvcheck)); 464 diff = PetscSqrtReal(PetscSqr(wvcheck - wv)); 465 PetscCheckFalse(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); 466 PetscCall(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; 473 474 isOdd = (PetscBool) ((k * (N - k)) & 1); 475 PetscCall(PetscMalloc3(Nk, &u, Nk, &starw, Nk, &starstarw)); 476 PetscCall(PetscDTAltVStar(N, k, 1, w, starw)); 477 PetscCall(PetscDTAltVStar(N, N-k, 1, starw, starstarw)); 478 if (verbose) { 479 PetscCall(PetscViewerASCIIPrintf(viewer, "star w:\n")); 480 PetscCall(PetscViewerASCIIPushTab(viewer)); 481 if (Nk) PetscCall(PetscRealView(Nk, starw, viewer)); 482 PetscCall(PetscViewerASCIIPopTab(viewer)); 483 484 PetscCall(PetscViewerASCIIPrintf(viewer, "star star w:\n")); 485 PetscCall(PetscViewerASCIIPushTab(viewer)); 486 if (Nk) PetscCall(PetscRealView(Nk, starstarw, viewer)); 487 PetscCall(PetscViewerASCIIPopTab(viewer)); 488 } 489 for (l = 0; l < Nk; l++) PetscCall(PetscRandomGetValueReal(rand,&u[l])); 490 PetscCall(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); 494 PetscCheckFalse(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); 495 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); 504 PetscCheckFalse(diff > PETSC_SMALL * norm,PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Hodge star check: star(star(w)) != (-1)^(N*(N-k)) w"); 505 PetscCall(PetscFree3(u, starw, starstarw)); 506 } 507 PetscCall(PetscFree(v)); 508 PetscCall(PetscFree(w)); 509 PetscCall(PetscViewerASCIIPopTab(viewer)); 510 } 511 PetscCall(PetscViewerASCIIPopTab(viewer)); 512 } 513 PetscCall(PetscRandomDestroy(&rand)); 514 PetscCall(PetscFinalize()); 515 return 0; 516 } 517 518 /*TEST 519 test: 520 suffix: 1234 521 args: -verbose 522 test: 523 suffix: 56 524 args: -N 5,6 525 TEST*/ 526