```/* nag_glopt_bnd_pso (e05sac) Example Program.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*/
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage05.h>
#include <nagx02.h>
#include <nagx04.h>

#ifdef __cplusplus
extern "C"
{
#endif
static void NAG_CALL objfun_schwefel(Integer *mode, Integer ndim,
const double x[], double *objf,
double vecout[], Integer nstate,
Nag_Comm *comm);
static void NAG_CALL monmod(Integer ndim, Integer npar, double x[],
const double xb[], double fb,
const double xbest[], const double fbest[],
const Integer itt[], Nag_Comm *comm,
Integer *inform);
#ifdef __cplusplus
}
#endif

static void display_result(Integer ndim, const double xb[], const double
x_target[], double fb, double f_target,
const Integer itt[], Integer inform);
static void display_option(const char *optstr, Nag_VariableType optype,
Integer ivalue, double rvalue, const char *cvalue);
static void get_known_solution(Integer ndim, double x_target[],
double *f_target);

/* Global constants - set the behaviour of the monitoring function.*/
static const Integer detail_level = 0;
static const Integer report_freq = 100;

int main(void)
{
/* This example program demonstrates how to use
* nag_glopt_bnd_pso (e05sac) in standard execution, and with a
* selection of coupled local minimizers.
* The non-default option 'Repeatability = On' is used here, giving
* repeatable results.
*/
/* Scalars */
Integer ndim = 2, npar = 5;
Integer exit_status = 0, lcvalue = 17;
Integer liopts = 100, lopts = 100;
double fb, f_target, rvalue;
Integer i, inform, ivalue;
/* Arrays */
static double ruser[2] = { -1.0, -1.0 };
char cvalue[17], optstr[81];
double *bl = 0, *bu = 0, opts[100], *xb = 0, *x_target = 0;
Integer iopts[100], itt[6];
/* Nag Types */
Nag_VariableType optype;
NagError fail;
Nag_Comm comm;

printf("nag_glopt_bnd_pso (e05sac) Example Program Results\n\n");

/* For communication with user-supplied functions: */
comm.user = ruser;

printf("Minimization of the Schwefel function.\n\n");

/* Allocate memory. */
if (!(bl = NAG_ALLOC(ndim, double)) || !(bu = NAG_ALLOC(ndim, double)) ||
!(xb = NAG_ALLOC(ndim, double))
|| !(x_target = NAG_ALLOC(ndim, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}

/* Store the known solution of the problem for a comparison. */
get_known_solution(ndim, x_target, &f_target);

/* Set box bounds for problem. */
for (i = 0; i < ndim; i++) {
bl[i] = -500.0;
bu[i] = 500.0;
}

/* Initialize fail structures. */
INIT_FAIL(fail);

/* Initialize the option arrays for nag_glopt_bnd_pso (e05sac)
* using nag_glopt_opt_set (e05zkc).
*/
nag_glopt_opt_set("Initialize = e05sac", iopts, liopts, opts, lopts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_glopt_opt_set (e05zkc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}

/* Query some default option values. */
printf("  Default Option Queries:\n\n");
/* nag_glopt_opt_get (e05zlc).
* Option getting routine for nag_glopt_bnd_pso (e05sac).
*/
nag_glopt_opt_get("Boundary", &ivalue, &rvalue, cvalue, lcvalue, &optype,
iopts, opts, &fail);
if (fail.code == NE_NOERROR) {
display_option("Boundary", optype, ivalue, rvalue, cvalue);
nag_glopt_opt_get("Maximum Iterations Completed", &ivalue, &rvalue,
cvalue, lcvalue, &optype, iopts, opts, &fail);
}
if (fail.code == NE_NOERROR) {
display_option("Maximum Iterations Completed", optype, ivalue, rvalue,
cvalue);
nag_glopt_opt_get("Distance Tolerance", &ivalue, &rvalue, cvalue, lcvalue,
&optype, iopts, opts, &fail);
}
if (fail.code != NE_NOERROR) {
printf("Error from nag_glopt_opt_set (e05zkc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
display_option("Distance Tolerance", optype, ivalue, rvalue, cvalue);

/* ------------------------------------------------------------------ */

printf("\n1. Solution without using coupled local minimizer.\n\n");

/* Set various options to non-default values if required. */
nag_glopt_opt_set("Repeatability = On", iopts, liopts, opts, lopts, &fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Verify Gradients = Off", iopts, liopts, opts, lopts,
&fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Boundary = Hyperspherical", iopts, liopts, opts, lopts,
&fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Maximum iterations static = 150", iopts, liopts, opts,
lopts, &fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Repulsion Initialize = 30", iopts, liopts, opts, lopts,
&fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Repulsion Finalize = 30", iopts, liopts, opts, lopts,
&fail);

if (fail.code != NE_NOERROR) {
printf("Error from nag_glopt_opt_set (e05zkc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}

/* nag_glopt_bnd_pso (e05sac).
* Global optimization using particle swarm algorithm (PSO),
* bound constraints only.
*/
nag_glopt_bnd_pso(ndim, npar, xb, &fb, bl, bu, objfun_schwefel, monmod,
iopts, opts, &comm, itt, &inform, &fail);
/* It is essential to test fail.code on exit. */
switch (fail.code) {
case NE_NOERROR:
case NW_FAST_SOLUTION:
case NW_SOLUTION_NOT_GUARANTEED:
/* No errors, best found solution at xb returned in fb. */
display_result(ndim, xb, x_target, fb, f_target, itt, inform);
break;
case NE_USER_STOP:
/* Exit flag set in objfun or monmod and returned in inform. */
display_result(ndim, xb, x_target, fb, f_target, itt, inform);
break;
default: /* An error was detected. */
exit_status = 1;
printf("Error from nag_glopt_bnd_pso (e05sac)\n%s\n", fail.message);
goto END;
}

/* ------------------------------------------------------------------ */

printf("2. Solution using coupled local minimizer "
"nag_opt_simplex_easy (e04cbc).\n\n");

/* Initialize fail structures. */
INIT_FAIL(fail);

/* Set an objective target. */
sprintf(optstr, "Target Objective Value = %32.16e", f_target);
nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, &fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Target Objective Tolerance = 1.0e-5", iopts, liopts,
opts, lopts, &fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Target Objective Safeguard = 1.0e-8", iopts, liopts,
opts, lopts, &fail);

/* Set the local minimizer to be nag_opt_simplex_easy (e04cbc)
* and set corresponding options.
*/
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Local Minimizer = e04cbc", iopts, liopts, opts, lopts,
&fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Local Interior Iterations = 10", iopts, liopts, opts,
lopts, &fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Local Exterior Iterations = 20", iopts, liopts, opts,
lopts, &fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Local Interior Tolerance = 1.0e-4", iopts, liopts,
opts, lopts, &fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Local Exterior Tolerance = 1.0e-4", iopts, liopts,
opts, lopts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_glopt_opt_set (e05zkc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}

/* Search for the global optimum. */
nag_glopt_bnd_pso(ndim, npar, xb, &fb, bl, bu, objfun_schwefel, monmod,
iopts, opts, &comm, itt, &inform, &fail);
/* It is essential to test fail.code on exit. */
switch (fail.code) {
case NE_NOERROR:
case NW_FAST_SOLUTION:
case NW_SOLUTION_NOT_GUARANTEED:
/* No errors, best found solution at xb returned in fb. */
display_result(ndim, xb, x_target, fb, f_target, itt, inform);
break;
case NE_USER_STOP:
/* Exit flag set in objfun or monmod and returned in inform. */
display_result(ndim, xb, x_target, fb, f_target, itt, inform);
break;
default: /* An error was detected. */
exit_status = 1;
printf("Error from nag_glopt_bnd_pso (e05sac)\n%s\n", fail.message);
goto END;
}

/* ----------------------------------------------------------------- */

printf("3. Solution using coupled local minimizer "

/* Initialize fail structures. */
INIT_FAIL(fail);

/* set the local minimizer to be nag_opt_conj_grad (e04dgc)
* and set corresponding options.
*/
nag_glopt_opt_set("Local Minimizer = e04dgc", iopts, liopts, opts, lopts,
&fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Local Interior Iterations = 5", iopts, liopts, opts,
lopts, &fail);
if (fail.code == NE_NOERROR)
nag_glopt_opt_set("Local Exterior Iterations = 20", iopts, liopts, opts,
lopts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_glopt_opt_set (e05zkc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}

/* Search for the global optimum. */
nag_glopt_bnd_pso(ndim, npar, xb, &fb, bl, bu, objfun_schwefel, monmod,
iopts, opts, &comm, itt, &inform, &fail);
/* It is essential to test fail.code on exit. */
switch (fail.code) {
case NE_NOERROR:
case NW_FAST_SOLUTION:
case NW_SOLUTION_NOT_GUARANTEED:
/* No errors, best found solution at xb returned in fb. */
display_result(ndim, xb, x_target, fb, f_target, itt, inform);
break;
case NE_USER_STOP:
/* Exit flag set in objfun or monmod and returned in inform. */
display_result(ndim, xb, x_target, fb, f_target, itt, inform);
break;
default: /* An error was detected. */
exit_status = 1;
printf("Error from nag_glopt_bnd_pso (e05sac)\n%s\n", fail.message);
goto END;
}

END:
/* Clean up memory. */
NAG_FREE(bl);
NAG_FREE(bu);
NAG_FREE(xb);
NAG_FREE(x_target);

return exit_status;
}

static void NAG_CALL objfun_schwefel(Integer *mode, Integer ndim,
const double x[], double *objf,
double vecout[], Integer nstate,
Nag_Comm *comm)
{
/* Objective function routine returning the schwefel function and
*/
Nag_Boolean evalobjf, evalobjg;
Integer i;
if (comm->user[0] == -1.0) {
printf("(User-supplied callback objfun_schwefel, first invocation.)\n");
comm->user[0] = 0.0;
}
/* Test nstate to indicate what stage of computation has been reached. */
switch (nstate) {
case 2:
/* objfun is called for the very first time. */
break;
case 1:
/* objfun is called on entry to a NAG local minimizer. */
break;
default: /* This will be the normal value of nstate. */
;
}
/* Test mode to determine whether to calculate objf and/or objgrd. */
evalobjf = Nag_FALSE;
evalobjg = Nag_FALSE;
switch (*mode) {
case 0:
case 5:
/* Only the value of the objective function is needed. */
evalobjf = Nag_TRUE;
break;
case 1:
case 6:
/* Only the values of the ndim gradients are required. */
evalobjg = Nag_TRUE;
break;
case 2:
case 7:
/* Both the objective function and the ndim gradients are required. */
evalobjf = Nag_TRUE;
evalobjg = Nag_TRUE;
}
if (evalobjf) {
/* Evaluate the objective function. */
*objf = 0.0;
for (i = 0; i < ndim; i++)
*objf += x[i] * sin(sqrt(fabs(x[i])));
}
if (evalobjg) {
/* Calculate the gradient of the objective function,
* and return the result in vecout.
*/
for (i = 0; i < ndim; i++) {
vecout[i] = sqrt(fabs(x[i]));
vecout[i] = sin(vecout[i]) + 0.5 * vecout[i] * cos(vecout[i]);
}
}
}

static void NAG_CALL monmod(Integer ndim, Integer npar, double x[],
const double xb[], double fb,
const double xbest[], const double fbest[],
const Integer itt[], Nag_Comm *comm,
Integer *inform)
{
Integer i, j;
#define X(J, I) x[(J-1)*ndim + (I-1)]
#define XBEST(J, I) xbest[(J-1)*ndim + (I-1)]
if (comm->user[1] == -1.0) {
printf("(User-supplied callback monmod, first invocation.)\n");
comm->user[1] = 0.0;
}
if (detail_level) {
/* Report on the first iteration, and every report_freq iterations. */
if (itt[0] == 1 || itt[0] % report_freq == 0) {
printf("* Locations of particles\n");
for (j = 1; j <= npar; j++) {
printf("  * Particle %2" NAG_IFMT "\n", j);
for (i = 1; i <= ndim; i++)
printf("    %2" NAG_IFMT "  %13.5f\n", i, X(j, i));
}
printf("* Cognitive memory\n");
for (j = 1; j <= npar; j++) {
printf("  * Particle %2" NAG_IFMT "\n", j);
printf("    * Best position\n");
for (i = 1; i <= ndim; i++)
printf("      %2" NAG_IFMT "  %13.5f\n", i, XBEST(j, i));
printf("    * Function value at best position\n");
printf("      %13.5f\n", fbest[j - 1]);
}
printf("* Current global optimum candidate\n");
for (i = 1; i <= ndim; i++)
printf("  %2" NAG_IFMT "  %13.5f\n", i, xb[i - 1]);
printf("* Current global optimum value\n");
printf("     %13.5f\n\n", fb);
}
}
/* If required set *inform<0 to force exit. */
*inform = 0;
#undef XBEST
#undef X
}

static void display_option(const char *optstr, Nag_VariableType optype,
Integer ivalue, double rvalue, const char *cvalue)
{
/* Subroutine to query optype and print the appropriate option values. */
switch (optype) {
case Nag_Integer:
printf("%-38s: %13" NAG_IFMT "\n", optstr, ivalue);
break;
case Nag_Real:
printf("%-38s: %13.4f\n", optstr, rvalue);
break;
case Nag_Character:
printf("%-38s: %13s\n", optstr, cvalue);
break;
printf("%-38s: %13" NAG_IFMT " %16s\n", optstr, ivalue, cvalue);
break;
printf("%-38s: %13.4f %16s\n", optstr, rvalue, cvalue);
break;
default:;
}
}

static void display_result(Integer ndim, const double xb[], const double
x_target[], double fb, double f_target,
const Integer itt[], Integer inform)
{
/* Display final results in comparison to known global optimum. */
Integer i;

/* Display final counters. */
printf(" Algorithm Statistics\n");
printf(" --------------------\n");
printf("%-38s: %13" NAG_IFMT "\n", "Total complete iterations", itt[0]);
printf("%-38s: %13" NAG_IFMT "\n", "Complete iterations since improvement",
itt[1]);
printf("%-38s: %13" NAG_IFMT "\n", "Total particles converged to xb",
itt[2]);
printf("%-38s: %13" NAG_IFMT "\n", "Total improvements to global optimum",
itt[3]);
printf("%-38s: %13" NAG_IFMT "\n", "Total function evaluations", itt[4]);
printf("%-38s: %13" NAG_IFMT "\n\n", "Total particles re-initialized",
itt[5]);
/* Display why finalization occurred. */
switch (inform) {
case 1:
printf("Solution Status : Target value achieved\n");
break;
case 2:
printf("Solution Status : Minimum swarm standard deviation obtained\n");
break;
case 3:
printf("Solution Status : Sufficient number of particles converged\n");
break;
case 4:
printf("Solution Status : Maximum static iterations attained\n");
break;
case 5:
printf("Solution Status : Maximum complete iterations attained\n");
break;
case 6:
printf("Solution Status : Maximum function evaluations exceeded\n");
break;
default:
if (inform < 0)
printf("Solution Status : User termination, inform = %16" NAG_IFMT "\n",
inform);
else
printf("Solution Status : Termination, an error has been detected\n");
break;
}
/* Display final objective value and location. */
printf("  Known objective optimum             : %13.5f\n", f_target);
printf("  Achieved objective value            : %13.5f\n\n", fb);

printf("  Comparison between the known optimum and the achieved solution.\n");
printf("          x_target            xb\n");
for (i = 0; i < ndim; i++)
printf("  %2" NAG_IFMT "  %12.2f  %12.2f\n", i + 1, x_target[i], xb[i]);

printf("\n");
}

static void get_known_solution(Integer ndim, double x_target[],
double *f_target)
{
/* Fill in the known solution of multidimensional Schwefel's function. */
Integer i;

if (f_target && x_target && ndim > 0) {
*f_target = -418.9828872724337 * ndim;
for (i = 0; i < ndim; i++)
x_target[i] = -420.9687463599820;
}
}
```