1 #include <petsc/private/petscimpl.h>
2 #include <petscblaslapack.h>
3
4 /*@
5 PetscLinearRegression - Gives the best least-squares linear fit to some x-y data points
6
7 Input Parameters:
8 + n - The number of points
9 . x - The x-values
10 - y - The y-values
11
12 Output Parameters:
13 + slope - The slope of the best-fit line
14 - intercept - The y-intercept of the best-fit line
15
16 Level: intermediate
17
18 .seealso: `PetscConvEstGetConvRate()`
19 @*/
PetscLinearRegression(PetscInt n,const PetscReal x[],const PetscReal y[],PetscReal * slope,PetscReal * intercept)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];
24
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];
49
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 }
72