/*
	minimize x1^2 + (1/3)x2^2
	subject to
	-x1 <= 0
	-x2 <= 0
	-x1 -x2 +1 = 0
	
	answer: (x1,x2) = (1/4,3/4)
		f(x1,x2)=1/4
					*/

#include <stdio.h>
#include <math.h>
#include <values.h>
#include "opt.h"

#define	Debug		0
#define	Static		static

Static dbl obj(int n, dbl *x)
{
	dbl	x1, x2;
	
	x1 = x[0];
	x2 = x[1];
	return x1*x1 + x2*x2/3.00;
}

Static dbl *gradobj(int n, dbl *x, dbl *d)
{
	dbl	x1, x2;
	
	x1 = x[0];
	x2 = x[1];
	d[0] = 2.00*x1;
	d[1] = 2.00*x2/3.00;
}

Static dbl *ineq(int n, dbl *x, int m, dbl *d)
{
	dbl	x1, x2;
	
	x1 = x[0];
	x2 = x[1];
	d[0] = -x1;
	d[1] = -x2;
}

Static dbl **gradineq(int n, dbl *x, int m, dbl **d)
{
	dbl	x1, x2;
	
	x1 = x[0];
	x2 = x[1];
	(d[0])[0] = -1.00;
	(d[0])[1] =  0.00;
	(d[1])[0] =  0.00;
	(d[1])[1] = -1.00;
}

Static dbl *cond(int n, dbl *x, int l, dbl *d)
{
	dbl	x1, x2;
	
	x1 = x[0];
	x2 = x[1];
	d[0] = - x1 - x2 + 1;
}

Static dbl **gradcond(int n, dbl *x, int l, dbl **d)
{
	dbl	x1, x2;
	
	x1 = x[0];
	x2 = x[1];
	(d[0])[0] = -1.00;
	(d[0])[1] = -1.00;
}

#define	NVAR	2
#define NINQ	2
#define	NCND	1

main()
{
	dbl	x[100], y;
	enum status	stat;
	
	x[0] = 100.00001;
	x[1] = 100.00001;
	matrixprintformat(" %l20.15f ");
	stat = MultiplierMethod(NVAR, obj, gradobj,
				NINQ, ineq, gradineq,
				NCND, cond, gradcond,
				x, 10, 1.0e-6);
	if (stat == success) {
		fprintf(stderr, "----- solution -----\n");
		vectorfprint(stderr, NVAR, x);
	} else {
		fprintf(stderr, "status = %d\n", stat);
	}
}