1 #include <petsc/private/petscimpl.h> 2 #include <petscblaslapack.h> 3 4 /*@C 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 @*/ 20 PetscErrorCode PetscLinearRegression(PetscInt n, const PetscReal x[], const PetscReal y[], PetscReal *slope, PetscReal *intercept) { 21 PetscScalar H[4]; 22 PetscReal *X, *Y, beta[2]; 23 24 PetscFunctionBegin; 25 if (n) { 26 PetscValidRealPointer(x, 2); 27 PetscValidRealPointer(y, 3); 28 } 29 PetscValidRealPointer(slope, 4); 30 PetscValidRealPointer(intercept, 5); 31 PetscCall(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]; 48 49 PetscCall(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 PetscCall(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 PetscCall(PetscFree2(X, Y)); 67 *intercept = beta[0]; 68 *slope = beta[1]; 69 PetscFunctionReturn(0); 70 } 71