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*/