/*
 
      Testexample for GBIT: Solution of a linear system
      using GBIT.
 
 *  Written by        L. Weimann 
 *  Purpose           Testexample for code GBIT
 *  Version           1.0
 *  Revision          June 2006
 *  Latest Change     June 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 "itlin.h"

void matvec(int n, double *x, double *b);
void precon(int n, double *x, double *b);

static double diag = 6.01;

int main()
{ struct ITLIN_OPT   *opt   = malloc(sizeof(struct ITLIN_OPT));
  struct ITLIN_INFO *info   = malloc(sizeof(struct ITLIN_INFO));
  int  j, rcode;
  double *x, *xsol, *b, *xscale, *res;
  double xh, resnorm, xdnorm;
  int n = 101, maxiter = 1000;  
  LOGICAL rescale = True;
  fprintf(stdout,"n=%i,  maxiter=%i,  rescale=%i,  diag=%e\n",
                 n,maxiter,rescale,diag);
  b      = (double*) malloc( n*sizeof(double) );
  x      = (double*) malloc( n*sizeof(double) );
  xsol   = (double*) malloc( n*sizeof(double) );
  xscale = (double*) malloc( n*sizeof(double) );
  res    = (double*) malloc( n*sizeof(double) );
  if ( !x || !b || !xsol || !xscale || !res )
    { fprintf(stdout,"\n memory allocation error!\n");
      exit(-10);
    };
  fprintf(stdout,"\n Test of GBIT\n\n");
  /* prescribed true solution */
  xh = (double) -(n-1)/2.0;
  for (j=0;j<n;j++)  xsol[j] = xh + j;
  /* compute right hand side b of linear system */
  matvec(n,xsol,b);
  /* initialize start vector for iteration */
  for (j=0;j<n;j++)  x[j] = 0.0;
  /* initialize solver options */
  opt->tol = 1.0e-10;
  opt->rho = 0.25;
  opt->i_max = 10;
  opt->rescale = rescale;
  opt->maxiter = maxiter;
  opt->errorlevel = Verbose;
  opt->monitorlevel = Verbose;
  opt->datalevel = Verbose;
  opt->errorfile = stdout;
  opt->monitorfile = stdout;
  opt->datafile = fopen("test_gbit.data","w");
  if (!opt->datafile) fprintf(stdout,
     "\n open of test_gbit.data failed!\n");
  opt->iterfile = fopen("test_gbit_iter.data","w");
  opt->resfile  = NULL;
  opt->miscfile = fopen("test_gbit_misc.data","w");
  /* call solver */
  gbit(n,x,&matvec,&precon,b,opt,info);
  /* close output files and print statistics */
  fclose(opt->datafile);
  fclose(opt->iterfile);
  fclose(opt->miscfile);
  fprintf(stdout,"\n Return code:         %6i\n",info->rcode);
  fprintf(stdout," Iterations:          %6i\n",info->iter);
  fprintf(stdout," Matvec calls:        %6i\n",info->nomatvec);
  fprintf(stdout," Precon calls:        %6i\n",info->noprecon);
  fprintf(stdout," Estimated Precision: %e\n",info->precision);
  /* compute and print residuum */
  matvec(n,x,res);
  for (j=0;j<n;j++) res[j] = b[j]-res[j];
  resnorm = zibnum_norm2(n,res);
  fprintf(stdout,"\n Norm of residuum: %e\n",resnorm);
  /* compute and print deviation from true solution */
  for (j=0;j<n;j++) xscale[j] = MAX(1.0,fabs(xsol[j]));
  for (j=0;j<n;j++) x[j] -= xsol[j];
  xdnorm = zibnum_scaled_norm2(n,x,xscale);
  fprintf(stdout," True precision:   %e\n",xdnorm);
  return 0;
}

void matvec(int n, double *x, double *b)
{ int j;
  b[0] = diag*x[0]-2*x[1]-x[2];
  b[1] = -2*x[0]+diag*x[1]-2*x[2]-x[3];
  for (j=2;j<n-2;j++) b[j] = -x[j-2]-2*x[j-1]+diag*x[j]-2*x[j+1]-x[j+2];
  b[n-2] = -x[n-4]-2*x[n-3]+diag*x[n-2]-2*x[n-1];
  b[n-1] = -x[n-3]-2*x[n-2]+diag*x[n-1];
}

void precon(int n, double *x, double *b)
{ int j; double div = ( diag != 0.0 ? diag : 1.0);
  for (j=0;j<n;j++) b[j] = x[j]/div;
}