iteration = solve_cg(k,l,q_off,eps_sq);
iteration = solve_cg(k,l,q_off,eps_sq, rel_prec);
int bicg(int k, int l, double q_off, double eps_sq) {    static double rho0,omega0,rho1,omega1,alpha,be,err,d1,d2;  int iteration;  if(ITER_MAX_BCG == 0) {    gamma5(DUM_SOLVER+1,l);    if(g_use_clover_flag == 1){      M_psi(DUM_SOLVER,k,q_off);     }    else {      Mtm_plus_psi(spinor_field[DUM_SOLVER], spinor_field[k]);    }    diff(spinor_field[DUM_SOLVER+1],spinor_field[DUM_SOLVER+1],spinor_field[DUM_SOLVER], VOLUME/2);    assign(spinor_field[DUM_SOLVER], spinor_field[DUM_SOLVER+1], VOLUME/2);    zero_spinor_field(spinor_field[DUM_SOLVER+2]);    zero_spinor_field(spinor_field[DUM_SOLVER+3]);    rho0=1.0; omega0=1.0; alpha=1.0;     /* main loop */    for(iteration=1;iteration<=ITER_MAX_BCG;iteration++) {      /*     if(g_proc_id == 0) printf("%d %e\n", iteration, rho1); */      if(rho1 == 0. && iteration > 1){	iteration = ITER_MAX_BCG+1;	break;      }      square_and_prod_r(&err,&rho1, spinor_field[DUM_SOLVER+1],  spinor_field[DUM_SOLVER], VOLUME/2);      /*     if(g_proc_id == 0) printf("%d %e %e\n", iteration, err, rho1); */      if(err <= eps_sq){	break;      }            be = (rho1/rho0)*(alpha/omega0);      /*  assign_add_mul_r(DUM_SOLVER+3,-omega0,DUM_SOLVER+2, VOLUME/2);	  assign_mul_add_r(DUM_SOLVER+3,be,DUM_SOLVER+1, VOLUME/2);  */      assign_mul_bra_add_mul_ket_add_r(spinor_field[DUM_SOLVER+3], spinor_field[DUM_SOLVER+2], spinor_field[DUM_SOLVER+1],				       -omega0, be, VOLUME/2);       if(g_use_clover_flag == 1){	M_psi(DUM_SOLVER+2, DUM_SOLVER+3, q_off);       }      else{	Mtm_plus_psi(spinor_field[DUM_SOLVER+2], spinor_field[DUM_SOLVER+3]);      }      alpha=rho1/scalar_prod_r(spinor_field[DUM_SOLVER], spinor_field[DUM_SOLVER+2], VOLUME/2);      assign_add_mul_r(spinor_field[DUM_SOLVER+1], spinor_field[DUM_SOLVER+2], -alpha, VOLUME/2);      if(g_use_clover_flag == 1){	M_psi(DUM_SOLVER+4, DUM_SOLVER+1, q_off);      }      else{	Mtm_plus_psi(spinor_field[DUM_SOLVER+4], spinor_field[DUM_SOLVER+1]);      }      square_and_prod_r(&d1,&d2,  spinor_field[DUM_SOLVER+4],  spinor_field[DUM_SOLVER+1], VOLUME/2);      omega1=d2/d1;      assign_add_mul_r_add_mul(spinor_field[k], spinor_field[DUM_SOLVER+3], spinor_field[DUM_SOLVER+1], alpha, omega1, VOLUME/2);       assign_add_mul_r(spinor_field[DUM_SOLVER+1], spinor_field[DUM_SOLVER+4], -omega1, VOLUME/2);      /*copy back */      rho0=rho1; omega0=omega1;    }    /* if bicg fails, redo with conjugate gradient */    if(g_proc_id==0) {      sout = fopen(solvout, "a");      fprintf(sout, "BiCGstab: iterations: %d mu: %f eps_sq: %e\n", iteration, g_mu, eps_sq);       fclose(sout);    }  }  else iteration = ITER_MAX_BCG;  if(iteration>=ITER_MAX_BCG){/*     zero_spinor_field(spinor_field[k]); */    iteration = solve_cg(k,l,q_off,eps_sq);    assign(spinor_field[DUM_DERI+6], spinor_field[k], VOLUME/2);    if(g_use_clover_flag == 1){      Q_psi(k,k,q_off);    }    else{      Qtm_minus_psi(spinor_field[k], spinor_field[k]);;    }  }  return iteration;}