1 #include <petsc/private/petscimpl.h> 2 #include <petsc/private/dtimpl.h> 3 4 PetscErrorCode PetscDTAltVApply(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wv) 5 { 6 PetscErrorCode ierr; 7 8 PetscFunctionBegin; 9 if (N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension"); 10 if (k < 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 11 if (N <= 3) { 12 if (!k) { 13 *wv = w[0]; 14 } else { 15 if (N == 1) {*wv = w[0] * v[0];} 16 else if (N == 2) { 17 if (k == 1) {*wv = w[0] * v[0] + w[1] * v[1];} 18 else {*wv = w[0] * (v[0] * v[3] - v[1] * v[2]);} 19 } else { 20 if (k == 1) {*wv = w[0] * v[0] + w[1] * v[1] + w[2] * v[2];} 21 else if (k == 2) { 22 *wv = w[0] * (v[0] * v[4] - v[1] * v[3]) + 23 w[1] * (v[0] * v[5] - v[2] * v[3]) + 24 w[2] * (v[1] * v[5] - v[2] * v[4]); 25 } else { 26 *wv = w[0] * (v[0] * (v[4] * v[8] - v[5] * v[7]) + 27 v[1] * (v[5] * v[6] - v[3] * v[8]) + 28 v[2] * (v[3] * v[7] - v[4] * v[6])); 29 } 30 } 31 } 32 } else { 33 PetscInt Nk, Nf; 34 PetscInt *subset, *work, *perm; 35 PetscInt i, j, l; 36 PetscReal sum = 0.; 37 38 ierr = PetscDTFactorialInt_Internal(k, &Nf);CHKERRQ(ierr); 39 ierr = PetscDTBinomial(N, k, &Nk);CHKERRQ(ierr); 40 ierr = PetscMalloc3(k, &subset, k, &work, k, &perm);CHKERRQ(ierr); 41 for (i = 0; i < Nk; i++) { 42 PetscReal subsum = 0.; 43 44 ierr = PetscDTEnumSubset(N, k, i, subset);CHKERRQ(ierr); 45 for (j = 0; j < Nf; j++) { 46 PetscBool permOdd; 47 PetscReal prod; 48 49 ierr = PetscDTEnumPerm(k, j, work, perm, &permOdd);CHKERRQ(ierr); 50 prod = permOdd ? -1. : 1.; 51 for (l = 0; l < k; l++) { 52 prod *= v[perm[l] * N + subset[l]]; 53 } 54 subsum += prod; 55 } 56 sum += w[i] * subsum; 57 } 58 ierr = PetscFree3(subset, work, perm);CHKERRQ(ierr); 59 *wv = sum; 60 } 61 PetscFunctionReturn(0); 62 } 63 64 PetscErrorCode PetscDTAltVWedge(PetscInt N, PetscInt j, PetscInt k, const PetscReal *a, const PetscReal *b, PetscReal *awedgeb) 65 { 66 PetscInt i; 67 PetscErrorCode ierr; 68 69 PetscFunctionBegin; 70 if (N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension"); 71 if (j < 0 || k < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "negative form degree"); 72 if (j + k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Wedge greater than dimension"); 73 if (N <= 3) { 74 PetscInt Njk; 75 76 ierr = PetscDTBinomial(N, j+k, &Njk);CHKERRQ(ierr); 77 if (!j) {for (i = 0; i < Njk; i++) {awedgeb[i] = a[0] * b[i];}} 78 else if (!k) {for (i = 0; i < Njk; i++) {awedgeb[i] = a[i] * b[0];}} 79 else { 80 if (N == 2) {awedgeb[0] = a[0] * b[1] - a[1] * b[0];} 81 else { 82 if (j+k == 2) { 83 awedgeb[0] = a[0] * b[1] - a[1] * b[0]; 84 awedgeb[1] = a[0] * b[2] - a[2] * b[0]; 85 awedgeb[2] = a[1] * b[2] - a[2] * b[1]; 86 } else { 87 awedgeb[0] = a[0] * b[2] - a[1] * b[1] + a[2] * b[0]; 88 } 89 } 90 } 91 } else { 92 PetscInt Njk; 93 PetscInt JKj; 94 PetscInt *subset, *subsetjk, *subsetj, *subsetk; 95 PetscInt i; 96 97 ierr = PetscDTBinomial(N, j+k, &Njk);CHKERRQ(ierr); 98 ierr = PetscDTBinomial(j+k, j, &JKj);CHKERRQ(ierr); 99 ierr = PetscMalloc4(j+k, &subset, j+k, &subsetjk, j, &subsetj, k, &subsetk);CHKERRQ(ierr); 100 for (i = 0; i < Njk; i++) { 101 PetscReal sum = 0.; 102 PetscInt l; 103 104 ierr = PetscDTEnumSubset(N, j+k, i, subset);CHKERRQ(ierr); 105 for (l = 0; l < JKj; l++) { 106 PetscBool jkOdd; 107 PetscInt m, jInd, kInd; 108 109 ierr = PetscDTEnumSplit(j+k, j, l, subsetjk, &jkOdd);CHKERRQ(ierr); 110 for (m = 0; m < j; m++) { 111 subsetj[m] = subset[subsetjk[m]]; 112 } 113 for (m = 0; m < k; m++) { 114 subsetk[m] = subset[subsetjk[j+m]]; 115 } 116 ierr = PetscDTSubsetIndex(N, j, subsetj, &jInd);CHKERRQ(ierr); 117 ierr = PetscDTSubsetIndex(N, k, subsetk, &kInd);CHKERRQ(ierr); 118 sum += jkOdd ? -(a[jInd] * b[kInd]) : (a[jInd] * b[kInd]); 119 } 120 awedgeb[i] = sum; 121 } 122 ierr = PetscFree4(subset, subsetjk, subsetj, subsetk);CHKERRQ(ierr); 123 } 124 PetscFunctionReturn(0); 125 } 126 127 PetscErrorCode PetscDTAltVWedgeMatrix(PetscInt N, PetscInt j, PetscInt k, const PetscReal *a, PetscReal *awedgeMat) 128 { 129 PetscInt i; 130 PetscErrorCode ierr; 131 132 PetscFunctionBegin; 133 if (N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimension"); 134 if (j < 0 || k < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "negative form degree"); 135 if (j + k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Wedge greater than dimension"); 136 if (N <= 3) { 137 PetscInt Njk; 138 139 ierr = PetscDTBinomial(N, j+k, &Njk);CHKERRQ(ierr); 140 if (!j) { 141 for (i = 0; i < Njk * Njk; i++) {awedgeMat[i] = 0.;} 142 for (i = 0; i < Njk; i++) {awedgeMat[i * (Njk + 1)] = a[0];} 143 } else if (!k) { 144 for (i = 0; i < Njk; i++) {awedgeMat[i] = a[i];} 145 } else { 146 if (N == 2) { 147 awedgeMat[0] = -a[1]; awedgeMat[1] = a[0]; 148 } else { 149 if (j+k == 2) { 150 awedgeMat[0] = -a[1]; awedgeMat[1] = a[0]; awedgeMat[2] = 0.; 151 awedgeMat[3] = -a[2]; awedgeMat[4] = 0.; awedgeMat[5] = a[0]; 152 awedgeMat[6] = 0.; awedgeMat[7] = -a[2]; awedgeMat[8] = a[1]; 153 } else { 154 awedgeMat[0] = a[2]; awedgeMat[1] = -a[1]; awedgeMat[2] = a[0]; 155 } 156 } 157 } 158 } else { 159 PetscInt Njk; 160 PetscInt Nk; 161 PetscInt JKj, i; 162 PetscInt *subset, *subsetjk, *subsetj, *subsetk; 163 164 ierr = PetscDTBinomial(N, k, &Nk);CHKERRQ(ierr); 165 ierr = PetscDTBinomial(N, j+k, &Njk);CHKERRQ(ierr); 166 ierr = PetscDTBinomial(j+k, j, &JKj);CHKERRQ(ierr); 167 ierr = PetscMalloc4(j+k, &subset, j+k, &subsetjk, j, &subsetj, k, &subsetk);CHKERRQ(ierr); 168 for (i = 0; i < Njk * Nk; i++) awedgeMat[i] = 0.; 169 for (i = 0; i < Njk; i++) { 170 PetscInt l; 171 172 ierr = PetscDTEnumSubset(N, j+k, i, subset);CHKERRQ(ierr); 173 for (l = 0; l < JKj; l++) { 174 PetscBool jkOdd; 175 PetscInt m, jInd, kInd; 176 177 ierr = PetscDTEnumSplit(j+k, j, l, subsetjk, &jkOdd);CHKERRQ(ierr); 178 for (m = 0; m < j; m++) { 179 subsetj[m] = subset[subsetjk[m]]; 180 } 181 for (m = 0; m < k; m++) { 182 subsetk[m] = subset[subsetjk[j+m]]; 183 } 184 ierr = PetscDTSubsetIndex(N, j, subsetj, &jInd);CHKERRQ(ierr); 185 ierr = PetscDTSubsetIndex(N, k, subsetk, &kInd);CHKERRQ(ierr); 186 awedgeMat[i * Nk + kInd] += jkOdd ? - a[jInd] : a[jInd]; 187 } 188 } 189 ierr = PetscFree4(subset, subsetjk, subsetj, subsetk);CHKERRQ(ierr); 190 } 191 PetscFunctionReturn(0); 192 } 193 194 /* L: V -> W [|W| by |V| array], L*: altW -> altV */ 195 PetscErrorCode PetscDTAltVPullback(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, const PetscReal *w, PetscReal *Lstarw) 196 { 197 PetscInt i, j, Nk, Mk; 198 PetscErrorCode ierr; 199 200 PetscFunctionBegin; 201 if (N < 0 || M < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimensions"); 202 if (PetscAbsInt(k) > N || PetscAbsInt(k) > M) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 203 if (N <= 3 && M <= 3) { 204 205 ierr = PetscDTBinomial(M, PetscAbsInt(k), &Mk);CHKERRQ(ierr); 206 ierr = PetscDTBinomial(N, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 207 if (!k) { 208 Lstarw[0] = w[0]; 209 } else if (k == 1) { 210 for (i = 0; i < Nk; i++) { 211 PetscReal sum = 0.; 212 213 for (j = 0; j < Mk; j++) {sum += L[j * Nk + i] * w[j];} 214 Lstarw[i] = sum; 215 } 216 } else if (k == -1) { 217 PetscReal mult[3] = {1., -1., 1.}; 218 219 for (i = 0; i < Nk; i++) { 220 PetscReal sum = 0.; 221 222 for (j = 0; j < Mk; j++) { 223 sum += L[(Mk - 1 - j) * Nk + (Nk - 1 - i)] * w[j] * mult[j]; 224 } 225 Lstarw[i] = mult[i] * sum; 226 } 227 } else if (k == 2) { 228 PetscInt pairs[3][2] = {{0,1},{0,2},{1,2}}; 229 230 for (i = 0; i < Nk; i++) { 231 PetscReal sum = 0.; 232 for (j = 0; j < Mk; j++) { 233 sum += (L[pairs[j][0] * N + pairs[i][0]] * L[pairs[j][1] * N + pairs[i][1]] - 234 L[pairs[j][1] * N + pairs[i][0]] * L[pairs[j][0] * N + pairs[i][1]]) * w[j]; 235 } 236 Lstarw[i] = sum; 237 } 238 } else if (k == -2) { 239 PetscInt pairs[3][2] = {{1,2},{2,0},{0,1}}; 240 PetscInt offi = (N == 2) ? 2 : 0; 241 PetscInt offj = (M == 2) ? 2 : 0; 242 243 for (i = 0; i < Nk; i++) { 244 PetscReal sum = 0.; 245 246 for (j = 0; j < Mk; j++) { 247 sum += (L[pairs[offj + j][0] * N + pairs[offi + i][0]] * 248 L[pairs[offj + j][1] * N + pairs[offi + i][1]] - 249 L[pairs[offj + j][1] * N + pairs[offi + i][0]] * 250 L[pairs[offj + j][0] * N + pairs[offi + i][1]]) * w[j]; 251 252 } 253 Lstarw[i] = sum; 254 } 255 } else { 256 PetscReal detL = L[0] * (L[4] * L[8] - L[5] * L[7]) + 257 L[1] * (L[5] * L[6] - L[3] * L[8]) + 258 L[2] * (L[3] * L[7] - L[4] * L[6]); 259 260 for (i = 0; i < Nk; i++) {Lstarw[i] = detL * w[i];} 261 } 262 } else { 263 PetscInt Nf, l, p; 264 PetscReal *Lw, *Lwv; 265 PetscInt *subsetw, *subsetv; 266 PetscInt *work, *perm; 267 PetscReal *walloc = NULL; 268 const PetscReal *ww = NULL; 269 PetscBool negative = PETSC_FALSE; 270 271 ierr = PetscDTBinomial(M, PetscAbsInt(k), &Mk);CHKERRQ(ierr); 272 ierr = PetscDTBinomial(N, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 273 ierr = PetscDTFactorialInt_Internal(PetscAbsInt(k), &Nf);CHKERRQ(ierr); 274 if (k < 0) { 275 negative = PETSC_TRUE; 276 k = -k; 277 ierr = PetscMalloc1(Mk, &walloc);CHKERRQ(ierr); 278 ierr = PetscDTAltVStar(M, M - k, 1, w, walloc);CHKERRQ(ierr); 279 ww = walloc; 280 } else { 281 ww = w; 282 } 283 ierr = PetscMalloc6(k, &subsetw, k, &subsetv, k, &work, k, &perm, N * k, &Lw, k * k, &Lwv);CHKERRQ(ierr); 284 for (i = 0; i < Nk; i++) Lstarw[i] = 0.; 285 for (i = 0; i < Mk; i++) { 286 ierr = PetscDTEnumSubset(M, k, i, subsetw);CHKERRQ(ierr); 287 for (j = 0; j < Nk; j++) { 288 ierr = PetscDTEnumSubset(N, k, j, subsetv);CHKERRQ(ierr); 289 for (p = 0; p < Nf; p++) { 290 PetscReal prod; 291 PetscBool isOdd; 292 293 ierr = PetscDTEnumPerm(k, p, work, perm, &isOdd);CHKERRQ(ierr); 294 prod = isOdd ? -ww[i] : ww[i]; 295 for (l = 0; l < k; l++) { 296 prod *= L[subsetw[perm[l]] * N + subsetv[l]]; 297 } 298 Lstarw[j] += prod; 299 } 300 } 301 } 302 if (negative) { 303 PetscReal *sLsw; 304 305 ierr = PetscMalloc1(Nk, &sLsw);CHKERRQ(ierr); 306 ierr = PetscDTAltVStar(N, N - k, -1, Lstarw, sLsw);CHKERRQ(ierr); 307 for (i = 0; i < Nk; i++) Lstarw[i] = sLsw[i]; 308 ierr = PetscFree(sLsw);CHKERRQ(ierr); 309 } 310 ierr = PetscFree6(subsetw, subsetv, work, perm, Lw, Lwv);CHKERRQ(ierr); 311 ierr = PetscFree(walloc);CHKERRQ(ierr); 312 } 313 PetscFunctionReturn(0); 314 } 315 316 PetscErrorCode PetscDTAltVPullbackMatrix(PetscInt N, PetscInt M, const PetscReal *L, PetscInt k, PetscReal *Lstar) 317 { 318 PetscInt Nk, Mk, Nf, i, j, l, p; 319 PetscReal *Lw, *Lwv; 320 PetscInt *subsetw, *subsetv; 321 PetscInt *work, *perm; 322 PetscBool negative = PETSC_FALSE; 323 PetscErrorCode ierr; 324 325 PetscFunctionBegin; 326 if (N < 0 || M < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid dimensions"); 327 if (PetscAbsInt(k) > N || PetscAbsInt(k) > M) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 328 if (N <= 3 && M <= 3) { 329 PetscReal mult[3] = {1., -1., 1.}; 330 331 ierr = PetscDTBinomial(M, PetscAbsInt(k), &Mk);CHKERRQ(ierr); 332 ierr = PetscDTBinomial(N, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 333 if (!k) { 334 Lstar[0] = 1.; 335 } else if (k == 1) { 336 for (i = 0; i < Nk; i++) {for (j = 0; j < Mk; j++) {Lstar[i * Mk + j] = L[j * Nk + i];}} 337 } else if (k == -1) { 338 for (i = 0; i < Nk; i++) { 339 for (j = 0; j < Mk; j++) { 340 Lstar[i * Mk + j] = L[(Mk - 1 - j) * Nk + (Nk - 1 - i)] * mult[i] * mult[j]; 341 } 342 } 343 } else if (k == 2) { 344 PetscInt pairs[3][2] = {{0,1},{0,2},{1,2}}; 345 346 for (i = 0; i < Nk; i++) { 347 for (j = 0; j < Mk; j++) { 348 Lstar[i * Mk + j] = L[pairs[j][0] * N + pairs[i][0]] * 349 L[pairs[j][1] * N + pairs[i][1]] - 350 L[pairs[j][1] * N + pairs[i][0]] * 351 L[pairs[j][0] * N + pairs[i][1]]; 352 } 353 } 354 } else if (k == -2) { 355 PetscInt pairs[3][2] = {{1,2},{2,0},{0,1}}; 356 PetscInt offi = (N == 2) ? 2 : 0; 357 PetscInt offj = (M == 2) ? 2 : 0; 358 359 for (i = 0; i < Nk; i++) { 360 for (j = 0; j < Mk; j++) { 361 Lstar[i * Mk + j] = L[pairs[offj + j][0] * N + pairs[offi + i][0]] * 362 L[pairs[offj + j][1] * N + pairs[offi + i][1]] - 363 L[pairs[offj + j][1] * N + pairs[offi + i][0]] * 364 L[pairs[offj + j][0] * N + pairs[offi + i][1]]; 365 } 366 } 367 } else { 368 PetscReal detL = L[0] * (L[4] * L[8] - L[5] * L[7]) + 369 L[1] * (L[5] * L[6] - L[3] * L[8]) + 370 L[2] * (L[3] * L[7] - L[4] * L[6]); 371 372 for (i = 0; i < Nk; i++) {Lstar[i] = detL;} 373 } 374 } else { 375 if (k < 0) { 376 negative = PETSC_TRUE; 377 k = -k; 378 } 379 ierr = PetscDTBinomial(M, PetscAbsInt(k), &Mk);CHKERRQ(ierr); 380 ierr = PetscDTBinomial(N, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 381 ierr = PetscDTFactorialInt_Internal(PetscAbsInt(k), &Nf);CHKERRQ(ierr); 382 ierr = PetscMalloc6(M, &subsetw, N, &subsetv, k, &work, k, &perm, N * k, &Lw, k * k, &Lwv);CHKERRQ(ierr); 383 for (i = 0; i < Nk * Mk; i++) Lstar[i] = 0.; 384 for (i = 0; i < Mk; i++) { 385 PetscBool iOdd; 386 PetscInt iidx, jidx; 387 388 ierr = PetscDTEnumSplit(M, k, i, subsetw, &iOdd);CHKERRQ(ierr); 389 iidx = negative ? Mk - 1 - i : i; 390 iOdd = negative ? iOdd ^ ((k * (M-k)) & 1) : PETSC_FALSE; 391 for (j = 0; j < Nk; j++) { 392 PetscBool jOdd; 393 394 ierr = PetscDTEnumSplit(N, k, j, subsetv, &jOdd);CHKERRQ(ierr); 395 jidx = negative ? Nk - 1 - j : j; 396 jOdd = negative ? iOdd ^ jOdd ^ ((k * (N-k)) & 1) : PETSC_FALSE; 397 for (p = 0; p < Nf; p++) { 398 PetscReal prod; 399 PetscBool isOdd; 400 401 ierr = PetscDTEnumPerm(k, p, work, perm, &isOdd);CHKERRQ(ierr); 402 isOdd ^= jOdd; 403 prod = isOdd ? -1. : 1.; 404 for (l = 0; l < k; l++) { 405 prod *= L[subsetw[perm[l]] * N + subsetv[l]]; 406 } 407 Lstar[jidx * Mk + iidx] += prod; 408 } 409 } 410 } 411 ierr = PetscFree6(subsetw, subsetv, work, perm, Lw, Lwv);CHKERRQ(ierr); 412 } 413 PetscFunctionReturn(0); 414 } 415 416 PetscErrorCode PetscDTAltVInterior(PetscInt N, PetscInt k, const PetscReal *w, const PetscReal *v, PetscReal *wIntv) 417 { 418 PetscInt i, Nk, Nkm; 419 PetscErrorCode ierr; 420 421 PetscFunctionBegin; 422 if (k <= 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 423 ierr = PetscDTBinomial(N, k, &Nk);CHKERRQ(ierr); 424 ierr = PetscDTBinomial(N, k-1, &Nkm);CHKERRQ(ierr); 425 if (N <= 3) { 426 if (k == 1) { 427 PetscReal sum = 0.; 428 429 for (i = 0; i < N; i++) { 430 sum += w[i] * v[i]; 431 } 432 wIntv[0] = sum; 433 } else if (k == N) { 434 PetscReal mult[3] = {1., -1., 1.}; 435 436 for (i = 0; i < N; i++) { 437 wIntv[N - 1 - i] = w[0] * v[i] * mult[i]; 438 } 439 } else { 440 wIntv[0] = - w[0]*v[1] - w[1]*v[2]; 441 wIntv[1] = w[0]*v[0] - w[2]*v[2]; 442 wIntv[2] = w[1]*v[0] + w[2]*v[1]; 443 } 444 } else { 445 PetscInt *subset, *work; 446 447 ierr = PetscMalloc2(k, &subset, k, &work);CHKERRQ(ierr); 448 for (i = 0; i < Nkm; i++) wIntv[i] = 0.; 449 for (i = 0; i < Nk; i++) { 450 PetscInt j, l, m; 451 452 ierr = PetscDTEnumSubset(N, k, i, subset);CHKERRQ(ierr); 453 for (j = 0; j < k; j++) { 454 PetscInt idx; 455 PetscBool flip = (j & 1); 456 457 for (l = 0, m = 0; l < k; l++) { 458 if (l != j) work[m++] = subset[l]; 459 } 460 ierr = PetscDTSubsetIndex(N, k - 1, work, &idx);CHKERRQ(ierr); 461 wIntv[idx] += flip ? -(w[i] * v[subset[j]]) : (w[i] * v[subset[j]]); 462 } 463 } 464 ierr = PetscFree2(subset, work);CHKERRQ(ierr); 465 } 466 PetscFunctionReturn(0); 467 } 468 469 PetscErrorCode PetscDTAltVInteriorMatrix(PetscInt N, PetscInt k, const PetscReal *v, PetscReal *intvMat) 470 { 471 PetscInt i, Nk, Nkm; 472 PetscErrorCode ierr; 473 474 PetscFunctionBegin; 475 if (k <= 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 476 ierr = PetscDTBinomial(N, k, &Nk);CHKERRQ(ierr); 477 ierr = PetscDTBinomial(N, k-1, &Nkm);CHKERRQ(ierr); 478 if (N <= 3) { 479 if (k == 1) { 480 for (i = 0; i < N; i++) intvMat[i] = v[i]; 481 } else if (k == N) { 482 PetscReal mult[3] = {1., -1., 1.}; 483 484 for (i = 0; i < N; i++) intvMat[N - 1 - i] = v[i] * mult[i]; 485 } else { 486 intvMat[0] = -v[1]; intvMat[1] = -v[2]; intvMat[2] = 0.; 487 intvMat[3] = v[0]; intvMat[4] = 0.; intvMat[5] = -v[2]; 488 intvMat[6] = 0.; intvMat[7] = v[0]; intvMat[8] = v[1]; 489 } 490 } else { 491 PetscInt *subset, *work; 492 493 ierr = PetscMalloc2(k, &subset, k, &work);CHKERRQ(ierr); 494 for (i = 0; i < Nk * Nkm; i++) intvMat[i] = 0.; 495 for (i = 0; i < Nk; i++) { 496 PetscInt j, l, m; 497 498 ierr = PetscDTEnumSubset(N, k, i, subset);CHKERRQ(ierr); 499 for (j = 0; j < k; j++) { 500 PetscInt idx; 501 PetscBool flip = (j & 1); 502 503 for (l = 0, m = 0; l < k; l++) { 504 if (l != j) work[m++] = subset[l]; 505 } 506 ierr = PetscDTSubsetIndex(N, k - 1, work, &idx);CHKERRQ(ierr); 507 intvMat[idx * Nk + i] += flip ? -v[subset[j]] : v[subset[j]]; 508 } 509 } 510 ierr = PetscFree2(subset, work);CHKERRQ(ierr); 511 } 512 PetscFunctionReturn(0); 513 } 514 515 PetscErrorCode PetscDTAltVInteriorPattern(PetscInt N, PetscInt k, PetscInt (*indices)[3]) 516 { 517 PetscInt i, Nk, Nkm; 518 PetscErrorCode ierr; 519 520 PetscFunctionBegin; 521 if (k <= 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 522 ierr = PetscDTBinomial(N, k, &Nk);CHKERRQ(ierr); 523 ierr = PetscDTBinomial(N, k-1, &Nkm);CHKERRQ(ierr); 524 if (N <= 3) { 525 if (k == 1) { 526 for (i = 0; i < N; i++) { 527 indices[i][0] = 0; 528 indices[i][1] = i; 529 indices[i][2] = i; 530 } 531 } else if (k == N) { 532 PetscInt val[3] = {0, -2, 2}; 533 534 for (i = 0; i < N; i++) { 535 indices[i][0] = N - 1 - i; 536 indices[i][1] = 0; 537 indices[i][2] = val[i]; 538 } 539 } else { 540 indices[0][0] = 0; indices[0][1] = 0; indices[0][2] = -(1 + 1); 541 indices[1][0] = 0; indices[1][1] = 1; indices[1][2] = -(2 + 1); 542 indices[2][0] = 1; indices[2][1] = 0; indices[2][2] = 0; 543 indices[3][0] = 1; indices[3][1] = 2; indices[3][2] = -(2 + 1); 544 indices[4][0] = 2; indices[4][1] = 1; indices[4][2] = 0; 545 indices[5][0] = 2; indices[5][1] = 2; indices[5][2] = 1; 546 } 547 } else { 548 PetscInt *subset, *work; 549 550 ierr = PetscMalloc2(k, &subset, k, &work);CHKERRQ(ierr); 551 for (i = 0; i < Nk; i++) { 552 PetscInt j, l, m; 553 554 ierr = PetscDTEnumSubset(N, k, i, subset);CHKERRQ(ierr); 555 for (j = 0; j < k; j++) { 556 PetscInt idx; 557 PetscBool flip = (j & 1); 558 559 for (l = 0, m = 0; l < k; l++) { 560 if (l != j) work[m++] = subset[l]; 561 } 562 ierr = PetscDTSubsetIndex(N, k - 1, work, &idx);CHKERRQ(ierr); 563 indices[i * k + j][0] = idx; 564 indices[i * k + j][1] = i; 565 indices[i * k + j][2] = flip ? -(subset[j] + 1) : subset[j]; 566 } 567 } 568 ierr = PetscFree2(subset, work);CHKERRQ(ierr); 569 } 570 PetscFunctionReturn(0); 571 } 572 573 PetscErrorCode PetscDTAltVStar(PetscInt N, PetscInt k, PetscInt pow, const PetscReal *w, PetscReal *starw) 574 { 575 PetscInt Nk, i; 576 PetscErrorCode ierr; 577 578 PetscFunctionBegin; 579 if (k < 0 || k > N) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "invalid form degree"); 580 ierr = PetscDTBinomial(N, k, &Nk);CHKERRQ(ierr); 581 pow = pow % 4; 582 pow = (pow + 4) % 4; /* make non-negative */ 583 /* pow is now 0, 1, 2, 3 */ 584 if (N <= 3) { 585 if (pow & 1) { 586 PetscReal mult[3] = {1., -1., 1.}; 587 588 for (i = 0; i < Nk; i++) starw[Nk - 1 - i] = w[i] * mult[i]; 589 } else { 590 for (i = 0; i < Nk; i++) starw[i] = w[i]; 591 } 592 if (pow > 1 && ((k * (N - k)) & 1)) { 593 for (i = 0; i < Nk; i++) starw[i] = -starw[i]; 594 } 595 } else { 596 PetscInt *subset; 597 598 ierr = PetscMalloc1(N, &subset);CHKERRQ(ierr); 599 if (pow % 2) { 600 PetscInt l = (pow == 1) ? k : N - k; 601 for (i = 0; i < Nk; i++) { 602 PetscBool sOdd; 603 PetscInt j, idx; 604 605 ierr = PetscDTEnumSplit(N, l, i, subset, &sOdd);CHKERRQ(ierr); 606 ierr = PetscDTSubsetIndex(N, l, subset, &idx);CHKERRQ(ierr); 607 ierr = PetscDTSubsetIndex(N, N-l, &subset[l], &j);CHKERRQ(ierr); 608 starw[j] = sOdd ? -w[idx] : w[idx]; 609 } 610 } else { 611 for (i = 0; i < Nk; i++) starw[i] = w[i]; 612 } 613 /* star^2 = -1^(k * (N - k)) */ 614 if (pow > 1 && (k * (N - k)) % 2) { 615 for (i = 0; i < Nk; i++) starw[i] = -starw[i]; 616 } 617 ierr = PetscFree(subset);CHKERRQ(ierr); 618 } 619 PetscFunctionReturn(0); 620 } 621