1 const char help[] = "Tests PetscDTPTrimmedEvalJet()"; 2 3 #include <petscdt.h> 4 #include <petscblaslapack.h> 5 #include <petscmat.h> 6 7 static PetscErrorCode constructTabulationAndMass(PetscInt dim, PetscInt deg, PetscInt form, PetscInt jetDegree, PetscInt npoints, const PetscReal *points, const PetscReal *weights, PetscInt *_Nb, PetscInt *_Nf, PetscInt *_Nk, PetscReal **B, PetscScalar **M) { 8 PetscInt Nf; // Number of form components 9 PetscInt Nbpt; // number of trimmed polynomials 10 PetscInt Nk; // jet size 11 PetscReal *p_trimmed; 12 13 PetscFunctionBegin; 14 PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(form), &Nf)); 15 PetscCall(PetscDTPTrimmedSize(dim, deg, form, &Nbpt)); 16 PetscCall(PetscDTBinomialInt(dim + jetDegree, dim, &Nk)); 17 PetscCall(PetscMalloc1(Nbpt * Nf * Nk * npoints, &p_trimmed)); 18 PetscCall(PetscDTPTrimmedEvalJet(dim, npoints, points, deg, form, jetDegree, p_trimmed)); 19 20 // compute the direct mass matrix 21 PetscScalar *M_trimmed; 22 PetscCall(PetscCalloc1(Nbpt * Nbpt, &M_trimmed)); 23 for (PetscInt i = 0; i < Nbpt; i++) { 24 for (PetscInt j = 0; j < Nbpt; j++) { 25 PetscReal v = 0.; 26 27 for (PetscInt f = 0; f < Nf; f++) { 28 const PetscReal *p_i = &p_trimmed[(i * Nf + f) * Nk * npoints]; 29 const PetscReal *p_j = &p_trimmed[(j * Nf + f) * Nk * npoints]; 30 31 for (PetscInt pt = 0; pt < npoints; pt++) { v += p_i[pt] * p_j[pt] * weights[pt]; } 32 } 33 M_trimmed[i * Nbpt + j] += v; 34 } 35 } 36 *_Nb = Nbpt; 37 *_Nf = Nf; 38 *_Nk = Nk; 39 *B = p_trimmed; 40 *M = M_trimmed; 41 PetscFunctionReturn(0); 42 } 43 44 static PetscErrorCode test(PetscInt dim, PetscInt deg, PetscInt form, PetscInt jetDegree, PetscBool cond) { 45 PetscQuadrature q; 46 PetscInt npoints; 47 const PetscReal *points; 48 const PetscReal *weights; 49 PetscInt Nf; // Number of form components 50 PetscInt Nk; // jet size 51 PetscInt Nbpt; // number of trimmed polynomials 52 PetscReal *p_trimmed; 53 PetscScalar *M_trimmed; 54 PetscReal *p_scalar; 55 PetscInt Nbp; // number of scalar polynomials 56 PetscScalar *Mcopy; 57 PetscScalar *M_moments; 58 PetscReal frob_err = 0.; 59 Mat mat_trimmed; 60 Mat mat_moments_T; 61 Mat AinvB; 62 PetscInt Nbm1; 63 Mat Mm1; 64 PetscReal *p_trimmed_copy; 65 PetscReal *M_moment_real; 66 67 PetscFunctionBegin; 68 // Construct an appropriate quadrature 69 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, deg + 2, -1., 1., &q)); 70 PetscCall(PetscQuadratureGetData(q, NULL, NULL, &npoints, &points, &weights)); 71 72 PetscCall(constructTabulationAndMass(dim, deg, form, jetDegree, npoints, points, weights, &Nbpt, &Nf, &Nk, &p_trimmed, &M_trimmed)); 73 74 PetscCall(PetscDTBinomialInt(dim + deg, dim, &Nbp)); 75 PetscCall(PetscMalloc1(Nbp * Nk * npoints, &p_scalar)); 76 PetscCall(PetscDTPKDEvalJet(dim, npoints, points, deg, jetDegree, p_scalar)); 77 78 PetscCall(PetscMalloc1(Nbpt * Nbpt, &Mcopy)); 79 // Print the condition numbers (useful for testing out different bases internally in PetscDTPTrimmedEvalJet()) 80 #if !defined(PETSC_USE_COMPLEX) 81 if (cond) { 82 PetscReal *S; 83 PetscScalar *work; 84 PetscBLASInt n = Nbpt; 85 PetscBLASInt lwork = 5 * Nbpt; 86 PetscBLASInt lierr; 87 88 PetscCall(PetscMalloc1(Nbpt, &S)); 89 PetscCall(PetscMalloc1(5 * Nbpt, &work)); 90 PetscCall(PetscArraycpy(Mcopy, M_trimmed, Nbpt * Nbpt)); 91 92 PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("N", "N", &n, &n, Mcopy, &n, S, NULL, &n, NULL, &n, work, &lwork, &lierr)); 93 PetscReal cond = S[0] / S[Nbpt - 1]; 94 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "dimension %" PetscInt_FMT ", degree %" PetscInt_FMT ", form %" PetscInt_FMT ": condition number %g\n", dim, deg, form, (double)cond)); 95 PetscCall(PetscFree(work)); 96 PetscCall(PetscFree(S)); 97 } 98 #endif 99 100 // compute the moments with the orthonormal polynomials 101 PetscCall(PetscCalloc1(Nbpt * Nbp * Nf, &M_moments)); 102 for (PetscInt i = 0; i < Nbp; i++) { 103 for (PetscInt j = 0; j < Nbpt; j++) { 104 for (PetscInt f = 0; f < Nf; f++) { 105 PetscReal v = 0.; 106 const PetscReal *p_i = &p_scalar[i * Nk * npoints]; 107 const PetscReal *p_j = &p_trimmed[(j * Nf + f) * Nk * npoints]; 108 109 for (PetscInt pt = 0; pt < npoints; pt++) { v += p_i[pt] * p_j[pt] * weights[pt]; } 110 M_moments[(i * Nf + f) * Nbpt + j] += v; 111 } 112 } 113 } 114 115 // subtract M_moments^T * M_moments from M_trimmed: because the trimmed polynomials should be contained in 116 // the full polynomials, the result should be zero 117 PetscCall(PetscArraycpy(Mcopy, M_trimmed, Nbpt * Nbpt)); 118 { 119 PetscBLASInt m = Nbpt; 120 PetscBLASInt n = Nbpt; 121 PetscBLASInt k = Nbp * Nf; 122 PetscScalar mone = -1.; 123 PetscScalar one = 1.; 124 125 PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &m, &n, &k, &mone, M_moments, &m, M_moments, &m, &one, Mcopy, &m)); 126 } 127 128 frob_err = 0.; 129 for (PetscInt i = 0; i < Nbpt * Nbpt; i++) frob_err += PetscRealPart(Mcopy[i]) * PetscRealPart(Mcopy[i]); 130 frob_err = PetscSqrtReal(frob_err); 131 132 PetscCheck(frob_err <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "dimension %" PetscInt_FMT ", degree %" PetscInt_FMT ", form %" PetscInt_FMT ": trimmed projection error %g", dim, deg, form, (double)frob_err); 133 134 // P trimmed is also supposed to contain the polynomials of one degree less: construction M_moment[0:sub,:] * M_trimmed^{-1} * M_moments[0:sub,:]^T should be the identity matrix 135 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, Nbpt, Nbpt, M_trimmed, &mat_trimmed)); 136 PetscCall(PetscDTBinomialInt(dim + deg - 1, dim, &Nbm1)); 137 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, Nbpt, Nbm1 * Nf, M_moments, &mat_moments_T)); 138 PetscCall(MatDuplicate(mat_moments_T, MAT_DO_NOT_COPY_VALUES, &AinvB)); 139 PetscCall(MatLUFactor(mat_trimmed, NULL, NULL, NULL)); 140 PetscCall(MatMatSolve(mat_trimmed, mat_moments_T, AinvB)); 141 PetscCall(MatTransposeMatMult(mat_moments_T, AinvB, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Mm1)); 142 PetscCall(MatShift(Mm1, -1.)); 143 PetscCall(MatNorm(Mm1, NORM_FROBENIUS, &frob_err)); 144 PetscCheck(frob_err <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "dimension %" PetscInt_FMT ", degree %" PetscInt_FMT ", form %" PetscInt_FMT ": trimmed reverse projection error %g", dim, deg, form, (double)frob_err); 145 PetscCall(MatDestroy(&Mm1)); 146 PetscCall(MatDestroy(&AinvB)); 147 PetscCall(MatDestroy(&mat_moments_T)); 148 149 // The Koszul differential applied to P trimmed (Lambda k+1) should be contained in P trimmed (Lambda k) 150 if (PetscAbsInt(form) < dim) { 151 PetscInt Nf1, Nbpt1, Nk1; 152 PetscReal *p_trimmed1; 153 PetscScalar *M_trimmed1; 154 PetscInt(*pattern)[3]; 155 PetscReal *p_koszul; 156 PetscScalar *M_koszul; 157 PetscScalar *M_k_moment; 158 Mat mat_koszul; 159 Mat mat_k_moment_T; 160 Mat AinvB; 161 Mat prod; 162 163 PetscCall(constructTabulationAndMass(dim, deg, form < 0 ? form - 1 : form + 1, 0, npoints, points, weights, &Nbpt1, &Nf1, &Nk1, &p_trimmed1, &M_trimmed1)); 164 165 PetscCall(PetscMalloc1(Nf1 * (PetscAbsInt(form) + 1), &pattern)); 166 PetscCall(PetscDTAltVInteriorPattern(dim, PetscAbsInt(form) + 1, pattern)); 167 168 // apply the Koszul operator 169 PetscCall(PetscCalloc1(Nbpt1 * Nf * npoints, &p_koszul)); 170 for (PetscInt b = 0; b < Nbpt1; b++) { 171 for (PetscInt a = 0; a < Nf1 * (PetscAbsInt(form) + 1); a++) { 172 PetscInt i, j, k; 173 PetscReal sign; 174 PetscReal *p_i; 175 const PetscReal *p_j; 176 177 i = pattern[a][0]; 178 if (form < 0) { i = Nf - 1 - i; } 179 j = pattern[a][1]; 180 if (form < 0) { j = Nf1 - 1 - j; } 181 k = pattern[a][2] < 0 ? -(pattern[a][2] + 1) : pattern[a][2]; 182 sign = pattern[a][2] < 0 ? -1 : 1; 183 if (form < 0 && (i & 1) ^ (j & 1)) { sign = -sign; } 184 185 p_i = &p_koszul[(b * Nf + i) * npoints]; 186 p_j = &p_trimmed1[(b * Nf1 + j) * npoints]; 187 for (PetscInt pt = 0; pt < npoints; pt++) { p_i[pt] += p_j[pt] * points[pt * dim + k] * sign; } 188 } 189 } 190 191 // mass matrix of the result 192 PetscCall(PetscMalloc1(Nbpt1 * Nbpt1, &M_koszul)); 193 for (PetscInt i = 0; i < Nbpt1; i++) { 194 for (PetscInt j = 0; j < Nbpt1; j++) { 195 PetscReal val = 0.; 196 197 for (PetscInt v = 0; v < Nf; v++) { 198 const PetscReal *p_i = &p_koszul[(i * Nf + v) * npoints]; 199 const PetscReal *p_j = &p_koszul[(j * Nf + v) * npoints]; 200 201 for (PetscInt pt = 0; pt < npoints; pt++) { val += p_i[pt] * p_j[pt] * weights[pt]; } 202 } 203 M_koszul[i * Nbpt1 + j] = val; 204 } 205 } 206 207 // moment matrix between the result and P trimmed 208 PetscCall(PetscMalloc1(Nbpt * Nbpt1, &M_k_moment)); 209 for (PetscInt i = 0; i < Nbpt1; i++) { 210 for (PetscInt j = 0; j < Nbpt; j++) { 211 PetscReal val = 0.; 212 213 for (PetscInt v = 0; v < Nf; v++) { 214 const PetscReal *p_i = &p_koszul[(i * Nf + v) * npoints]; 215 const PetscReal *p_j = &p_trimmed[(j * Nf + v) * Nk * npoints]; 216 217 for (PetscInt pt = 0; pt < npoints; pt++) { val += p_i[pt] * p_j[pt] * weights[pt]; } 218 } 219 M_k_moment[i * Nbpt + j] = val; 220 } 221 } 222 223 // M_k_moment M_trimmed^{-1} M_k_moment^T == M_koszul 224 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, Nbpt1, Nbpt1, M_koszul, &mat_koszul)); 225 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, Nbpt, Nbpt1, M_k_moment, &mat_k_moment_T)); 226 PetscCall(MatDuplicate(mat_k_moment_T, MAT_DO_NOT_COPY_VALUES, &AinvB)); 227 PetscCall(MatMatSolve(mat_trimmed, mat_k_moment_T, AinvB)); 228 PetscCall(MatTransposeMatMult(mat_k_moment_T, AinvB, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &prod)); 229 PetscCall(MatAXPY(prod, -1., mat_koszul, SAME_NONZERO_PATTERN)); 230 PetscCall(MatNorm(prod, NORM_FROBENIUS, &frob_err)); 231 if (frob_err > PETSC_SMALL) { 232 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "dimension %" PetscInt_FMT ", degree %" PetscInt_FMT ", forms (%" PetscInt_FMT ", %" PetscInt_FMT "): koszul projection error %g", dim, deg, form, form < 0 ? (form - 1) : (form + 1), (double)frob_err); 233 } 234 235 PetscCall(MatDestroy(&prod)); 236 PetscCall(MatDestroy(&AinvB)); 237 PetscCall(MatDestroy(&mat_k_moment_T)); 238 PetscCall(MatDestroy(&mat_koszul)); 239 PetscCall(PetscFree(M_k_moment)); 240 PetscCall(PetscFree(M_koszul)); 241 PetscCall(PetscFree(p_koszul)); 242 PetscCall(PetscFree(pattern)); 243 PetscCall(PetscFree(p_trimmed1)); 244 PetscCall(PetscFree(M_trimmed1)); 245 } 246 247 // M_moments has shape [Nbp][Nf][Nbpt] 248 // p_scalar has shape [Nbp][Nk][npoints] 249 // contracting on [Nbp] should be the same shape as 250 // p_trimmed, which is [Nbpt][Nf][Nk][npoints] 251 PetscCall(PetscCalloc1(Nbpt * Nf * Nk * npoints, &p_trimmed_copy)); 252 PetscCall(PetscMalloc1(Nbp * Nf * Nbpt, &M_moment_real)); 253 for (PetscInt i = 0; i < Nbp * Nf * Nbpt; i++) { M_moment_real[i] = PetscRealPart(M_moments[i]); } 254 for (PetscInt f = 0; f < Nf; f++) { 255 PetscBLASInt m = Nk * npoints; 256 PetscBLASInt n = Nbpt; 257 PetscBLASInt k = Nbp; 258 PetscBLASInt lda = Nk * npoints; 259 PetscBLASInt ldb = Nf * Nbpt; 260 PetscBLASInt ldc = Nf * Nk * npoints; 261 PetscReal alpha = 1.0; 262 PetscReal beta = 1.0; 263 264 PetscCallBLAS("BLASREALgemm", BLASREALgemm_("N", "T", &m, &n, &k, &alpha, p_scalar, &lda, &M_moment_real[f * Nbpt], &ldb, &beta, &p_trimmed_copy[f * Nk * npoints], &ldc)); 265 } 266 frob_err = 0.; 267 for (PetscInt i = 0; i < Nbpt * Nf * Nk * npoints; i++) { frob_err += (p_trimmed_copy[i] - p_trimmed[i]) * (p_trimmed_copy[i] - p_trimmed[i]); } 268 frob_err = PetscSqrtReal(frob_err); 269 270 PetscCheck(frob_err < 10 * PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "dimension %" PetscInt_FMT ", degree %" PetscInt_FMT ", form %" PetscInt_FMT ": jet error %g", dim, deg, form, (double)frob_err); 271 272 PetscCall(PetscFree(M_moment_real)); 273 PetscCall(PetscFree(p_trimmed_copy)); 274 PetscCall(MatDestroy(&mat_trimmed)); 275 PetscCall(PetscFree(Mcopy)); 276 PetscCall(PetscFree(M_moments)); 277 PetscCall(PetscFree(M_trimmed)); 278 PetscCall(PetscFree(p_trimmed)); 279 PetscCall(PetscFree(p_scalar)); 280 PetscCall(PetscQuadratureDestroy(&q)); 281 PetscFunctionReturn(0); 282 } 283 284 int main(int argc, char **argv) { 285 PetscInt max_dim = 3; 286 PetscInt max_deg = 4; 287 PetscInt k = 3; 288 PetscBool cond = PETSC_FALSE; 289 290 PetscFunctionBeginUser; 291 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 292 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for PetscDTPTrimmedEvalJet() tests", "none"); 293 PetscCall(PetscOptionsInt("-max_dim", "Maximum dimension of the simplex", __FILE__, max_dim, &max_dim, NULL)); 294 PetscCall(PetscOptionsInt("-max_degree", "Maximum degree of the trimmed polynomial space", __FILE__, max_deg, &max_deg, NULL)); 295 PetscCall(PetscOptionsInt("-max_jet", "The number of derivatives to test", __FILE__, k, &k, NULL)); 296 PetscCall(PetscOptionsBool("-cond", "Compute the condition numbers of the mass matrices of the bases", __FILE__, cond, &cond, NULL)); 297 PetscOptionsEnd(); 298 for (PetscInt dim = 2; dim <= max_dim; dim++) { 299 for (PetscInt deg = 1; deg <= max_deg; deg++) { 300 for (PetscInt form = -dim + 1; form <= dim; form++) { PetscCall(test(dim, deg, form, PetscMax(1, k), cond)); } 301 } 302 } 303 PetscCall(PetscFinalize()); 304 return 0; 305 } 306 307 /*TEST 308 309 test: 310 requires: !single 311 args: 312 313 TEST*/ 314