/**
# Helper macro to invert (linear) spatial operators

The macro below can be used to easily invert linear systems described
by stencils i.e.
$$
\mathcal{L}(a) = b
$$
For example, let us consider the Poisson equation
$$
\nabla^2 a = b
$$
where $a$ is unknown and $b$ is given. This can be discretised as

~~~literatec
(a[1] + a[-1] + a[0,1] + a[0,-1] - 4.*a[])/sq(Delta) = b[];
~~~

This can be solved using the macro below and a [multigrid
solver](poisson.h#mg_solve) with

~~~literatec
solve (a, (a[1] + a[-1] + a[0,1] + a[0,-1] - 4.*a[])/sq(Delta), b);
~~~

The macro can take the same optional arguments as
[mg_solve()](poisson.h#mg_solve) to tune the multigrid solver.

The macro returns [multigrid statistics](poisson.h#mgstats). */

#include "poisson.h"

/**
## Implementation

Note that the macro below is a slightly simplified version of
the [mg_solve()](poisson.h#mg_solve) and
[mg_cycle()](poisson.h#mg_cycle) functions where more
documentation can be found. */

macro
mgstats solve (scalar a, double func, double rhs,
	       int nrelax = 4,
	       int minlevel = 0,
	       double tolerance = TOLERANCE)
{{
  mgstats _s = (mgstats){0};
  scalar _res[], _da[];
  scalar_clone (_da, a);
  for (int b = 0; b < nboundary; b++)
    _da.boundary[b] = _da.boundary_homogeneous[b];
  _s.nrelax = nrelax;
  double _resb;
  {
    double maxres = 0.;
    foreach (reduction(max:maxres)) {
      _res[] = rhs - func;
      if (fabs (_res[]) > maxres)
	maxres = fabs (_res[]);
    }
    _resb = _s.resb = _s.resa = maxres;
  }
  for (_s.i = 0; _s.i < NITERMAX && (_s.i < NITERMIN || _s.resa > tolerance); _s.i++) {
    {
      restriction ({_res});
      int _maxlevel = grid->maxdepth;
      int _minlevel = min (minlevel, _maxlevel);
      for (int l = _minlevel; l <= _maxlevel; l++) {
	if (l == _minlevel)
	  foreach_level_or_leaf (l)
	    foreach_blockf (_da)
	      _da[] = 0.;
	else
	  foreach_level (l)
	    foreach_blockf (_da)
	    _da[] = bilinear (point, _da);
	boundary_level ({_da}, l);
	for (int i = 0; i < _s.nrelax; i++) {
	  scalar a = _da;
	  foreach_level_or_leaf (l) {
	    a[] = 0.;
	    double _n = _res[] - func, _d;
	    diagonalize(a)
	      _d = func;
	    a[] = _n/_d;
	  }
	  boundary_level ({_da}, l);
	}
      }      
      foreach()
	foreach_blockf (a)
	  a[] += _da[];
    }
    {
      double maxres = 0.;
      foreach (reduction(max:maxres)) {
	_res[] = rhs - func;
	if (fabs (_res[]) > maxres)
	  maxres = fabs (_res[]);
      }
      _s.resa = maxres;
    }
    if (_s.resa > tolerance) {
      if (_resb/_s.resa < 1.2 && _s.nrelax < 100)
	_s.nrelax++;
      else if (_resb/_s.resa > 10 && _s.nrelax > 2)
	_s.nrelax--;
    }
    _resb = _s.resa;
  }
  _s.minlevel = minlevel;
  if (_s.resa > tolerance)
    fprintf (stderr,
	     "src/solve.h:%d: warning: convergence for %s not reached after %d iterations\n"
	     "  res: %g nrelax: %d\n", LINENO, a.name,
	     _s.i, _s.resa, _s.nrelax),
      fflush (ferr);
  return _s;
}}