Actual source code: dtprob.c
1: #include <petscdt.h>
3: #include <petscvec.h>
4: #include <petscdraw.h>
5: #include <petsc/private/petscimpl.h>
7: const char *const DTProbDensityTypes[] = {"constant", "gaussian", "maxwell_boltzmann", "DTProb Density", "DTPROB_DENSITY_", NULL};
9: /*@
10: PetscPDFMaxwellBoltzmann1D - PDF for the Maxwell-Boltzmann distribution in 1D
12: Not collective
14: Input Parameter:
15: . x - Speed in [0, \infty]
17: Output Parameter:
18: . p - The probability density at x
20: Level: beginner
22: .seealso: `PetscCDFMaxwellBoltzmann1D`
23: @*/
24: PetscErrorCode PetscPDFMaxwellBoltzmann1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
25: {
26: p[0] = PetscSqrtReal(2. / PETSC_PI) * PetscExpReal(-0.5 * PetscSqr(x[0]));
27: return 0;
28: }
30: /*@
31: PetscCDFMaxwellBoltzmann1D - CDF for the Maxwell-Boltzmann distribution in 1D
33: Not collective
35: Input Parameter:
36: . x - Speed in [0, \infty]
38: Output Parameter:
39: . p - The probability density at x
41: Level: beginner
43: .seealso: `PetscPDFMaxwellBoltzmann1D`
44: @*/
45: PetscErrorCode PetscCDFMaxwellBoltzmann1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
46: {
47: p[0] = PetscErfReal(x[0] / PETSC_SQRT2);
48: return 0;
49: }
51: /*@
52: PetscPDFMaxwellBoltzmann2D - PDF for the Maxwell-Boltzmann distribution in 2D
54: Not collective
56: Input Parameter:
57: . x - Speed in [0, \infty]
59: Output Parameter:
60: . p - The probability density at x
62: Level: beginner
64: .seealso: `PetscCDFMaxwellBoltzmann2D`
65: @*/
66: PetscErrorCode PetscPDFMaxwellBoltzmann2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
67: {
68: p[0] = x[0] * PetscExpReal(-0.5 * PetscSqr(x[0]));
69: return 0;
70: }
72: /*@
73: PetscCDFMaxwellBoltzmann2D - CDF for the Maxwell-Boltzmann distribution in 2D
75: Not collective
77: Input Parameter:
78: . x - Speed in [0, \infty]
80: Output Parameter:
81: . p - The probability density at x
83: Level: beginner
85: .seealso: `PetscPDFMaxwellBoltzmann2D`
86: @*/
87: PetscErrorCode PetscCDFMaxwellBoltzmann2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
88: {
89: p[0] = 1. - PetscExpReal(-0.5 * PetscSqr(x[0]));
90: return 0;
91: }
93: /*@
94: PetscPDFMaxwellBoltzmann3D - PDF for the Maxwell-Boltzmann distribution in 3D
96: Not collective
98: Input Parameter:
99: . x - Speed in [0, \infty]
101: Output Parameter:
102: . p - The probability density at x
104: Level: beginner
106: .seealso: `PetscCDFMaxwellBoltzmann3D`
107: @*/
108: PetscErrorCode PetscPDFMaxwellBoltzmann3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
109: {
110: p[0] = PetscSqrtReal(2. / PETSC_PI) * PetscSqr(x[0]) * PetscExpReal(-0.5 * PetscSqr(x[0]));
111: return 0;
112: }
114: /*@
115: PetscCDFMaxwellBoltzmann3D - CDF for the Maxwell-Boltzmann distribution in 3D
117: Not collective
119: Input Parameter:
120: . x - Speed in [0, \infty]
122: Output Parameter:
123: . p - The probability density at x
125: Level: beginner
127: .seealso: `PetscPDFMaxwellBoltzmann3D`
128: @*/
129: PetscErrorCode PetscCDFMaxwellBoltzmann3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
130: {
131: p[0] = PetscErfReal(x[0] / PETSC_SQRT2) - PetscSqrtReal(2. / PETSC_PI) * x[0] * PetscExpReal(-0.5 * PetscSqr(x[0]));
132: return 0;
133: }
135: /*@
136: PetscPDFGaussian1D - PDF for the Gaussian distribution in 1D
138: Not collective
140: Input Parameter:
141: . x - Coordinate in [-\infty, \infty]
143: Output Parameter:
144: . p - The probability density at x
146: Level: beginner
148: .seealso: `PetscPDFMaxwellBoltzmann3D`
149: @*/
150: PetscErrorCode PetscPDFGaussian1D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
151: {
152: const PetscReal sigma = scale ? scale[0] : 1.;
153: p[0] = PetscSqrtReal(1. / (2. * PETSC_PI)) * PetscExpReal(-0.5 * PetscSqr(x[0] / sigma)) / sigma;
154: return 0;
155: }
157: PetscErrorCode PetscCDFGaussian1D(const PetscReal x[], const PetscReal scale[], PetscReal p[])
158: {
159: const PetscReal sigma = scale ? scale[0] : 1.;
160: p[0] = 0.5 * (1. + PetscErfReal(x[0] / PETSC_SQRT2 / sigma));
161: return 0;
162: }
164: /*@
165: PetscPDFSampleGaussian1D - Sample uniformly from a Gaussian distribution in 1D
167: Not collective
169: Input Parameters:
170: + p - A uniform variable on [0, 1]
171: - dummy - ignored
173: Output Parameter:
174: . x - Coordinate in [-\infty, \infty]
176: Level: beginner
178: Note: http://www.mimirgames.com/articles/programming/approximations-of-the-inverse-error-function/
179: https://stackoverflow.com/questions/27229371/inverse-error-function-in-c
180: @*/
181: PetscErrorCode PetscPDFSampleGaussian1D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
182: {
183: const PetscReal q = 2 * p[0] - 1.;
184: const PetscInt maxIter = 100;
185: PetscReal ck[100], r = 0.;
186: PetscInt k, m;
189: /* Transform input to [-1, 1] since the code below computes the inverse error function */
190: for (k = 0; k < maxIter; ++k) ck[k] = 0.;
191: ck[0] = 1;
192: r = ck[0] * (PetscSqrtReal(PETSC_PI) / 2.) * q;
193: for (k = 1; k < maxIter; ++k) {
194: const PetscReal temp = 2. * k + 1.;
196: for (m = 0; m <= k - 1; ++m) {
197: PetscReal denom = (m + 1.) * (2. * m + 1.);
199: ck[k] += (ck[m] * ck[k - 1 - m]) / denom;
200: }
201: r += (ck[k] / temp) * PetscPowReal((PetscSqrtReal(PETSC_PI) / 2.) * q, 2. * k + 1.);
202: }
203: /* Scale erfinv() by \sqrt{\pi/2} */
204: x[0] = PetscSqrtReal(PETSC_PI * 0.5) * r;
205: return 0;
206: }
208: /*@
209: PetscPDFGaussian2D - PDF for the Gaussian distribution in 2D
211: Not collective
213: Input Parameters:
214: + x - Coordinate in [-\infty, \infty]^2
215: - dummy - ignored
217: Output Parameter:
218: . p - The probability density at x
220: Level: beginner
222: .seealso: `PetscPDFSampleGaussian2D()`, `PetscPDFMaxwellBoltzmann3D()`
223: @*/
224: PetscErrorCode PetscPDFGaussian2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
225: {
226: p[0] = (1. / PETSC_PI) * PetscExpReal(-0.5 * (PetscSqr(x[0]) + PetscSqr(x[1])));
227: return 0;
228: }
230: /*@
231: PetscPDFSampleGaussian2D - Sample uniformly from a Gaussian distribution in 2D
233: Not collective
235: Input Parameters:
236: + p - A uniform variable on [0, 1]^2
237: - dummy - ignored
239: Output Parameter:
240: . x - Coordinate in [-\infty, \infty]^2
242: Level: beginner
244: Note: https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
246: .seealso: `PetscPDFGaussian2D()`, `PetscPDFMaxwellBoltzmann3D()`
247: @*/
248: PetscErrorCode PetscPDFSampleGaussian2D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
249: {
250: const PetscReal mag = PetscSqrtReal(-2.0 * PetscLogReal(p[0]));
251: x[0] = mag * PetscCosReal(2.0 * PETSC_PI * p[1]);
252: x[1] = mag * PetscSinReal(2.0 * PETSC_PI * p[1]);
253: return 0;
254: }
256: /*@
257: PetscPDFGaussian3D - PDF for the Gaussian distribution in 3D
259: Not collective
261: Input Parameters:
262: + x - Coordinate in [-\infty, \infty]^3
263: - dummy - ignored
265: Output Parameter:
266: . p - The probability density at x
268: Level: beginner
270: .seealso: `PetscPDFSampleGaussian3D()`, `PetscPDFMaxwellBoltzmann3D()`
271: @*/
272: PetscErrorCode PetscPDFGaussian3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
273: {
274: p[0] = (1. / PETSC_PI * PetscSqrtReal(PETSC_PI)) * PetscExpReal(-0.5 * (PetscSqr(x[0]) + PetscSqr(x[1]) + PetscSqr(x[2])));
275: return 0;
276: }
278: /*@
279: PetscPDFSampleGaussian3D - Sample uniformly from a Gaussian distribution in 3D
281: Not collective
283: Input Parameters:
284: + p - A uniform variable on [0, 1]^3
285: - dummy - ignored
287: Output Parameter:
288: . x - Coordinate in [-\infty, \infty]^3
290: Level: beginner
292: Note: https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
294: .seealso: `PetscPDFGaussian3D()`, `PetscPDFMaxwellBoltzmann3D()`
295: @*/
296: PetscErrorCode PetscPDFSampleGaussian3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
297: {
298: PetscPDFSampleGaussian1D(p, dummy, x);
299: PetscPDFSampleGaussian2D(&p[1], dummy, &x[1]);
300: return 0;
301: }
303: /*@
304: PetscPDFConstant1D - PDF for the uniform distribution in 1D
306: Not collective
308: Input Parameter:
309: . x - Coordinate in [-1, 1]
311: Output Parameter:
312: . p - The probability density at x
314: Level: beginner
316: .seealso: `PetscCDFConstant1D()`, `PetscPDFSampleConstant1D()`, `PetscPDFConstant2D()`, `PetscPDFConstant3D()`
317: @*/
318: PetscErrorCode PetscPDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
319: {
320: p[0] = x[0] > -1. && x[0] <= 1. ? 0.5 : 0.;
321: return 0;
322: }
324: /*@
325: PetscCDFConstant1D - CDF for the uniform distribution in 1D
327: Not collective
329: Input Parameter:
330: . x - Coordinate in [-1, 1]
332: Output Parameter:
333: . p - The cumulative probability at x
335: Level: beginner
337: .seealso: `PetscPDFConstant1D()`, `PetscPDFSampleConstant1D()`, `PetscCDFConstant2D()`, `PetscCDFConstant3D()`
338: @*/
339: PetscErrorCode PetscCDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
340: {
341: p[0] = x[0] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5 * (x[0] + 1.));
342: return 0;
343: }
345: /*@
346: PetscPDFSampleConstant1D - Sample uniformly from a uniform distribution on [-1, 1] in 1D
348: Not collective
350: Input Parameter:
351: . p - A uniform variable on [0, 1]
353: Output Parameter:
354: . x - Coordinate in [-1, 1]
356: Level: beginner
358: .seealso: `PetscPDFConstant1D()`, `PetscCDFConstant1D()`, `PetscPDFSampleConstant2D()`, `PetscPDFSampleConstant3D()`
359: @*/
360: PetscErrorCode PetscPDFSampleConstant1D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
361: {
362: x[0] = 2. * p[0] - 1.;
363: return 0;
364: }
366: /*@
367: PetscPDFConstant2D - PDF for the uniform distribution in 2D
369: Not collective
371: Input Parameter:
372: . x - Coordinate in [-1, 1] x [-1, 1]
374: Output Parameter:
375: . p - The probability density at x
377: Level: beginner
379: .seealso: `PetscCDFConstant2D()`, `PetscPDFSampleConstant2D()`, `PetscPDFConstant1D()`, `PetscPDFConstant3D()`
380: @*/
381: PetscErrorCode PetscPDFConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
382: {
383: p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. ? 0.25 : 0.;
384: return 0;
385: }
387: /*@
388: PetscCDFConstant2D - CDF for the uniform distribution in 2D
390: Not collective
392: Input Parameter:
393: . x - Coordinate in [-1, 1] x [-1, 1]
395: Output Parameter:
396: . p - The cumulative probability at x
398: Level: beginner
400: .seealso: `PetscPDFConstant2D()`, `PetscPDFSampleConstant2D()`, `PetscCDFConstant1D()`, `PetscCDFConstant3D()`
401: @*/
402: PetscErrorCode PetscCDFConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
403: {
404: p[0] = x[0] <= -1. || x[1] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5 * (x[0] + 1.)) * (x[1] > 1. ? 1. : 0.5 * (x[1] + 1.));
405: return 0;
406: }
408: /*@
409: PetscPDFSampleConstant2D - Sample uniformly from a uniform distribution on [-1, 1] x [-1, 1] in 2D
411: Not collective
413: Input Parameter:
414: . p - Two uniform variables on [0, 1]
416: Output Parameter:
417: . x - Coordinate in [-1, 1] x [-1, 1]
419: Level: beginner
421: .seealso: `PetscPDFConstant2D()`, `PetscCDFConstant2D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant3D()`
422: @*/
423: PetscErrorCode PetscPDFSampleConstant2D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
424: {
425: x[0] = 2. * p[0] - 1.;
426: x[1] = 2. * p[1] - 1.;
427: return 0;
428: }
430: /*@
431: PetscPDFConstant3D - PDF for the uniform distribution in 3D
433: Not collective
435: Input Parameter:
436: . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1]
438: Output Parameter:
439: . p - The probability density at x
441: Level: beginner
443: .seealso: `PetscCDFConstant3D()`, `PetscPDFSampleConstant3D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant2D()`
444: @*/
445: PetscErrorCode PetscPDFConstant3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
446: {
447: p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. && x[2] > -1. && x[2] <= 1. ? 0.125 : 0.;
448: return 0;
449: }
451: /*@
452: PetscCDFConstant3D - CDF for the uniform distribution in 3D
454: Not collective
456: Input Parameter:
457: . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1]
459: Output Parameter:
460: . p - The cumulative probability at x
462: Level: beginner
464: .seealso: `PetscPDFConstant3D()`, `PetscPDFSampleConstant3D()`, `PetscCDFConstant1D()`, `PetscCDFConstant2D()`
465: @*/
466: PetscErrorCode PetscCDFConstant3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
467: {
468: p[0] = x[0] <= -1. || x[1] <= -1. || x[2] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5 * (x[0] + 1.)) * (x[1] > 1. ? 1. : 0.5 * (x[1] + 1.)) * (x[2] > 1. ? 1. : 0.5 * (x[2] + 1.));
469: return 0;
470: }
472: /*@
473: PetscPDFSampleConstant3D - Sample uniformly from a uniform distribution on [-1, 1] x [-1, 1] in 3D
475: Not collective
477: Input Parameter:
478: . p - Three uniform variables on [0, 1]
480: Output Parameter:
481: . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1]
483: Level: beginner
485: .seealso: `PetscPDFConstant3D()`, `PetscCDFConstant3D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant2D()`
486: @*/
487: PetscErrorCode PetscPDFSampleConstant3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
488: {
489: x[0] = 2. * p[0] - 1.;
490: x[1] = 2. * p[1] - 1.;
491: x[2] = 2. * p[2] - 1.;
492: return 0;
493: }
495: /*@C
496: PetscProbCreateFromOptions - Return the probability distribution specified by the argumetns and options
498: Not collective
500: Input Parameters:
501: + dim - The dimension of sample points
502: . prefix - The options prefix, or NULL
503: - name - The option name for the probility distribution type
505: Output Parameters:
506: + pdf - The PDF of this type
507: . cdf - The CDF of this type
508: - sampler - The PDF sampler of this type
510: Level: intermediate
512: .seealso: `PetscPDFMaxwellBoltzmann1D()`, `PetscPDFGaussian1D()`, `PetscPDFConstant1D()`
513: @*/
514: PetscErrorCode PetscProbCreateFromOptions(PetscInt dim, const char prefix[], const char name[], PetscProbFunc *pdf, PetscProbFunc *cdf, PetscProbFunc *sampler)
515: {
516: DTProbDensityType den = DTPROB_DENSITY_GAUSSIAN;
518: PetscOptionsBegin(PETSC_COMM_SELF, prefix, "PetscProb Options", "DT");
519: PetscOptionsEnum(name, "Method to compute PDF <constant, gaussian>", "", DTProbDensityTypes, (PetscEnum)den, (PetscEnum *)&den, NULL);
520: PetscOptionsEnd();
522: if (pdf) {
524: *pdf = NULL;
525: }
526: if (cdf) {
528: *cdf = NULL;
529: }
530: if (sampler) {
532: *sampler = NULL;
533: }
534: switch (den) {
535: case DTPROB_DENSITY_CONSTANT:
536: switch (dim) {
537: case 1:
538: if (pdf) *pdf = PetscPDFConstant1D;
539: if (cdf) *cdf = PetscCDFConstant1D;
540: if (sampler) *sampler = PetscPDFSampleConstant1D;
541: break;
542: case 2:
543: if (pdf) *pdf = PetscPDFConstant2D;
544: if (cdf) *cdf = PetscCDFConstant2D;
545: if (sampler) *sampler = PetscPDFSampleConstant2D;
546: break;
547: case 3:
548: if (pdf) *pdf = PetscPDFConstant3D;
549: if (cdf) *cdf = PetscCDFConstant3D;
550: if (sampler) *sampler = PetscPDFSampleConstant3D;
551: break;
552: default:
553: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
554: }
555: break;
556: case DTPROB_DENSITY_GAUSSIAN:
557: switch (dim) {
558: case 1:
559: if (pdf) *pdf = PetscPDFGaussian1D;
560: if (cdf) *cdf = PetscCDFGaussian1D;
561: if (sampler) *sampler = PetscPDFSampleGaussian1D;
562: break;
563: case 2:
564: if (pdf) *pdf = PetscPDFGaussian2D;
565: if (sampler) *sampler = PetscPDFSampleGaussian2D;
566: break;
567: case 3:
568: if (pdf) *pdf = PetscPDFGaussian3D;
569: if (sampler) *sampler = PetscPDFSampleGaussian3D;
570: break;
571: default:
572: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
573: }
574: break;
575: case DTPROB_DENSITY_MAXWELL_BOLTZMANN:
576: switch (dim) {
577: case 1:
578: if (pdf) *pdf = PetscPDFMaxwellBoltzmann1D;
579: if (cdf) *cdf = PetscCDFMaxwellBoltzmann1D;
580: break;
581: case 2:
582: if (pdf) *pdf = PetscPDFMaxwellBoltzmann2D;
583: if (cdf) *cdf = PetscCDFMaxwellBoltzmann2D;
584: break;
585: case 3:
586: if (pdf) *pdf = PetscPDFMaxwellBoltzmann3D;
587: if (cdf) *cdf = PetscCDFMaxwellBoltzmann3D;
588: break;
589: default:
590: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
591: }
592: break;
593: default:
594: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Density type %s is not supported", DTProbDensityTypes[PetscMax(0, PetscMin(den, DTPROB_NUM_DENSITY))]);
595: }
596: return 0;
597: }
599: #ifdef PETSC_HAVE_KS
600: EXTERN_C_BEGIN
601: #include <KolmogorovSmirnovDist.h>
602: EXTERN_C_END
603: #endif
605: /*@C
606: PetscProbComputeKSStatistic - Compute the Kolmogorov-Smirnov statistic for the empirical distribution for an input vector, compared to an analytic CDF.
608: Collective on v
610: Input Parameters:
611: + v - The data vector, blocksize is the sample dimension
612: - cdf - The analytic CDF
614: Output Parameter:
615: . alpha - The KS statisic
617: Level: advanced
619: Note: The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is
620: $
621: $ D_n = \sup_x \left| F_n(x) - F(x) \right|
622: $
623: where $\sup_x$ is the supremum of the set of distances, and the empirical distribution
624: function $F_n(x)$ is discrete, and given by
625: $
626: $ F_n = # of samples <= x / n
627: $
628: The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep''
629: cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes
630: the largest absolute difference between the two distribution functions across all $x$ values.
632: .seealso: `PetscProbFunc`
633: @*/
634: PetscErrorCode PetscProbComputeKSStatistic(Vec v, PetscProbFunc cdf, PetscReal *alpha)
635: {
636: #if !defined(PETSC_HAVE_KS)
637: SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "No support for Kolmogorov-Smirnov test.\nReconfigure using --download-ks");
638: #else
639: PetscViewer viewer = NULL;
640: PetscViewerFormat format;
641: PetscDraw draw;
642: PetscDrawLG lg;
643: PetscDrawAxis axis;
644: const PetscScalar *a;
645: PetscReal *speed, Dn = PETSC_MIN_REAL;
646: PetscBool isascii = PETSC_FALSE, isdraw = PETSC_FALSE, flg;
647: PetscInt n, p, dim, d;
648: PetscMPIInt size;
649: const char *names[2] = {"Analytic", "Empirical"};
650: char title[PETSC_MAX_PATH_LEN];
651: PetscOptions options;
652: const char *prefix;
653: MPI_Comm comm;
655: PetscObjectGetComm((PetscObject)v, &comm);
656: PetscObjectGetOptionsPrefix((PetscObject)v, &prefix);
657: PetscObjectGetOptions((PetscObject)v, &options);
658: PetscOptionsGetViewer(comm, options, prefix, "-ks_monitor", &viewer, &format, &flg);
659: if (flg) {
660: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
661: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
662: }
663: if (isascii) {
664: PetscViewerPushFormat(viewer, format);
665: } else if (isdraw) {
666: PetscViewerPushFormat(viewer, format);
667: PetscViewerDrawGetDraw(viewer, 0, &draw);
668: PetscDrawLGCreate(draw, 2, &lg);
669: PetscDrawLGSetLegend(lg, names);
670: }
672: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
673: VecGetLocalSize(v, &n);
674: VecGetBlockSize(v, &dim);
675: n /= dim;
676: /* TODO Parallel algorithm is harder */
677: if (size == 1) {
678: PetscMalloc1(n, &speed);
679: VecGetArrayRead(v, &a);
680: for (p = 0; p < n; ++p) {
681: PetscReal mag = 0.;
683: for (d = 0; d < dim; ++d) mag += PetscSqr(PetscRealPart(a[p * dim + d]));
684: speed[p] = PetscSqrtReal(mag);
685: }
686: PetscSortReal(n, speed);
687: VecRestoreArrayRead(v, &a);
688: for (p = 0; p < n; ++p) {
689: const PetscReal x = speed[p], Fn = ((PetscReal)p) / n;
690: PetscReal F, vals[2];
692: cdf(&x, NULL, &F);
693: Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
694: if (isascii) PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn);
695: if (isdraw) {
696: vals[0] = F;
697: vals[1] = Fn;
698: PetscDrawLGAddCommonPoint(lg, x, vals);
699: }
700: }
701: if (speed[n - 1] < 6.) {
702: const PetscReal k = (PetscInt)(6. - speed[n - 1]) + 1, dx = (6. - speed[n - 1]) / k;
703: for (p = 0; p <= k; ++p) {
704: const PetscReal x = speed[n - 1] + p * dx, Fn = 1.0;
705: PetscReal F, vals[2];
707: cdf(&x, NULL, &F);
708: Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
709: if (isascii) PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn);
710: if (isdraw) {
711: vals[0] = F;
712: vals[1] = Fn;
713: PetscDrawLGAddCommonPoint(lg, x, vals);
714: }
715: }
716: }
717: PetscFree(speed);
718: }
719: if (isdraw) {
720: PetscDrawLGGetAxis(lg, &axis);
721: PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Kolmogorov-Smirnov Test (Dn %.2g)", Dn);
722: PetscDrawAxisSetLabels(axis, title, "x", "CDF(x)");
723: PetscDrawLGDraw(lg);
724: PetscDrawLGDestroy(&lg);
725: }
726: if (viewer) {
727: PetscViewerFlush(viewer);
728: PetscViewerPopFormat(viewer);
729: PetscViewerDestroy(&viewer);
730: }
731: *alpha = KSfbar((int)n, (double)Dn);
732: return 0;
733: #endif
734: }