|
Packit |
67cb25 |
#define penalty2_N 8 /* 2*p */
|
|
Packit |
67cb25 |
#define penalty2_P 4
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double penalty2_x0[penalty2_P] = { 0.5, 0.5, 0.5, 0.5 };
|
|
Packit |
67cb25 |
static double penalty2_epsrel = 1.0e-12;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double penalty2_J[penalty2_N * penalty2_P];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
penalty2_checksol(const double x[], const double sumsq,
|
|
Packit |
67cb25 |
const double epsrel, const char *sname,
|
|
Packit |
67cb25 |
const char *pname)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double sumsq_exact = 9.37629300735544219e-06;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(sumsq, sumsq_exact, epsrel, "%s/%s sumsq",
|
|
Packit |
67cb25 |
sname, pname);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
(void)x; /* avoid unused parameter warning */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
penalty2_f (const gsl_vector * x, void *params, gsl_vector * f)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double alpha = 1.0e-5;
|
|
Packit |
67cb25 |
const double sqrt_alpha = sqrt(alpha);
|
|
Packit |
67cb25 |
double x1 = gsl_vector_get(x, 0);
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
double sum = penalty2_P * x1 * x1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set(f, 0, x1 - 0.2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* rows [2:p] */
|
|
Packit |
67cb25 |
for (i = 1; i < penalty2_P; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = gsl_vector_get(x, i);
|
|
Packit |
67cb25 |
double xim1 = gsl_vector_get(x, i - 1);
|
|
Packit |
67cb25 |
double yi = exp(0.1*(i + 1.0)) + exp(0.1*i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set(f, i, sqrt_alpha*(exp(0.1*xi) + exp(0.1*xim1) - yi));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
sum += (penalty2_P - i) * xi * xi;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* rows [p+1:2p-1] */
|
|
Packit |
67cb25 |
for (i = penalty2_P; i < penalty2_N - 1; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = gsl_vector_get(x, i - penalty2_P + 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set(f, i, sqrt_alpha*(exp(0.1*xi) - exp(-0.1)));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* row 2p */
|
|
Packit |
67cb25 |
gsl_vector_set(f, penalty2_N - 1, sum - 1.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
(void)params; /* avoid unused parameter warning */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
penalty2_df (CBLAS_TRANSPOSE_t TransJ, const gsl_vector * x,
|
|
Packit |
67cb25 |
const gsl_vector * u, void * params, gsl_vector * v,
|
|
Packit |
67cb25 |
gsl_matrix * JTJ)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view J = gsl_matrix_view_array(penalty2_J, penalty2_N, penalty2_P);
|
|
Packit |
67cb25 |
const double alpha = 1.0e-5;
|
|
Packit |
67cb25 |
const double sqrt_alpha = sqrt(alpha);
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < penalty2_P; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xj = gsl_vector_get(x, j);
|
|
Packit |
67cb25 |
double delta1j = (j == 0) ? 1.0 : 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* first and last rows */
|
|
Packit |
67cb25 |
gsl_matrix_set(&J.matrix, 0, j, delta1j);
|
|
Packit |
67cb25 |
gsl_matrix_set(&J.matrix, penalty2_N - 1, j, 2.0 * (penalty2_P - j) * xj);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* rows [2:p] */
|
|
Packit |
67cb25 |
for (i = 1; i < penalty2_P; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = gsl_vector_get(x, i);
|
|
Packit |
67cb25 |
double xim1 = gsl_vector_get(x, i - 1);
|
|
Packit |
67cb25 |
double Jij;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i == j)
|
|
Packit |
67cb25 |
Jij = exp(0.1 * xi);
|
|
Packit |
67cb25 |
else if (i - 1 == j)
|
|
Packit |
67cb25 |
Jij = exp(0.1 * xim1);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
Jij = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Jij *= 0.1 * sqrt_alpha;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set(&J.matrix, i, j, Jij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* rows [p+1:2p-1] */
|
|
Packit |
67cb25 |
for (i = penalty2_P; i < penalty2_N - 1; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = gsl_vector_get(x, i - penalty2_P + 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i - penalty2_P + 1 == j)
|
|
Packit |
67cb25 |
gsl_matrix_set(&J.matrix, i, j, 0.1 * sqrt_alpha * exp(0.1*xi));
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
gsl_matrix_set(&J.matrix, i, j, 0.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (v)
|
|
Packit |
67cb25 |
gsl_blas_dgemv(TransJ, 1.0, &J.matrix, u, 0.0, v);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (JTJ)
|
|
Packit |
67cb25 |
gsl_blas_dsyrk(CblasLower, CblasTrans, 1.0, &J.matrix, 0.0, JTJ);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
(void)params; /* avoid unused parameter warning */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
penalty2_fvv (const gsl_vector * x, const gsl_vector * v,
|
|
Packit |
67cb25 |
void *params, gsl_vector * fvv)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double alpha = 1.0e-5;
|
|
Packit |
67cb25 |
const double sqrt_alpha = sqrt(alpha);
|
|
Packit |
67cb25 |
double v1 = gsl_vector_get(v, 0);
|
|
Packit |
67cb25 |
double sum = penalty2_P * v1 * v1;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* first row */
|
|
Packit |
67cb25 |
gsl_vector_set(fvv, 0, 0.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* rows [2:p] */
|
|
Packit |
67cb25 |
for (i = 1; i < penalty2_P; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = gsl_vector_get(x, i);
|
|
Packit |
67cb25 |
double xim1 = gsl_vector_get(x, i - 1);
|
|
Packit |
67cb25 |
double vi = gsl_vector_get(v, i);
|
|
Packit |
67cb25 |
double vim1 = gsl_vector_get(v, i - 1);
|
|
Packit |
67cb25 |
double term1 = exp(xi / 10.0);
|
|
Packit |
67cb25 |
double term2 = exp(xim1 / 10.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set(fvv, i,
|
|
Packit |
67cb25 |
sqrt_alpha / 100.0 * (term1 * vi * vi + term2 * vim1 * vim1));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
sum += (penalty2_P - i) * vi * vi;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* last row */
|
|
Packit |
67cb25 |
gsl_vector_set(fvv, penalty2_N - 1, 2.0 * sum);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* rows [p+1:2p-1] */
|
|
Packit |
67cb25 |
for (i = penalty2_P; i < penalty2_N - 1; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = gsl_vector_get(x, i - penalty2_P + 1);
|
|
Packit |
67cb25 |
double vi = gsl_vector_get(v, i - penalty2_P + 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set(fvv, i, sqrt_alpha / 100.0 * exp(xi / 10.0) * vi * vi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
(void)params; /* avoid unused parameter warning */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static gsl_multilarge_nlinear_fdf penalty2_func =
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
penalty2_f,
|
|
Packit |
67cb25 |
penalty2_df,
|
|
Packit |
67cb25 |
penalty2_fvv,
|
|
Packit |
67cb25 |
penalty2_N,
|
|
Packit |
67cb25 |
penalty2_P,
|
|
Packit |
67cb25 |
NULL,
|
|
Packit |
67cb25 |
0,
|
|
Packit |
67cb25 |
0,
|
|
Packit |
67cb25 |
0,
|
|
Packit |
67cb25 |
0
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static test_fdf_problem penalty2_problem =
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
"penalty2",
|
|
Packit |
67cb25 |
penalty2_x0,
|
|
Packit |
67cb25 |
NULL,
|
|
Packit |
67cb25 |
&penalty2_epsrel,
|
|
Packit |
67cb25 |
&penalty2_checksol,
|
|
Packit |
67cb25 |
&penalty2_func
|
|
Packit |
67cb25 |
};
|