Actual source code: mathfit.c

  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:   if (n) {
 28:   }
 31:   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];

 49:     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:     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:   PetscFree2(X, Y);
 67:   *intercept = beta[0];
 68:   *slope     = beta[1];
 69:   return 0;
 70: }