_SO(if(g_proc_id == g_stdio_proc){printf("%d\t%g\n", restart*m+j, alpha[j+1].re*alpha[j+1].re); fflush(stdout);});
if(g_proc_id == g_stdio_proc){ printf("%d\t%g residue\n", restart*m+j, alpha[j+1].re*alpha[j+1].re); fflush(stdout); }
int gmres(spinor * const P,spinor * const Q, 	   const int m, const int max_restarts,	   const double eps_sq, matrix_mult f){  int restart, i, j, k;  double beta, eps;  complex tmp1, tmp2;  int N=VOLUME/2;/*   init_solver_field(3); */  eps=sqrt(eps_sq);  init_gmres(m, VOLUMEPLUSRAND);  assign(spinor_field[DUM_SOLVER+2], P, N);  for(restart = 0; restart < max_restarts; restart++){    /* r_0=Q-AP  (b=Q, x+0=P) */    f(spinor_field[DUM_SOLVER], spinor_field[DUM_SOLVER+2]);    diff(spinor_field[DUM_SOLVER], Q, spinor_field[DUM_SOLVER], N);    /* v_0=r_0/||r_0|| */    alpha[0].re=sqrt(square_norm(spinor_field[DUM_SOLVER], N));    if(alpha[0].re==0.){      assign(P, spinor_field[DUM_SOLVER+2], N);      return(restart*m);    }    mul_r(V[0], 1./alpha[0].re, spinor_field[DUM_SOLVER], N);    for(j = 0; j < m; j++){      /* spinor_field[DUM_SOLVER]=A*v_j */      f(spinor_field[DUM_SOLVER], V[j]);      /* Set h_ij and omega_j */      /* spinor_field[DUM_SOLVER+1] <- omega_j */      assign(spinor_field[DUM_SOLVER+1], spinor_field[DUM_SOLVER], N);      for(i = 0; i <= j; i++){	H[i][j] = scalar_prod(V[i], spinor_field[DUM_SOLVER+1], N);	assign_diff_mul(spinor_field[DUM_SOLVER+1], V[i], H[i][j], N);      }      _complex_set(H[j+1][j], sqrt(square_norm(spinor_field[DUM_SOLVER+1], N)), 0.);      for(i = 0; i < j; i++){	tmp1 = H[i][j];	tmp2 = H[i+1][j];	_mult_real(H[i][j], tmp2, s[i]);	_add_assign_complex_conj(H[i][j], c[i], tmp1);	_mult_real(H[i+1][j], tmp1, s[i]);	_diff_assign_complex(H[i+1][j], c[i], tmp2);      }      /* Set beta, s, c, alpha[j],[j+1] */      beta = sqrt(_complex_square_norm(H[j][j]) + _complex_square_norm(H[j+1][j]));      s[j] = H[j+1][j].re / beta;      _mult_real(c[j], H[j][j], 1./beta);      _complex_set(H[j][j], beta, 0.);      _mult_real(alpha[j+1], alpha[j], s[j]);      tmp1 = alpha[j];      _mult_assign_complex_conj(alpha[j], c[j], tmp1);      /* precision reached? */      _SO(if(g_proc_id == g_stdio_proc){printf("%d\t%g\n", restart*m+j, alpha[j+1].re*alpha[j+1].re); fflush(stdout);});      if(alpha[j+1].re <= eps){	_mult_real(alpha[j], alpha[j], 1./H[j][j].re);	assign_add_mul(spinor_field[DUM_SOLVER+2], V[j], alpha[j], N);	for(i = j-1; i >= 0; i--){	  for(k = i+1; k <= j; k++){ 	    _mult_assign_complex(tmp1, H[i][k], alpha[k]); 	    _diff_complex(alpha[i], tmp1);	  }	  _mult_real(alpha[i], alpha[i], 1./H[i][i].re);	  assign_add_mul(spinor_field[DUM_SOLVER+2], V[i], alpha[i], N);	}	for(i = 0; i < m; i++){	  alpha[i].im = 0.;	}	assign(P, spinor_field[DUM_SOLVER+2], N);	return(restart*m+j);      }      /* if not */      else{	if(j != m-1){	  mul_r(V[(j+1)], 1./H[j+1][j].re, spinor_field[DUM_SOLVER+1], N);	}      }    }    j=m-1;    /* prepare for restart */    _mult_real(alpha[j], alpha[j], 1./H[j][j].re);    assign_add_mul(spinor_field[DUM_SOLVER+2], V[j], alpha[j], N);    for(i = j-1; i >= 0; i--){      for(k = i+1; k <= j; k++){	_mult_assign_complex(tmp1, H[i][k], alpha[k]);	_diff_complex(alpha[i], tmp1);      }      _mult_real(alpha[i], alpha[i], 1./H[i][i].re);      assign_add_mul(spinor_field[DUM_SOLVER+2], V[i], alpha[i], N);    }    for(i = 0; i < m; i++){      alpha[i].im = 0.;    }  }  /* If maximal number of restarts is reached */  assign(P, spinor_field[DUM_SOLVER+2], N);  return(-1);}