|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <stdio.h>
|
|
Packit |
67cb25 |
#include <sys/time.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_multilarge_nlinear.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_spblas.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_spmatrix.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* parameters for functions */
|
|
Packit |
67cb25 |
struct model_params
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double alpha;
|
|
Packit |
67cb25 |
gsl_spmatrix *J;
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* penalty function */
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
penalty_f (const gsl_vector * x, void *params, gsl_vector * f)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
struct model_params *par = (struct model_params *) params;
|
|
Packit |
67cb25 |
const double sqrt_alpha = sqrt(par->alpha);
|
|
Packit |
67cb25 |
const size_t p = x->size;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
double sum = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < p; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = gsl_vector_get(x, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set(f, i, sqrt_alpha*(xi - 1.0));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
sum += xi * xi;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set(f, p, sum - 0.25);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
penalty_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 |
struct model_params *par = (struct model_params *) params;
|
|
Packit |
67cb25 |
const size_t p = x->size;
|
|
Packit |
67cb25 |
size_t j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* store 2*x in last row of J */
|
|
Packit |
67cb25 |
for (j = 0; j < p; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xj = gsl_vector_get(x, j);
|
|
Packit |
67cb25 |
gsl_spmatrix_set(par->J, p, j, 2.0 * xj);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute v = op(J) u */
|
|
Packit |
67cb25 |
if (v)
|
|
Packit |
67cb25 |
gsl_spblas_dgemv(TransJ, 1.0, par->J, u, 0.0, v);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (JTJ)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view diag = gsl_matrix_diagonal(JTJ);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute J^T J = [ alpha*I_p + 4 x x^T ] */
|
|
Packit |
67cb25 |
gsl_matrix_set_zero(JTJ);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* store 4 x x^T in lower half of JTJ */
|
|
Packit |
67cb25 |
gsl_blas_dsyr(CblasLower, 4.0, x, JTJ);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* add alpha to diag(JTJ) */
|
|
Packit |
67cb25 |
gsl_vector_add_constant(&diag.vector, par->alpha);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
penalty_fvv (const gsl_vector * x, const gsl_vector * v,
|
|
Packit |
67cb25 |
void *params, gsl_vector * fvv)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t p = x->size;
|
|
Packit |
67cb25 |
double normv = gsl_blas_dnrm2(v);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set_zero(fvv);
|
|
Packit |
67cb25 |
gsl_vector_set(fvv, p, 2.0 * normv * normv);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
(void)params; /* avoid unused parameter warning */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
solve_system(const gsl_vector *x0, gsl_multilarge_nlinear_fdf *fdf,
|
|
Packit |
67cb25 |
gsl_multilarge_nlinear_parameters *params)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const gsl_multilarge_nlinear_type *T = gsl_multilarge_nlinear_trust;
|
|
Packit |
67cb25 |
const size_t max_iter = 200;
|
|
Packit |
67cb25 |
const double xtol = 1.0e-8;
|
|
Packit |
67cb25 |
const double gtol = 1.0e-8;
|
|
Packit |
67cb25 |
const double ftol = 1.0e-8;
|
|
Packit |
67cb25 |
const size_t n = fdf->n;
|
|
Packit |
67cb25 |
const size_t p = fdf->p;
|
|
Packit |
67cb25 |
gsl_multilarge_nlinear_workspace *work =
|
|
Packit |
67cb25 |
gsl_multilarge_nlinear_alloc(T, params, n, p);
|
|
Packit |
67cb25 |
gsl_vector * f = gsl_multilarge_nlinear_residual(work);
|
|
Packit |
67cb25 |
gsl_vector * x = gsl_multilarge_nlinear_position(work);
|
|
Packit |
67cb25 |
int info;
|
|
Packit |
67cb25 |
double chisq0, chisq, rcond, xsq;
|
|
Packit |
67cb25 |
struct timeval tv0, tv1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gettimeofday(&tv0, NULL);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* initialize solver */
|
|
Packit |
67cb25 |
gsl_multilarge_nlinear_init(x0, fdf, work);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* store initial cost */
|
|
Packit |
67cb25 |
gsl_blas_ddot(f, f, &chisq0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* iterate until convergence */
|
|
Packit |
67cb25 |
gsl_multilarge_nlinear_driver(max_iter, xtol, gtol, ftol,
|
|
Packit |
67cb25 |
NULL, NULL, &info, work);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gettimeofday(&tv1, NULL);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* store final cost */
|
|
Packit |
67cb25 |
gsl_blas_ddot(f, f, &chisq);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute final ||x||^2 */
|
|
Packit |
67cb25 |
gsl_blas_ddot(x, x, &xsq);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* store cond(J(x)) */
|
|
Packit |
67cb25 |
gsl_multilarge_nlinear_rcond(&rcond, work);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* print summary */
|
|
Packit |
67cb25 |
fprintf(stderr, "%-25s %-5zu %-4zu %-5zu %-6zu %-4zu %-10.4e %-10.4e %-7.2f %-11.4e %.2f\n",
|
|
Packit |
67cb25 |
gsl_multilarge_nlinear_trs_name(work),
|
|
Packit |
67cb25 |
gsl_multilarge_nlinear_niter(work),
|
|
Packit |
67cb25 |
fdf->nevalf,
|
|
Packit |
67cb25 |
fdf->nevaldfu,
|
|
Packit |
67cb25 |
fdf->nevaldf2,
|
|
Packit |
67cb25 |
fdf->nevalfvv,
|
|
Packit |
67cb25 |
chisq0,
|
|
Packit |
67cb25 |
chisq,
|
|
Packit |
67cb25 |
1.0 / rcond,
|
|
Packit |
67cb25 |
xsq,
|
|
Packit |
67cb25 |
(tv1.tv_sec - tv0.tv_sec) + 1.0e-6 * (tv1.tv_usec - tv0.tv_usec));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_multilarge_nlinear_free(work);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
main (void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t p = 2000;
|
|
Packit |
67cb25 |
const size_t n = p + 1;
|
|
Packit |
67cb25 |
gsl_vector *f = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
gsl_vector *x = gsl_vector_alloc(p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* allocate sparse Jacobian matrix with 2*p non-zero elements in triplet format */
|
|
Packit |
67cb25 |
gsl_spmatrix *J = gsl_spmatrix_alloc_nzmax(n, p, 2 * p, GSL_SPMATRIX_TRIPLET);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_multilarge_nlinear_fdf fdf;
|
|
Packit |
67cb25 |
gsl_multilarge_nlinear_parameters fdf_params =
|
|
Packit |
67cb25 |
gsl_multilarge_nlinear_default_parameters();
|
|
Packit |
67cb25 |
struct model_params params;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
params.alpha = 1.0e-5;
|
|
Packit |
67cb25 |
params.J = J;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* define function to be minimized */
|
|
Packit |
67cb25 |
fdf.f = penalty_f;
|
|
Packit |
67cb25 |
fdf.df = penalty_df;
|
|
Packit |
67cb25 |
fdf.fvv = penalty_fvv;
|
|
Packit |
67cb25 |
fdf.n = n;
|
|
Packit |
67cb25 |
fdf.p = p;
|
|
Packit |
67cb25 |
fdf.params = ¶m;;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < p; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* starting point */
|
|
Packit |
67cb25 |
gsl_vector_set(x, i, i + 1.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* store sqrt(alpha)*I_p in upper p-by-p block of J */
|
|
Packit |
67cb25 |
gsl_spmatrix_set(J, i, i, sqrt(params.alpha));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fprintf(stderr, "%-25s %-4s %-4s %-5s %-6s %-4s %-10s %-10s %-7s %-11s %-10s\n",
|
|
Packit |
67cb25 |
"Method", "NITER", "NFEV", "NJUEV", "NJTJEV", "NAEV", "Init Cost",
|
|
Packit |
67cb25 |
"Final cost", "cond(J)", "Final |x|^2", "Time (s)");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fdf_params.scale = gsl_multilarge_nlinear_scale_levenberg;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fdf_params.trs = gsl_multilarge_nlinear_trs_lm;
|
|
Packit |
67cb25 |
solve_system(x, &fdf, &fdf_params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fdf_params.trs = gsl_multilarge_nlinear_trs_lmaccel;
|
|
Packit |
67cb25 |
solve_system(x, &fdf, &fdf_params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fdf_params.trs = gsl_multilarge_nlinear_trs_dogleg;
|
|
Packit |
67cb25 |
solve_system(x, &fdf, &fdf_params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fdf_params.trs = gsl_multilarge_nlinear_trs_ddogleg;
|
|
Packit |
67cb25 |
solve_system(x, &fdf, &fdf_params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fdf_params.trs = gsl_multilarge_nlinear_trs_subspace2D;
|
|
Packit |
67cb25 |
solve_system(x, &fdf, &fdf_params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fdf_params.trs = gsl_multilarge_nlinear_trs_cgst;
|
|
Packit |
67cb25 |
solve_system(x, &fdf, &fdf_params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free(f);
|
|
Packit |
67cb25 |
gsl_vector_free(x);
|
|
Packit |
67cb25 |
gsl_spmatrix_free(J);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return 0;
|
|
Packit |
67cb25 |
}
|