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, lwork, lierr; 87 88 PetscCall(PetscBLASIntCast(Nbpt, &n)); 89 PetscCall(PetscBLASIntCast(5 * Nbpt, &lwork)); 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, n, k; 122 PetscScalar mone = -1.; 123 PetscScalar one = 1.; 124 125 PetscCall(PetscBLASIntCast(Nbpt, &m)); 126 PetscCall(PetscBLASIntCast(Nbpt, &n)); 127 PetscCall(PetscBLASIntCast(Nbp * Nf, &k)); 128 PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &m, &n, &k, &mone, M_moments, &m, M_moments, &m, &one, Mcopy, &m)); 129 } 130 131 frob_err = 0.; 132 for (PetscInt i = 0; i < Nbpt * Nbpt; i++) frob_err += PetscRealPart(Mcopy[i]) * PetscRealPart(Mcopy[i]); 133 frob_err = PetscSqrtReal(frob_err); 134 135 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); 136 137 // 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 138 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, Nbpt, Nbpt, M_trimmed, &mat_trimmed)); 139 PetscCall(PetscDTBinomialInt(dim + deg - 1, dim, &Nbm1)); 140 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, Nbpt, Nbm1 * Nf, M_moments, &mat_moments_T)); 141 PetscCall(MatDuplicate(mat_moments_T, MAT_DO_NOT_COPY_VALUES, &AinvB)); 142 PetscCall(MatLUFactor(mat_trimmed, NULL, NULL, NULL)); 143 PetscCall(MatMatSolve(mat_trimmed, mat_moments_T, AinvB)); 144 PetscCall(MatTransposeMatMult(mat_moments_T, AinvB, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Mm1)); 145 PetscCall(MatShift(Mm1, -1.)); 146 PetscCall(MatNorm(Mm1, NORM_FROBENIUS, &frob_err)); 147 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); 148 PetscCall(MatDestroy(&Mm1)); 149 PetscCall(MatDestroy(&AinvB)); 150 PetscCall(MatDestroy(&mat_moments_T)); 151 152 // The Koszul differential applied to P trimmed (Lambda k+1) should be contained in P trimmed (Lambda k) 153 if (PetscAbsInt(form) < dim) { 154 PetscInt Nf1, Nbpt1, Nk1; 155 PetscReal *p_trimmed1; 156 PetscScalar *M_trimmed1; 157 PetscInt (*pattern)[3]; 158 PetscReal *p_koszul; 159 PetscScalar *M_koszul; 160 PetscScalar *M_k_moment; 161 Mat mat_koszul; 162 Mat mat_k_moment_T; 163 Mat AinvB; 164 Mat prod; 165 166 PetscCall(constructTabulationAndMass(dim, deg, form < 0 ? form - 1 : form + 1, 0, npoints, points, weights, &Nbpt1, &Nf1, &Nk1, &p_trimmed1, &M_trimmed1)); 167 168 PetscCall(PetscMalloc1(Nf1 * (PetscAbsInt(form) + 1), &pattern)); 169 PetscCall(PetscDTAltVInteriorPattern(dim, PetscAbsInt(form) + 1, pattern)); 170 171 // apply the Koszul operator 172 PetscCall(PetscCalloc1(Nbpt1 * Nf * npoints, &p_koszul)); 173 for (PetscInt b = 0; b < Nbpt1; b++) { 174 for (PetscInt a = 0; a < Nf1 * (PetscAbsInt(form) + 1); a++) { 175 PetscInt i, j, k; 176 PetscReal sign; 177 PetscReal *p_i; 178 const PetscReal *p_j; 179 180 i = pattern[a][0]; 181 if (form < 0) i = Nf - 1 - i; 182 j = pattern[a][1]; 183 if (form < 0) j = Nf1 - 1 - j; 184 k = pattern[a][2] < 0 ? -(pattern[a][2] + 1) : pattern[a][2]; 185 sign = pattern[a][2] < 0 ? -1 : 1; 186 if (form < 0 && (i & 1) ^ (j & 1)) sign = -sign; 187 188 p_i = &p_koszul[(b * Nf + i) * npoints]; 189 p_j = &p_trimmed1[(b * Nf1 + j) * npoints]; 190 for (PetscInt pt = 0; pt < npoints; pt++) p_i[pt] += p_j[pt] * points[pt * dim + k] * sign; 191 } 192 } 193 194 // mass matrix of the result 195 PetscCall(PetscMalloc1(Nbpt1 * Nbpt1, &M_koszul)); 196 for (PetscInt i = 0; i < Nbpt1; i++) { 197 for (PetscInt j = 0; j < Nbpt1; j++) { 198 PetscReal val = 0.; 199 200 for (PetscInt v = 0; v < Nf; v++) { 201 const PetscReal *p_i = &p_koszul[(i * Nf + v) * npoints]; 202 const PetscReal *p_j = &p_koszul[(j * Nf + v) * npoints]; 203 204 for (PetscInt pt = 0; pt < npoints; pt++) val += p_i[pt] * p_j[pt] * weights[pt]; 205 } 206 M_koszul[i * Nbpt1 + j] = val; 207 } 208 } 209 210 // moment matrix between the result and P trimmed 211 PetscCall(PetscMalloc1(Nbpt * Nbpt1, &M_k_moment)); 212 for (PetscInt i = 0; i < Nbpt1; i++) { 213 for (PetscInt j = 0; j < Nbpt; j++) { 214 PetscReal val = 0.; 215 216 for (PetscInt v = 0; v < Nf; v++) { 217 const PetscReal *p_i = &p_koszul[(i * Nf + v) * npoints]; 218 const PetscReal *p_j = &p_trimmed[(j * Nf + v) * Nk * npoints]; 219 220 for (PetscInt pt = 0; pt < npoints; pt++) val += p_i[pt] * p_j[pt] * weights[pt]; 221 } 222 M_k_moment[i * Nbpt + j] = val; 223 } 224 } 225 226 // M_k_moment M_trimmed^{-1} M_k_moment^T == M_koszul 227 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, Nbpt1, Nbpt1, M_koszul, &mat_koszul)); 228 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, Nbpt, Nbpt1, M_k_moment, &mat_k_moment_T)); 229 PetscCall(MatDuplicate(mat_k_moment_T, MAT_DO_NOT_COPY_VALUES, &AinvB)); 230 PetscCall(MatMatSolve(mat_trimmed, mat_k_moment_T, AinvB)); 231 PetscCall(MatTransposeMatMult(mat_k_moment_T, AinvB, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &prod)); 232 PetscCall(MatAXPY(prod, -1., mat_koszul, SAME_NONZERO_PATTERN)); 233 PetscCall(MatNorm(prod, NORM_FROBENIUS, &frob_err)); 234 PetscCheck(frob_err <= PETSC_SMALL, 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 PetscCall(MatDestroy(&prod)); 237 PetscCall(MatDestroy(&AinvB)); 238 PetscCall(MatDestroy(&mat_k_moment_T)); 239 PetscCall(MatDestroy(&mat_koszul)); 240 PetscCall(PetscFree(M_k_moment)); 241 PetscCall(PetscFree(M_koszul)); 242 PetscCall(PetscFree(p_koszul)); 243 PetscCall(PetscFree(pattern)); 244 PetscCall(PetscFree(p_trimmed1)); 245 PetscCall(PetscFree(M_trimmed1)); 246 } 247 248 // M_moments has shape [Nbp][Nf][Nbpt] 249 // p_scalar has shape [Nbp][Nk][npoints] 250 // contracting on [Nbp] should be the same shape as 251 // p_trimmed, which is [Nbpt][Nf][Nk][npoints] 252 PetscCall(PetscCalloc1(Nbpt * Nf * Nk * npoints, &p_trimmed_copy)); 253 PetscCall(PetscMalloc1(Nbp * Nf * Nbpt, &M_moment_real)); 254 for (PetscInt i = 0; i < Nbp * Nf * Nbpt; i++) M_moment_real[i] = PetscRealPart(M_moments[i]); 255 for (PetscInt f = 0; f < Nf; f++) { 256 PetscBLASInt m, n, k, lda, ldb, ldc; 257 PetscReal alpha = 1.0; 258 PetscReal beta = 1.0; 259 260 PetscCall(PetscBLASIntCast(Nk * npoints, &m)); 261 PetscCall(PetscBLASIntCast(Nbpt, &n)); 262 PetscCall(PetscBLASIntCast(Nbp, &k)); 263 PetscCall(PetscBLASIntCast(Nk * npoints, &lda)); 264 PetscCall(PetscBLASIntCast(Nf * Nbpt, &ldb)); 265 PetscCall(PetscBLASIntCast(Nf * Nk * npoints, &ldc)); 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 output_file: output/empty.out 316 317 TEST*/ 318