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