Actual source code: dt.c
petsc-3.11.4 2019-09-28
1: /* Discretization tools */
3: #include <petscdt.h>
4: #include <petscblaslapack.h>
5: #include <petsc/private/petscimpl.h>
6: #include <petsc/private/dtimpl.h>
7: #include <petscviewer.h>
8: #include <petscdmplex.h>
9: #include <petscdmshell.h>
11: #if defined(PETSC_HAVE_MPFR)
12: #include <mpfr.h>
13: #endif
15: static PetscBool GaussCite = PETSC_FALSE;
16: const char GaussCitation[] = "@article{GolubWelsch1969,\n"
17: " author = {Golub and Welsch},\n"
18: " title = {Calculation of Quadrature Rules},\n"
19: " journal = {Math. Comp.},\n"
20: " volume = {23},\n"
21: " number = {106},\n"
22: " pages = {221--230},\n"
23: " year = {1969}\n}\n";
25: /*@
26: PetscQuadratureCreate - Create a PetscQuadrature object
28: Collective on MPI_Comm
30: Input Parameter:
31: . comm - The communicator for the PetscQuadrature object
33: Output Parameter:
34: . q - The PetscQuadrature object
36: Level: beginner
38: .keywords: PetscQuadrature, quadrature, create
39: .seealso: PetscQuadratureDestroy(), PetscQuadratureGetData()
40: @*/
41: PetscErrorCode PetscQuadratureCreate(MPI_Comm comm, PetscQuadrature *q)
42: {
47: PetscSysInitializePackage();
48: PetscHeaderCreate(*q,PETSC_OBJECT_CLASSID,"PetscQuadrature","Quadrature","DT",comm,PetscQuadratureDestroy,PetscQuadratureView);
49: (*q)->dim = -1;
50: (*q)->Nc = 1;
51: (*q)->order = -1;
52: (*q)->numPoints = 0;
53: (*q)->points = NULL;
54: (*q)->weights = NULL;
55: return(0);
56: }
58: /*@
59: PetscQuadratureDuplicate - Create a deep copy of the PetscQuadrature object
61: Collective on PetscQuadrature
63: Input Parameter:
64: . q - The PetscQuadrature object
66: Output Parameter:
67: . r - The new PetscQuadrature object
69: Level: beginner
71: .keywords: PetscQuadrature, quadrature, clone
72: .seealso: PetscQuadratureCreate(), PetscQuadratureDestroy(), PetscQuadratureGetData()
73: @*/
74: PetscErrorCode PetscQuadratureDuplicate(PetscQuadrature q, PetscQuadrature *r)
75: {
76: PetscInt order, dim, Nc, Nq;
77: const PetscReal *points, *weights;
78: PetscReal *p, *w;
79: PetscErrorCode ierr;
83: PetscQuadratureCreate(PetscObjectComm((PetscObject) q), r);
84: PetscQuadratureGetOrder(q, &order);
85: PetscQuadratureSetOrder(*r, order);
86: PetscQuadratureGetData(q, &dim, &Nc, &Nq, &points, &weights);
87: PetscMalloc1(Nq*dim, &p);
88: PetscMalloc1(Nq*Nc, &w);
89: PetscMemcpy(p, points, Nq*dim * sizeof(PetscReal));
90: PetscMemcpy(w, weights, Nc * Nq * sizeof(PetscReal));
91: PetscQuadratureSetData(*r, dim, Nc, Nq, p, w);
92: return(0);
93: }
95: /*@
96: PetscQuadratureDestroy - Destroys a PetscQuadrature object
98: Collective on PetscQuadrature
100: Input Parameter:
101: . q - The PetscQuadrature object
103: Level: beginner
105: .keywords: PetscQuadrature, quadrature, destroy
106: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
107: @*/
108: PetscErrorCode PetscQuadratureDestroy(PetscQuadrature *q)
109: {
113: if (!*q) return(0);
115: if (--((PetscObject)(*q))->refct > 0) {
116: *q = NULL;
117: return(0);
118: }
119: PetscFree((*q)->points);
120: PetscFree((*q)->weights);
121: PetscHeaderDestroy(q);
122: return(0);
123: }
125: /*@
126: PetscQuadratureGetOrder - Return the order of the method
128: Not collective
130: Input Parameter:
131: . q - The PetscQuadrature object
133: Output Parameter:
134: . order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
136: Level: intermediate
138: .seealso: PetscQuadratureSetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
139: @*/
140: PetscErrorCode PetscQuadratureGetOrder(PetscQuadrature q, PetscInt *order)
141: {
145: *order = q->order;
146: return(0);
147: }
149: /*@
150: PetscQuadratureSetOrder - Return the order of the method
152: Not collective
154: Input Parameters:
155: + q - The PetscQuadrature object
156: - order - The order of the quadrature, i.e. the highest degree polynomial that is exactly integrated
158: Level: intermediate
160: .seealso: PetscQuadratureGetOrder(), PetscQuadratureGetData(), PetscQuadratureSetData()
161: @*/
162: PetscErrorCode PetscQuadratureSetOrder(PetscQuadrature q, PetscInt order)
163: {
166: q->order = order;
167: return(0);
168: }
170: /*@
171: PetscQuadratureGetNumComponents - Return the number of components for functions to be integrated
173: Not collective
175: Input Parameter:
176: . q - The PetscQuadrature object
178: Output Parameter:
179: . Nc - The number of components
181: Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
183: Level: intermediate
185: .seealso: PetscQuadratureSetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
186: @*/
187: PetscErrorCode PetscQuadratureGetNumComponents(PetscQuadrature q, PetscInt *Nc)
188: {
192: *Nc = q->Nc;
193: return(0);
194: }
196: /*@
197: PetscQuadratureSetNumComponents - Return the number of components for functions to be integrated
199: Not collective
201: Input Parameters:
202: + q - The PetscQuadrature object
203: - Nc - The number of components
205: Note: We are performing an integral int f(x) . w(x) dx, where both f and w (the weight) have Nc components.
207: Level: intermediate
209: .seealso: PetscQuadratureGetNumComponents(), PetscQuadratureGetData(), PetscQuadratureSetData()
210: @*/
211: PetscErrorCode PetscQuadratureSetNumComponents(PetscQuadrature q, PetscInt Nc)
212: {
215: q->Nc = Nc;
216: return(0);
217: }
219: /*@C
220: PetscQuadratureGetData - Returns the data defining the quadrature
222: Not collective
224: Input Parameter:
225: . q - The PetscQuadrature object
227: Output Parameters:
228: + dim - The spatial dimension
229: . Nc - The number of components
230: . npoints - The number of quadrature points
231: . points - The coordinates of each quadrature point
232: - weights - The weight of each quadrature point
234: Level: intermediate
236: Fortran Notes:
237: From Fortran you must call PetscQuadratureRestoreData() when you are done with the data
239: .keywords: PetscQuadrature, quadrature
240: .seealso: PetscQuadratureCreate(), PetscQuadratureSetData()
241: @*/
242: PetscErrorCode PetscQuadratureGetData(PetscQuadrature q, PetscInt *dim, PetscInt *Nc, PetscInt *npoints, const PetscReal *points[], const PetscReal *weights[])
243: {
246: if (dim) {
248: *dim = q->dim;
249: }
250: if (Nc) {
252: *Nc = q->Nc;
253: }
254: if (npoints) {
256: *npoints = q->numPoints;
257: }
258: if (points) {
260: *points = q->points;
261: }
262: if (weights) {
264: *weights = q->weights;
265: }
266: return(0);
267: }
269: /*@C
270: PetscQuadratureSetData - Sets the data defining the quadrature
272: Not collective
274: Input Parameters:
275: + q - The PetscQuadrature object
276: . dim - The spatial dimension
277: . Nc - The number of components
278: . npoints - The number of quadrature points
279: . points - The coordinates of each quadrature point
280: - weights - The weight of each quadrature point
282: Note: This routine owns the references to points and weights, so they must be allocated using PetscMalloc() and the user should not free them.
284: Level: intermediate
286: .keywords: PetscQuadrature, quadrature
287: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
288: @*/
289: PetscErrorCode PetscQuadratureSetData(PetscQuadrature q, PetscInt dim, PetscInt Nc, PetscInt npoints, const PetscReal points[], const PetscReal weights[])
290: {
293: if (dim >= 0) q->dim = dim;
294: if (Nc >= 0) q->Nc = Nc;
295: if (npoints >= 0) q->numPoints = npoints;
296: if (points) {
298: q->points = points;
299: }
300: if (weights) {
302: q->weights = weights;
303: }
304: return(0);
305: }
307: static PetscErrorCode PetscQuadratureView_Ascii(PetscQuadrature quad, PetscViewer v)
308: {
309: PetscInt q, d, c;
310: PetscViewerFormat format;
311: PetscErrorCode ierr;
314: if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "Quadrature on %D points with %D components of order %D\n", quad->numPoints, quad->Nc, quad->order);}
315: else {PetscViewerASCIIPrintf(v, "Quadrature on %D points of order %D\n", quad->numPoints, quad->order);}
316: PetscViewerGetFormat(v, &format);
317: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
318: for (q = 0; q < quad->numPoints; ++q) {
319: PetscViewerASCIIPrintf(v, "(");
320: PetscViewerASCIIUseTabs(v, PETSC_FALSE);
321: for (d = 0; d < quad->dim; ++d) {
322: if (d) {PetscViewerASCIIPrintf(v, ", ");}
323: PetscViewerASCIIPrintf(v, "%+g", (double)quad->points[q*quad->dim+d]);
324: }
325: PetscViewerASCIIPrintf(v, ") ");
326: if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, "(");}
327: for (c = 0; c < quad->Nc; ++c) {
328: if (c) {PetscViewerASCIIPrintf(v, ", ");}
329: PetscViewerASCIIPrintf(v, "%g", (double)quad->weights[q*quad->Nc+c]);
330: }
331: if (quad->Nc > 1) {PetscViewerASCIIPrintf(v, ")");}
332: PetscViewerASCIIPrintf(v, "\n");
333: PetscViewerASCIIUseTabs(v, PETSC_TRUE);
334: }
335: return(0);
336: }
338: /*@C
339: PetscQuadratureView - Views a PetscQuadrature object
341: Collective on PetscQuadrature
343: Input Parameters:
344: + quad - The PetscQuadrature object
345: - viewer - The PetscViewer object
347: Level: beginner
349: .keywords: PetscQuadrature, quadrature, view
350: .seealso: PetscQuadratureCreate(), PetscQuadratureGetData()
351: @*/
352: PetscErrorCode PetscQuadratureView(PetscQuadrature quad, PetscViewer viewer)
353: {
354: PetscBool iascii;
360: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) quad), &viewer);}
361: PetscObjectPrintClassNamePrefixType((PetscObject)quad, viewer);
362: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
363: PetscViewerASCIIPushTab(viewer);
364: if (iascii) {PetscQuadratureView_Ascii(quad, viewer);}
365: PetscViewerASCIIPopTab(viewer);
366: return(0);
367: }
369: /*@C
370: PetscQuadratureExpandComposite - Return a quadrature over the composite element, which has the original quadrature in each subelement
372: Not collective
374: Input Parameter:
375: + q - The original PetscQuadrature
376: . numSubelements - The number of subelements the original element is divided into
377: . v0 - An array of the initial points for each subelement
378: - jac - An array of the Jacobian mappings from the reference to each subelement
380: Output Parameters:
381: . dim - The dimension
383: Note: Together v0 and jac define an affine mapping from the original reference element to each subelement
385: Not available from Fortran
387: Level: intermediate
389: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
390: @*/
391: PetscErrorCode PetscQuadratureExpandComposite(PetscQuadrature q, PetscInt numSubelements, const PetscReal v0[], const PetscReal jac[], PetscQuadrature *qref)
392: {
393: const PetscReal *points, *weights;
394: PetscReal *pointsRef, *weightsRef;
395: PetscInt dim, Nc, order, npoints, npointsRef, c, p, cp, d, e;
396: PetscErrorCode ierr;
403: PetscQuadratureCreate(PETSC_COMM_SELF, qref);
404: PetscQuadratureGetOrder(q, &order);
405: PetscQuadratureGetData(q, &dim, &Nc, &npoints, &points, &weights);
406: npointsRef = npoints*numSubelements;
407: PetscMalloc1(npointsRef*dim,&pointsRef);
408: PetscMalloc1(npointsRef*Nc, &weightsRef);
409: for (c = 0; c < numSubelements; ++c) {
410: for (p = 0; p < npoints; ++p) {
411: for (d = 0; d < dim; ++d) {
412: pointsRef[(c*npoints + p)*dim+d] = v0[c*dim+d];
413: for (e = 0; e < dim; ++e) {
414: pointsRef[(c*npoints + p)*dim+d] += jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
415: }
416: }
417: /* Could also use detJ here */
418: for (cp = 0; cp < Nc; ++cp) weightsRef[(c*npoints+p)*Nc+cp] = weights[p*Nc+cp]/numSubelements;
419: }
420: }
421: PetscQuadratureSetOrder(*qref, order);
422: PetscQuadratureSetData(*qref, dim, Nc, npointsRef, pointsRef, weightsRef);
423: return(0);
424: }
426: /*@
427: PetscDTLegendreEval - evaluate Legendre polynomial at points
429: Not Collective
431: Input Arguments:
432: + npoints - number of spatial points to evaluate at
433: . points - array of locations to evaluate at
434: . ndegree - number of basis degrees to evaluate
435: - degrees - sorted array of degrees to evaluate
437: Output Arguments:
438: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
439: . D - row-oriented derivative evaluation matrix (or NULL)
440: - D2 - row-oriented second derivative evaluation matrix (or NULL)
442: Level: intermediate
444: .seealso: PetscDTGaussQuadrature()
445: @*/
446: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
447: {
448: PetscInt i,maxdegree;
451: if (!npoints || !ndegree) return(0);
452: maxdegree = degrees[ndegree-1];
453: for (i=0; i<npoints; i++) {
454: PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
455: PetscInt j,k;
456: x = points[i];
457: pm2 = 0;
458: pm1 = 1;
459: pd2 = 0;
460: pd1 = 0;
461: pdd2 = 0;
462: pdd1 = 0;
463: k = 0;
464: if (degrees[k] == 0) {
465: if (B) B[i*ndegree+k] = pm1;
466: if (D) D[i*ndegree+k] = pd1;
467: if (D2) D2[i*ndegree+k] = pdd1;
468: k++;
469: }
470: for (j=1; j<=maxdegree; j++,k++) {
471: PetscReal p,d,dd;
472: p = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
473: d = pd2 + (2*j-1)*pm1;
474: dd = pdd2 + (2*j-1)*pd1;
475: pm2 = pm1;
476: pm1 = p;
477: pd2 = pd1;
478: pd1 = d;
479: pdd2 = pdd1;
480: pdd1 = dd;
481: if (degrees[k] == j) {
482: if (B) B[i*ndegree+k] = p;
483: if (D) D[i*ndegree+k] = d;
484: if (D2) D2[i*ndegree+k] = dd;
485: }
486: }
487: }
488: return(0);
489: }
491: /*@
492: PetscDTGaussQuadrature - create Gauss quadrature
494: Not Collective
496: Input Arguments:
497: + npoints - number of points
498: . a - left end of interval (often-1)
499: - b - right end of interval (often +1)
501: Output Arguments:
502: + x - quadrature points
503: - w - quadrature weights
505: Level: intermediate
507: References:
508: . 1. - Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 1969.
510: .seealso: PetscDTLegendreEval()
511: @*/
512: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
513: {
515: PetscInt i;
516: PetscReal *work;
517: PetscScalar *Z;
518: PetscBLASInt N,LDZ,info;
521: PetscCitationsRegister(GaussCitation, &GaussCite);
522: /* Set up the Golub-Welsch system */
523: for (i=0; i<npoints; i++) {
524: x[i] = 0; /* diagonal is 0 */
525: if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
526: }
527: PetscMalloc2(npoints*npoints,&Z,PetscMax(1,2*npoints-2),&work);
528: PetscBLASIntCast(npoints,&N);
529: LDZ = N;
530: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
531: PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
532: PetscFPTrapPop();
533: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");
535: for (i=0; i<(npoints+1)/2; i++) {
536: PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
537: x[i] = (a+b)/2 - y*(b-a)/2;
538: if (x[i] == -0.0) x[i] = 0.0;
539: x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;
541: w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
542: }
543: PetscFree2(Z,work);
544: return(0);
545: }
547: /*@
548: PetscDTGaussTensorQuadrature - creates a tensor-product Gauss quadrature
550: Not Collective
552: Input Arguments:
553: + dim - The spatial dimension
554: . Nc - The number of components
555: . npoints - number of points in one dimension
556: . a - left end of interval (often-1)
557: - b - right end of interval (often +1)
559: Output Argument:
560: . q - A PetscQuadrature object
562: Level: intermediate
564: .seealso: PetscDTGaussQuadrature(), PetscDTLegendreEval()
565: @*/
566: PetscErrorCode PetscDTGaussTensorQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
567: {
568: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints, i, j, k, c;
569: PetscReal *x, *w, *xw, *ww;
573: PetscMalloc1(totpoints*dim,&x);
574: PetscMalloc1(totpoints*Nc,&w);
575: /* Set up the Golub-Welsch system */
576: switch (dim) {
577: case 0:
578: PetscFree(x);
579: PetscFree(w);
580: PetscMalloc1(1, &x);
581: PetscMalloc1(Nc, &w);
582: x[0] = 0.0;
583: for (c = 0; c < Nc; ++c) w[c] = 1.0;
584: break;
585: case 1:
586: PetscMalloc1(npoints,&ww);
587: PetscDTGaussQuadrature(npoints, a, b, x, ww);
588: for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = ww[i];
589: PetscFree(ww);
590: break;
591: case 2:
592: PetscMalloc2(npoints,&xw,npoints,&ww);
593: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
594: for (i = 0; i < npoints; ++i) {
595: for (j = 0; j < npoints; ++j) {
596: x[(i*npoints+j)*dim+0] = xw[i];
597: x[(i*npoints+j)*dim+1] = xw[j];
598: for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = ww[i] * ww[j];
599: }
600: }
601: PetscFree2(xw,ww);
602: break;
603: case 3:
604: PetscMalloc2(npoints,&xw,npoints,&ww);
605: PetscDTGaussQuadrature(npoints, a, b, xw, ww);
606: for (i = 0; i < npoints; ++i) {
607: for (j = 0; j < npoints; ++j) {
608: for (k = 0; k < npoints; ++k) {
609: x[((i*npoints+j)*npoints+k)*dim+0] = xw[i];
610: x[((i*npoints+j)*npoints+k)*dim+1] = xw[j];
611: x[((i*npoints+j)*npoints+k)*dim+2] = xw[k];
612: for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = ww[i] * ww[j] * ww[k];
613: }
614: }
615: }
616: PetscFree2(xw,ww);
617: break;
618: default:
619: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
620: }
621: PetscQuadratureCreate(PETSC_COMM_SELF, q);
622: PetscQuadratureSetOrder(*q, 2*npoints-1);
623: PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
624: PetscObjectChangeTypeName((PetscObject)*q,"GaussTensor");
625: return(0);
626: }
628: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
629: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
630: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
631: {
632: PetscReal f = 1.0;
633: PetscInt i;
636: for (i = 1; i < n+1; ++i) f *= i;
637: *factorial = f;
638: return(0);
639: }
641: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
642: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
643: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
644: {
645: PetscReal apb, pn1, pn2;
646: PetscInt k;
649: if (!n) {*P = 1.0; return(0);}
650: if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
651: apb = a + b;
652: pn2 = 1.0;
653: pn1 = 0.5 * (a - b + (apb + 2.0) * x);
654: *P = 0.0;
655: for (k = 2; k < n+1; ++k) {
656: PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
657: PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
658: PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
659: PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);
661: a2 = a2 / a1;
662: a3 = a3 / a1;
663: a4 = a4 / a1;
664: *P = (a2 + a3 * x) * pn1 - a4 * pn2;
665: pn2 = pn1;
666: pn1 = *P;
667: }
668: return(0);
669: }
671: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
672: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
673: {
674: PetscReal nP;
678: if (!n) {*P = 0.0; return(0);}
679: PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
680: *P = 0.5 * (a + b + n + 1) * nP;
681: return(0);
682: }
684: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
685: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
686: {
688: *xi = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
689: *eta = y;
690: return(0);
691: }
693: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
694: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
695: {
697: *xi = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
698: *eta = 0.5 * (1.0 + y) * (1.0 - z) - 1.0;
699: *zeta = z;
700: return(0);
701: }
703: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
704: {
705: PetscInt maxIter = 100;
706: PetscReal eps = 1.0e-8;
707: PetscReal a1, a2, a3, a4, a5, a6;
708: PetscInt k;
713: a1 = PetscPowReal(2.0, a+b+1);
714: #if defined(PETSC_HAVE_TGAMMA)
715: a2 = PetscTGamma(a + npoints + 1);
716: a3 = PetscTGamma(b + npoints + 1);
717: a4 = PetscTGamma(a + b + npoints + 1);
718: #else
719: {
720: PetscInt ia, ib;
722: ia = (PetscInt) a;
723: ib = (PetscInt) b;
724: 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 */
725: PetscDTFactorial_Internal(ia + npoints, &a2);
726: PetscDTFactorial_Internal(ib + npoints, &a3);
727: PetscDTFactorial_Internal(ia + ib + npoints, &a4);
728: } else {
729: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
730: }
731: }
732: #endif
734: PetscDTFactorial_Internal(npoints, &a5);
735: a6 = a1 * a2 * a3 / a4 / a5;
736: /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
737: Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
738: for (k = 0; k < npoints; ++k) {
739: PetscReal r = -PetscCosReal((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP;
740: PetscInt j;
742: if (k > 0) r = 0.5 * (r + x[k-1]);
743: for (j = 0; j < maxIter; ++j) {
744: PetscReal s = 0.0, delta, f, fp;
745: PetscInt i;
747: for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
748: PetscDTComputeJacobi(a, b, npoints, r, &f);
749: PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
750: delta = f / (fp - f * s);
751: r = r - delta;
752: if (PetscAbsReal(delta) < eps) break;
753: }
754: x[k] = r;
755: PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
756: w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
757: }
758: return(0);
759: }
761: /*@
762: PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex
764: Not Collective
766: Input Arguments:
767: + dim - The simplex dimension
768: . Nc - The number of components
769: . npoints - The number of points in one dimension
770: . a - left end of interval (often-1)
771: - b - right end of interval (often +1)
773: Output Argument:
774: . q - A PetscQuadrature object
776: Level: intermediate
778: References:
779: . 1. - Karniadakis and Sherwin. FIAT
781: .seealso: PetscDTGaussTensorQuadrature(), PetscDTGaussQuadrature()
782: @*/
783: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt Nc, PetscInt npoints, PetscReal a, PetscReal b, PetscQuadrature *q)
784: {
785: PetscInt totpoints = dim > 1 ? dim > 2 ? npoints*PetscSqr(npoints) : PetscSqr(npoints) : npoints;
786: PetscReal *px, *wx, *py, *wy, *pz, *wz, *x, *w;
787: PetscInt i, j, k, c;
791: if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
792: PetscMalloc1(totpoints*dim, &x);
793: PetscMalloc1(totpoints*Nc, &w);
794: switch (dim) {
795: case 0:
796: PetscFree(x);
797: PetscFree(w);
798: PetscMalloc1(1, &x);
799: PetscMalloc1(Nc, &w);
800: x[0] = 0.0;
801: for (c = 0; c < Nc; ++c) w[c] = 1.0;
802: break;
803: case 1:
804: PetscMalloc1(npoints,&wx);
805: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, x, wx);
806: for (i = 0; i < npoints; ++i) for (c = 0; c < Nc; ++c) w[i*Nc+c] = wx[i];
807: PetscFree(wx);
808: break;
809: case 2:
810: PetscMalloc4(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy);
811: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
812: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
813: for (i = 0; i < npoints; ++i) {
814: for (j = 0; j < npoints; ++j) {
815: PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);
816: for (c = 0; c < Nc; ++c) w[(i*npoints+j)*Nc+c] = 0.5 * wx[i] * wy[j];
817: }
818: }
819: PetscFree4(px,wx,py,wy);
820: break;
821: case 3:
822: PetscMalloc6(npoints,&px,npoints,&wx,npoints,&py,npoints,&wy,npoints,&pz,npoints,&wz);
823: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
824: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
825: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 2.0, 0.0, pz, wz);
826: for (i = 0; i < npoints; ++i) {
827: for (j = 0; j < npoints; ++j) {
828: for (k = 0; k < npoints; ++k) {
829: 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]);
830: for (c = 0; c < Nc; ++c) w[((i*npoints+j)*npoints+k)*Nc+c] = 0.125 * wx[i] * wy[j] * wz[k];
831: }
832: }
833: }
834: PetscFree6(px,wx,py,wy,pz,wz);
835: break;
836: default:
837: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
838: }
839: PetscQuadratureCreate(PETSC_COMM_SELF, q);
840: PetscQuadratureSetOrder(*q, 2*npoints-1);
841: PetscQuadratureSetData(*q, dim, Nc, totpoints, x, w);
842: PetscObjectChangeTypeName((PetscObject)*q,"GaussJacobi");
843: return(0);
844: }
846: /*@
847: PetscDTTanhSinhTensorQuadrature - create tanh-sinh quadrature for a tensor product cell
849: Not Collective
851: Input Arguments:
852: + dim - The cell dimension
853: . level - The number of points in one dimension, 2^l
854: . a - left end of interval (often-1)
855: - b - right end of interval (often +1)
857: Output Argument:
858: . q - A PetscQuadrature object
860: Level: intermediate
862: .seealso: PetscDTGaussTensorQuadrature()
863: @*/
864: PetscErrorCode PetscDTTanhSinhTensorQuadrature(PetscInt dim, PetscInt level, PetscReal a, PetscReal b, PetscQuadrature *q)
865: {
866: const PetscInt p = 16; /* Digits of precision in the evaluation */
867: const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
868: const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
869: const PetscReal h = PetscPowReal(2.0, -level); /* Step size, length between x_k */
870: PetscReal xk; /* Quadrature point x_k on reference domain [-1, 1] */
871: PetscReal wk = 0.5*PETSC_PI; /* Quadrature weight at x_k */
872: PetscReal *x, *w;
873: PetscInt K, k, npoints;
874: PetscErrorCode ierr;
877: if (dim > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Dimension %d not yet implemented", dim);
878: if (!level) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a number of significant digits");
879: /* Find K such that the weights are < 32 digits of precision */
880: for (K = 1; PetscAbsReal(PetscLog10Real(wk)) < 2*p; ++K) {
881: wk = 0.5*h*PETSC_PI*PetscCoshReal(K*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(K*h)));
882: }
883: PetscQuadratureCreate(PETSC_COMM_SELF, q);
884: PetscQuadratureSetOrder(*q, 2*K+1);
885: npoints = 2*K-1;
886: PetscMalloc1(npoints*dim, &x);
887: PetscMalloc1(npoints, &w);
888: /* Center term */
889: x[0] = beta;
890: w[0] = 0.5*alpha*PETSC_PI;
891: for (k = 1; k < K; ++k) {
892: wk = 0.5*alpha*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
893: xk = PetscTanhReal(0.5*PETSC_PI*PetscSinhReal(k*h));
894: x[2*k-1] = -alpha*xk+beta;
895: w[2*k-1] = wk;
896: x[2*k+0] = alpha*xk+beta;
897: w[2*k+0] = wk;
898: }
899: PetscQuadratureSetData(*q, dim, 1, npoints, x, w);
900: return(0);
901: }
903: PetscErrorCode PetscDTTanhSinhIntegrate(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
904: {
905: const PetscInt p = 16; /* Digits of precision in the evaluation */
906: const PetscReal alpha = (b-a)/2.; /* Half-width of the integration interval */
907: const PetscReal beta = (b+a)/2.; /* Center of the integration interval */
908: PetscReal h = 1.0; /* Step size, length between x_k */
909: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
910: PetscReal osum = 0.0; /* Integral on last level */
911: PetscReal psum = 0.0; /* Integral on the level before the last level */
912: PetscReal sum; /* Integral on current level */
913: PetscReal yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
914: PetscReal lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
915: PetscReal wk; /* Quadrature weight at x_k */
916: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
917: PetscInt d; /* Digits of precision in the integral */
920: if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
921: /* Center term */
922: func(beta, &lval);
923: sum = 0.5*alpha*PETSC_PI*lval;
924: /* */
925: do {
926: PetscReal lterm, rterm, maxTerm = 0.0, d1, d2, d3, d4;
927: PetscInt k = 1;
929: ++l;
930: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
931: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
932: psum = osum;
933: osum = sum;
934: h *= 0.5;
935: sum *= 0.5;
936: do {
937: wk = 0.5*h*PETSC_PI*PetscCoshReal(k*h)/PetscSqr(PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
938: yk = 1.0/(PetscExpReal(0.5*PETSC_PI*PetscSinhReal(k*h)) * PetscCoshReal(0.5*PETSC_PI*PetscSinhReal(k*h)));
939: lx = -alpha*(1.0 - yk)+beta;
940: rx = alpha*(1.0 - yk)+beta;
941: func(lx, &lval);
942: func(rx, &rval);
943: lterm = alpha*wk*lval;
944: maxTerm = PetscMax(PetscAbsReal(lterm), maxTerm);
945: sum += lterm;
946: rterm = alpha*wk*rval;
947: maxTerm = PetscMax(PetscAbsReal(rterm), maxTerm);
948: sum += rterm;
949: ++k;
950: /* Only need to evaluate every other point on refined levels */
951: if (l != 1) ++k;
952: } while (PetscAbsReal(PetscLog10Real(wk)) < p); /* Only need to evaluate sum until weights are < 32 digits of precision */
954: d1 = PetscLog10Real(PetscAbsReal(sum - osum));
955: d2 = PetscLog10Real(PetscAbsReal(sum - psum));
956: d3 = PetscLog10Real(maxTerm) - p;
957: if (PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)) == 0.0) d4 = 0.0;
958: else d4 = PetscLog10Real(PetscMax(PetscAbsReal(lterm), PetscAbsReal(rterm)));
959: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
960: } while (d < digits && l < 12);
961: *sol = sum;
963: return(0);
964: }
966: #if defined(PETSC_HAVE_MPFR)
967: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
968: {
969: const PetscInt safetyFactor = 2; /* Calculate abcissa until 2*p digits */
970: PetscInt l = 0; /* Level of refinement, h = 2^{-l} */
971: mpfr_t alpha; /* Half-width of the integration interval */
972: mpfr_t beta; /* Center of the integration interval */
973: mpfr_t h; /* Step size, length between x_k */
974: mpfr_t osum; /* Integral on last level */
975: mpfr_t psum; /* Integral on the level before the last level */
976: mpfr_t sum; /* Integral on current level */
977: mpfr_t yk; /* Quadrature point 1 - x_k on reference domain [-1, 1] */
978: mpfr_t lx, rx; /* Quadrature points to the left and right of 0 on the real domain [a, b] */
979: mpfr_t wk; /* Quadrature weight at x_k */
980: PetscReal lval, rval; /* Terms in the quadature sum to the left and right of 0 */
981: PetscInt d; /* Digits of precision in the integral */
982: mpfr_t pi2, kh, msinh, mcosh, maxTerm, curTerm, tmp;
985: if (digits <= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must give a positive number of significant digits");
986: /* Create high precision storage */
987: 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);
988: /* Initialization */
989: mpfr_set_d(alpha, 0.5*(b-a), MPFR_RNDN);
990: mpfr_set_d(beta, 0.5*(b+a), MPFR_RNDN);
991: mpfr_set_d(osum, 0.0, MPFR_RNDN);
992: mpfr_set_d(psum, 0.0, MPFR_RNDN);
993: mpfr_set_d(h, 1.0, MPFR_RNDN);
994: mpfr_const_pi(pi2, MPFR_RNDN);
995: mpfr_mul_d(pi2, pi2, 0.5, MPFR_RNDN);
996: /* Center term */
997: func(0.5*(b+a), &lval);
998: mpfr_set(sum, pi2, MPFR_RNDN);
999: mpfr_mul(sum, sum, alpha, MPFR_RNDN);
1000: mpfr_mul_d(sum, sum, lval, MPFR_RNDN);
1001: /* */
1002: do {
1003: PetscReal d1, d2, d3, d4;
1004: PetscInt k = 1;
1006: ++l;
1007: mpfr_set_d(maxTerm, 0.0, MPFR_RNDN);
1008: /* PetscPrintf(PETSC_COMM_SELF, "LEVEL %D sum: %15.15f\n", l, sum); */
1009: /* At each level of refinement, h --> h/2 and sum --> sum/2 */
1010: mpfr_set(psum, osum, MPFR_RNDN);
1011: mpfr_set(osum, sum, MPFR_RNDN);
1012: mpfr_mul_d(h, h, 0.5, MPFR_RNDN);
1013: mpfr_mul_d(sum, sum, 0.5, MPFR_RNDN);
1014: do {
1015: mpfr_set_si(kh, k, MPFR_RNDN);
1016: mpfr_mul(kh, kh, h, MPFR_RNDN);
1017: /* Weight */
1018: mpfr_set(wk, h, MPFR_RNDN);
1019: mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
1020: mpfr_mul(msinh, msinh, pi2, MPFR_RNDN);
1021: mpfr_mul(mcosh, mcosh, pi2, MPFR_RNDN);
1022: mpfr_cosh(tmp, msinh, MPFR_RNDN);
1023: mpfr_sqr(tmp, tmp, MPFR_RNDN);
1024: mpfr_mul(wk, wk, mcosh, MPFR_RNDN);
1025: mpfr_div(wk, wk, tmp, MPFR_RNDN);
1026: /* Abscissa */
1027: mpfr_set_d(yk, 1.0, MPFR_RNDZ);
1028: mpfr_cosh(tmp, msinh, MPFR_RNDN);
1029: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1030: mpfr_exp(tmp, msinh, MPFR_RNDN);
1031: mpfr_div(yk, yk, tmp, MPFR_RNDZ);
1032: /* Quadrature points */
1033: mpfr_sub_d(lx, yk, 1.0, MPFR_RNDZ);
1034: mpfr_mul(lx, lx, alpha, MPFR_RNDU);
1035: mpfr_add(lx, lx, beta, MPFR_RNDU);
1036: mpfr_d_sub(rx, 1.0, yk, MPFR_RNDZ);
1037: mpfr_mul(rx, rx, alpha, MPFR_RNDD);
1038: mpfr_add(rx, rx, beta, MPFR_RNDD);
1039: /* Evaluation */
1040: func(mpfr_get_d(lx, MPFR_RNDU), &lval);
1041: func(mpfr_get_d(rx, MPFR_RNDD), &rval);
1042: /* Update */
1043: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1044: mpfr_mul_d(tmp, tmp, lval, MPFR_RNDN);
1045: mpfr_add(sum, sum, tmp, MPFR_RNDN);
1046: mpfr_abs(tmp, tmp, MPFR_RNDN);
1047: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1048: mpfr_set(curTerm, tmp, MPFR_RNDN);
1049: mpfr_mul(tmp, wk, alpha, MPFR_RNDN);
1050: mpfr_mul_d(tmp, tmp, rval, MPFR_RNDN);
1051: mpfr_add(sum, sum, tmp, MPFR_RNDN);
1052: mpfr_abs(tmp, tmp, MPFR_RNDN);
1053: mpfr_max(maxTerm, maxTerm, tmp, MPFR_RNDN);
1054: mpfr_max(curTerm, curTerm, tmp, MPFR_RNDN);
1055: ++k;
1056: /* Only need to evaluate every other point on refined levels */
1057: if (l != 1) ++k;
1058: mpfr_log10(tmp, wk, MPFR_RNDN);
1059: mpfr_abs(tmp, tmp, MPFR_RNDN);
1060: } while (mpfr_get_d(tmp, MPFR_RNDN) < safetyFactor*digits); /* Only need to evaluate sum until weights are < 32 digits of precision */
1061: mpfr_sub(tmp, sum, osum, MPFR_RNDN);
1062: mpfr_abs(tmp, tmp, MPFR_RNDN);
1063: mpfr_log10(tmp, tmp, MPFR_RNDN);
1064: d1 = mpfr_get_d(tmp, MPFR_RNDN);
1065: mpfr_sub(tmp, sum, psum, MPFR_RNDN);
1066: mpfr_abs(tmp, tmp, MPFR_RNDN);
1067: mpfr_log10(tmp, tmp, MPFR_RNDN);
1068: d2 = mpfr_get_d(tmp, MPFR_RNDN);
1069: mpfr_log10(tmp, maxTerm, MPFR_RNDN);
1070: d3 = mpfr_get_d(tmp, MPFR_RNDN) - digits;
1071: mpfr_log10(tmp, curTerm, MPFR_RNDN);
1072: d4 = mpfr_get_d(tmp, MPFR_RNDN);
1073: d = PetscAbsInt(PetscMin(0, PetscMax(PetscMax(PetscMax(PetscSqr(d1)/d2, 2*d1), d3), d4)));
1074: } while (d < digits && l < 8);
1075: *sol = mpfr_get_d(sum, MPFR_RNDN);
1076: /* Cleanup */
1077: mpfr_clears(alpha, beta, h, sum, osum, psum, yk, wk, lx, rx, tmp, maxTerm, curTerm, pi2, kh, msinh, mcosh, NULL);
1078: return(0);
1079: }
1080: #else
1082: PetscErrorCode PetscDTTanhSinhIntegrateMPFR(void (*func)(PetscReal, PetscReal *), PetscReal a, PetscReal b, PetscInt digits, PetscReal *sol)
1083: {
1084: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "This method will not work without MPFR. Reconfigure using --download-mpfr --download-gmp");
1085: }
1086: #endif
1088: /* Overwrites A. Can only handle full-rank problems with m>=n
1089: * A in column-major format
1090: * Ainv in row-major format
1091: * tau has length m
1092: * worksize must be >= max(1,n)
1093: */
1094: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1095: {
1097: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
1098: PetscScalar *A,*Ainv,*R,*Q,Alpha;
1101: #if defined(PETSC_USE_COMPLEX)
1102: {
1103: PetscInt i,j;
1104: PetscMalloc2(m*n,&A,m*n,&Ainv);
1105: for (j=0; j<n; j++) {
1106: for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
1107: }
1108: mstride = m;
1109: }
1110: #else
1111: A = A_in;
1112: Ainv = Ainv_out;
1113: #endif
1115: PetscBLASIntCast(m,&M);
1116: PetscBLASIntCast(n,&N);
1117: PetscBLASIntCast(mstride,&lda);
1118: PetscBLASIntCast(worksize,&ldwork);
1119: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1120: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info));
1121: PetscFPTrapPop();
1122: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1123: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
1125: /* Extract an explicit representation of Q */
1126: Q = Ainv;
1127: PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1128: K = N; /* full rank */
1129: PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info));
1130: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
1132: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1133: Alpha = 1.0;
1134: ldb = lda;
1135: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb));
1136: /* Ainv is Q, overwritten with inverse */
1138: #if defined(PETSC_USE_COMPLEX)
1139: {
1140: PetscInt i;
1141: for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
1142: PetscFree2(A,Ainv);
1143: }
1144: #endif
1145: return(0);
1146: }
1148: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
1149: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
1150: {
1152: PetscReal *Bv;
1153: PetscInt i,j;
1156: PetscMalloc1((ninterval+1)*ndegree,&Bv);
1157: /* Point evaluation of L_p on all the source vertices */
1158: PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
1159: /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
1160: for (i=0; i<ninterval; i++) {
1161: for (j=0; j<ndegree; j++) {
1162: if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1163: else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
1164: }
1165: }
1166: PetscFree(Bv);
1167: return(0);
1168: }
1170: /*@
1171: PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
1173: Not Collective
1175: Input Arguments:
1176: + degree - degree of reconstruction polynomial
1177: . nsource - number of source intervals
1178: . sourcex - sorted coordinates of source cell boundaries (length nsource+1)
1179: . ntarget - number of target intervals
1180: - targetx - sorted coordinates of target cell boundaries (length ntarget+1)
1182: Output Arguments:
1183: . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
1185: Level: advanced
1187: .seealso: PetscDTLegendreEval()
1188: @*/
1189: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
1190: {
1192: PetscInt i,j,k,*bdegrees,worksize;
1193: PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
1194: PetscScalar *tau,*work;
1200: 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);
1201: #if defined(PETSC_USE_DEBUG)
1202: for (i=0; i<nsource; i++) {
1203: 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]);
1204: }
1205: for (i=0; i<ntarget; i++) {
1206: 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]);
1207: }
1208: #endif
1209: xmin = PetscMin(sourcex[0],targetx[0]);
1210: xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
1211: center = (xmin + xmax)/2;
1212: hscale = (xmax - xmin)/2;
1213: worksize = nsource;
1214: PetscMalloc4(degree+1,&bdegrees,nsource+1,&sourcey,nsource*(degree+1),&Bsource,worksize,&work);
1215: PetscMalloc4(nsource,&tau,nsource*(degree+1),&Bsinv,ntarget+1,&targety,ntarget*(degree+1),&Btarget);
1216: for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
1217: for (i=0; i<=degree; i++) bdegrees[i] = i+1;
1218: PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
1219: PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
1220: for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
1221: PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
1222: for (i=0; i<ntarget; i++) {
1223: PetscReal rowsum = 0;
1224: for (j=0; j<nsource; j++) {
1225: PetscReal sum = 0;
1226: for (k=0; k<degree+1; k++) {
1227: sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
1228: }
1229: R[i*nsource+j] = sum;
1230: rowsum += sum;
1231: }
1232: for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
1233: }
1234: PetscFree4(bdegrees,sourcey,Bsource,work);
1235: PetscFree4(tau,Bsinv,targety,Btarget);
1236: return(0);
1237: }