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 = Nbpt;
87: PetscBLASInt lwork = 5 * Nbpt;
88: PetscBLASInt lierr;
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 = Nbpt;
122: PetscBLASInt n = Nbpt;
123: PetscBLASInt k = Nbp * Nf;
124: PetscScalar mone = -1.;
125: PetscScalar one = 1.;
127: PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &m, &n, &k, &mone, M_moments, &m, M_moments, &m, &one, Mcopy, &m));
128: }
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);
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);
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_DEFAULT, &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));
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;
165: PetscCall(constructTabulationAndMass(dim, deg, form < 0 ? form - 1 : form + 1, 0, npoints, points, weights, &Nbpt1, &Nf1, &Nk1, &p_trimmed1, &M_trimmed1));
167: PetscCall(PetscMalloc1(Nf1 * (PetscAbsInt(form) + 1), &pattern));
168: PetscCall(PetscDTAltVInteriorPattern(dim, PetscAbsInt(form) + 1, pattern));
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;
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;
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: }
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.;
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];
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: }
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.;
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];
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: }
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_DEFAULT, &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: }
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: }
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;
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);
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);
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: }
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;
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: }
310: /*TEST
312: test:
313: requires: !single
314: args:
316: TEST*/