Actual source code: mathfit.c

      1: #include <petsc/private/petscimpl.h>
      2: #include <petscblaslapack.h>

      4: /*@
      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: }