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