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, 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); 14 15 PetscFunctionBegin; 16 k = PetscAbsInt(k); 17 PetscCall(PetscDTBinomialInt(N, k, &Nk)); 18 PetscCall(PetscDTBinomialInt(M, k, &Mk)); 19 if (negative) { 20 PetscCall(PetscMalloc1(Mk, &walloc)); 21 PetscCall(PetscDTAltVStar(M, M - k, 1, w, walloc)); 22 ww = walloc; 23 } else { 24 ww = w; 25 } 26 PetscCall(PetscMalloc2(Nk, &Lstarw, M * k, &Lx)); 27 PetscCall(PetscMalloc2(Nk * Mk, &Lstar, Nk, &Lstarwcheck)); 28 PetscCall(PetscDTAltVPullback(N, M, L, negative ? -k : k, w, Lstarw)); 29 PetscCall(PetscDTAltVPullbackMatrix(N, M, L, negative ? -k : k, Lstar)); 30 if (negative) { 31 PetscReal *sLsw; 32 33 PetscCall(PetscMalloc1(Nk, &sLsw)); 34 PetscCall(PetscDTAltVStar(N, N - k, 1, Lstarw, sLsw)); 35 PetscCall(PetscDTAltVApply(N, k, sLsw, x, &Lstarwx)); 36 PetscCall(PetscFree(sLsw)); 37 } else { 38 PetscCall(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.; 43 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++) sum += Lstar[i * Mk + j] * w[j]; 53 Lstarwcheck[i] = sum; 54 diffMat += PetscSqr(PetscAbsReal(Lstarwcheck[i] - Lstarw[i])); 55 normMat += PetscSqr(Lstarwcheck[i]) + PetscSqr(Lstarw[i]); 56 } 57 diffMat = PetscSqrtReal(diffMat); 58 normMat = PetscSqrtReal(normMat); 59 if (verbose) { 60 PetscCall(PetscViewerASCIIPrintf(viewer, "L:\n")); 61 PetscCall(PetscViewerASCIIPushTab(viewer)); 62 if (M * N > 0) PetscCall(PetscRealView(M * N, L, viewer)); 63 PetscCall(PetscViewerASCIIPopTab(viewer)); 64 65 PetscCall(PetscViewerASCIIPrintf(viewer, "L*:\n")); 66 PetscCall(PetscViewerASCIIPushTab(viewer)); 67 if (Nk * Mk > 0) PetscCall(PetscRealView(Nk * Mk, Lstar, viewer)); 68 PetscCall(PetscViewerASCIIPopTab(viewer)); 69 70 PetscCall(PetscViewerASCIIPrintf(viewer, "L*w:\n")); 71 PetscCall(PetscViewerASCIIPushTab(viewer)); 72 if (Nk > 0) PetscCall(PetscRealView(Nk, Lstarw, viewer)); 73 PetscCall(PetscViewerASCIIPopTab(viewer)); 74 } 75 PetscCall(PetscDTAltVApply(M, k, ww, Lx, &wLx)); 76 diff = PetscAbsReal(wLx - Lstarwx); 77 PetscCheck(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)", (double)wLx, (double)Lstarwx); 78 PetscCheck(diffMat <= PETSC_SMALL * normMat, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "pullback check: pullback matrix does match matrix-free result"); 79 PetscCall(PetscFree2(Lstar, Lstarwcheck)); 80 PetscCall(PetscFree2(Lstarw, Lx)); 81 PetscCall(PetscFree(walloc)); 82 PetscFunctionReturn(PETSC_SUCCESS); 83 } 84 85 int main(int argc, char **argv) 86 { 87 PetscInt i, numTests = 5, n[5] = {0, 1, 2, 3, 4}; 88 PetscBool verbose = PETSC_FALSE; 89 PetscRandom rand; 90 PetscViewer viewer; 91 92 PetscFunctionBeginUser; 93 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 94 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for exterior algebra tests", "none"); 95 PetscCall(PetscOptionsIntArray("-N", "Up to 5 vector space dimensions to test", "ex7.c", n, &numTests, NULL)); 96 PetscCall(PetscOptionsBool("-verbose", "Verbose test output", "ex7.c", verbose, &verbose, NULL)); 97 PetscOptionsEnd(); 98 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand)); 99 PetscCall(PetscRandomSetInterval(rand, -1., 1.)); 100 PetscCall(PetscRandomSetFromOptions(rand)); 101 if (!numTests) numTests = 5; 102 viewer = PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD); 103 for (i = 0; i < numTests; i++) { 104 PetscInt k, N = n[i]; 105 106 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "N = %" PetscInt_FMT ":\n", N)); 107 PetscCall(PetscViewerASCIIPushTab(viewer)); 108 109 if (verbose) { 110 PetscInt *perm; 111 PetscInt fac = 1; 112 113 PetscCall(PetscMalloc1(N, &perm)); 114 115 for (k = 1; k <= N; k++) fac *= k; 116 PetscCall(PetscViewerASCIIPrintf(viewer, "Permutations of %" PetscInt_FMT ":\n", N)); 117 PetscCall(PetscViewerASCIIPushTab(viewer)); 118 for (k = 0; k < fac; k++) { 119 PetscBool isOdd, isOddCheck; 120 PetscInt j, kCheck; 121 122 PetscCall(PetscDTEnumPerm(N, k, perm, &isOdd)); 123 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ":", k)); 124 for (j = 0; j < N; j++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, perm[j])); 125 PetscCall(PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even")); 126 PetscCall(PetscDTPermIndex(N, perm, &kCheck, &isOddCheck)); 127 PetscCheck(kCheck == k && isOddCheck == isOdd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PetscDTEnumPerm / PetscDTPermIndex mismatch for (%" PetscInt_FMT ", %" PetscInt_FMT ")", N, k); 128 } 129 PetscCall(PetscViewerASCIIPopTab(viewer)); 130 PetscCall(PetscFree(perm)); 131 } 132 for (k = 0; k <= N; k++) { 133 PetscInt j, Nk, M; 134 PetscReal *w, *v, wv; 135 PetscInt *subset; 136 137 PetscCall(PetscDTBinomialInt(N, k, &Nk)); 138 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "k = %" PetscInt_FMT ":\n", k)); 139 PetscCall(PetscViewerASCIIPushTab(viewer)); 140 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT " choose %" PetscInt_FMT "): %" PetscInt_FMT "\n", N, k, Nk)); 141 142 /* Test subset and complement enumeration */ 143 PetscCall(PetscMalloc1(N, &subset)); 144 PetscCall(PetscViewerASCIIPushTab(viewer)); 145 for (j = 0; j < Nk; j++) { 146 PetscBool isOdd, isOddCheck; 147 PetscInt jCheck, kCheck; 148 149 PetscCall(PetscDTEnumSplit(N, k, j, subset, &isOdd)); 150 PetscCall(PetscDTPermIndex(N, subset, &kCheck, &isOddCheck)); 151 PetscCheck(isOddCheck == isOdd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PetscDTEnumSplit sign does not mmatch PetscDTPermIndex sign"); 152 if (verbose) { 153 PetscInt l; 154 155 PetscCall(PetscViewerASCIIPrintf(viewer, "subset %" PetscInt_FMT ":", j)); 156 for (l = 0; l < k; l++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, subset[l])); 157 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " |")); 158 for (l = k; l < N; l++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, subset[l])); 159 PetscCall(PetscPrintf(PETSC_COMM_WORLD, ", %s\n", isOdd ? "odd" : "even")); 160 } 161 PetscCall(PetscDTSubsetIndex(N, k, subset, &jCheck)); 162 PetscCheck(jCheck == j, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "jCheck (%" PetscInt_FMT ") != j (%" PetscInt_FMT ")", jCheck, j); 163 } 164 PetscCall(PetscViewerASCIIPopTab(viewer)); 165 PetscCall(PetscFree(subset)); 166 167 /* Make a random k form */ 168 PetscCall(PetscMalloc1(Nk, &w)); 169 for (j = 0; j < Nk; j++) PetscCall(PetscRandomGetValueReal(rand, &w[j])); 170 /* Make a set of random vectors */ 171 PetscCall(PetscMalloc1(N * k, &v)); 172 for (j = 0; j < N * k; j++) PetscCall(PetscRandomGetValueReal(rand, &v[j])); 173 174 PetscCall(PetscDTAltVApply(N, k, w, v, &wv)); 175 176 if (verbose) { 177 PetscCall(PetscViewerASCIIPrintf(viewer, "w:\n")); 178 PetscCall(PetscViewerASCIIPushTab(viewer)); 179 if (Nk) PetscCall(PetscRealView(Nk, w, viewer)); 180 PetscCall(PetscViewerASCIIPopTab(viewer)); 181 PetscCall(PetscViewerASCIIPrintf(viewer, "v:\n")); 182 PetscCall(PetscViewerASCIIPushTab(viewer)); 183 if (N * k > 0) PetscCall(PetscRealView(N * k, v, viewer)); 184 PetscCall(PetscViewerASCIIPopTab(viewer)); 185 PetscCall(PetscViewerASCIIPrintf(viewer, "w(v): %g\n", (double)wv)); 186 } 187 188 /* sanity checks */ 189 if (k == 1) { /* 1-forms are functionals (dot products) */ 190 PetscInt l; 191 PetscReal wvcheck = 0.; 192 PetscReal diff; 193 194 for (l = 0; l < N; l++) wvcheck += w[l] * v[l]; 195 diff = PetscSqrtReal(PetscSqr(wvcheck - wv)); 196 PetscCheck(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); 197 } 198 if (k == N && N < 5) { /* n-forms are scaled determinants */ 199 PetscReal det, wvcheck, diff; 200 201 switch (k) { 202 case 0: 203 det = 1.; 204 break; 205 case 1: 206 det = v[0]; 207 break; 208 case 2: 209 det = v[0] * v[3] - v[1] * v[2]; 210 break; 211 case 3: 212 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]); 213 break; 214 case 4: 215 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])); 216 break; 217 default: 218 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "invalid k"); 219 } 220 wvcheck = det * w[0]; 221 diff = PetscSqrtReal(PetscSqr(wvcheck - wv)); 222 PetscCheck(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); 223 } 224 if (k > 0) { /* k-forms are linear in each component */ 225 PetscReal alpha; 226 PetscReal *x, *axv, wx, waxv, waxvcheck; 227 PetscReal diff; 228 PetscReal rj; 229 PetscInt l; 230 231 PetscCall(PetscMalloc2(N * k, &x, N * k, &axv)); 232 PetscCall(PetscRandomGetValueReal(rand, &alpha)); 233 PetscCall(PetscRandomSetInterval(rand, 0, k)); 234 PetscCall(PetscRandomGetValueReal(rand, &rj)); 235 j = (PetscInt)rj; 236 PetscCall(PetscRandomSetInterval(rand, -1., 1.)); 237 for (l = 0; l < N * k; l++) x[l] = v[l]; 238 for (l = 0; l < N * k; l++) axv[l] = v[l]; 239 for (l = 0; l < N; l++) { 240 PetscReal val; 241 242 PetscCall(PetscRandomGetValueReal(rand, &val)); 243 x[j * N + l] = val; 244 axv[j * N + l] += alpha * val; 245 } 246 PetscCall(PetscDTAltVApply(N, k, w, x, &wx)); 247 PetscCall(PetscDTAltVApply(N, k, w, axv, &waxv)); 248 waxvcheck = alpha * wx + wv; 249 diff = waxv - waxvcheck; 250 PetscCheck(PetscAbsReal(diff) <= 10. * PETSC_SMALL * (PetscAbsReal(waxv) + PetscAbsReal(waxvcheck)), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "linearity check: component %" PetscInt_FMT ", waxvcheck (%g) != waxv (%g)", j, (double)waxvcheck, (double)waxv); 251 PetscCall(PetscFree2(x, axv)); 252 } 253 if (k > 1) { /* k-forms are antisymmetric */ 254 PetscReal rj, rl, *swapv, wswapv, diff; 255 PetscInt l, m; 256 257 PetscCall(PetscRandomSetInterval(rand, 0, k)); 258 PetscCall(PetscRandomGetValueReal(rand, &rj)); 259 j = (PetscInt)rj; 260 l = j; 261 while (l == j) { 262 PetscCall(PetscRandomGetValueReal(rand, &rl)); 263 l = (PetscInt)rl; 264 } 265 PetscCall(PetscRandomSetInterval(rand, -1., 1.)); 266 PetscCall(PetscMalloc1(N * k, &swapv)); 267 for (m = 0; m < N * k; m++) swapv[m] = v[m]; 268 for (m = 0; m < N; m++) { 269 swapv[j * N + m] = v[l * N + m]; 270 swapv[l * N + m] = v[j * N + m]; 271 } 272 PetscCall(PetscDTAltVApply(N, k, w, swapv, &wswapv)); 273 diff = PetscAbsReal(wswapv + wv); 274 PetscCheck(diff <= PETSC_SMALL * (PetscAbsReal(wswapv) + PetscAbsReal(wv)), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "antisymmetry check: components %" PetscInt_FMT " & %" PetscInt_FMT ", wswapv (%g) != -wv (%g)", j, l, (double)wswapv, (double)wv); 275 PetscCall(PetscFree(swapv)); 276 } 277 for (j = 0; j <= k && j + k <= N; j++) { /* wedge product */ 278 PetscInt Nj, Njk, l, JKj; 279 PetscReal *u, *uWw, *uWwcheck, *uWwmat, *x, *xsplit, uWwx, uWwxcheck, diff, norm; 280 PetscInt *split; 281 282 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "wedge j = %" PetscInt_FMT ":\n", j)); 283 PetscCall(PetscViewerASCIIPushTab(viewer)); 284 PetscCall(PetscDTBinomialInt(N, j, &Nj)); 285 PetscCall(PetscDTBinomialInt(N, j + k, &Njk)); 286 PetscCall(PetscMalloc4(Nj, &u, Njk, &uWw, N * (j + k), &x, N * (j + k), &xsplit)); 287 PetscCall(PetscMalloc1(j + k, &split)); 288 for (l = 0; l < Nj; l++) PetscCall(PetscRandomGetValueReal(rand, &u[l])); 289 for (l = 0; l < N * (j + k); l++) PetscCall(PetscRandomGetValueReal(rand, &x[l])); 290 PetscCall(PetscDTAltVWedge(N, j, k, u, w, uWw)); 291 PetscCall(PetscDTAltVApply(N, j + k, uWw, x, &uWwx)); 292 if (verbose) { 293 PetscCall(PetscViewerASCIIPrintf(viewer, "u:\n")); 294 PetscCall(PetscViewerASCIIPushTab(viewer)); 295 if (Nj) PetscCall(PetscRealView(Nj, u, viewer)); 296 PetscCall(PetscViewerASCIIPopTab(viewer)); 297 PetscCall(PetscViewerASCIIPrintf(viewer, "u wedge w:\n")); 298 PetscCall(PetscViewerASCIIPushTab(viewer)); 299 if (Njk) PetscCall(PetscRealView(Njk, uWw, viewer)); 300 PetscCall(PetscViewerASCIIPopTab(viewer)); 301 PetscCall(PetscViewerASCIIPrintf(viewer, "x:\n")); 302 PetscCall(PetscViewerASCIIPushTab(viewer)); 303 if (N * (j + k) > 0) PetscCall(PetscRealView(N * (j + k), x, viewer)); 304 PetscCall(PetscViewerASCIIPopTab(viewer)); 305 PetscCall(PetscViewerASCIIPrintf(viewer, "u wedge w(x): %g\n", (double)uWwx)); 306 } 307 /* verify wedge formula */ 308 uWwxcheck = 0.; 309 PetscCall(PetscDTBinomialInt(j + k, j, &JKj)); 310 for (l = 0; l < JKj; l++) { 311 PetscBool isOdd; 312 PetscReal ux, wx; 313 PetscInt m, p; 314 315 PetscCall(PetscDTEnumSplit(j + k, j, l, split, &isOdd)); 316 for (m = 0; m < j + k; m++) { 317 for (p = 0; p < N; p++) xsplit[m * N + p] = x[split[m] * N + p]; 318 } 319 PetscCall(PetscDTAltVApply(N, j, u, xsplit, &ux)); 320 PetscCall(PetscDTAltVApply(N, k, w, PetscSafePointerPlusOffset(xsplit, j * N), &wx)); 321 uWwxcheck += isOdd ? -(ux * wx) : (ux * wx); 322 } 323 diff = PetscAbsReal(uWwx - uWwxcheck); 324 PetscCheck(diff <= 10. * PETSC_SMALL * (PetscAbsReal(uWwx) + PetscAbsReal(uWwxcheck)), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "wedge check: forms %" PetscInt_FMT " & %" PetscInt_FMT ", uWwxcheck (%g) != uWwx (%g)", j, k, (double)uWwxcheck, (double)uWwx); 325 PetscCall(PetscFree(split)); 326 PetscCall(PetscMalloc2(Nk * Njk, &uWwmat, Njk, &uWwcheck)); 327 PetscCall(PetscDTAltVWedgeMatrix(N, j, k, u, uWwmat)); 328 if (verbose) { 329 PetscCall(PetscViewerASCIIPrintf(viewer, "(u wedge):\n")); 330 PetscCall(PetscViewerASCIIPushTab(viewer)); 331 if ((Nk * Njk) > 0) PetscCall(PetscRealView(Nk * Njk, uWwmat, viewer)); 332 PetscCall(PetscViewerASCIIPopTab(viewer)); 333 } 334 diff = 0.; 335 norm = 0.; 336 for (l = 0; l < Njk; l++) { 337 PetscInt m; 338 PetscReal sum = 0.; 339 340 for (m = 0; m < Nk; m++) sum += uWwmat[l * Nk + m] * w[m]; 341 uWwcheck[l] = sum; 342 diff += PetscSqr(uWwcheck[l] - uWw[l]); 343 norm += PetscSqr(uWwcheck[l]) + PetscSqr(uWw[l]); 344 } 345 diff = PetscSqrtReal(diff); 346 norm = PetscSqrtReal(norm); 347 PetscCheck(diff <= PETSC_SMALL * norm, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "wedge matrix check: wedge matrix application does not match wedge direct application"); 348 PetscCall(PetscFree2(uWwmat, uWwcheck)); 349 PetscCall(PetscFree4(u, uWw, x, xsplit)); 350 PetscCall(PetscViewerASCIIPopTab(viewer)); 351 } 352 for (M = PetscMax(1, k); M <= N; M++) { /* pullback */ 353 PetscReal *L, *u, *x; 354 PetscInt Mk, l; 355 356 PetscCall(PetscDTBinomialInt(M, k, &Mk)); 357 PetscCall(PetscMalloc3(M * N, &L, Mk, &u, M * k, &x)); 358 for (l = 0; l < M * N; l++) PetscCall(PetscRandomGetValueReal(rand, &L[l])); 359 for (l = 0; l < Mk; l++) PetscCall(PetscRandomGetValueReal(rand, &u[l])); 360 for (l = 0; l < M * k; l++) PetscCall(PetscRandomGetValueReal(rand, &x[l])); 361 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "pullback M = %" PetscInt_FMT ":\n", M)); 362 PetscCall(PetscViewerASCIIPushTab(viewer)); 363 PetscCall(CheckPullback(M, N, L, k, w, x, verbose, viewer)); 364 if (M != N) PetscCall(CheckPullback(N, M, L, k, u, v, PETSC_FALSE, viewer)); 365 PetscCall(PetscViewerASCIIPopTab(viewer)); 366 if ((k % N) && (N > 1)) { 367 if (verbose) PetscCall(PetscViewerASCIIPrintf(viewer, "negative pullback M = %" PetscInt_FMT ":\n", M)); 368 PetscCall(PetscViewerASCIIPushTab(viewer)); 369 PetscCall(CheckPullback(M, N, L, -k, w, x, verbose, viewer)); 370 if (M != N) PetscCall(CheckPullback(N, M, L, -k, u, v, PETSC_FALSE, viewer)); 371 PetscCall(PetscViewerASCIIPopTab(viewer)); 372 } 373 PetscCall(PetscFree3(L, u, x)); 374 } 375 if (k > 0) { /* Interior */ 376 PetscInt Nkm, l, m; 377 PetscReal *wIntv0, *wIntv0check, wvcheck, diff, diffMat, normMat; 378 PetscReal *intv0mat, *matcheck; 379 PetscInt(*indices)[3]; 380 381 PetscCall(PetscDTBinomialInt(N, k - 1, &Nkm)); 382 PetscCall(PetscMalloc5(Nkm, &wIntv0, Nkm, &wIntv0check, Nk * Nkm, &intv0mat, Nk * Nkm, &matcheck, Nk * k, &indices)); 383 PetscCall(PetscDTAltVInterior(N, k, w, v, wIntv0)); 384 PetscCall(PetscDTAltVInteriorMatrix(N, k, v, intv0mat)); 385 PetscCall(PetscDTAltVInteriorPattern(N, k, indices)); 386 if (verbose) { 387 PetscCall(PetscViewerASCIIPrintf(viewer, "interior product matrix pattern:\n")); 388 PetscCall(PetscViewerASCIIPushTab(viewer)); 389 for (l = 0; l < Nk * k; l++) { 390 PetscInt row = indices[l][0]; 391 PetscInt col = indices[l][1]; 392 PetscInt x = indices[l][2]; 393 394 PetscCall(PetscViewerASCIIPrintf(viewer, "intV[%" PetscInt_FMT ",%" PetscInt_FMT "] = %sV[%" PetscInt_FMT "]\n", row, col, x < 0 ? "-" : " ", x < 0 ? -(x + 1) : x)); 395 } 396 PetscCall(PetscViewerASCIIPopTab(viewer)); 397 } 398 for (l = 0; l < Nkm * Nk; l++) matcheck[l] = 0.; 399 for (l = 0; l < Nk * k; l++) { 400 PetscInt row = indices[l][0]; 401 PetscInt col = indices[l][1]; 402 PetscInt x = indices[l][2]; 403 404 if (x < 0) { 405 matcheck[row * Nk + col] = -v[-(x + 1)]; 406 } else { 407 matcheck[row * Nk + col] = v[x]; 408 } 409 } 410 diffMat = 0.; 411 normMat = 0.; 412 for (l = 0; l < Nkm * Nk; l++) { 413 diffMat += PetscSqr(PetscAbsReal(matcheck[l] - intv0mat[l])); 414 normMat += PetscSqr(matcheck[l]) + PetscSqr(intv0mat[l]); 415 } 416 diffMat = PetscSqrtReal(diffMat); 417 normMat = PetscSqrtReal(normMat); 418 PetscCheck(diffMat <= PETSC_SMALL * normMat, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior product check: matrix pattern does not match matrix"); 419 diffMat = 0.; 420 normMat = 0.; 421 for (l = 0; l < Nkm; l++) { 422 PetscReal sum = 0.; 423 424 for (m = 0; m < Nk; m++) sum += intv0mat[l * Nk + m] * w[m]; 425 wIntv0check[l] = sum; 426 427 diffMat += PetscSqr(PetscAbsReal(wIntv0check[l] - wIntv0[l])); 428 normMat += PetscSqr(wIntv0check[l]) + PetscSqr(wIntv0[l]); 429 } 430 diffMat = PetscSqrtReal(diffMat); 431 normMat = PetscSqrtReal(normMat); 432 PetscCheck(diffMat <= PETSC_SMALL * normMat, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Interior product check: application does not match matrix"); 433 if (verbose) { 434 PetscCall(PetscViewerASCIIPrintf(viewer, "(w int v_0):\n")); 435 PetscCall(PetscViewerASCIIPushTab(viewer)); 436 if (Nkm) PetscCall(PetscRealView(Nkm, wIntv0, viewer)); 437 PetscCall(PetscViewerASCIIPopTab(viewer)); 438 439 PetscCall(PetscViewerASCIIPrintf(viewer, "(int v_0):\n")); 440 PetscCall(PetscViewerASCIIPushTab(viewer)); 441 if (Nk * Nkm > 0) PetscCall(PetscRealView(Nk * Nkm, intv0mat, viewer)); 442 PetscCall(PetscViewerASCIIPopTab(viewer)); 443 } 444 PetscCall(PetscDTAltVApply(N, k - 1, wIntv0, &v[N], &wvcheck)); 445 diff = PetscSqrtReal(PetscSqr(wvcheck - wv)); 446 PetscCheck(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); 447 PetscCall(PetscFree5(wIntv0, wIntv0check, intv0mat, matcheck, indices)); 448 } 449 if (k >= N - k) { /* Hodge star */ 450 PetscReal *u, *starw, *starstarw, wu, starwdotu; 451 PetscReal diff, norm; 452 PetscBool isOdd; 453 PetscInt l; 454 455 isOdd = (PetscBool)((k * (N - k)) & 1); 456 PetscCall(PetscMalloc3(Nk, &u, Nk, &starw, Nk, &starstarw)); 457 PetscCall(PetscDTAltVStar(N, k, 1, w, starw)); 458 PetscCall(PetscDTAltVStar(N, N - k, 1, starw, starstarw)); 459 if (verbose) { 460 PetscCall(PetscViewerASCIIPrintf(viewer, "star w:\n")); 461 PetscCall(PetscViewerASCIIPushTab(viewer)); 462 if (Nk) PetscCall(PetscRealView(Nk, starw, viewer)); 463 PetscCall(PetscViewerASCIIPopTab(viewer)); 464 465 PetscCall(PetscViewerASCIIPrintf(viewer, "star star w:\n")); 466 PetscCall(PetscViewerASCIIPushTab(viewer)); 467 if (Nk) PetscCall(PetscRealView(Nk, starstarw, viewer)); 468 PetscCall(PetscViewerASCIIPopTab(viewer)); 469 } 470 for (l = 0; l < Nk; l++) PetscCall(PetscRandomGetValueReal(rand, &u[l])); 471 PetscCall(PetscDTAltVWedge(N, k, N - k, w, u, &wu)); 472 starwdotu = 0.; 473 for (l = 0; l < Nk; l++) starwdotu += starw[l] * u[l]; 474 diff = PetscAbsReal(wu - starwdotu); 475 PetscCheck(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); 476 477 diff = 0.; 478 norm = 0.; 479 for (l = 0; l < Nk; l++) { 480 diff += PetscSqr(w[l] - (isOdd ? -starstarw[l] : starstarw[l])); 481 norm += PetscSqr(w[l]) + PetscSqr(starstarw[l]); 482 } 483 diff = PetscSqrtReal(diff); 484 norm = PetscSqrtReal(norm); 485 PetscCheck(diff <= PETSC_SMALL * norm, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Hodge star check: star(star(w)) != (-1)^(N*(N-k)) w"); 486 PetscCall(PetscFree3(u, starw, starstarw)); 487 } 488 PetscCall(PetscFree(v)); 489 PetscCall(PetscFree(w)); 490 PetscCall(PetscViewerASCIIPopTab(viewer)); 491 } 492 PetscCall(PetscViewerASCIIPopTab(viewer)); 493 } 494 PetscCall(PetscRandomDestroy(&rand)); 495 PetscCall(PetscFinalize()); 496 return 0; 497 } 498 499 /*TEST 500 test: 501 suffix: 1234 502 args: -verbose 503 test: 504 suffix: 56 505 args: -N 5,6 506 TEST*/ 507