/*
 
      Testexample for GMRES: Solution of a linear system
      using GMRES.
 
 *  Written by        L. Weimann 
 *  Purpose           Testexample for code GMRES
 *  Version           1.0
 *  Revision          December 2006
 *  Latest Change     December 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 preconr(int n, double *x, double *b);
double gmres_norm2(int n, double *v);

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, *w, *res;
  double xh, resnorm, xdnorm;
  int n = 101, maxiter = 1000;  
  TERM_CHECK termcheck = CheckEachIter;
  fprintf(stdout,"n=%i,  maxiter=%i,  termcheck=%i,  diag=%e\n",
                 n,maxiter,termcheck,diag);
  b      = (double*) malloc( n*sizeof(double) );
  x      = (double*) malloc( n*sizeof(double) );
  xsol   = (double*) malloc( n*sizeof(double) );
  w      = (double*) malloc( n*sizeof(double) );
  res    = (double*) malloc( n*sizeof(double) );
  if ( !x || !b || !xsol || !w || !res )
    { fprintf(stdout,"\n memory allocation error!\n");
      exit(-10);
    };
  fprintf(stdout,"\n Test of GMRES\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] = b[j];
  /* initialize solver options */
  opt->tol = 1.0e-10;
  opt->i_max = 10;
  opt->termcheck = termcheck;
  opt->maxiter = maxiter;
  opt->errorlevel = Verbose;
  opt->monitorlevel = Verbose;
  opt->datalevel = Verbose;
  opt->errorfile = stdout;
  opt->monitorfile = stdout;
  opt->datafile = fopen("test_gmres.data","w");
  if (!opt->datafile) fprintf(stdout,
     "\n open of test_gmres.data failed!\n");
  opt->iterfile = fopen("test_gmres_iter.data","w");
  opt->resfile  = fopen("test_gmres_res.data","w");
  opt->miscfile = fopen("test_gmres_misc.data","w");
  /* call solver */
  gmres(n,x,&matvec,&preconr,NULL,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," Right Precon calls:           %6i\n",info->noprecr);
  fprintf(stdout," Left Precon calls:            %6i\n",info->noprecl);
  fprintf(stdout," Estimated residual reduction: %e\n",info->precision);
  /* compute and print residuum */
  preconr(n,x,w);
  matvec(n,w,res);
  for (j=0;j<n;j++) res[j] = b[j]-res[j];
  resnorm = gmres_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++) w[j] -= xsol[j];
  xdnorm = zibnum_norm2(n,w);
  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 preconr(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;
}