/*
 
      Testexample for NLEQ_RES: Computation of the Tchebychev
      Polynom's evaluation points for quadrature-formulas.
 
 *  Written by        L. Weimann 
 *  Purpose           Testexample for code NLEQ_RES
 *  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"

void printvec(char vname[], int n, double *v);

void f(int *n, double *x, double *fx, int *fail);
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, rcode, nn=9; int n;  double x[9];
  double *fscale;

  if (!problem || !opt || !info) 
    { fprintf(stdout,"\n Memory allocation error!\n"); exit(-10);};
  problem->fun = f;
  problem->jac = jac;
  opt->tol = 1.0e-6;
  opt->maxiter = 50;
  opt->errorlevel = Verbose;
  opt->monitorlevel = Verbose;
  opt->datalevel = Verbose;
  opt->errorfile = fopen("test_nleq_res.out","w");
  if (!opt->errorfile) fprintf(stdout,
     "\n open of test_nleq_res.out failed!\n");
  opt->monitorfile = opt->errorfile;
  opt->datafile = fopen("test_nleq_res.dat","w");
  if (!opt->datafile) fprintf(stdout,
     "\n open of test_nleq_res.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);
  fscale = opt->scale;
  if (opt->monitorfile) out=opt->monitorfile;
  else                  out=stdout;
  if ( opt->monitorfile && opt->datafile )
    fprintf(stdout," monitor: test_nleq_res.out , data: test_nleq_res.dat\n");

  for (n=2;n<=nn;n++)
    { for (i=0;i<n;i++) x[i]=  (double)(i+1.0) / (double)(n+1.0);
      for (i=0;i<n;i++) fscale[i]=  1.0;
      fprintf(out,"\n Tchebyquad problem, n=%i\n\n",n);
      printvec("start vector",n,x);
      printvec("scale vector",n,fscale);
      nleq_res(*problem,n,x,opt,info);
      if ( info->rcode == 0 )
        { fprintf(out,"\n");
          printvec("solution",n,x);
          fprintf(out,"\n residuum=%e\n",info->precision);
        }
      else
         fprintf(out,"\n NLEQ_RES 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 printvec(char vname[], int n, double *v)
{  int i;
   fprintf(out," %s: (",vname);
   for (i=0;i<n-1;i++) fprintf(out,"%12e,",v[i]);
   fprintf(out,"%12e)\n",v[n-1]);
}

void f(int *argn, double *x, double *fx, int *iflag)
{
   int i, i1, l, n=*argn;
   double ti2, ti1, ti, factt;

   for (i=1;i<n;i=i+2)
     { fx[i-1] = 0.0; i1=i+1; 
       fx[i] = (double)(n)/(double)(i1*i1-1.0); };
   if( n%2 == 1) fx[n-1]=0.0;
   for (l=0;l<n;l++)
     { factt = 4.0*x[l]-2.0;
       ti2 = 1.0;  ti1 = 0.5*factt;
       fx[0]=ti1+fx[0];
       for (i=1;i<n;i++)
         { ti = factt*ti1-ti2;
           fx[i]=ti+fx[i];
           ti2 = ti1; ti1 = ti;
         };
      };
    *iflag=0; return;
}

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

void jac(int *argn, int *lddfx, double *x,
         JACOBIAN *dfx, int *iflag)
{
      int i,j, n=*argn;
      double ti2,ti1,ti,factt,tabli2,tabli1,tabli;


      for (j=0;j<n;j++)
       { factt = 4.0*x[j]-2.0;
         ti2 = 1.0;     ti1 = 0.5*factt;
         tabli2 = 0.0;  tabli1 = 2.0;
         DFX(0,j)=tabli1;
         for (i=1;i<n;i++)
          { ti = factt*ti1-ti2;
            tabli = 4.0*ti1+factt*tabli1-tabli2;
            DFX(i,j)=tabli;
            ti2 = ti1;        ti1 = ti;
            tabli2 = tabli1;  tabli1 = tabli;
          };
       };
      *iflag=0; return;
}