Actual source code: mathfit.c
1: #include <petsc/private/petscimpl.h>
2: #include <petscblaslapack.h>
4: /*@
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: PetscFunctionBegin;
26: if (n) {
27: PetscAssertPointer(x, 2);
28: PetscAssertPointer(y, 3);
29: }
30: PetscAssertPointer(slope, 4);
31: PetscAssertPointer(intercept, 5);
32: PetscCall(PetscMalloc2(n * 2, &X, n * 2, &Y));
33: for (PetscInt k = 0; k < n; ++k) {
34: /* X[n,2] = [1, x] */
35: X[k * 2 + 0] = 1.0;
36: X[k * 2 + 1] = x[k];
37: }
38: /* H = X^T X */
39: for (PetscInt i = 0; i < 2; ++i) {
40: for (PetscInt j = 0; j < 2; ++j) {
41: H[i * 2 + j] = 0.0;
42: for (PetscInt k = 0; k < n; ++k) H[i * 2 + j] += X[k * 2 + i] * X[k * 2 + j];
43: }
44: }
45: /* H = (X^T X)^{-1} */
46: {
47: PetscBLASInt two = 2, ipiv[2], info;
48: PetscScalar work[2];
50: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
51: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&two, &two, H, &two, ipiv, &info));
52: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&two, H, &two, ipiv, work, &two, &info));
53: PetscCall(PetscFPTrapPop());
54: }
55: /* Y = H X^T */
56: for (PetscInt i = 0; i < 2; ++i) {
57: for (PetscInt k = 0; k < n; ++k) {
58: Y[i * n + k] = 0.0;
59: for (PetscInt j = 0; j < 2; ++j) Y[i * n + k] += PetscRealPart(H[i * 2 + j]) * X[k * 2 + j];
60: }
61: }
62: /* beta = Y error = [y-intercept, slope] */
63: for (PetscInt i = 0; i < 2; ++i) {
64: beta[i] = 0.0;
65: for (PetscInt k = 0; k < n; ++k) beta[i] += Y[i * n + k] * y[k];
66: }
67: PetscCall(PetscFree2(X, Y));
68: *intercept = beta[0];
69: *slope = beta[1];
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }