/*********************************************************
main.cpp
cpp example setup and call for nag routines
03/11/2004
gem@cs.brown.edu
**********************************************************/

#include <stdio.h>
#include <string.h>
#include <fstream.h>
#include <iostream.h>


typedef double* Matr;
Matr alloc_matr(int l, int c)
  {
  return new double [l*c];
  }

#include "e04uef.h"
#include "e04ucf.h"



void CONFUN(int* mode, int* ncnln, int* n, int* ldcj, int* needc,double* x, double* c, Matr cjac, int* nstate, int* iuser, double* user){

   int i,j;
   if(*nstate == 1)
    for (i = 0; i < 100; ++i) 
       cjac[i]= 0.0;

   if(needc[0] > 0) {
     if((*mode == 0) || (*mode == 2))
        c[0] = x[0]*x[0] + x[1]*x[1] + x[2]*x[2] +x[3]*x[3]; 
     if((*mode == 1) || (*mode == 2)){
        cjac[0*10 + 0] = 2*x[0];
        cjac[1*10 + 0] = 2*x[1];
        cjac[2*10 + 0] = 2*x[2];
        cjac[3*10 + 0] = 2*x[3]; //???
     };
    };
   if(needc[1] > 0) {
     if((*mode == 0) || (*mode == 2))
        c[1] = x[0]*x[1]*x[2]*x[3]; 
     if((*mode == 1) || (*mode == 2)){
        cjac[0*10+ 1] = x[1]*x[2]*x[3];
        cjac[1*10+ 1] = x[0]*x[2]*x[3];
        cjac[2*10+ 1] = x[0]*x[1]*x[3];
        cjac[3*10+ 1] = x[0]*x[1]*x[2]; //????
     };
    };

   
};

void OBJFUN(int* mode, int* n, double* x, double* objf, double* objgrd, int* nstate, int* iuser, double* user){

  if((*mode == 0) ||(*mode==2)) *objf = x[0]*x[3]*(x[0]+x[1]+x[2]) + x[2];
  if((*mode == 1) ||(*mode==2)){
      objgrd[0] = x[3] * (2*x[0]+x[1]+x[2]);
      objgrd[1] = x[0] * x[3];
      objgrd[2] = x[0] * x[3] + 1.0;
      objgrd[3] = x[0] * (x[0]+x[1]+x[2]);
     };
};





int main(){
  int i,j;
  int n,ncnln,nclin;
  n = 4;
  ncnln = 2;
  nclin = 1;
  int IFAIL;
  int lda = 10, ldcj = 10, ldr = 10;
  double OBJF;

  Matr A = alloc_matr(nclin,n);
  for(i=0;i<4;i++)
      A[i] = 1.0;
    
  /*
  for(i=0;i<nclin;i++)
    for(j=0;j<n;j++)
       A[i][j] = 1.0;
  */
  double* BL = new double [n + ncnln + nclin];

  double* BU = new double [n + ncnln + nclin];

  double* X = new double [n];

  double* C = new double [ncnln];

  Matr CJAC = alloc_matr(ldcj,10);
  double* CLAMDA = new double [ n + ncnln + nclin];
  double* OBJGRD = new double [10];
  Matr R = alloc_matr(ldr, 10);
  double USER;
  
  double* WORK = new double [1000];
 
  int* ISTATE = new int [n + ncnln + nclin];

  int ITER;
  int IUSER;
  int* IWORK = new int [100];
  int LIWORK = 100;
  int LWORK = 1000;
  for(i=0;i<4;i++) BL[i] = 1.0;

  BL[4] = -1.0e+25;
  BL[5] = -1.0e+25;
  BL[6] = 25.0;


  for(i=0;i<4;i++) BU[i] = 5.0;

  BU[4] = 20.0;
  BU[5] = 40.0;
  BU[6] = 1.0e+25;

  X[0] = 1.0;
  X[1] = 5.0;
  X[2] = 5.0;
  X[3] = 1.0;
 
  IFAIL = -1;

  e04ucf(&n, &nclin, &ncnln, &lda, &ldcj,&ldr,A,BL,BU, CONFUN, OBJFUN, &ITER,ISTATE,C, CJAC,CLAMDA, &OBJF,OBJGRD,R,X, IWORK, &LIWORK,WORK, &LWORK, &IUSER, &USER,&IFAIL);
  
  return 1;

};