Actual source code: mathfit.c
1: #include <petsc/private/petscimpl.h>
2: #include <petscblaslapack.h>
4: /*@C
5: PetscLinearRegression - Gives the best least-squares linear fit to some x-y data points
7: Input Parameters:
8: + n - The number of points
9: . x - The x-values
10: - y - The y-values
12: Output Parameters:
13: + slope - The slope of the best-fit line
14: - intercept - The y-intercept of the best-fit line
16: Level: intermediate
18: .seealso: `PetscConvEstGetConvRate()`
19: @*/
20: PetscErrorCode PetscLinearRegression(PetscInt n, const PetscReal x[], const PetscReal y[], PetscReal *slope, PetscReal *intercept)
21: {
22: PetscScalar H[4];
23: PetscReal *X, *Y, beta[2];
25: if (n) {
28: }
31: PetscMalloc2(n * 2, &X, n * 2, &Y);
32: for (PetscInt k = 0; k < n; ++k) {
33: /* X[n,2] = [1, x] */
34: X[k * 2 + 0] = 1.0;
35: X[k * 2 + 1] = x[k];
36: }
37: /* H = X^T X */
38: for (PetscInt i = 0; i < 2; ++i) {
39: for (PetscInt j = 0; j < 2; ++j) {
40: H[i * 2 + j] = 0.0;
41: for (PetscInt k = 0; k < n; ++k) H[i * 2 + j] += X[k * 2 + i] * X[k * 2 + j];
42: }
43: }
44: /* H = (X^T X)^{-1} */
45: {
46: PetscBLASInt two = 2, ipiv[2], info;
47: PetscScalar work[2];
49: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
50: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&two, &two, H, &two, ipiv, &info));
51: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&two, H, &two, ipiv, work, &two, &info));
52: PetscFPTrapPop();
53: }
54: /* Y = H X^T */
55: for (PetscInt i = 0; i < 2; ++i) {
56: for (PetscInt k = 0; k < n; ++k) {
57: Y[i * n + k] = 0.0;
58: for (PetscInt j = 0; j < 2; ++j) Y[i * n + k] += PetscRealPart(H[i * 2 + j]) * X[k * 2 + j];
59: }
60: }
61: /* beta = Y error = [y-intercept, slope] */
62: for (PetscInt i = 0; i < 2; ++i) {
63: beta[i] = 0.0;
64: for (PetscInt k = 0; k < n; ++k) beta[i] += Y[i * n + k] * y[k];
65: }
66: PetscFree2(X, Y);
67: *intercept = beta[0];
68: *slope = beta[1];
69: return 0;
70: }