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 @*/ 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