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