1: #include <petscsys.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];
 24:   PetscInt       i, j, k;

 28:   *slope = *intercept = 0.0;
 29:   PetscMalloc2(n*2, &X, n*2, &Y);
 30:   for (k = 0; k < n; ++k) {
 31:     /* X[n,2] = [1, x] */
 32:     X[k*2+0] = 1.0;
 33:     X[k*2+1] = x[k];
 34:   }
 35:   /* H = X^T X */
 36:   for (i = 0; i < 2; ++i) {
 37:     for (j = 0; j < 2; ++j) {
 38:       H[i*2+j] = 0.0;
 39:       for (k = 0; k < n; ++k) {
 40:         H[i*2+j] += X[k*2+i] * X[k*2+j];
 41:       }
 42:     }
 43:   }
 44:   /* H = (X^T X)^{-1} */
 45:   {
 46:     PetscBLASInt two = 2, ipiv[2], info;
 47:     PetscScalar  work[2];

 49:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 50:     PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&two, &two, H, &two, ipiv, &info));
 51:     PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&two, H, &two, ipiv, work, &two, &info));
 52:     PetscFPTrapPop();
 53:   }
 54:     /* Y = H X^T */
 55:   for (i = 0; i < 2; ++i) {
 56:     for (k = 0; k < n; ++k) {
 57:       Y[i*n+k] = 0.0;
 58:       for (j = 0; j < 2; ++j) {
 59:         Y[i*n+k] += PetscRealPart(H[i*2+j]) * X[k*2+j];
 60:       }
 61:     }
 62:   }
 63:   /* beta = Y error = [y-intercept, slope] */
 64:   for (i = 0; i < 2; ++i) {
 65:     beta[i] = 0.0;
 66:     for (k = 0; k < n; ++k) {
 67:       beta[i] += Y[i*n+k] * y[k];
 68:     }
 69:   }
 70:   PetscFree2(X, Y);
 71:   *intercept = beta[0];
 72:   *slope     = beta[1];
 73:   return(0);
 74: }