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