1 #include <petscsys.h> 2 #include <petscblaslapack.h> 3 4 /*@C 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 PetscInt i, j, k; 25 PetscErrorCode ierr; 26 27 PetscFunctionBegin; 28 *slope = *intercept = 0.0; 29 ierr = PetscMalloc2(n*2, &X, n*2, &Y);CHKERRQ(ierr); 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]; 48 49 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 50 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&two, &two, H, &two, ipiv, &info)); 51 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&two, H, &two, ipiv, work, &two, &info)); 52 ierr = PetscFPTrapPop();CHKERRQ(ierr); 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 ierr = PetscFree2(X, Y);CHKERRQ(ierr); 71 *intercept = beta[0]; 72 *slope = beta[1]; 73 PetscFunctionReturn(0); 74 } 75