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:   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: }