|
@@ -25,6 +25,11 @@
|
|
|
#include <gsl/gsl_multifit.h>
|
|
|
#endif //TESTGSL
|
|
|
|
|
|
+// From additional-lapack-func
|
|
|
+#ifdef MLR_MODEL
|
|
|
+#include "mindgels.h"
|
|
|
+#endif //MLR_MODEL
|
|
|
+
|
|
|
typedef struct { int h, w; double *x;} matrix_t, *matrix;
|
|
|
|
|
|
static long count_file_lines(FILE *f)
|
|
@@ -133,9 +138,90 @@ void gsl_multiple_reg_coeff(double *mpar, double *my, long n, unsigned ncoeff, u
|
|
|
}
|
|
|
#endif //TESTGSL
|
|
|
|
|
|
+#ifdef MLR_MODEL
|
|
|
void dgels_multiple_reg_coeff(double *mpar, double *my, long nn, unsigned ncoeff, unsigned nparameters, double *coeff, unsigned **combinations)
|
|
|
{
|
|
|
-
|
|
|
+ /* Arguments */
|
|
|
+/* ========= */
|
|
|
+
|
|
|
+/* TRANS (input) CHARACTER*1 */
|
|
|
+/* = 'N': the linear system involves A; */
|
|
|
+/* = 'T': the linear system involves A**T. */
|
|
|
+
|
|
|
+/* M (input) INTEGER */
|
|
|
+/* The number of rows of the matrix A. M >= 0. */
|
|
|
+
|
|
|
+/* N (input) INTEGER */
|
|
|
+/* The number of columns of the matrix A. N >= 0. */
|
|
|
+
|
|
|
+/* NRHS (input) INTEGER */
|
|
|
+/* The number of right hand sides, i.e., the number of */
|
|
|
+/* columns of the matrices B and X. NRHS >=0. */
|
|
|
+
|
|
|
+/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
|
|
|
+/* On entry, the M-by-N matrix A. */
|
|
|
+/* On exit, */
|
|
|
+/* if M >= N, A is overwritten by details of its QR */
|
|
|
+/* factorization as returned by DGEQRF; */
|
|
|
+/* if M < N, A is overwritten by details of its LQ */
|
|
|
+/* factorization as returned by DGELQF. */
|
|
|
+
|
|
|
+/* LDA (input) INTEGER */
|
|
|
+/* The leading dimension of the array A. LDA >= max(1,M). */
|
|
|
+
|
|
|
+/* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) */
|
|
|
+/* On entry, the matrix B of right hand side vectors, stored */
|
|
|
+/* columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS */
|
|
|
+/* if TRANS = 'T'. */
|
|
|
+/* On exit, if INFO = 0, B is overwritten by the solution */
|
|
|
+/* vectors, stored columnwise: */
|
|
|
+/* if TRANS = 'N' and m >= n, rows 1 to n of B contain the least */
|
|
|
+/* squares solution vectors; the residual sum of squares for the */
|
|
|
+/* solution in each column is given by the sum of squares of */
|
|
|
+/* elements N+1 to M in that column; */
|
|
|
+/* if TRANS = 'N' and m < n, rows 1 to N of B contain the */
|
|
|
+/* minimum norm solution vectors; */
|
|
|
+/* if TRANS = 'T' and m >= n, rows 1 to M of B contain the */
|
|
|
+/* minimum norm solution vectors; */
|
|
|
+/* if TRANS = 'T' and m < n, rows 1 to M of B contain the */
|
|
|
+/* least squares solution vectors; the residual sum of squares */
|
|
|
+/* for the solution in each column is given by the sum of */
|
|
|
+/* squares of elements M+1 to N in that column. */
|
|
|
+
|
|
|
+/* LDB (input) INTEGER */
|
|
|
+/* The leading dimension of the array B. LDB >= MAX(1,M,N). */
|
|
|
+
|
|
|
+/* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
|
|
|
+/* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
|
|
|
+
|
|
|
+/* LWORK (input) INTEGER */
|
|
|
+/* The dimension of the array WORK. */
|
|
|
+/* LWORK >= max( 1, MN + max( MN, NRHS ) ). */
|
|
|
+/* For optimal performance, */
|
|
|
+/* LWORK >= max( 1, MN + max( MN, NRHS )*NB ). */
|
|
|
+/* where MN = min(M,N) and NB is the optimum block size. */
|
|
|
+
|
|
|
+/* If LWORK = -1, then a workspace query is assumed; the routine */
|
|
|
+/* only calculates the optimal size of the WORK array, returns */
|
|
|
+/* this value as the first entry of the WORK array, and no error */
|
|
|
+/* message related to LWORK is issued by XERBLA. */
|
|
|
+
|
|
|
+/* INFO (output) INTEGER */
|
|
|
+/* = 0: successful exit */
|
|
|
+/* < 0: if INFO = -i, the i-th argument had an illegal value */
|
|
|
+/* > 0: if INFO = i, the i-th diagonal element of the */
|
|
|
+/* triangular factor of A is zero, so that A does not have */
|
|
|
+/* full rank; the least squares solution could not be */
|
|
|
+/* computed. */
|
|
|
+
|
|
|
+/* ===================================================================== */
|
|
|
+
|
|
|
+ if(nn <= ncoeff)
|
|
|
+ {
|
|
|
+ printf("\nERROR: This function is not intended for the use when number of parameters is larger than the number of observations. Check how your matrices A and B were allocated or simply add more benchmarks.\n");
|
|
|
+ exit(1);
|
|
|
+ }
|
|
|
+
|
|
|
char trans = 'N';
|
|
|
integer m = nn;
|
|
|
integer n = ncoeff;
|
|
@@ -149,12 +235,14 @@ void dgels_multiple_reg_coeff(double *mpar, double *my, long nn, unsigned ncoeff
|
|
|
Y[i] = my[i];
|
|
|
X[i*n] = 1.;
|
|
|
for (int j=1; j < n; j++)
|
|
|
+ {
|
|
|
coefficient = 1.;
|
|
|
for(int k=0; k < nparameters; k++)
|
|
|
{
|
|
|
coefficient *= pow(mpar[i*nparameters+k],combinations[j-1][k]);
|
|
|
}
|
|
|
X[i*n+j] = coefficient;
|
|
|
+ }
|
|
|
}
|
|
|
|
|
|
integer lda = m;
|
|
@@ -170,7 +258,7 @@ void dgels_multiple_reg_coeff(double *mpar, double *my, long nn, unsigned ncoeff
|
|
|
/* Check for the full rank */
|
|
|
if( info != 0 )
|
|
|
{
|
|
|
- printf( "Problems with DGELS; info=%i\n");
|
|
|
+ printf( "Problems with DGELS; info=%ld\n", info);
|
|
|
exit(1);
|
|
|
}
|
|
|
|
|
@@ -181,6 +269,7 @@ void dgels_multiple_reg_coeff(double *mpar, double *my, long nn, unsigned ncoeff
|
|
|
free(Y);
|
|
|
free(work);
|
|
|
}
|
|
|
+#endif //MLR_MODEL
|
|
|
|
|
|
int _starpu_multiple_regression(struct starpu_perfmodel_history_list *ptr, double *coeff, unsigned ncoeff, unsigned nparameters, unsigned **combinations, char *codelet_name)
|
|
|
{
|
|
@@ -221,9 +310,10 @@ int _starpu_multiple_regression(struct starpu_perfmodel_history_list *ptr, doubl
|
|
|
// Computing coefficients using multiple linear regression
|
|
|
#ifdef TESTGSL
|
|
|
gsl_multiple_reg_coeff(mpar, my, n, ncoeff, nparameters, coeff, combinations);
|
|
|
-#elseif
|
|
|
- dgels_multiple_reg_coeff(mpar, my, n, ncoeff, nparameters, coeff, combinations);
|
|
|
#endif
|
|
|
+#ifdef MLR_MODEL
|
|
|
+ dgels_multiple_reg_coeff(mpar, my, n, ncoeff, nparameters, coeff, combinations);
|
|
|
+#endif //MLR_MODEL
|
|
|
|
|
|
// Preparing new output calibration file
|
|
|
if (calibrate==2)
|