/* E05UC_P0W_F C++ Header Example Program.
*
* Copyright 2019 Numerical Algorithms Group.
* Mark 27, 2019.
*/
#include <nag.h>
#include <nagad.h>
#include <stdio.h>
#include <math.h>
#include <iostream>
#include <string>
using namespace std;
extern "C"
{
static void NAG_CALL objfun(void* &ad_handle,
Integer &mode,
const Integer &n,
const double x[],
double &objf,
double objgrd[],
const Integer &nstate,
Integer iuser[],
double ruser[]);
static void NAG_CALL confun(void* &ad_handle,
Integer &mode,
const Integer &ncnln,
const Integer &n,
const Integer &ldcj,
const Integer needc[],
const double x[],
double c[],
double cjac[],
const Integer &nstate,
Integer iuser[],
double ruser[]);
static void NAG_CALL mystart(void* &ad_handle,
const Integer &npts,
double quas[],
const Integer &n,
const logical &repeat,
const double bl[],
const double bu[],
Integer iuser[],
double ruser[],
Integer &mode);
}
int main(void)
{
// Scalars
int exit_status = 0;
cout << "E05UC_P0W_F C++ Header Example Program Results\n\n";
Integer ifail = 0;
void *ad_handle = 0;
// Skip first line of data file
string mystr;
getline (cin, mystr);
// problem sizes
const Integer n=2, nclin=1, ncnln=2;
Integer nb, npts;
cin >> nb >> npts;
logical repeat = true;
const Integer liopts=740, lopts=485;
Integer lda = nclin, sda = n, ldc = ncnln, ldcj = ncnln, ldr = n;
Integer lb = n + nclin + ncnln, listat = n + nclin +ncnln;
double *a=0, *bl=0, *bu=0, *c=0, *cjac=0, *objf=0;
double *objgrd=0, *clamda=0, *r=0, *x=0, *work=0, *opts=0;
Integer *iopts=0, *info=0, *iter=0, *istate=0;
a = new double [lda*sda];
bl = new double [lb];
bu = new double [lb];
c = new double [ldc*nb];
cjac = new double [ldcj*n*nb];
clamda = new double [lb*nb];
r = new double [ldr*n*nb];
x = new double [n*nb];
objgrd = new double [n*nb];
work = new double [nclin];
opts = new double [lopts];
objf = new double [nb];
istate = new Integer [listat*nb];
iopts = new Integer [liopts];
info = new Integer [nb];
iter = new Integer [nb];
bl[0] = -500.0;
bl[1] = -500.0;
bl[2] = -10000.0;
bl[3] = -1.0;
bl[4] = -0.9;
bu[0] = 500.0;
bu[1] = 500.0;
bu[2] = 10.0;
bu[3] = 500000.0;
bu[4] = 0.9;
a[0] = 3.0;
a[1] = -2.0;
double ruser[6];
ruser[0] = 1.0;
ruser[1] = 1.0;
ruser[2] = 1.0;
ruser[3] = 3.0;
ruser[4] = 0.005;
ruser[5] = 0.01;
// Initialize the solver.
ifail = 0;
e05zk_p0w_f_(ad_handle,"Initialize = E05UCF",iopts,liopts,opts,lopts,
ifail,19);
// Solve the problem
Integer iuser[1];
ifail = -1;
e05uc_p0w_f_(ad_handle,n,nclin,ncnln,a,lda,bl,bu,confun,
objfun,npts,x,n,mystart,repeat,nb,objf,
objgrd,n,iter,c,ldc,cjac,ldcj,n,r,ldr,n,
clamda,lb,istate,listat,iopts,opts,iuser,ruser,
info,ifail);
// Primal results
cout.setf(ios::scientific,ios::floatfield);
cout.precision(4);
Integer l;
if (ifail==0) {
l = nb;
} else if (ifail==8) {
l = info[nb-1];
cout.width(16); cout << iter[nb-1] << " starting points converged" << endl;
} else {
goto END;
}
for (int i=0; i<l; i++) {
cout << "\n Solution number " << i+1 << endl;
cout << "\n Local minimization exited with code " << info[i] << endl;
cout << "\n Varbl Istate Value Lagr Mult" << endl;
cout << endl;
for (int j=0; j<n; j++) {
cout << " V ";
cout.width(4); cout << j+1;
cout.width(4); cout << istate[j+i*lb];
cout.width(12); cout << x[j+i*n];
cout.width(12); cout << clamda[j+i*lb] << endl;
}
if (nclin>0) {
cout << "\n L con Istate Value Lagr Mult" << endl;
const double alpha = 1.0;
const double beta = 0.0;
ifail = 0;
f06pa_p0w_f_(ad_handle,"N",nclin,n,alpha,a,lda,&x[n*i],1,beta,work,1,
ifail,1);
cout << endl;
for (int k=n; k<n+nclin; k++) {
int j = k - n;
cout << " L ";
cout.width(4); cout << j+1;
cout.width(4); cout << istate[k+i*lb];
cout.width(12); cout << work[j];
cout.width(12); cout << clamda[k+i*lb] << endl;
}
}
if (ncnln>0) {
cout << "\n N con Istate Value Lagr Mult" << endl;
cout << endl;
for (int k=n+nclin; k<n+nclin+ncnln; k++) {
int j = k - n - nclin;
cout << " N ";
cout.width(4); cout << j+1;
cout.width(4); cout << istate[k+i*lb];
cout.width(12); cout << c[j+i*ncnln];
cout.width(12); cout << clamda[k+i*lb] << endl;
}
}
cout << "\n Final objective value = ";
cout.width(12); cout << objf[i];
cout << "\n QP multipliers" << endl;
for (int k=0; k<n+nclin+ncnln; k++) {
cout.width(12); cout << clamda[k+i*lb] << endl;
}
cout << endl;
if (l>1) {
cout << "\n ---------------------------------------------------------\n";
}
}
END:
delete [] a;
delete [] bl;
delete [] bu;
delete [] c;
delete [] cjac;
delete [] clamda;
delete [] r;
delete [] x;
delete [] objgrd;
delete [] work;
delete [] opts;
delete [] objf;
delete [] istate;
delete [] iopts;
delete [] info;
delete [] iter;
return exit_status;
}
static void NAG_CALL objfun(void* &ad_handle,
Integer &mode,
const Integer &n,
const double x[],
double &objf,
double objgrd[],
const Integer &nstate,
Integer iuser[],
double ruser[])
{
if (mode==0 || mode==2) {
objf = 0.0;
for (int i=0; i<n; i++) {
if (x[i] >= 0.0) {
objf += x[i]*sin(ruser[0]*sqrt(x[i]));
} else {
objf += x[i]*sin(ruser[0]*sqrt(-x[i]));
}
}
}
if (mode==1 || mode==2) {
for (int i=0; i<n; i++) {
if (x[i] >= 0.0) {
objgrd[i] = sin(ruser[0]*sqrt(x[i])) +
0.5*ruser[0]*sqrt(x[i])*cos(ruser[0]*sqrt(x[i]));
} else {
objgrd[i] = sin(ruser[0]*sqrt(-x[i])) +
0.5*ruser[0]*sqrt(-x[i])*cos(ruser[0]*sqrt(-x[i]));
}
}
}
return;
}
static void NAG_CALL confun(void* &ad_handle,
Integer &mode,
const Integer &ncnln,
const Integer &n,
const Integer &ldcj,
const Integer needc[],
const double x[],
double c[],
double cjac[],
const Integer &nstate,
Integer iuser[],
double ruser[])
{
for (int k = 0; k<ncnln; ++k) {
if (mode==0 || mode==2) {
if (k==0) {
c[k] = ruser[1]*x[0]*x[0] - ruser[2]*x[1]*x[1] + ruser[3]*x[0]*x[1];
} else {
c[k] = cos(x[0]*x[0]*ruser[4]*ruser[4] + x[1]*ruser[5]);
}
}
if (mode==1 || mode==2) {
if (k==0) {
cjac[0] = 2.0*ruser[1]*x[0] + ruser[3]*x[1];
cjac[ncnln] = -2.0*ruser[2]*x[1] + ruser[3]*x[0];
} else {
double theta;
theta = x[0]*x[0]*ruser[4]*ruser[4] + x[1]*ruser[5];
cjac[1] = -sin(theta)*2.0*x[0]*ruser[4]*ruser[4];
cjac[ncnln+1] = -sin(theta)*ruser[5];
}
}
}
return;
}
static void NAG_CALL mystart(void* &ad_handle,
const Integer &npts,
double quas[],
const Integer &n,
const logical &repeat,
const double bl[],
const double bu[],
Integer iuser[],
double ruser[],
Integer &mode)
{
if (repeat) {
// Generate a uniform spread of points between bl and bu.
for (int i = 0; i<npts; i++) {
double rq = ((double) (i-1))/((double)(npts-1));
for (int j=0; j<n; j++) {
quas[j+i*n] = bl[j] + rq*(bu[j]-bl[j]);
}
}
} else {
// Generate a non-repeatable spread of points between bl and bu.
const Integer genid = 2, subid = 53;
Integer lstate = -1, sdum[1];
Integer ifail = 0;
g05kgf_(genid,subid,sdum,lstate,ifail);
Integer *state=0;
double *rquas=0;
state = new Integer [lstate];
rquas = new double[n];
ifail = 0;
g05kgf_(genid,subid,state,lstate,ifail);
for (int i=0; i<npts;i++) {
ifail = 0;
g05saf_(n,state,rquas,ifail);
for (int j=0; j<n;j++) {
quas[j+n*i] = bl[j] + (bu[j]-bl[j])*rquas[j];
}
}
delete [] state;
delete [] rquas;
}
// Set mode negative to terminate execution for any reason.
mode = 0;
return;
}