Actual source code: ex13.c
1: const char help[] = "Tests PetscDTPTrimmedEvalJet()";
3: #include <petscdt.h>
4: #include <petscblaslapack.h>
5: #include <petscmat.h>
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;
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));
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.;
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];
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: }
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;
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));
74: PetscCall(constructTabulationAndMass(dim, deg, form, jetDegree, npoints, points, weights, &Nbpt, &Nf, &Nk, &p_trimmed, &M_trimmed));
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));
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;
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));
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
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];
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: }
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.;
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: }
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);
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);
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));
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;
166: PetscCall(constructTabulationAndMass(dim, deg, form < 0 ? form - 1 : form + 1, 0, npoints, points, weights, &Nbpt1, &Nf1, &Nk1, &p_trimmed1, &M_trimmed1));
168: PetscCall(PetscMalloc1(Nf1 * (PetscAbsInt(form) + 1), &pattern));
169: PetscCall(PetscDTAltVInteriorPattern(dim, PetscAbsInt(form) + 1, pattern));
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;
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;
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: }
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.;
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];
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: }
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.;
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];
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: }
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: if (frob_err > PETSC_SMALL) {
235: 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);
236: }
238: PetscCall(MatDestroy(&prod));
239: PetscCall(MatDestroy(&AinvB));
240: PetscCall(MatDestroy(&mat_k_moment_T));
241: PetscCall(MatDestroy(&mat_koszul));
242: PetscCall(PetscFree(M_k_moment));
243: PetscCall(PetscFree(M_koszul));
244: PetscCall(PetscFree(p_koszul));
245: PetscCall(PetscFree(pattern));
246: PetscCall(PetscFree(p_trimmed1));
247: PetscCall(PetscFree(M_trimmed1));
248: }
250: // M_moments has shape [Nbp][Nf][Nbpt]
251: // p_scalar has shape [Nbp][Nk][npoints]
252: // contracting on [Nbp] should be the same shape as
253: // p_trimmed, which is [Nbpt][Nf][Nk][npoints]
254: PetscCall(PetscCalloc1(Nbpt * Nf * Nk * npoints, &p_trimmed_copy));
255: PetscCall(PetscMalloc1(Nbp * Nf * Nbpt, &M_moment_real));
256: for (PetscInt i = 0; i < Nbp * Nf * Nbpt; i++) M_moment_real[i] = PetscRealPart(M_moments[i]);
257: for (PetscInt f = 0; f < Nf; f++) {
258: PetscBLASInt m, n, k, lda, ldb, ldc;
259: PetscReal alpha = 1.0;
260: PetscReal beta = 1.0;
262: PetscCall(PetscBLASIntCast(Nk * npoints, &m));
263: PetscCall(PetscBLASIntCast(Nbpt, &n));
264: PetscCall(PetscBLASIntCast(Nbp, &k));
265: PetscCall(PetscBLASIntCast(Nk * npoints, &lda));
266: PetscCall(PetscBLASIntCast(Nf * Nbpt, &ldb));
267: PetscCall(PetscBLASIntCast(Nf * Nk * npoints, &ldc));
268: 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));
269: }
270: frob_err = 0.;
271: 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]);
272: frob_err = PetscSqrtReal(frob_err);
274: 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);
276: PetscCall(PetscFree(M_moment_real));
277: PetscCall(PetscFree(p_trimmed_copy));
278: PetscCall(MatDestroy(&mat_trimmed));
279: PetscCall(PetscFree(Mcopy));
280: PetscCall(PetscFree(M_moments));
281: PetscCall(PetscFree(M_trimmed));
282: PetscCall(PetscFree(p_trimmed));
283: PetscCall(PetscFree(p_scalar));
284: PetscCall(PetscQuadratureDestroy(&q));
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: int main(int argc, char **argv)
289: {
290: PetscInt max_dim = 3;
291: PetscInt max_deg = 4;
292: PetscInt k = 3;
293: PetscBool cond = PETSC_FALSE;
295: PetscFunctionBeginUser;
296: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
297: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for PetscDTPTrimmedEvalJet() tests", "none");
298: PetscCall(PetscOptionsInt("-max_dim", "Maximum dimension of the simplex", __FILE__, max_dim, &max_dim, NULL));
299: PetscCall(PetscOptionsInt("-max_degree", "Maximum degree of the trimmed polynomial space", __FILE__, max_deg, &max_deg, NULL));
300: PetscCall(PetscOptionsInt("-max_jet", "The number of derivatives to test", __FILE__, k, &k, NULL));
301: PetscCall(PetscOptionsBool("-cond", "Compute the condition numbers of the mass matrices of the bases", __FILE__, cond, &cond, NULL));
302: PetscOptionsEnd();
303: for (PetscInt dim = 2; dim <= max_dim; dim++) {
304: for (PetscInt deg = 1; deg <= max_deg; deg++) {
305: for (PetscInt form = -dim + 1; form <= dim; form++) PetscCall(test(dim, deg, form, PetscMax(1, k), cond));
306: }
307: }
308: PetscCall(PetscFinalize());
309: return 0;
310: }
312: /*TEST
314: test:
315: requires: !single
316: args:
318: TEST*/