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