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