/*
 
      Testexample for QNRES: Solution of the
      Rosenbrock problem (modified start vector).
 
 *  Written by        L. Weimann 
 *  Purpose           Testexample for code QNRES
 *  Version           1.1
 *  Revision          May 2006
 *  Latest Change     May 2006
 *  Library           NewtonLib
 *  Code              C, Double Precision
 *  Environment       Standard C environment on PC's,
                      workstations and hosts.
 
*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "nleq.h"

extern void f(int *n, double *x, double *fx, int *fail);
extern void jac(int *argn, int *lddfx, double *x,
         JACOBIAN *dfx, int *iflag);

static FILE *out;

int main()
{ struct NLEQ_FUN *problem = malloc(sizeof(struct NLEQ_FUN));
  struct NLEQ_OPT *opt     = malloc(sizeof(struct NLEQ_OPT));
  struct NLEQ_INFO *info   = malloc(sizeof(struct NLEQ_INFO));
  int  i, n=2, rcode;  double x[2];
  double *fscale;

  if (!problem || !opt || !info) 
    { fprintf(stdout,"\n Memory allocation error!\n"); exit(-10);};
  problem->fun = f;
  opt->tol = 1.0e-10;
  opt->maxiter = 50;
  opt->errorlevel = Verbose;
  opt->monitorlevel = Verbose;
  opt->datalevel = Verbose;
  opt->errorfile = fopen("test_qnres.out","w");
  if (!opt->errorfile) fprintf(stdout,
     "\n open of test_qnres.out failed!\n");
  opt->monitorfile = opt->errorfile;
  opt->datafile = fopen("test_qnres.dat","w");
  if (!opt->datafile) fprintf(stdout,
     "\n open of test_qnres.dat failed!\n");
  opt->iterfile = NULL;
  opt->resfile  = NULL;
  opt->miscfile = NULL;
  opt->nonlin = Highly_Nonlinear;
  opt->restricted = False;
  opt->scaleopt = StandardScale;
  rcode = nleq_fwalloc(n,&(opt->scale),"opt->scale");
  if ( rcode != 0 ) exit(-10);
  opt->nleqcalled = False;
  fscale = opt->scale;
  if (opt->monitorfile) out=opt->monitorfile;
  else                  out=stdout;
  if ( opt->monitorfile && opt->datafile )
    fprintf(stdout," monitor: test_qnres.out , data: test_qnres.dat\n");
  fprintf(out,"\n Rosenbrock problem\n");
  for (i=1;i<=2;i++)
    {
      if ( i%2 == 1 )  problem->jac = jac;
      else             problem->jac = NULL;
      x[0]= 0.0 ;  x[1]= 20.0;
      fscale[0] = 1.0; fscale[1] = 1.0;
      fprintf(out,"\n start vector: (%12e,%12e)\n",x[0],x[1]);
      fprintf(out," scale vector: (%12e,%12e)\n",
                     fscale[0],fscale[1]);
      qnres(*problem,n,x,opt,info);
      if ( info->rcode == 0 )
        { fprintf(out,"\n solution: (%12e,%12e)\n",x[0],x[1]);
          fprintf(out,"\n residuum=%e\n",info->precision);
        }
      else
         fprintf(out,"\n QNRES failed - no solution found\n");
      fprintf(out," iter    = %i\n",info->iter);
      fprintf(out," rcode   = %i\n",info->rcode);
      fprintf(out," subcode = %i\n",info->subcode);
      fprintf(out," nfun    = %i\n",info->nofunevals);
      fprintf(out," njac    = %i\n",info->nojacevals);
    };
  exit(0);
}

void f(int *argn, double *x, double *fx, int *iflag)
{
    fx[0] = 1.0 - x[0];
    fx[1] = 10.0*(x[1] - x[0]*x[0]);
    *iflag=0; return;
}

#define DFX(I,J) dfx[(*lddfx)*I+J]

void jac(int *argn, int *lddfx, double *x,
         JACOBIAN *dfx, int *iflag)
{
      DFX(0,0) = -1.0;
      DFX(0,1) =  0.0;
      DFX(1,0) = -20.0*x[0];
      DFX(1,1) =  10.0;
      *iflag=0; return;
}