Actual source code: dt.c

petsc-master 2017-12-15
Report Typos and Errors
  1: /* Discretization tools */

  3: #include <petscconf.h>
  4: #if defined(PETSC_HAVE_MATHIMF_H)
  5: #include <mathimf.h>           /* this needs to be included before math.h */
  6: #endif
  7: #ifdef PETSC_HAVE_MPFR
  8: #include <mpfr.h>
  9: #endif

 11:  #include <petscdt.h>
 12:  #include <petscblaslapack.h>
 13:  #include <petsc/private/petscimpl.h>
 14:  #include <petsc/private/dtimpl.h>
 15:  #include <petscviewer.h>
 16:  #include <petscdmplex.h>
 17:  #include <petscdmshell.h>

 19: static PetscBool GaussCite       = PETSC_FALSE;
 20: const char       GaussCitation[] = "@article{GolubWelsch1969,\n"
 21:                                    "  author  = {Golub and Welsch},\n"
 22:                                    "  title   = {Calculation of Quadrature Rules},\n"
 23:                                    "  journal = {Math. Comp.},\n"
 24:                                    "  volume  = {23},\n"
 25:                                    "  number  = {106},\n"
 26:                                    "  pages   = {221--230},\n"
 27:                                    "  year    = {1969}\n}\n";

 29: /*@
 30:   PetscQuadratureCreate - Create a PetscQuadrature object

 32:   Collective on MPI_Comm

 34:   Input Parameter:
 35: . comm - The communicator for the PetscQuadrature object

 37:   Output Parameter:
 38: . q  - The PetscQuadrature object

 40:   Level: beginner

 42: .keywords: PetscQuadrature, quadrature, create
 43: .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
 44: @*/
 45: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
 46: {

 51:   PetscSysInitializePackage();
 52:   PetscHeaderCreate(*q,PETSC_OBJECT_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
 53:   (*q)->dim       = -1;
 54:   (*q)->Nc        =  1;
 55:   (*q)->order     = -1;
 56:   (*q)->numPoints = 0;
 57:   (*q)->points    = NULL;
 58:   (*q)->weights   = NULL;
 59:   return(0);
 60: }

 62: /*@
 63:   PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object

 65:   Collective on PetscQuadrature

 67:   Input Parameter:
 68: . q  - The PetscQuadrature object

 70:   Output Parameter:
 71: . r  - The new PetscQuadrature object

 73:   Level: beginner

 75: .keywords: PetscQuadrature, quadrature, clone
 76: .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy(), PetscQuadratureGetData()
 77: @*/
 78: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
 79: {
 80:   PetscInt         order, dim, Nc, Nq;
 81:   const PetscReal *points, *weights;
 82:   PetscReal       *p, *w;
 83:   PetscErrorCode   ierr;

 87:   PetscQuadratureCreate(PetscObjectComm((PetscObject) q), r);
 88:   PetscQuadratureGetOrder(q, &order);
 89:   PetscQuadratureSetOrder(*r, order);
 90:   PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights);
 91:   PetscMalloc1(Nq*dim, &p);
 92:   PetscMalloc1(Nq*Nc, &w);
 93:   PetscMemcpy(p, points, Nq*dim * sizeof(PetscReal));
 94:   PetscMemcpy(w, weights, Nc * Nq * sizeof(PetscReal));
 95:   PetscQuadratureSetData(*r, dim, Nc, Nq, p, w);
 96:   return(0);
 97: }

 99: /*@
100:   PetscQuadratureDestroy - Destroys a PetscQuadrature object

102:   Collective on PetscQuadrature

104:   Input Parameter:
105: . q  - The PetscQuadrature object

107:   Level: beginner

109: .keywords: PetscQuadrature, quadrature, destroy
110: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
111: @*/
112: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
113: {

117:   if (!*q) return(0);
119:   if (--((PetscObject)(*q))->refct > 0) {
120:     *q = NULL;
121:     return(0);
122:   }
123:   PetscFree((*q)->points);
124:   PetscFree((*q)->weights);
125:   PetscHeaderDestroy(q);
126:   return(0);
127: }

129: /*@
130:   PetscQuadratureGetOrder - Return the order of the method

132:   Not collective

134:   Input Parameter:
135: . q - The PetscQuadrature object

137:   Output Parameter:
138: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated

140:   Level: intermediate

142: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
143: @*/
144: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
145: {
149:   *order = q->order;
150:   return(0);
151: }

153: /*@
154:   PetscQuadratureSetOrder - Return the order of the method

156:   Not collective

158:   Input Parameters:
159: + q - The PetscQuadrature object
160: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated

162:   Level: intermediate

164: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
165: @*/
166: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
167: {
170:   q->order = order;
171:   return(0);
172: }

174: /*@
175:   PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated

177:   Not collective

179:   Input Parameter:
180: . q - The PetscQuadrature object

182:   Output Parameter:
183: . Nc - The number of components

185:   Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.

187:   Level: intermediate

189: .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
190: @*/
191: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
192: {
196:   *Nc = q->Nc;
197:   return(0);
198: }

200: /*@
201:   PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated

203:   Not collective

205:   Input Parameters:
206: + q  - The PetscQuadrature object
207: - Nc - The number of components

209:   Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.

211:   Level: intermediate

213: .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
214: @*/
215: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
216: {
219:   q->Nc = Nc;
220:   return(0);
221: }

223: /*@C
224:   PetscQuadratureGetData - Returns the data defining the quadrature

226:   Not collective

228:   Input Parameter:
229: . q  - The PetscQuadrature object

231:   Output Parameters:
232: + dim - The spatial dimension
233: , Nc - The number of components
234: . npoints - The number of quadrature points
235: . points - The coordinates of each quadrature point
236: - weights - The weight of each quadrature point

238:   Level: intermediate

240: .keywords: PetscQuadrature, quadrature
241: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
242: @*/
243: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
244: {
247:   if (dim) {
249:     *dim = q->dim;
250:   }
251:   if (Nc) {
253:     *Nc = q->Nc;
254:   }
255:   if (npoints) {
257:     *npoints = q->numPoints;
258:   }
259:   if (points) {
261:     *points = q->points;
262:   }
263:   if (weights) {
265:     *weights = q->weights;
266:   }
267:   return(0);
268: }

270: /*@C
271:   PetscQuadratureSetData - Sets the data defining the quadrature

273:   Not collective

275:   Input Parameters:
276: + q  - The PetscQuadrature object
277: . dim - The spatial dimension
278: , Nc - The number of components
279: . npoints - The number of quadrature points
280: . points - The coordinates of each quadrature point
281: - weights - The weight of each quadrature point

283:   Level: intermediate

285: .keywords: PetscQuadrature, quadrature
286: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
287: @*/
288: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
289: {
292:   if (dim >= 0)     q->dim       = dim;
293:   if (Nc >= 0)      q->Nc        = Nc;
294:   if (npoints >= 0) q->numPoints = npoints;
295:   if (points) {
297:     q->points = points;
298:   }
299:   if (weights) {
301:     q->weights = weights;
302:   }
303:   return(0);
304: }

306: /*@C
307:   PetscQuadratureView - Views a PetscQuadrature object

309:   Collective on PetscQuadrature

311:   Input Parameters:
312: + q  - The PetscQuadrature object
313: - viewer - The PetscViewer object

315:   Level: beginner

317: .keywords: PetscQuadrature, quadrature, view
318: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
319: @*/
320: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
321: {
322:   PetscInt       q, d, c;

326:   PetscObjectPrintClassNamePrefixType((PetscObject)quad,viewer);
327:   if (quad->Nc > 1) {PetscViewerASCIIPrintf(viewer, "Quadrature on %D points with %D components\n  (", quad->numPoints, quad->Nc);}
328:   else              {PetscViewerASCIIPrintf(viewer, "Quadrature on %D points\n  (", quad->numPoints);}
329:   for (q = 0; q < quad->numPoints; ++q) {
330:     for (d = 0; d < quad->dim; ++d) {
331:       if (d) PetscViewerASCIIPrintf(viewer, ", ");
332:       PetscViewerASCIIPrintf(viewer, "%g\n", (double)quad->points[q*quad->dim+d]);
333:     }
334:     if (quad->Nc > 1) {
335:       PetscViewerASCIIPrintf(viewer, ") (");
336:       for (c = 0; c < quad->Nc; ++c) {
337:         if (c) PetscViewerASCIIPrintf(viewer, ", ");
338:         PetscViewerASCIIPrintf(viewer, "%g", (double)quad->weights[q*quad->Nc+c]);
339:       }
340:       PetscViewerASCIIPrintf(viewer, ")\n");
341:     } else {
342:       PetscViewerASCIIPrintf(viewer, ") %g\n", (double)quad->weights[q]);
343:     }
344:   }
345:   return(0);
346: }

348: /*@C
349:   PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement

351:   Not collective

353:   Input Parameter:
354: + q - The original PetscQuadrature
355: . numSubelements - The number of subelements the original element is divided into
356: . v0 - An array of the initial points for each subelement
357: - jac - An array of the Jacobian mappings from the reference to each subelement

359:   Output Parameters:
360: . dim - The dimension

362:   Note: Together v0 and jac define an affine mapping from the original reference element to each subelement

364:   Level: intermediate

366: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
367: @*/
368: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
369: {
370:   const PetscReal *points,    *weights;
371:   PetscReal       *pointsRef, *weightsRef;
372:   PetscInt         dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
373:   PetscErrorCode   ierr;

380:   PetscQuadratureCreate(PETSC_COMM_SELF, qref);
381:   PetscQuadratureGetOrder(q, &order);
382:   PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
383:   npointsRef = npoints*numSubelements;
384:   PetscMalloc1(npointsRef*dim,&pointsRef);
385:   PetscMalloc1(npointsRef*Nc, &weightsRef);
386:   for (c = 0; c < numSubelements; ++c) {
387:     for (p = 0; p < npoints; ++p) {
388:       for (d = 0; d < dim; ++d) {
389:         pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
390:         for (e = 0; e < dim; ++e) {
391:           pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
392:         }
393:       }
394:       /* Could also use detJ here */
395:       for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements;
396:     }
397:   }
398:   PetscQuadratureSetOrder(*qref, order);
399:   PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);
400:   return(0);
401: }

403: /*@
404:    PetscDTLegendreEval - evaluate Legendre polynomial at points

406:    Not Collective

408:    Input Arguments:
409: +  npoints - number of spatial points to evaluate at
410: .  points - array of locations to evaluate at
411: .  ndegree - number of basis degrees to evaluate
412: -  degrees - sorted array of degrees to evaluate

414:    Output Arguments:
415: +  B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
416: .  D - row-oriented derivative evaluation matrix (or NULL)
417: -  D2 - row-oriented second derivative evaluation matrix (or NULL)

419:    Level: intermediate

421: .seealso: PetscDTGaussQuadrature()
422: @*/
423: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
424: {
425:   PetscInt i,maxdegree;

428:   if (!npoints || !ndegree) return(0);
429:   maxdegree = degrees[ndegree-1];
430:   for (i=0; i<npoints; i++) {
431:     PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
432:     PetscInt  j,k;
433:     x    = points[i];
434:     pm2  = 0;
435:     pm1  = 1;
436:     pd2  = 0;
437:     pd1  = 0;
438:     pdd2 = 0;
439:     pdd1 = 0;
440:     k    = 0;
441:     if (degrees[k] == 0) {
442:       if (B) B[i*ndegree+k] = pm1;
443:       if (D) D[i*ndegree+k] = pd1;
444:       if (D2) D2[i*ndegree+k] = pdd1;
445:       k++;
446:     }
447:     for (j=1; j<=maxdegree; j++,k++) {
448:       PetscReal p,d,dd;
449:       p    = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
450:       d    = pd2 + (2*j-1)*pm1;
451:       dd   = pdd2 + (2*j-1)*pd1;
452:       pm2  = pm1;
453:       pm1  = p;
454:       pd2  = pd1;
455:       pd1  = d;
456:       pdd2 = pdd1;
457:       pdd1 = dd;
458:       if (degrees[k] == j) {
459:         if (B) B[i*ndegree+k] = p;
460:         if (D) D[i*ndegree+k] = d;
461:         if (D2) D2[i*ndegree+k] = dd;
462:       }
463:     }
464:   }
465:   return(0);
466: }

468: /*@
469:    PetscDTGaussQuadrature - create Gauss quadrature

471:    Not Collective

473:    Input Arguments:
474: +  npoints - number of points
475: .  a - left end of interval (often-1)
476: -  b - right end of interval (often +1)

478:    Output Arguments:
479: +  x - quadrature points
480: -  w - quadrature weights

482:    Level: intermediate

484:    References:
485: .   1. - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969.

487: .seealso: PetscDTLegendreEval()
488: @*/
489: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
490: {
492:   PetscInt       i;
493:   PetscReal      *work;
494:   PetscScalar    *Z;
495:   PetscBLASInt   N,LDZ,info;

498:   PetscCitationsRegister(GaussCitation, &GaussCite);
499:   /* Set up the Golub-Welsch system */
500:   for (i=0; i<npoints; i++) {
501:     x[i] = 0;                   /* diagonal is 0 */
502:     if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
503:   }
504:   PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
505:   PetscBLASIntCast(npoints,&N);
506:   LDZ  = N;
507:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
508:   PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
509:   PetscFPTrapPop();
510:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");

512:   for (i=0; i<(npoints+1)/2; i++) {
513:     PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
514:     x[i]           = (a+b)/2 - y*(b-a)/2;
515:     if (x[i] == -0.0) x[i] = 0.0;
516:     x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;

518:     w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
519:   }
520:   PetscFree2(Z,work);
521:   return(0);
522: }

524: /*@
525:   PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature

527:   Not Collective

529:   Input Arguments:
530: + dim     - The spatial dimension
531: . Nc      - The number of components
532: . npoints - number of points in one dimension
533: . a       - left end of interval (often-1)
534: - b       - right end of interval (often +1)

536:   Output Argument:
537: . q - A PetscQuadrature object

539:   Level: intermediate

541: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
542: @*/
543: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
544: {
545:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
546:   PetscReal     *x, *w, *xw, *ww;

550:   PetscMalloc1(totpoints*dim,&x);
551:   PetscMalloc1(totpoints*Nc,&w);
552:   /* Set up the Golub-Welsch system */
553:   switch (dim) {
554:   case 0:
555:     PetscFree(x);
556:     PetscFree(w);
557:     PetscMalloc1(1, &x);
558:     PetscMalloc1(Nc, &w);
559:     x[0] = 0.0;
560:     for (c = 0; c < Nc; ++c) w[c] = 1.0;
561:     break;
562:   case 1:
563:     PetscMalloc1(npoints,&ww);
564:     PetscDTGaussQuadrature(npoints, a, b, x, ww);
565:     for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i];
566:     PetscFree(ww);
567:     break;
568:   case 2:
569:     PetscMalloc2(npoints,&xw,npoints,&ww);
570:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
571:     for (i = 0; i < npoints; ++i) {
572:       for (j = 0; j < npoints; ++j) {
573:         x[(i*npoints+j)*dim+0] = xw[i];
574:         x[(i*npoints+j)*dim+1] = xw[j];
575:         for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j];
576:       }
577:     }
578:     PetscFree2(xw,ww);
579:     break;
580:   case 3:
581:     PetscMalloc2(npoints,&xw,npoints,&ww);
582:     PetscDTGaussQuadrature(npoints, a, b, xw, ww);
583:     for (i = 0; i < npoints; ++i) {
584:       for (j = 0; j < npoints; ++j) {
585:         for (k = 0; k < npoints; ++k) {
586:           x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
587:           x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
588:           x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
589:           for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k];
590:         }
591:       }
592:     }
593:     PetscFree2(xw,ww);
594:     break;
595:   default:
596:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
597:   }
598:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
599:   PetscQuadratureSetOrder(*q, npoints-1);
600:   PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
601:   return(0);
602: }

604: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
605:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
606: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
607: {
608:   PetscReal f = 1.0;
609:   PetscInt  i;

612:   for (i = 1; i < n+1; ++i) f *= i;
613:   *factorial = f;
614:   return(0);
615: }

617: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
618:    Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
619: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
620: {
621:   PetscReal apb, pn1, pn2;
622:   PetscInt  k;

625:   if (!n) {*P = 1.0; return(0);}
626:   if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
627:   apb = a + b;
628:   pn2 = 1.0;
629:   pn1 = 0.5 * (a - b + (apb + 2.0) * x);
630:   *P  = 0.0;
631:   for (k = 2; k < n+1; ++k) {
632:     PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
633:     PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
634:     PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
635:     PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);

637:     a2  = a2 / a1;
638:     a3  = a3 / a1;
639:     a4  = a4 / a1;
640:     *P  = (a2 + a3 * x) * pn1 - a4 * pn2;
641:     pn2 = pn1;
642:     pn1 = *P;
643:   }
644:   return(0);
645: }

647: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
648: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
649: {
650:   PetscReal      nP;

654:   if (!n) {*P = 0.0; return(0);}
655:   PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
656:   *P   = 0.5 * (a + b + n + 1) * nP;
657:   return(0);
658: }

660: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
661: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
662: {
664:   *xi  = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
665:   *eta = y;
666:   return(0);
667: }

669: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
670: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
671: {
673:   *xi   = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
674:   *eta  = 0.5  * (1.0 + y) * (1.0 - z) - 1.0;
675:   *zeta = z;
676:   return(0);
677: }

679: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
680: {
681:   PetscInt       maxIter = 100;
682:   PetscReal      eps     = 1.0e-8;
683:   PetscReal      a1, a2, a3, a4, a5, a6;
684:   PetscInt       k;


689:   a1      = PetscPowReal(2.0, a+b+1);
690: #if defined(PETSC_HAVE_TGAMMA)
691:   a2      = PetscTGamma(a + npoints + 1);
692:   a3      = PetscTGamma(b + npoints + 1);
693:   a4      = PetscTGamma(a + b + npoints + 1);
694: #else
695:   {
696:     PetscInt ia, ib;

698:     ia = (PetscInt) a;
699:     ib = (PetscInt) b;
700:     if (ia == a && ib == b && ia + npoints + 1 > 0 && ib + npoints + 1 > 0 && ia + ib + npoints + 1 > 0) { /* All gamma(x) terms are (x-1)! terms */
701:       PetscDTFactorial_Internal(ia + npoints, &a2);
702:       PetscDTFactorial_Internal(ib + npoints, &a3);
703:       PetscDTFactorial_Internal(ia + ib + npoints, &a4);
704:     } else {
705:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
706:     }
707:   }
708: #endif

710:   PetscDTFactorial_Internal(npoints, &a5);
711:   a6   = a1 * a2 * a3 / a4 / a5;
712:   /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
713:    Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
714:   for (k = 0; k < npoints; ++k) {
715:     PetscReal r = -PetscCosReal((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP;
716:     PetscInt  j;

718:     if (k > 0) r = 0.5 * (r + x[k-1]);
719:     for (j = 0; j < maxIter; ++j) {
720:       PetscReal s = 0.0, delta, f, fp;
721:       PetscInt  i;

723:       for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
724:       PetscDTComputeJacobi(a, b, npoints, r, &f);
725:       PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
726:       delta = f / (fp - f * s);
727:       r     = r - delta;
728:       if (PetscAbsReal(delta) < eps) break;
729:     }
730:     x[k] = r;
731:     PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
732:     w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
733:   }
734:   return(0);
735: }

737: /*@C
738:   PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex

740:   Not Collective

742:   Input Arguments:
743: + dim     - The simplex dimension
744: . Nc      - The number of components
745: . npoints - The number of points in one dimension
746: . a       - left end of interval (often-1)
747: - b       - right end of interval (often +1)

749:   Output Argument:
750: . q - A PetscQuadrature object

752:   Level: intermediate

754:   References:
755: .  1. - Karniadakis and Sherwin.  FIAT

757: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
758: @*/
759: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
760: {
761:   PetscInt       totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints;
762:   PetscReal     *px, *wx, *py, *wy, *pz, *wz, *x, *w;
763:   PetscInt       i, j, k, c;

767:   if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
768:   PetscMalloc1(totpoints*dim, &x);
769:   PetscMalloc1(totpoints*Nc, &w);
770:   switch (dim) {
771:   case 0:
772:     PetscFree(x);
773:     PetscFree(w);
774:     PetscMalloc1(1, &x);
775:     PetscMalloc1(Nc, &w);
776:     x[0] = 0.0;
777:     for (c = 0; c < Nc; ++c) w[c] = 1.0;
778:     break;
779:   case 1:
780:     PetscMalloc1(npoints,&wx);
781:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, x, wx);
782:     for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = wx[i];
783:     PetscFree(wx);
784:     break;
785:   case 2:
786:     PetscMalloc4(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy);
787:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
788:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
789:     for (i = 0; i < npoints; ++i) {
790:       for (j = 0; j < npoints; ++j) {
791:         PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);
792:         for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = 0.5 * wx[i] * wy[j];
793:       }
794:     }
795:     PetscFree4(px,wx,py,wy);
796:     break;
797:   case 3:
798:     PetscMalloc6(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy,npoints,&pz,npoints,&wz);
799:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
800:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
801:     PetscDTGaussJacobiQuadrature1D_Internal(npoints, 2.0, 0.0, pz, wz);
802:     for (i = 0; i < npoints; ++i) {
803:       for (j = 0; j < npoints; ++j) {
804:         for (k = 0; k < npoints; ++k) {
805:           PetscDTMapCubeToTetrahedron_Internal(px[i], py[j], pz[k], &x[((i*npoints+j)*npoints+k)*3+0], &x[((i*npoints+j)*npoints+k)*3+1], &x[((i*npoints+j)*npoints+k)*3+2]);
806:           for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = 0.125 * wx[i] * wy[j] * wz[k];
807:         }
808:       }
809:     }
810:     PetscFree6(px,wx,py,wy,pz,wz);
811:     break;
812:   default:
813:     SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
814:   }
815:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
816:   PetscQuadratureSetOrder(*q, npoints-1);
817:   PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
818:   return(0);
819: }

821: /*@C
822:   PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell

824:   Not Collective

826:   Input Arguments:
827: + dim   - The cell dimension
828: . level - The number of points in one dimension, 2^l
829: . a     - left end of interval (often-1)
830: - b     - right end of interval (often +1)

832:   Output Argument:
833: . q - A PetscQuadrature object

835:   Level: intermediate

837: .seealso: PetscDTGaussTensorQuadrature()
838: @*/
839: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
840: {
841:   const PetscInt  p     = 16;                        /* Digits of precision in the evaluation */
842:   const PetscReal alpha = (b-a)/2.;                  /* Half-width of the integration interval */
843:   const PetscReal beta  = (b+a)/2.;                  /* Center of the integration interval */
844:   const PetscReal h     = PetscPowReal(2.0, -level); /* Step size, length between x_k */
845:   PetscReal       xk;                                /* Quadrature point x_k on reference domain [-1, 1] */
846:   PetscReal       wk    = 0.5*PETSC_PI;              /* Quadrature weight at x_k */
847:   PetscReal      *x, *w;
848:   PetscInt        K, k, npoints;
849:   PetscErrorCode  ierr;

852:   if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
853:   if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
854:   /* Find K such that the weights are < 32 digits of precision */
855:   for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
856:     wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
857:   }
858:   PetscQuadratureCreate(PETSC_COMM_SELF, q);
859:   PetscQuadratureSetOrder(*q, 2*K+1);
860:   npoints = 2*K-1;
861:   PetscMalloc1(npoints*dim, &x);
862:   PetscMalloc1(npoints, &w);
863:   /* Center term */
864:   x[0] = beta;
865:   w[0] = 0.5*alpha*PETSC_PI;
866:   for (k = 1; k < K; ++k) {
867:     wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
868:     xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h));
869:     x[2*k-1] = -alpha*xk+beta;
870:     w[2*k-1] = wk;
871:     x[2*k+0] =  alpha*xk+beta;
872:     w[2*k+0] = wk;
873:   }
874:   PetscQuadratureSetData(*q, dim, 1, npoints, x, w);
875:   return(0);
876: }

878: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
879: {
880:   const PetscInt  p     = 16;        /* Digits of precision in the evaluation */
881:   const PetscReal alpha = (b-a)/2.;  /* Half-width of the integration interval */
882:   const PetscReal beta  = (b+a)/2.;  /* Center of the integration interval */
883:   PetscReal       h     = 1.0;       /* Step size, length between x_k */
884:   PetscInt        l     = 0;         /* Level of refinement, h = 2^{-l} */
885:   PetscReal       osum  = 0.0;       /* Integral on last level */
886:   PetscReal       psum  = 0.0;       /* Integral on the level before the last level */
887:   PetscReal       sum;               /* Integral on current level */
888:   PetscReal       yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
889:   PetscReal       lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
890:   PetscReal       wk;                /* Quadrature weight at x_k */
891:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
892:   PetscInt        d;                 /* Digits of precision in the integral */

895:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
896:   /* Center term */
897:   func(beta, &lval);
898:   sum = 0.5*alpha*PETSC_PI*lval;
899:   /* */
900:   do {
901:     PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
902:     PetscInt  k = 1;

904:     ++l;
905:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
906:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
907:     psum = osum;
908:     osum = sum;
909:     h   *= 0.5;
910:     sum *= 0.5;
911:     do {
912:       wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
913:       yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
914:       lx = -alpha*(1.0 - yk)+beta;
915:       rx =  alpha*(1.0 - yk)+beta;
916:       func(lx, &lval);
917:       func(rx, &rval);
918:       lterm   = alpha*wk*lval;
919:       maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
920:       sum    += lterm;
921:       rterm   = alpha*wk*rval;
922:       maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
923:       sum    += rterm;
924:       ++k;
925:       /* Only need to evaluate every other point on refined levels */
926:       if (l != 1) ++k;
927:     } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */

929:     d1 = PetscLog10Real(PetscAbsReal(sum - osum));
930:     d2 = PetscLog10Real(PetscAbsReal(sum - psum));
931:     d3 = PetscLog10Real(maxTerm) - p;
932:     if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
933:     else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
934:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
935:   } while (d < digits && l < 12);
936:   *sol = sum;

938:   return(0);
939: }

941: #ifdef PETSC_HAVE_MPFR
942: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
943: {
944:   const PetscInt  safetyFactor = 2;  /* Calculate abcissa until 2*p digits */
945:   PetscInt        l            = 0;  /* Level of refinement, h = 2^{-l} */
946:   mpfr_t          alpha;             /* Half-width of the integration interval */
947:   mpfr_t          beta;              /* Center of the integration interval */
948:   mpfr_t          h;                 /* Step size, length between x_k */
949:   mpfr_t          osum;              /* Integral on last level */
950:   mpfr_t          psum;              /* Integral on the level before the last level */
951:   mpfr_t          sum;               /* Integral on current level */
952:   mpfr_t          yk;                /* Quadrature point 1 - x_k on reference domain [-1, 1] */
953:   mpfr_t          lx, rx;            /* Quadrature points to the left and right of 0 on the real domain [a, b] */
954:   mpfr_t          wk;                /* Quadrature weight at x_k */
955:   PetscReal       lval, rval;        /* Terms in the quadature sum to the left and right of 0 */
956:   PetscInt        d;                 /* Digits of precision in the integral */
957:   mpfr_t          pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;

960:   if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
961:   /* Create high precision storage */
962:   mpfr_inits2(PetscCeilReal(safetyFactor*digits*PetscLogReal(10.)/PetscLogReal(2.)), alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
963:   /* Initialization */
964:   mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
965:   mpfr_set_d(beta,  0.5*(b+a), MPFR_RNDN);
966:   mpfr_set_d(osum,  0.0,       MPFR_RNDN);
967:   mpfr_set_d(psum,  0.0,       MPFR_RNDN);
968:   mpfr_set_d(h,     1.0,       MPFR_RNDN);
969:   mpfr_const_pi(pi2, MPFR_RNDN);
970:   mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
971:   /* Center term */
972:   func(0.5*(b+a), &lval);
973:   mpfr_set(sum, pi2, MPFR_RNDN);
974:   mpfr_mul(sum, sum, alpha, MPFR_RNDN);
975:   mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
976:   /* */
977:   do {
978:     PetscReal d1, d2, d3, d4;
979:     PetscInt  k = 1;

981:     ++l;
982:     mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
983:     /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
984:     /* At each level of refinement, h --> h/2 and sum --> sum/2 */
985:     mpfr_set(psum, osum, MPFR_RNDN);
986:     mpfr_set(osum,  sum, MPFR_RNDN);
987:     mpfr_mul_d(h,   h,   0.5, MPFR_RNDN);
988:     mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
989:     do {
990:       mpfr_set_si(kh, k, MPFR_RNDN);
991:       mpfr_mul(kh, kh, h, MPFR_RNDN);
992:       /* Weight */
993:       mpfr_set(wk, h, MPFR_RNDN);
994:       mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
995:       mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
996:       mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
997:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
998:       mpfr_sqr(tmp, tmp, MPFR_RNDN);
999:       mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
1000:       mpfr_div(wk, wk, tmp, MPFR_RNDN);
1001:       /* Abscissa */
1002:       mpfr_set_d(yk, 1.0, MPFR_RNDZ);
1003:       mpfr_cosh(tmp, msinh, MPFR_RNDN);
1004:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1005:       mpfr_exp(tmp, msinh, MPFR_RNDN);
1006:       mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1007:       /* Quadrature points */
1008:       mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
1009:       mpfr_mul(lx, lx, alpha, MPFR_RNDU);
1010:       mpfr_add(lx, lx, beta, MPFR_RNDU);
1011:       mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
1012:       mpfr_mul(rx, rx, alpha, MPFR_RNDD);
1013:       mpfr_add(rx, rx, beta, MPFR_RNDD);
1014:       /* Evaluation */
1015:       func(mpfr_get_d(lx, MPFR_RNDU), &lval);
1016:       func(mpfr_get_d(rx, MPFR_RNDD), &rval);
1017:       /* Update */
1018:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1019:       mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
1020:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
1021:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1022:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1023:       mpfr_set(curTerm, tmp, MPFR_RNDN);
1024:       mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1025:       mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
1026:       mpfr_add(sum, sum, tmp, MPFR_RNDN);
1027:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1028:       mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1029:       mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
1030:       ++k;
1031:       /* Only need to evaluate every other point on refined levels */
1032:       if (l != 1) ++k;
1033:       mpfr_log10(tmp, wk, MPFR_RNDN);
1034:       mpfr_abs(tmp, tmp, MPFR_RNDN);
1035:     } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
1036:     mpfr_sub(tmp, sum, osum, MPFR_RNDN);
1037:     mpfr_abs(tmp, tmp, MPFR_RNDN);
1038:     mpfr_log10(tmp, tmp, MPFR_RNDN);
1039:     d1 = mpfr_get_d(tmp, MPFR_RNDN);
1040:     mpfr_sub(tmp, sum, psum, MPFR_RNDN);
1041:     mpfr_abs(tmp, tmp, MPFR_RNDN);
1042:     mpfr_log10(tmp, tmp, MPFR_RNDN);
1043:     d2 = mpfr_get_d(tmp, MPFR_RNDN);
1044:     mpfr_log10(tmp, maxTerm, MPFR_RNDN);
1045:     d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
1046:     mpfr_log10(tmp, curTerm, MPFR_RNDN);
1047:     d4 = mpfr_get_d(tmp, MPFR_RNDN);
1048:     d  = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1049:   } while (d < digits && l < 8);
1050:   *sol = mpfr_get_d(sum, MPFR_RNDN);
1051:   /* Cleanup */
1052:   mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
1053:   return(0);
1054: }
1055: #else

1057: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1058: {
1059:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
1060: }
1061: #endif

1063: /* Overwrites A. Can only handle full-rank problems with m>=n
1064:  * A in column-major format
1065:  * Ainv in row-major format
1066:  * tau has length m
1067:  * worksize must be >= max(1,n)
1068:  */
1069: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1070: {
1072:   PetscBLASInt   M,N,K,lda,ldb,ldwork,info;
1073:   PetscScalar    *A,*Ainv,*R,*Q,Alpha;

1076: #if defined(PETSC_USE_COMPLEX)
1077:   {
1078:     PetscInt i,j;
1079:     PetscMalloc2(m*n,&A,m*n,&Ainv);
1080:     for (j=0; j<n; j++) {
1081:       for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
1082:     }
1083:     mstride = m;
1084:   }
1085: #else
1086:   A = A_in;
1087:   Ainv = Ainv_out;
1088: #endif

1090:   PetscBLASIntCast(m,&M);
1091:   PetscBLASIntCast(n,&N);
1092:   PetscBLASIntCast(mstride,&lda);
1093:   PetscBLASIntCast(worksize,&ldwork);
1094:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1095:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1096:   PetscFPTrapPop();
1097:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1098:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1100:   /* Extract an explicit representation of Q */
1101:   Q = Ainv;
1102:   PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1103:   K = N;                        /* full rank */
1104:   PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1105:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

1107:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1108:   Alpha = 1.0;
1109:   ldb = lda;
1110:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
1111:   /* Ainv is Q, overwritten with inverse */

1113: #if defined(PETSC_USE_COMPLEX)
1114:   {
1115:     PetscInt i;
1116:     for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
1117:     PetscFree2(A,Ainv);
1118:   }
1119: #endif
1120:   return(0);
1121: }

1123: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
1124: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
1125: {
1127:   PetscReal      *Bv;
1128:   PetscInt       i,j;

1131:   PetscMalloc1((ninterval+1)*ndegree,&Bv);
1132:   /* Point evaluation of L_p on all the source vertices */
1133:   PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
1134:   /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
1135:   for (i=0; i<ninterval; i++) {
1136:     for (j=0; j<ndegree; j++) {
1137:       if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1138:       else           B[i*ndegree+j]   = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1139:     }
1140:   }
1141:   PetscFree(Bv);
1142:   return(0);
1143: }

1145: /*@
1146:    PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals

1148:    Not Collective

1150:    Input Arguments:
1151: +  degree - degree of reconstruction polynomial
1152: .  nsource - number of source intervals
1153: .  sourcex - sorted coordinates of source cell boundaries (length nsource+1)
1154: .  ntarget - number of target intervals
1155: -  targetx - sorted coordinates of target cell boundaries (length ntarget+1)

1157:    Output Arguments:
1158: .  R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]

1160:    Level: advanced

1162: .seealso: PetscDTLegendreEval()
1163: @*/
1164: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
1165: {
1167:   PetscInt       i,j,k,*bdegrees,worksize;
1168:   PetscReal      xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
1169:   PetscScalar    *tau,*work;

1175:   if (degree >= nsource) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Reconstruction degree %D must be less than number of source intervals %D",degree,nsource);
1176: #if defined(PETSC_USE_DEBUG)
1177:   for (i=0; i<nsource; i++) {
1178:     if (sourcex[i] >= sourcex[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Source interval %D has negative orientation (%g,%g)",i,(double)sourcex[i],(double)sourcex[i+1]);
1179:   }
1180:   for (i=0; i<ntarget; i++) {
1181:     if (targetx[i] >= targetx[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Target interval %D has negative orientation (%g,%g)",i,(double)targetx[i],(double)targetx[i+1]);
1182:   }
1183: #endif
1184:   xmin = PetscMin(sourcex[0],targetx[0]);
1185:   xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
1186:   center = (xmin + xmax)/2;
1187:   hscale = (xmax - xmin)/2;
1188:   worksize = nsource;
1189:   PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
1190:   PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
1191:   for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
1192:   for (i=0; i<=degree; i++) bdegrees[i] = i+1;
1193:   PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
1194:   PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
1195:   for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
1196:   PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
1197:   for (i=0; i<ntarget; i++) {
1198:     PetscReal rowsum = 0;
1199:     for (j=0; j<nsource; j++) {
1200:       PetscReal sum = 0;
1201:       for (k=0; k<degree+1; k++) {
1202:         sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
1203:       }
1204:       R[i*nsource+j] = sum;
1205:       rowsum += sum;
1206:     }
1207:     for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
1208:   }
1209:   PetscFree4(bdegrees,sourcey,Bsource,work);
1210:   PetscFree4(tau,Bsinv,targety,Btarget);
1211:   return(0);
1212: }