set_params.cpp

Go to the documentation of this file.

#include <algorithm>
#include <command_line_params.h>
#include <host_utils.h>
#include "misc.h"
 
void setGaugeParam(QudaGaugeParam &gauge_param)
{
  gauge_param.type = QUDA_SU3_LINKS;
 
  gauge_param.X[0] = xdim;
  gauge_param.X[1] = ydim;
  gauge_param.X[2] = zdim;
  gauge_param.X[3] = tdim;
 
  gauge_param.cpu_prec = cpu_prec;
  gauge_param.cuda_prec = cuda_prec;
  gauge_param.cuda_prec_sloppy = cuda_prec;
  gauge_param.cuda_prec_precondition = cuda_prec;
  gauge_param.cuda_prec_eigensolver = cuda_prec;
 
  gauge_param.reconstruct = link_recon;
  gauge_param.reconstruct_sloppy = link_recon;
  gauge_param.reconstruct_precondition = link_recon;
  gauge_param.reconstruct_eigensolver = link_recon;
  gauge_param.reconstruct_refinement_sloppy = link_recon;
 
  gauge_param.anisotropy = 1.0;
  gauge_param.tadpole_coeff = 1.0;
 
  gauge_param.ga_pad = 0;
  gauge_param.mom_ga_pad = 0;
  gauge_param.gauge_fix = QUDA_GAUGE_FIXED_NO;
 
  gauge_param.struct_size = sizeof(gauge_param);
}
 
void setWilsonGaugeParam(QudaGaugeParam &gauge_param)
{
  setGaugeParam(gauge_param);
  gauge_param.anisotropy = anisotropy;
  gauge_param.type = QUDA_WILSON_LINKS;
  gauge_param.gauge_order = QUDA_QDP_GAUGE_ORDER;
  gauge_param.t_boundary = fermion_t_boundary;
 
  gauge_param.cuda_prec_sloppy = cuda_prec_sloppy;
  gauge_param.cuda_prec_precondition = cuda_prec_precondition;
  gauge_param.cuda_prec_eigensolver = cuda_prec_eigensolver;
  gauge_param.cuda_prec_refinement_sloppy = cuda_prec_refinement_sloppy;
 
  gauge_param.reconstruct_sloppy = link_recon_sloppy;
  gauge_param.reconstruct_precondition = link_recon_precondition;
  gauge_param.reconstruct_eigensolver = link_recon_eigensolver;
  gauge_param.reconstruct_refinement_sloppy = link_recon_sloppy;
 
  int pad_size = 0;
  // For multi-GPU, ga_pad must be large enough to store a time-slice
#ifdef MULTI_GPU
  int x_face_size = gauge_param.X[1] * gauge_param.X[2] * gauge_param.X[3] / 2;
  int y_face_size = gauge_param.X[0] * gauge_param.X[2] * gauge_param.X[3] / 2;
  int z_face_size = gauge_param.X[0] * gauge_param.X[1] * gauge_param.X[3] / 2;
  int t_face_size = gauge_param.X[0] * gauge_param.X[1] * gauge_param.X[2] / 2;
  pad_size = std::max({x_face_size, y_face_size, z_face_size, t_face_size});
#endif
  gauge_param.ga_pad = pad_size;
 
  gauge_param.struct_size = sizeof(gauge_param);
}
 
void setStaggeredGaugeParam(QudaGaugeParam &gauge_param)
{
  setGaugeParam(gauge_param);
 
  gauge_param.cuda_prec_sloppy = prec_sloppy;
  gauge_param.cuda_prec_refinement_sloppy = prec_refinement_sloppy;
  gauge_param.cuda_prec_precondition = prec_precondition;
  gauge_param.cuda_prec_eigensolver = prec_eigensolver;
  gauge_param.reconstruct_sloppy = link_recon_sloppy;
  gauge_param.reconstruct_precondition = link_recon_precondition;
  gauge_param.reconstruct_eigensolver = link_recon_eigensolver;
  gauge_param.reconstruct_refinement_sloppy = link_recon_sloppy;
 
  // For HISQ, this must always be set to 1.0, since the tadpole
  // correction is baked into the coefficients for the first fattening.
  // The tadpole doesn't mean anything for the second fattening
  // since the input fields are unitarized.
  gauge_param.tadpole_coeff = 1.0;
 
  if (dslash_type == QUDA_ASQTAD_DSLASH) {
    gauge_param.scale = -1.0 / 24.0;
    if (eps_naik != 0) { gauge_param.scale *= (1.0 + eps_naik); }
  } else {
    gauge_param.scale = 1.0;
  }
 
  gauge_param.gauge_order = QUDA_MILC_GAUGE_ORDER;
  gauge_param.t_boundary = fermion_t_boundary;
  gauge_param.staggered_phase_type = QUDA_STAGGERED_PHASE_MILC;
  gauge_param.type = QUDA_WILSON_LINKS;
 
  int pad_size = 0;
#ifdef MULTI_GPU
  int x_face_size = gauge_param.X[1] * gauge_param.X[2] * gauge_param.X[3] / 2;
  int y_face_size = gauge_param.X[0] * gauge_param.X[2] * gauge_param.X[3] / 2;
  int z_face_size = gauge_param.X[0] * gauge_param.X[1] * gauge_param.X[3] / 2;
  int t_face_size = gauge_param.X[0] * gauge_param.X[1] * gauge_param.X[2] / 2;
  pad_size = std::max({x_face_size, y_face_size, z_face_size, t_face_size});
#endif
  gauge_param.ga_pad = pad_size;
 
  gauge_param.struct_size = sizeof(gauge_param);
}
 
void setInvertParam(QudaInvertParam &inv_param)
{
  // Set dslash type
  inv_param.dslash_type = dslash_type;
 
  // Use kappa or mass normalisation
  if (kappa == -1.0) {
    inv_param.mass = mass;
    inv_param.kappa = 1.0 / (2.0 * (1 + 3 / anisotropy + mass));
    if (dslash_type == QUDA_LAPLACE_DSLASH) inv_param.kappa = 1.0 / (8 + mass);
  } else {
    inv_param.kappa = kappa;
    inv_param.mass = 0.5 / kappa - (1.0 + 3.0 / anisotropy);
    if (dslash_type == QUDA_LAPLACE_DSLASH) inv_param.mass = 1.0 / kappa - 8.0;
  }
  printfQuda("Kappa = %.8f Mass = %.8f\n", inv_param.kappa, inv_param.mass);
 
  // Use 3D or 4D laplace
  inv_param.laplace3D = laplace3D;
 
  // Some fermion specific parameters
  if (dslash_type == QUDA_TWISTED_MASS_DSLASH || dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
    inv_param.mu = mu;
    inv_param.epsilon = epsilon;
    inv_param.twist_flavor = twist_flavor;
    inv_param.Ls = (inv_param.twist_flavor == QUDA_TWIST_NONDEG_DOUBLET) ? 2 : 1;
  } else if (dslash_type == QUDA_DOMAIN_WALL_DSLASH || dslash_type == QUDA_DOMAIN_WALL_4D_DSLASH
             || dslash_type == QUDA_MOBIUS_DWF_DSLASH || dslash_type == QUDA_MOBIUS_DWF_EOFA_DSLASH) {
    inv_param.m5 = m5;
    kappa5 = 0.5 / (5 + inv_param.m5);
    inv_param.Ls = Lsdim;
    for (int k = 0; k < Lsdim; k++) { // for mobius only
      // b5[k], c[k] values are chosen for arbitrary values,
      // but the difference of them are same as 1.0
      inv_param.b_5[k] = b5;
      inv_param.c_5[k] = c5;
    }
    inv_param.eofa_pm = eofa_pm;
    inv_param.eofa_shift = eofa_shift;
    inv_param.mq1 = eofa_mq1;
    inv_param.mq2 = eofa_mq2;
    inv_param.mq3 = eofa_mq3;
  } else {
    inv_param.Ls = 1;
  }
 
  // Set clover specific parameters
  if (dslash_type == QUDA_CLOVER_WILSON_DSLASH || dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
    inv_param.clover_cpu_prec = cpu_prec;
    inv_param.clover_cuda_prec = cuda_prec;
    inv_param.clover_cuda_prec_sloppy = cuda_prec_sloppy;
    inv_param.clover_cuda_prec_precondition = cuda_prec_precondition;
    inv_param.clover_cuda_prec_eigensolver = cuda_prec_eigensolver;
    inv_param.clover_cuda_prec_refinement_sloppy = cuda_prec_sloppy;
    inv_param.clover_order = QUDA_PACKED_CLOVER_ORDER;
    // Use kappa * csw or supplied clover_coeff
    inv_param.clover_csw = clover_csw;
    if (clover_coeff == 0.0) {
      inv_param.clover_coeff = clover_csw * inv_param.kappa;
    } else {
      inv_param.clover_coeff = clover_coeff;
    }
  }
 
  // General parameter setup
  inv_param.inv_type = inv_type;
  inv_param.solution_type = solution_type;
  inv_param.solve_type = solve_type;
  inv_param.matpc_type = matpc_type;
  inv_param.dagger = QUDA_DAG_NO;
  inv_param.mass_normalization = normalization;
  inv_param.solver_normalization = QUDA_DEFAULT_NORMALIZATION;
  inv_param.pipeline = pipeline;
  inv_param.Nsteps = 2;
  inv_param.gcrNkrylov = gcrNkrylov;
  inv_param.ca_basis = ca_basis;
  inv_param.ca_lambda_min = ca_lambda_min;
  inv_param.ca_lambda_max = ca_lambda_max;
  inv_param.tol = tol;
  inv_param.tol_restart = tol_restart;
  if (tol_hq == 0 && tol == 0) {
    errorQuda("qudaInvert: requesting zero residual\n");
    exit(1);
  }
 
  // require both L2 relative and heavy quark residual to determine convergence
  inv_param.residual_type = static_cast<QudaResidualType_s>(0);
  inv_param.residual_type = (tol != 0) ?
    static_cast<QudaResidualType_s>(inv_param.residual_type | QUDA_L2_RELATIVE_RESIDUAL) :
    inv_param.residual_type;
  inv_param.residual_type = (tol_hq != 0) ?
    static_cast<QudaResidualType_s>(inv_param.residual_type | QUDA_HEAVY_QUARK_RESIDUAL) :
    inv_param.residual_type;
 
  inv_param.tol_hq = tol_hq; // specify a tolerance for the residual for heavy quark residual
 
  // Offsets used only by multi-shift solver
  // These should be set in the application code. We set the them here by way of
  // example
  inv_param.num_offset = multishift;
  for (int i = 0; i < inv_param.num_offset; i++) inv_param.offset[i] = 0.06 + i * i * 0.1;
  // these can be set individually
  for (int i = 0; i < inv_param.num_offset; i++) {
    inv_param.tol_offset[i] = inv_param.tol;
    inv_param.tol_hq_offset[i] = inv_param.tol_hq;
  }
  inv_param.maxiter = niter;
  inv_param.reliable_delta = reliable_delta;
  inv_param.use_alternative_reliable = alternative_reliable;
  inv_param.use_sloppy_partial_accumulator = 0;
  inv_param.solution_accumulator_pipeline = solution_accumulator_pipeline;
  inv_param.max_res_increase = 1;
 
  // domain decomposition preconditioner parameters
  inv_param.inv_type_precondition = precon_type;
 
  inv_param.schwarz_type = precon_schwarz_type;
  inv_param.precondition_cycle = precon_schwarz_cycle;
  inv_param.tol_precondition = tol_precondition;
  inv_param.maxiter_precondition = maxiter_precondition;
  inv_param.verbosity_precondition = mg_verbosity[0];
  inv_param.cuda_prec_precondition = cuda_prec_precondition;
  inv_param.cuda_prec_eigensolver = cuda_prec_eigensolver;
  inv_param.omega = 1.0;
 
  inv_param.cpu_prec = cpu_prec;
  inv_param.cuda_prec = cuda_prec;
  inv_param.cuda_prec_sloppy = cuda_prec_sloppy;
  inv_param.cuda_prec_refinement_sloppy = cuda_prec_refinement_sloppy;
  inv_param.preserve_source = QUDA_PRESERVE_SOURCE_YES;
  inv_param.gamma_basis = QUDA_DEGRAND_ROSSI_GAMMA_BASIS;
  inv_param.dirac_order = QUDA_DIRAC_ORDER;
 
  inv_param.input_location = QUDA_CPU_FIELD_LOCATION;
  inv_param.output_location = QUDA_CPU_FIELD_LOCATION;
 
  inv_param.sp_pad = 0;
  inv_param.cl_pad = 0;
 
  inv_param.verbosity = verbosity;
 
  inv_param.extlib_type = solver_ext_lib;
 
  // Whether or not to use native BLAS LAPACK
  inv_param.native_blas_lapack = (native_blas_lapack ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE);
 
  inv_param.struct_size = sizeof(inv_param);
}
 
// Parameters defining the eigensolver
void setEigParam(QudaEigParam &eig_param)
{
  eig_param.eig_type = eig_type;
  eig_param.spectrum = eig_spectrum;
  if ((eig_type == QUDA_EIG_TR_LANCZOS || eig_type == QUDA_EIG_BLK_TR_LANCZOS)
      && !(eig_spectrum == QUDA_SPECTRUM_LR_EIG || eig_spectrum == QUDA_SPECTRUM_SR_EIG)) {
    errorQuda("Only real spectrum type (LR or SR) can be passed to Lanczos type solver.");
  }
 
  // The solver will exit when n_conv extremal eigenpairs have converged
  if (eig_n_conv < 0) {
    eig_param.n_conv = eig_n_ev;
    eig_n_conv = eig_n_ev;
  } else {
    eig_param.n_conv = eig_n_conv;
  }
 
  // Inverters will deflate only this number of vectors.
  if (eig_n_ev_deflate < 0) {
    eig_param.n_ev_deflate = eig_n_conv;
    eig_n_ev_deflate = eig_n_conv;
  } else {
    if (eig_n_ev_deflate > eig_n_conv) errorQuda("Can not deflate more that eig_n_conv eigenvectors.");
    eig_param.n_ev_deflate = eig_n_ev_deflate;
  }
 
  eig_param.block_size
    = (eig_param.eig_type == QUDA_EIG_TR_LANCZOS || eig_param.eig_type == QUDA_EIG_IR_ARNOLDI) ? 1 : eig_block_size;
  eig_param.n_ev = eig_n_ev;
  eig_param.n_kr = eig_n_kr;
  eig_param.tol = eig_tol;
  eig_param.qr_tol = eig_qr_tol;
  eig_param.batched_rotate = eig_batched_rotate;
  eig_param.require_convergence = eig_require_convergence ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  eig_param.check_interval = eig_check_interval;
  eig_param.max_restarts = eig_max_restarts;
 
  eig_param.use_norm_op = eig_use_normop ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  eig_param.use_dagger = eig_use_dagger ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  eig_param.compute_gamma5 = eig_compute_gamma5 ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  eig_param.compute_svd = eig_compute_svd ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  if (eig_compute_svd) {
    eig_param.use_dagger = QUDA_BOOLEAN_FALSE;
    eig_param.use_norm_op = QUDA_BOOLEAN_TRUE;
  }
 
  eig_param.use_eigen_qr = eig_use_eigen_qr ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  eig_param.use_poly_acc = eig_use_poly_acc ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  eig_param.poly_deg = eig_poly_deg;
  eig_param.a_min = eig_amin;
  eig_param.a_max = eig_amax;
 
  eig_param.arpack_check = eig_arpack_check ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  strcpy(eig_param.arpack_logfile, eig_arpack_logfile);
 
  strcpy(eig_param.vec_infile, eig_vec_infile);
  strcpy(eig_param.vec_outfile, eig_vec_outfile);
  eig_param.save_prec = eig_save_prec;
  eig_param.io_parity_inflate = eig_io_parity_inflate ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  eig_param.struct_size = sizeof(eig_param);
}
 
void setMultigridParam(QudaMultigridParam &mg_param)
{
  QudaInvertParam &inv_param = *mg_param.invert_param; // this will be used to setup SolverParam parent in MGParam class
 
  // Whether or not to use native BLAS LAPACK
  inv_param.native_blas_lapack = (native_blas_lapack ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE);
 
  inv_param.Ls = 1;
 
  inv_param.sp_pad = 0;
  inv_param.cl_pad = 0;
 
  inv_param.cpu_prec = cpu_prec;
  inv_param.cuda_prec = cuda_prec;
  inv_param.cuda_prec_sloppy = cuda_prec_sloppy;
  inv_param.cuda_prec_precondition = cuda_prec_precondition;
  inv_param.cuda_prec_eigensolver = cuda_prec_eigensolver;
  inv_param.preserve_source = QUDA_PRESERVE_SOURCE_NO;
  inv_param.gamma_basis = QUDA_DEGRAND_ROSSI_GAMMA_BASIS;
  inv_param.dirac_order = QUDA_DIRAC_ORDER;
 
  if (kappa == -1.0) {
    inv_param.mass = mass;
    inv_param.kappa = 1.0 / (2.0 * (1 + 3 / anisotropy + mass));
  } else {
    inv_param.kappa = kappa;
    inv_param.mass = 0.5 / kappa - (1 + 3 / anisotropy);
  }
 
  if (dslash_type == QUDA_CLOVER_WILSON_DSLASH || dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
    inv_param.clover_cpu_prec = cpu_prec;
    inv_param.clover_cuda_prec = cuda_prec;
    inv_param.clover_cuda_prec_sloppy = cuda_prec_sloppy;
    inv_param.clover_cuda_prec_precondition = cuda_prec_precondition;
    inv_param.clover_cuda_prec_eigensolver = cuda_prec_eigensolver;
    inv_param.clover_cuda_prec_refinement_sloppy = cuda_prec_sloppy;
    inv_param.clover_order = QUDA_PACKED_CLOVER_ORDER;
    // Use kappa * csw or supplied clover_coeff
    inv_param.clover_csw = clover_csw;
    if (clover_coeff == 0.0) {
      inv_param.clover_coeff = clover_csw * inv_param.kappa;
    } else {
      inv_param.clover_coeff = clover_coeff;
    }
  }
 
  inv_param.input_location = QUDA_CPU_FIELD_LOCATION;
  inv_param.output_location = QUDA_CPU_FIELD_LOCATION;
 
  inv_param.dslash_type = dslash_type;
 
  if (dslash_type == QUDA_TWISTED_MASS_DSLASH || dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
    inv_param.mu = mu;
    inv_param.epsilon = epsilon;
    inv_param.twist_flavor = twist_flavor;
    inv_param.Ls = (inv_param.twist_flavor == QUDA_TWIST_NONDEG_DOUBLET) ? 2 : 1;
 
    if (twist_flavor == QUDA_TWIST_NONDEG_DOUBLET) {
      printfQuda("Twisted-mass doublet non supported (yet)\n");
      exit(0);
    }
  }
 
  inv_param.dagger = QUDA_DAG_NO;
  inv_param.mass_normalization = QUDA_KAPPA_NORMALIZATION;
 
  inv_param.matpc_type = matpc_type;
  inv_param.solution_type = QUDA_MAT_SOLUTION;
 
  inv_param.solve_type = QUDA_DIRECT_SOLVE;
 
  mg_param.invert_param = &inv_param;
  mg_param.n_level = mg_levels;
  for (int i = 0; i < mg_param.n_level; i++) {
    for (int j = 0; j < 4; j++) {
      // if not defined use 4
      mg_param.geo_block_size[i][j] = geo_block_size[i][j] ? geo_block_size[i][j] : 4;
    }
    for (int j = 4; j < QUDA_MAX_DIM; j++) mg_param.geo_block_size[i][j] = 1;
    mg_param.use_eig_solver[i] = mg_eig[i] ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
    mg_param.verbosity[i] = mg_verbosity[i];
    mg_param.setup_inv_type[i] = setup_inv[i];
    mg_param.num_setup_iter[i] = num_setup_iter[i];
    mg_param.setup_tol[i] = setup_tol[i];
    mg_param.setup_maxiter[i] = setup_maxiter[i];
    mg_param.setup_maxiter_refresh[i] = setup_maxiter_refresh[i];
 
    // Basis to use for CA-CGN(E/R) setup
    mg_param.setup_ca_basis[i] = setup_ca_basis[i];
 
    // Basis size for CACG setup
    mg_param.setup_ca_basis_size[i] = setup_ca_basis_size[i];
 
    // Minimum and maximum eigenvalue for Chebyshev CA basis setup
    mg_param.setup_ca_lambda_min[i] = setup_ca_lambda_min[i];
    mg_param.setup_ca_lambda_max[i] = setup_ca_lambda_max[i];
 
    mg_param.spin_block_size[i] = 1;
    mg_param.n_vec[i] = nvec[i] == 0 ? 24 : nvec[i];          // default to 24 vectors if not set
    mg_param.n_block_ortho[i] = n_block_ortho[i];             // number of times to Gram-Schmidt
    mg_param.precision_null[i] = prec_null;                   // precision to store the null-space basis
    mg_param.smoother_halo_precision[i] = smoother_halo_prec; // precision of the halo exchange in the smoother
    mg_param.nu_pre[i] = nu_pre[i];
    mg_param.nu_post[i] = nu_post[i];
    mg_param.mu_factor[i] = mu_factor[i];
 
    mg_param.cycle_type[i] = QUDA_MG_CYCLE_RECURSIVE;
 
    // Is not a staggered solve, always aggregate
    mg_param.transfer_type[i] = QUDA_TRANSFER_AGGREGATE;
 
    // set the coarse solver wrappers including bottom solver
    mg_param.coarse_solver[i] = coarse_solver[i];
    mg_param.coarse_solver_tol[i] = coarse_solver_tol[i];
    mg_param.coarse_solver_maxiter[i] = coarse_solver_maxiter[i];
 
    // Basis to use for CA-CGN(E/R) coarse solver
    mg_param.coarse_solver_ca_basis[i] = coarse_solver_ca_basis[i];
 
    // Basis size for CACG coarse solver/
    mg_param.coarse_solver_ca_basis_size[i] = coarse_solver_ca_basis_size[i];
 
    // Minimum and maximum eigenvalue for Chebyshev CA basis
    mg_param.coarse_solver_ca_lambda_min[i] = coarse_solver_ca_lambda_min[i];
    mg_param.coarse_solver_ca_lambda_max[i] = coarse_solver_ca_lambda_max[i];
 
    mg_param.smoother[i] = smoother_type[i];
 
    // set the smoother / bottom solver tolerance (for MR smoothing this will be ignored)
    mg_param.smoother_tol[i] = smoother_tol[i];
 
    // set to QUDA_DIRECT_SOLVE for no even/odd preconditioning on the smoother
    // set to QUDA_DIRECT_PC_SOLVE for to enable even/odd preconditioning on the smoother
    mg_param.smoother_solve_type[i] = smoother_solve_type[i];
 
    // set to QUDA_ADDITIVE_SCHWARZ for Additive Schwarz precondioned smoother (presently only impelemented for MR)
    mg_param.smoother_schwarz_type[i] = mg_schwarz_type[i];
 
    // if using Schwarz preconditioning then use local reductions only
    mg_param.global_reduction[i] = (mg_schwarz_type[i] == QUDA_INVALID_SCHWARZ) ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
    // set number of Schwarz cycles to apply
    mg_param.smoother_schwarz_cycle[i] = mg_schwarz_cycle[i];
 
    // Set set coarse_grid_solution_type: this defines which linear
    // system we are solving on a given level
    // * QUDA_MAT_SOLUTION - we are solving the full system and inject
    //   a full field into coarse grid
    // * QUDA_MATPC_SOLUTION - we are solving the e/o-preconditioned
    //   system, and only inject single parity field into coarse grid
    //
    // Multiple possible scenarios here
    //
    // 1. **Direct outer solver and direct smoother**: here we use
    // full-field residual coarsening, and everything involves the
    // full system so coarse_grid_solution_type = QUDA_MAT_SOLUTION
    //
    // 2. **Direct outer solver and preconditioned smoother**: here,
    // only the smoothing uses e/o preconditioning, so
    // coarse_grid_solution_type = QUDA_MAT_SOLUTION_TYPE.
    // We reconstruct the full residual prior to coarsening after the
    // pre-smoother, and then need to project the solution for post
    // smoothing.
    //
    // 3. **Preconditioned outer solver and preconditioned smoother**:
    // here we use single-parity residual coarsening throughout, so
    // coarse_grid_solution_type = QUDA_MATPC_SOLUTION.  This is a bit
    // questionable from a theoretical point of view, since we don't
    // coarsen the preconditioned operator directly, rather we coarsen
    // the full operator and preconditioned that, but it just works.
    // This is the optimal combination in general for Wilson-type
    // operators: although there is an occasional increase in
    // iteration or two), by working completely in the preconditioned
    // space, we save the cost of reconstructing the full residual
    // from the preconditioned smoother, and re-projecting for the
    // subsequent smoother, as well as reducing the cost of the
    // ancillary blas operations in the coarse-grid solve.
    //
    // Note, we cannot use preconditioned outer solve with direct
    // smoother
    //
    // Finally, we have to treat the top level carefully: for all
    // other levels the entry into and out of the grid will be a
    // full-field, which we can then work in Schur complement space or
    // not (e.g., freedom to choose coarse_grid_solution_type).  For
    // the top level, if the outer solver is for the preconditioned
    // system, then we must use preconditoning, e.g., option 3.) above.
 
    if (i == 0) { // top-level treatment
      if (coarse_solve_type[0] != solve_type)
        errorQuda("Mismatch between top-level MG solve type %d and outer solve type %d", coarse_solve_type[0],
                  solve_type);
 
      if (solve_type == QUDA_DIRECT_SOLVE) {
        mg_param.coarse_grid_solution_type[i] = QUDA_MAT_SOLUTION;
      } else if (solve_type == QUDA_DIRECT_PC_SOLVE) {
        mg_param.coarse_grid_solution_type[i] = QUDA_MATPC_SOLUTION;
      } else {
        errorQuda("Unexpected solve_type = %d\n", solve_type);
      }
 
    } else {
 
      if (coarse_solve_type[i] == QUDA_DIRECT_SOLVE) {
        mg_param.coarse_grid_solution_type[i] = QUDA_MAT_SOLUTION;
      } else if (coarse_solve_type[i] == QUDA_DIRECT_PC_SOLVE) {
        mg_param.coarse_grid_solution_type[i] = QUDA_MATPC_SOLUTION;
      } else {
        errorQuda("Unexpected solve_type = %d\n", coarse_solve_type[i]);
      }
    }
 
    mg_param.omega[i] = omega; // over/under relaxation factor
 
    // CPU setup disabled in release/1.1.x
    mg_param.location[i] = QUDA_CUDA_FIELD_LOCATION;
    mg_param.setup_location[i] = QUDA_CUDA_FIELD_LOCATION;
  }
 
  // whether to run GPU setup but putting temporaries into mapped (slow CPU) memory
  mg_param.setup_minimize_memory = QUDA_BOOLEAN_FALSE;
 
  // only coarsen the spin on the first restriction
  mg_param.spin_block_size[0] = 2;
 
  mg_param.setup_type = setup_type;
  mg_param.pre_orthonormalize = pre_orthonormalize ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  mg_param.post_orthonormalize = post_orthonormalize ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  mg_param.compute_null_vector = generate_nullspace ? QUDA_COMPUTE_NULL_VECTOR_YES : QUDA_COMPUTE_NULL_VECTOR_NO;
 
  mg_param.generate_all_levels = generate_all_levels ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  mg_param.run_verify = verify_results ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  mg_param.run_low_mode_check = low_mode_check ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  mg_param.run_oblique_proj_check = oblique_proj_check ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  mg_param.use_mma = mg_use_mma ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  // Whether or not to use thin restarts in the evolve tests
  mg_param.thin_update_only = mg_evolve_thin_updates ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  // set file i/o parameters
  for (int i = 0; i < mg_param.n_level; i++) {
    strcpy(mg_param.vec_infile[i], mg_vec_infile[i]);
    strcpy(mg_param.vec_outfile[i], mg_vec_outfile[i]);
    if (strcmp(mg_param.vec_infile[i], "") != 0) mg_param.vec_load[i] = QUDA_BOOLEAN_TRUE;
    if (strcmp(mg_param.vec_outfile[i], "") != 0) mg_param.vec_store[i] = QUDA_BOOLEAN_TRUE;
  }
 
  mg_param.coarse_guess = mg_eig_coarse_guess ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  mg_param.struct_size = sizeof(mg_param);
 
  // these need to tbe set for now but are actually ignored by the MG setup
  // needed to make it pass the initialization test
  inv_param.inv_type = QUDA_GCR_INVERTER;
  inv_param.tol = 1e-10;
  inv_param.maxiter = 1000;
  inv_param.reliable_delta = reliable_delta;
  inv_param.gcrNkrylov = 10;
 
  inv_param.verbosity = verbosity;
  inv_param.verbosity_precondition = verbosity;
 
  // Use kappa * csw or supplied clover_coeff
  inv_param.clover_csw = clover_csw;
  if (clover_coeff == 0.0) {
    inv_param.clover_coeff = clover_csw * inv_param.kappa;
  } else {
    inv_param.clover_coeff = clover_coeff;
  }  
 
}
 
void setMultigridInvertParam(QudaInvertParam &inv_param)
{
  inv_param.Ls = 1;
 
  inv_param.sp_pad = 0;
  inv_param.cl_pad = 0;
 
  inv_param.cpu_prec = cpu_prec;
  inv_param.cuda_prec = cuda_prec;
  inv_param.cuda_prec_sloppy = cuda_prec_sloppy;
 
  inv_param.cuda_prec_precondition = cuda_prec_precondition;
  inv_param.cuda_prec_eigensolver = cuda_prec_eigensolver;
  inv_param.preserve_source = QUDA_PRESERVE_SOURCE_NO;
  inv_param.gamma_basis = QUDA_DEGRAND_ROSSI_GAMMA_BASIS;
  inv_param.dirac_order = QUDA_DIRAC_ORDER;
 
  if (kappa == -1.0) {
    inv_param.mass = mass;
    inv_param.kappa = 1.0 / (2.0 * (1 + 3 / anisotropy + mass));
  } else {
    inv_param.kappa = kappa;
    inv_param.mass = 0.5 / kappa - (1 + 3 / anisotropy);
  }
  
  if (dslash_type == QUDA_CLOVER_WILSON_DSLASH || dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
    inv_param.clover_cpu_prec = cpu_prec;
    inv_param.clover_cuda_prec = cuda_prec;
    inv_param.clover_cuda_prec_sloppy = cuda_prec_sloppy;
    inv_param.clover_cuda_prec_precondition = cuda_prec_precondition;
    inv_param.clover_cuda_prec_eigensolver = cuda_prec_eigensolver;
    inv_param.clover_cuda_prec_refinement_sloppy = cuda_prec_sloppy;
    inv_param.clover_order = QUDA_PACKED_CLOVER_ORDER;
 
    // Use kappa * csw or supplied clover_coeff
    inv_param.clover_csw = clover_csw;
    if (clover_coeff == 0.0) {
      inv_param.clover_coeff = clover_csw * inv_param.kappa;
    } else {
      inv_param.clover_coeff = clover_coeff;
    }
  }
 
  inv_param.input_location = QUDA_CPU_FIELD_LOCATION;
  inv_param.output_location = QUDA_CPU_FIELD_LOCATION;
 
  inv_param.dslash_type = dslash_type;
 
  if (dslash_type == QUDA_TWISTED_MASS_DSLASH || dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
    inv_param.mu = mu;
    inv_param.epsilon = epsilon;
    inv_param.twist_flavor = twist_flavor;
    inv_param.Ls = (inv_param.twist_flavor == QUDA_TWIST_NONDEG_DOUBLET) ? 2 : 1;
 
    if (twist_flavor == QUDA_TWIST_NONDEG_DOUBLET) {
      printfQuda("Twisted-mass doublet non supported (yet)\n");
      exit(0);
    }
  }
 
  inv_param.clover_coeff = clover_coeff;
 
  inv_param.dagger = QUDA_DAG_NO;
  inv_param.mass_normalization = QUDA_KAPPA_NORMALIZATION;
 
  // do we want full solution or single-parity solution
  inv_param.solution_type = QUDA_MAT_SOLUTION;
 
  // do we want to use an even-odd preconditioned solve or not
  inv_param.solve_type = solve_type;
  inv_param.matpc_type = matpc_type;
 
  inv_param.inv_type = QUDA_GCR_INVERTER;
 
  inv_param.verbosity = verbosity;
  inv_param.verbosity_precondition = mg_verbosity[0];
 
  inv_param.inv_type_precondition = QUDA_MG_INVERTER;
  inv_param.pipeline = pipeline;
  inv_param.gcrNkrylov = gcrNkrylov;
  inv_param.tol = tol;
 
  // require both L2 relative and heavy quark residual to determine convergence
  inv_param.residual_type = static_cast<QudaResidualType>(QUDA_L2_RELATIVE_RESIDUAL);
  inv_param.tol_hq = tol_hq; // specify a tolerance for the residual for heavy quark residual
 
  // Offsets used only by multi-shift solver
  // should be set in application
  inv_param.num_offset = multishift;
  for (int i = 0; i < inv_param.num_offset; i++) inv_param.offset[i] = 0.06 + i * i * 0.1;
  // these can be set individually
  for (int i = 0; i < inv_param.num_offset; i++) {
    inv_param.tol_offset[i] = inv_param.tol;
    inv_param.tol_hq_offset[i] = inv_param.tol_hq;
  }
  inv_param.maxiter = niter;
  inv_param.reliable_delta = reliable_delta;
 
  // domain decomposition preconditioner is disabled when using MG
  inv_param.schwarz_type = QUDA_INVALID_SCHWARZ;
  inv_param.precondition_cycle = 1;
  inv_param.tol_precondition = 1e-1;
  inv_param.maxiter_precondition = 1;
  inv_param.omega = 1.0;
 
  // Whether or not to use native BLAS LAPACK
  inv_param.native_blas_lapack = (native_blas_lapack ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE);
 
  inv_param.struct_size = sizeof(inv_param);
}
 
// Parameters defining the eigensolver
void setMultigridEigParam(QudaEigParam &mg_eig_param, int level)
{
  mg_eig_param.eig_type = mg_eig_type[level];
  mg_eig_param.spectrum = mg_eig_spectrum[level];
  if ((mg_eig_type[level] == QUDA_EIG_TR_LANCZOS || mg_eig_type[level] == QUDA_EIG_BLK_TR_LANCZOS)
      && !(mg_eig_spectrum[level] == QUDA_SPECTRUM_LR_EIG || mg_eig_spectrum[level] == QUDA_SPECTRUM_SR_EIG)) {
    errorQuda("Only real spectrum type (LR or SR) can be passed to the a Lanczos type solver");
  }
 
  mg_eig_param.block_size
    = (mg_eig_param.eig_type == QUDA_EIG_TR_LANCZOS || mg_eig_param.eig_type == QUDA_EIG_IR_ARNOLDI) ?
    1 :
    mg_eig_block_size[level];
  mg_eig_param.n_ev = mg_eig_n_ev[level];
  mg_eig_param.n_kr = mg_eig_n_kr[level];
  mg_eig_param.n_conv = nvec[level];
 
  // Inverters will deflate only this number of vectors.
  if (mg_eig_n_ev_deflate[level] > 0 && mg_eig_n_ev_deflate[level] < mg_eig_param.n_conv)
    mg_eig_param.n_ev_deflate = mg_eig_n_ev_deflate[level];
  else if (mg_eig_n_ev_deflate[level] > mg_eig_param.n_conv)
    errorQuda("Can not deflate more than nvec[%d] eigenvectors.", mg_eig_param.n_conv);
  else {
    mg_eig_param.n_ev_deflate = mg_eig_param.n_conv;
    mg_eig_n_ev_deflate[level] = mg_eig_param.n_conv;
  }
 
  mg_eig_param.batched_rotate = mg_eig_batched_rotate[level];
  mg_eig_param.require_convergence = mg_eig_require_convergence[level] ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  mg_eig_param.tol = mg_eig_tol[level];
  mg_eig_param.qr_tol = mg_eig_qr_tol[level];
  mg_eig_param.check_interval = mg_eig_check_interval[level];
  mg_eig_param.max_restarts = mg_eig_max_restarts[level];
 
  mg_eig_param.compute_svd = QUDA_BOOLEAN_FALSE;
  mg_eig_param.compute_gamma5 = QUDA_BOOLEAN_FALSE;
  mg_eig_param.use_norm_op = mg_eig_use_normop[level] ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  mg_eig_param.use_dagger = mg_eig_use_dagger[level] ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  mg_eig_param.use_eigen_qr = mg_eig_use_eigen_qr[level] ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  mg_eig_param.use_poly_acc = mg_eig_use_poly_acc[level] ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  mg_eig_param.poly_deg = mg_eig_poly_deg[level];
  mg_eig_param.a_min = mg_eig_amin[level];
  mg_eig_param.a_max = mg_eig_amax[level];
 
  // set file i/o parameters
  // Give empty strings, Multigrid will handle IO.
  strcpy(mg_eig_param.vec_infile, "");
  strcpy(mg_eig_param.vec_outfile, "");
  mg_eig_param.save_prec = mg_eig_save_prec[level];
  mg_eig_param.io_parity_inflate = QUDA_BOOLEAN_FALSE;
 
  mg_eig_param.struct_size = sizeof(mg_eig_param);
}
 
void setContractInvertParam(QudaInvertParam &inv_param)
{
  inv_param.Ls = 1;
  inv_param.sp_pad = 0;
  inv_param.cl_pad = 0;
 
  inv_param.cpu_prec = cpu_prec;
  inv_param.cuda_prec = cuda_prec;
  inv_param.cuda_prec_sloppy = cuda_prec_sloppy;
  inv_param.cuda_prec_precondition = cuda_prec_precondition;
  inv_param.cuda_prec_eigensolver = cuda_prec_eigensolver;
 
  inv_param.preserve_source = QUDA_PRESERVE_SOURCE_NO;
  inv_param.dirac_order = QUDA_DIRAC_ORDER;
  // Quda performs contractions in Degrand-Rossi gamma basis,
  // but the user may suppy vectors in any supported order.
  inv_param.gamma_basis = QUDA_DEGRAND_ROSSI_GAMMA_BASIS;
 
  inv_param.input_location = QUDA_CPU_FIELD_LOCATION;
  inv_param.output_location = QUDA_CPU_FIELD_LOCATION;
 
  inv_param.struct_size = sizeof(inv_param);
}
 
void setStaggeredMGInvertParam(QudaInvertParam &inv_param)
{
  // Solver params
  inv_param.verbosity = QUDA_VERBOSE;
  inv_param.mass = mass;
 
  // outer solver parameters
  inv_param.inv_type = QUDA_GCR_INVERTER;
  inv_param.tol = tol;
  inv_param.maxiter = niter;
  inv_param.reliable_delta = 1e-4;
  inv_param.pipeline = pipeline;
 
  inv_param.Ls = 1;
 
  if (tol_hq == 0 && tol == 0) {
    errorQuda("qudaInvert: requesting zero residual\n");
    exit(1);
  }
 
  // require both L2 relative and heavy quark residual to determine convergence
  inv_param.residual_type = static_cast<QudaResidualType>(QUDA_L2_RELATIVE_RESIDUAL);
  inv_param.tol_hq = tol_hq; // specify a tolerance for the residual for heavy quark residual
 
  /* ESW HACK: comment this out to do a non-MG solve. */
  inv_param.inv_type_precondition = QUDA_MG_INVERTER;
  inv_param.verbosity_precondition = mg_verbosity[0];
  inv_param.cuda_prec_precondition = cuda_prec_precondition;
  inv_param.cuda_prec_eigensolver = cuda_prec_eigensolver;
 
  // Specify Krylov sub-size for GCR, BICGSTAB(L)
  inv_param.gcrNkrylov = gcrNkrylov;
 
  // do we want full solution or single-parity solution
  inv_param.solution_type = QUDA_MAT_SOLUTION;
 
  // do we want to use an even-odd preconditioned solve or not
  inv_param.solve_type = solve_type;
  inv_param.matpc_type = matpc_type;
  inv_param.dagger = QUDA_DAG_NO;
  inv_param.mass_normalization = QUDA_MASS_NORMALIZATION;
 
  inv_param.cpu_prec = cpu_prec;
  inv_param.cuda_prec = cuda_prec;
  inv_param.cuda_prec_sloppy = cuda_prec_sloppy;
  inv_param.preserve_source = QUDA_PRESERVE_SOURCE_YES;
  inv_param.gamma_basis = QUDA_DEGRAND_ROSSI_GAMMA_BASIS;
  inv_param.dirac_order = QUDA_DIRAC_ORDER;
 
  inv_param.dslash_type = dslash_type;
 
  inv_param.input_location = QUDA_CPU_FIELD_LOCATION;
  inv_param.output_location = QUDA_CPU_FIELD_LOCATION;
 
  inv_param.sp_pad = 0;
  inv_param.cl_pad = 0;
 
  // these can be set individually
  for (int i = 0; i < inv_param.num_offset; i++) {
    inv_param.tol_offset[i] = inv_param.tol;
    inv_param.tol_hq_offset[i] = inv_param.tol_hq;
  }
 
  // domain decomposition preconditioner is disabled when using MG
  inv_param.schwarz_type = QUDA_INVALID_SCHWARZ;
  inv_param.precondition_cycle = 1;
  inv_param.tol_precondition = 1e-1;
  inv_param.maxiter_precondition = 1;
  inv_param.omega = 1.0;
 
  // Whether or not to use native BLAS LAPACK
  inv_param.native_blas_lapack = (native_blas_lapack ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE);
 
  inv_param.struct_size = sizeof(inv_param);
}
 
void setStaggeredInvertParam(QudaInvertParam &inv_param)
{
  // Solver params
  inv_param.verbosity = verbosity;
  inv_param.mass = mass;
  inv_param.kappa = kappa = 1.0 / (8.0 + mass); // for Laplace operator
  inv_param.laplace3D = laplace3D;              // for Laplace operator
 
  // outer solver parameters
  inv_param.inv_type = inv_type;
  inv_param.tol = tol;
  inv_param.tol_restart = tol_restart;
  inv_param.maxiter = niter;
  inv_param.reliable_delta = reliable_delta;
  inv_param.use_alternative_reliable = alternative_reliable;
  inv_param.use_sloppy_partial_accumulator = false;
  inv_param.solution_accumulator_pipeline = solution_accumulator_pipeline;
  inv_param.pipeline = pipeline;
 
  inv_param.Ls = 1; // Nsrc
 
  if (tol_hq == 0 && tol == 0) {
    errorQuda("qudaInvert: requesting zero residual\n");
    exit(1);
  }
  // require both L2 relative and heavy quark residual to determine convergence
  inv_param.residual_type = static_cast<QudaResidualType_s>(0);
  inv_param.residual_type = (tol != 0) ?
    static_cast<QudaResidualType_s>(inv_param.residual_type | QUDA_L2_RELATIVE_RESIDUAL) :
    inv_param.residual_type;
  inv_param.residual_type = (tol_hq != 0) ?
    static_cast<QudaResidualType_s>(inv_param.residual_type | QUDA_HEAVY_QUARK_RESIDUAL) :
    inv_param.residual_type;
  inv_param.heavy_quark_check = (inv_param.residual_type & QUDA_HEAVY_QUARK_RESIDUAL ? 5 : 0);
 
  inv_param.tol_hq = tol_hq; // specify a tolerance for the residual for heavy quark residual
 
  inv_param.Nsteps = 2;
 
  // domain decomposition preconditioner parameters
  inv_param.inv_type_precondition = precon_type;
  inv_param.tol_precondition = tol_precondition;
  inv_param.maxiter_precondition = maxiter_precondition;
  inv_param.verbosity_precondition = QUDA_SILENT;
  inv_param.cuda_prec_precondition = prec_precondition;
  inv_param.cuda_prec_eigensolver = prec_eigensolver;
 
  // Specify Krylov sub-size for GCR, BICGSTAB(L), basis size for CA-CG, CA-GCR
  inv_param.gcrNkrylov = gcrNkrylov;
 
  // Specify basis for CA-CG, lambda min/max for Chebyshev basis
  //   lambda_max < lambda_max . use power iters to generate
  inv_param.ca_basis = ca_basis;
  inv_param.ca_lambda_min = ca_lambda_min;
  inv_param.ca_lambda_max = ca_lambda_max;
 
  inv_param.solution_type = solution_type;
  inv_param.solve_type = solve_type;
  inv_param.matpc_type = matpc_type;
  inv_param.dagger = QUDA_DAG_NO;
  inv_param.mass_normalization = QUDA_MASS_NORMALIZATION;
 
  inv_param.cpu_prec = cpu_prec;
  inv_param.cuda_prec = prec;
  inv_param.cuda_prec_sloppy = prec_sloppy;
  inv_param.cuda_prec_refinement_sloppy = prec_refinement_sloppy;
  inv_param.preserve_source = QUDA_PRESERVE_SOURCE_YES;
  inv_param.gamma_basis = QUDA_DEGRAND_ROSSI_GAMMA_BASIS; // this is meaningless, but must be thus set
  inv_param.dirac_order = QUDA_DIRAC_ORDER;
 
  inv_param.dslash_type = dslash_type;
 
  inv_param.input_location = QUDA_CPU_FIELD_LOCATION;
  inv_param.output_location = QUDA_CPU_FIELD_LOCATION;
 
  inv_param.sp_pad = 0;
 
  // Whether or not to use native BLAS LAPACK
  inv_param.native_blas_lapack = (native_blas_lapack ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE);
 
  inv_param.struct_size = sizeof(inv_param);
}
 
void setStaggeredMultigridParam(QudaMultigridParam &mg_param)
{
  QudaInvertParam &inv_param = *mg_param.invert_param; // this will be used to setup SolverParam parent in MGParam class
 
  // Whether or not to use native BLAS LAPACK
  inv_param.native_blas_lapack = (native_blas_lapack ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE);
 
  inv_param.Ls = 1;
 
  inv_param.sp_pad = 0;
  inv_param.cl_pad = 0;
 
  inv_param.cpu_prec = cpu_prec;
  inv_param.cuda_prec = cuda_prec;
  inv_param.cuda_prec_sloppy = cuda_prec_sloppy;
  inv_param.cuda_prec_precondition = cuda_prec_precondition;
  inv_param.cuda_prec_eigensolver = cuda_prec_eigensolver;
  inv_param.preserve_source = QUDA_PRESERVE_SOURCE_NO;
  inv_param.gamma_basis = QUDA_DEGRAND_ROSSI_GAMMA_BASIS;
  inv_param.dirac_order = QUDA_DIRAC_ORDER;
 
  inv_param.input_location = QUDA_CPU_FIELD_LOCATION;
  inv_param.output_location = QUDA_CPU_FIELD_LOCATION;
 
  inv_param.dslash_type = dslash_type;
 
  inv_param.mass = mass;
  inv_param.kappa = 1.0 / (2.0 * (4.0 + inv_param.mass));
 
  inv_param.dagger = QUDA_DAG_NO;
  inv_param.mass_normalization = QUDA_MASS_NORMALIZATION;
 
  inv_param.matpc_type = matpc_type;
  inv_param.solution_type = QUDA_MAT_SOLUTION;
 
  inv_param.solve_type = QUDA_DIRECT_SOLVE;
 
  mg_param.use_mma = mg_use_mma ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  mg_param.invert_param = &inv_param;
  mg_param.n_level = mg_levels;
  for (int i = 0; i < mg_param.n_level; i++) {
    for (int j = 0; j < 4; j++) {
      // if not defined use 4
      mg_param.geo_block_size[i][j] = geo_block_size[i][j] ? geo_block_size[i][j] : 4;
    }
    mg_param.use_eig_solver[i] = mg_eig[i] ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
    mg_param.verbosity[i] = mg_verbosity[i];
    mg_param.setup_inv_type[i] = setup_inv[i];
    mg_param.num_setup_iter[i] = num_setup_iter[i];
    mg_param.setup_tol[i] = setup_tol[i];
    mg_param.setup_maxiter[i] = setup_maxiter[i];
 
    // Basis to use for CA-CGN(E/R) setup
    mg_param.setup_ca_basis[i] = setup_ca_basis[i];
 
    // Basis size for CACG setup
    mg_param.setup_ca_basis_size[i] = setup_ca_basis_size[i];
 
    // Minimum and maximum eigenvalue for Chebyshev CA basis setup
    mg_param.setup_ca_lambda_min[i] = setup_ca_lambda_min[i];
    mg_param.setup_ca_lambda_max[i] = setup_ca_lambda_max[i];
 
    mg_param.spin_block_size[i] = 1;
    mg_param.n_vec[i] = nvec[i] == 0 ? 64 : nvec[i];          // default to 64 vectors if not set
    mg_param.n_block_ortho[i] = n_block_ortho[i];             // number of times to Gram-Schmidt
    mg_param.precision_null[i] = prec_null;                   // precision to store the null-space basis
    mg_param.smoother_halo_precision[i] = smoother_halo_prec; // precision of the halo exchange in the smoother
    mg_param.nu_pre[i] = nu_pre[i];
    mg_param.nu_post[i] = nu_post[i];
    mg_param.mu_factor[i] = mu_factor[i];
 
    mg_param.transfer_type[i] = (i == 0) ? staggered_transfer_type : QUDA_TRANSFER_AGGREGATE;
 
    mg_param.cycle_type[i] = QUDA_MG_CYCLE_RECURSIVE;
 
    // set the coarse solver wrappers including bottom solver
    mg_param.coarse_solver[i] = coarse_solver[i];
    mg_param.coarse_solver_tol[i] = coarse_solver_tol[i];
    mg_param.coarse_solver_maxiter[i] = coarse_solver_maxiter[i];
 
    // Basis to use for CA-CGN(E/R) coarse solver
    mg_param.coarse_solver_ca_basis[i] = coarse_solver_ca_basis[i];
 
    // Basis size for CACG coarse solver/
    mg_param.coarse_solver_ca_basis_size[i] = coarse_solver_ca_basis_size[i];
 
    // Minimum and maximum eigenvalue for Chebyshev CA basis
    mg_param.coarse_solver_ca_lambda_min[i] = coarse_solver_ca_lambda_min[i];
    mg_param.coarse_solver_ca_lambda_max[i] = coarse_solver_ca_lambda_max[i];
 
    mg_param.smoother[i] = smoother_type[i];
 
    // set the smoother / bottom solver tolerance (for MR smoothing this will be ignored)
    mg_param.smoother_tol[i] = smoother_tol[i];
 
    // set to QUDA_DIRECT_SOLVE for no even/odd preconditioning on the smoother
    // set to QUDA_DIRECT_PC_SOLVE for to enable even/odd preconditioning on the smoother
    mg_param.smoother_solve_type[i] = smoother_solve_type[i];
 
    // set to QUDA_ADDITIVE_SCHWARZ for Additive Schwarz precondioned smoother (presently only impelemented for MR)
    mg_param.smoother_schwarz_type[i] = mg_schwarz_type[i];
 
    // if using Schwarz preconditioning then use local reductions only
    mg_param.global_reduction[i] = (mg_schwarz_type[i] == QUDA_INVALID_SCHWARZ) ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
    // set number of Schwarz cycles to apply
    mg_param.smoother_schwarz_cycle[i] = mg_schwarz_cycle[i];
 
    // Set set coarse_grid_solution_type: this defines which linear
    // system we are solving on a given level
    // * QUDA_MAT_SOLUTION - we are solving the full system and inject
    //   a full field into coarse grid
    // * QUDA_MATPC_SOLUTION - we are solving the e/o-preconditioned
    //   system, and only inject single parity field into coarse grid
    //
    // Multiple possible scenarios here
    //
    // 1. **Direct outer solver and direct smoother**: here we use
    // full-field residual coarsening, and everything involves the
    // full system so coarse_grid_solution_type = QUDA_MAT_SOLUTION
    //
    // 2. **Direct outer solver and preconditioned smoother**: here,
    // only the smoothing uses e/o preconditioning, so
    // coarse_grid_solution_type = QUDA_MAT_SOLUTION_TYPE.
    // We reconstruct the full residual prior to coarsening after the
    // pre-smoother, and then need to project the solution for post
    // smoothing.
    //
    // 3. **Preconditioned outer solver and preconditioned smoother**:
    // here we use single-parity residual coarsening throughout, so
    // coarse_grid_solution_type = QUDA_MATPC_SOLUTION.  This is a bit
    // questionable from a theoretical point of view, since we don't
    // coarsen the preconditioned operator directly, rather we coarsen
    // the full operator and preconditioned that, but it just works.
    // This is the optimal combination in general for Wilson-type
    // operators: although there is an occasional increase in
    // iteration or two), by working completely in the preconditioned
    // space, we save the cost of reconstructing the full residual
    // from the preconditioned smoother, and re-projecting for the
    // subsequent smoother, as well as reducing the cost of the
    // ancillary blas operations in the coarse-grid solve.
    //
    // Note, we cannot use preconditioned outer solve with direct
    // smoother
    //
    // Finally, we have to treat the top level carefully: for all
    // other levels the entry into and out of the grid will be a
    // full-field, which we can then work in Schur complement space or
    // not (e.g., freedom to choose coarse_grid_solution_type).  For
    // the top level, if the outer solver is for the preconditioned
    // system, then we must use preconditoning, e.g., option 3.) above.
 
    if (i == 0) { // top-level treatment
      if (coarse_solve_type[0] != solve_type)
        errorQuda("Mismatch between top-level MG solve type %s and outer solve type %s",
                  get_solve_str(coarse_solve_type[0]), get_solve_str(solve_type));
 
      if (solve_type == QUDA_DIRECT_SOLVE) {
        mg_param.coarse_grid_solution_type[i] = QUDA_MAT_SOLUTION;
      } else if (solve_type == QUDA_DIRECT_PC_SOLVE) {
        mg_param.coarse_grid_solution_type[i] = QUDA_MATPC_SOLUTION;
      } else {
        errorQuda("Unexpected solve_type = %s\n", get_solve_str(solve_type));
      }
 
    } else {
 
      if (coarse_solve_type[i] == QUDA_DIRECT_SOLVE) {
        mg_param.coarse_grid_solution_type[i] = QUDA_MAT_SOLUTION;
      } else if (coarse_solve_type[i] == QUDA_DIRECT_PC_SOLVE) {
        mg_param.coarse_grid_solution_type[i] = QUDA_MATPC_SOLUTION;
      } else {
        errorQuda("Unexpected solve_type = %s\n", get_solve_str(coarse_solve_type[i]));
      }
    }
 
    mg_param.omega[i] = omega; // over/under relaxation factor
 
    // CPU setup disabled in release/1.1.x
    mg_param.location[i] = QUDA_CUDA_FIELD_LOCATION;
    mg_param.setup_location[i] = QUDA_CUDA_FIELD_LOCATION;
    nu_pre[i] = 2;
    nu_post[i] = 2;
  }
 
  // whether to run GPU setup but putting temporaries into mapped (slow CPU) memory
  mg_param.setup_minimize_memory = QUDA_BOOLEAN_FALSE;
 
  // coarsening the spin on the first restriction is undefined for staggered fields.
  mg_param.spin_block_size[0] = 0;
 
  if (staggered_transfer_type == QUDA_TRANSFER_OPTIMIZED_KD) {
    mg_param.spin_block_size[1] = 0; // we're coarsening the optimized KD op
  }
 
  mg_param.setup_type = setup_type;
  mg_param.pre_orthonormalize = pre_orthonormalize ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  mg_param.post_orthonormalize = post_orthonormalize ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  mg_param.compute_null_vector = generate_nullspace ? QUDA_COMPUTE_NULL_VECTOR_YES : QUDA_COMPUTE_NULL_VECTOR_NO;
 
  mg_param.generate_all_levels = generate_all_levels ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  mg_param.run_verify = verify_results ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  mg_param.run_low_mode_check = low_mode_check ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
  mg_param.run_oblique_proj_check = oblique_proj_check ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  // set file i/o parameters
  for (int i = 0; i < mg_param.n_level; i++) {
    strcpy(mg_param.vec_infile[i], mg_vec_infile[i]);
    strcpy(mg_param.vec_outfile[i], mg_vec_outfile[i]);
    if (strcmp(mg_param.vec_infile[i], "") != 0) mg_param.vec_load[i] = QUDA_BOOLEAN_TRUE;
    if (strcmp(mg_param.vec_outfile[i], "") != 0) mg_param.vec_store[i] = QUDA_BOOLEAN_TRUE;
  }
 
  mg_param.coarse_guess = mg_eig_coarse_guess ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
 
  // these need to tbe set for now but are actually ignored by the MG setup
  // needed to make it pass the initialization test
  inv_param.inv_type = QUDA_GCR_INVERTER;
  inv_param.tol = 1e-10;
  inv_param.maxiter = 1000;
  inv_param.reliable_delta = reliable_delta;
  inv_param.gcrNkrylov = 10;
 
  inv_param.verbosity = verbosity;
  inv_param.verbosity_precondition = verbosity;
 
  inv_param.struct_size = sizeof(inv_param);
}
 
void setDeflatedInvertParam(QudaInvertParam &inv_param)
{
  inv_param.Ls = 1;
 
  inv_param.sp_pad = 0;
  inv_param.cl_pad = 0;
 
  inv_param.cpu_prec = cpu_prec;
  inv_param.cuda_prec = cuda_prec;
  inv_param.cuda_prec_sloppy = cuda_prec_sloppy;
  inv_param.cuda_prec_refinement_sloppy = cuda_prec_refinement_sloppy;
 
  inv_param.cuda_prec_precondition = cuda_prec_precondition;
  inv_param.cuda_prec_eigensolver = cuda_prec_eigensolver;
  inv_param.preserve_source = QUDA_PRESERVE_SOURCE_NO;
  inv_param.gamma_basis = QUDA_DEGRAND_ROSSI_GAMMA_BASIS;
  inv_param.dirac_order = QUDA_DIRAC_ORDER;
 
  if (dslash_type == QUDA_CLOVER_WILSON_DSLASH || dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
    inv_param.clover_cpu_prec = cpu_prec;
    inv_param.clover_cuda_prec = cuda_prec;
    inv_param.clover_cuda_prec_sloppy = cuda_prec_sloppy;
    inv_param.clover_cuda_prec_precondition = cuda_prec_precondition;
    inv_param.clover_order = QUDA_PACKED_CLOVER_ORDER;
  }
 
  inv_param.input_location = QUDA_CPU_FIELD_LOCATION;
  inv_param.output_location = QUDA_CPU_FIELD_LOCATION;
 
  inv_param.dslash_type = dslash_type;
 
  if (kappa == -1.0) {
    inv_param.mass = mass;
    inv_param.kappa = 1.0 / (2.0 * (1 + 3 / anisotropy + mass));
  } else {
    inv_param.kappa = kappa;
    inv_param.mass = 0.5 / kappa - (1 + 3 / anisotropy);
  }
 
  if (dslash_type == QUDA_TWISTED_MASS_DSLASH || dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
    inv_param.mu = mu;
    inv_param.twist_flavor = twist_flavor;
    inv_param.Ls = (inv_param.twist_flavor == QUDA_TWIST_NONDEG_DOUBLET) ? 2 : 1;
 
    if (twist_flavor == QUDA_TWIST_NONDEG_DOUBLET) {
      printfQuda("Twisted-mass doublet non supported (yet)\n");
      exit(0);
    }
  }
 
  // Use kappa * csw or supplied clover_coeff
  inv_param.clover_csw = clover_csw;
  if (clover_coeff == 0.0) {
    inv_param.clover_coeff = clover_csw * inv_param.kappa;
  } else {
    inv_param.clover_coeff = clover_coeff;
  }
 
  inv_param.dagger = QUDA_DAG_NO;
  inv_param.mass_normalization = normalization;
 
  // do we want full solution or single-parity solution
  inv_param.solution_type = QUDA_MAT_SOLUTION;
  // inv_param.solution_type = QUDA_MATPC_SOLUTION;
 
  // do we want to use an even-odd preconditioned solve or not
  inv_param.solve_type = solve_type;
  inv_param.matpc_type = matpc_type;
 
  if (inv_type != QUDA_EIGCG_INVERTER && inv_type != QUDA_INC_EIGCG_INVERTER && inv_type != QUDA_GMRESDR_INVERTER)
    errorQuda("Requested solver %s is not a deflated solver type", get_solver_str(inv_type));
 
  inv_param.inv_type = inv_type;
  inv_param.tol = tol;
  inv_param.tol_hq = tol_hq; // specify a tolerance for the residual for heavy quark residual
 
  inv_param.rhs_idx = 0;
 
  inv_param.n_ev = n_ev;
  inv_param.max_search_dim = max_search_dim;
  inv_param.deflation_grid = deflation_grid;
  inv_param.tol_restart = tol_restart;
  inv_param.eigcg_max_restarts = eigcg_max_restarts;
  inv_param.max_restart_num = max_restart_num;
  inv_param.inc_tol = inc_tol;
  inv_param.eigenval_tol = eigenval_tol;
 
  if (inv_param.inv_type == QUDA_EIGCG_INVERTER || inv_param.inv_type == QUDA_INC_EIGCG_INVERTER) {
    inv_param.solve_type = QUDA_NORMOP_PC_SOLVE;
  } else if (inv_param.inv_type == QUDA_GMRESDR_INVERTER) {
    inv_param.solve_type = QUDA_DIRECT_PC_SOLVE;
    inv_param.tol_restart = 0.0; // restart is not requested...
  }
 
  inv_param.verbosity = verbosity;
  inv_param.verbosity_precondition = verbosity;
 
  inv_param.inv_type_precondition = precon_type;
  inv_param.gcrNkrylov = 6;
 
  // require both L2 relative and heavy quark residual to determine convergence
  inv_param.residual_type = static_cast<QudaResidualType>(QUDA_L2_RELATIVE_RESIDUAL);
  // Offsets used only by multi-shift solver
  // should be set in application
  inv_param.num_offset = multishift;
  for (int i = 0; i < inv_param.num_offset; i++) inv_param.offset[i] = 0.06 + i * i * 0.1;
  // these can be set individually
  for (int i = 0; i < inv_param.num_offset; i++) {
    inv_param.tol_offset[i] = inv_param.tol;
    inv_param.tol_hq_offset[i] = inv_param.tol_hq;
  }
  inv_param.maxiter = niter;
  inv_param.reliable_delta = 1e-1;
 
  // domain decomposition preconditioner parameters
  inv_param.schwarz_type = precon_schwarz_type;
  inv_param.precondition_cycle = precon_schwarz_cycle;
  inv_param.tol_precondition = tol_precondition;
  inv_param.maxiter_precondition = maxiter_precondition;
  inv_param.omega = 1.0;
 
  inv_param.extlib_type = solver_ext_lib;
 
  // Whether or not to use native BLAS LAPACK
  inv_param.native_blas_lapack = (native_blas_lapack ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE);
 
  inv_param.struct_size = sizeof(inv_param);
}
 
void setDeflationParam(QudaEigParam &df_param)
{
  df_param.import_vectors = QUDA_BOOLEAN_FALSE;
  df_param.run_verify = QUDA_BOOLEAN_FALSE;
 
  df_param.nk = df_param.invert_param->n_ev;
  df_param.np = df_param.invert_param->n_ev * df_param.invert_param->deflation_grid;
  df_param.extlib_type = deflation_ext_lib;
 
  df_param.cuda_prec_ritz = prec_ritz;
  df_param.location = location_ritz;
  df_param.mem_type_ritz = mem_type_ritz;
 
  // set file i/o parameters
  strcpy(df_param.vec_infile, eig_vec_infile);
  strcpy(df_param.vec_outfile, eig_vec_outfile);
  df_param.io_parity_inflate = eig_io_parity_inflate ? QUDA_BOOLEAN_TRUE : QUDA_BOOLEAN_FALSE;
}
 
void setQudaStaggeredInvTestParams()
{
  if (dslash_type == QUDA_LAPLACE_DSLASH) {
    if (test_type != 0) { errorQuda("Test type %d is not supported for the Laplace operator.\n", test_type); }
 
    solve_type = QUDA_DIRECT_SOLVE;
    solution_type = QUDA_MAT_SOLUTION;
    matpc_type = QUDA_MATPC_EVEN_EVEN; // doesn't matter
 
  } else {
 
    if (test_type == 0 && (inv_type == QUDA_CG_INVERTER || inv_type == QUDA_PCG_INVERTER)
        && solve_type != QUDA_NORMOP_SOLVE && solve_type != QUDA_DIRECT_PC_SOLVE) {
      warningQuda("The full spinor staggered operator (test 0) can't be inverted with (P)CG. Switching to BiCGstab.\n");
      inv_type = QUDA_BICGSTAB_INVERTER;
    }
 
    if (solve_type == QUDA_INVALID_SOLVE) {
      if (test_type == 0) {
        solve_type = QUDA_DIRECT_SOLVE;
      } else {
        solve_type = QUDA_DIRECT_PC_SOLVE;
      }
    }
 
    if (test_type == 1 || test_type == 3 || test_type == 5) {
      matpc_type = QUDA_MATPC_EVEN_EVEN;
    } else if (test_type == 2 || test_type == 4 || test_type == 6) {
      matpc_type = QUDA_MATPC_ODD_ODD;
    } else if (test_type == 0) {
      matpc_type = QUDA_MATPC_EVEN_EVEN; // it doesn't matter
    }
 
    if (test_type == 0 || test_type == 1 || test_type == 2) {
      solution_type = QUDA_MAT_SOLUTION;
    } else {
      solution_type = QUDA_MATPC_SOLUTION;
    }
  }
 
  if (prec_sloppy == QUDA_INVALID_PRECISION) { prec_sloppy = prec; }
 
  if (prec_refinement_sloppy == QUDA_INVALID_PRECISION) { prec_refinement_sloppy = prec_sloppy; }
  if (link_recon_sloppy == QUDA_RECONSTRUCT_INVALID) { link_recon_sloppy = link_recon; }
 
  if (inv_type != QUDA_CG_INVERTER && (test_type == 5 || test_type == 6)) {
    errorQuda("Preconditioning is currently not supported in multi-shift solver solvers");
  }
 
  // Set n_naiks to 2 if eps_naik != 0.0
  if (dslash_type == QUDA_ASQTAD_DSLASH) {
    if (eps_naik != 0.0) {
      if (compute_fatlong) {
        n_naiks = 2;
        printfQuda("Note: epsilon-naik != 0, testing epsilon correction links.\n");
      } else {
        eps_naik = 0.0;
        printfQuda("Not computing fat-long, ignoring epsilon correction.\n");
      }
    } else {
      printfQuda("Note: epsilon-naik = 0, testing original HISQ links.\n");
    }
  }
}
 
void setQudaStaggeredEigTestParams()
{
  if (dslash_type == QUDA_LAPLACE_DSLASH) {
    // LAPLACE operator path, only DIRECT solves feasible.
    if (test_type != 0) { errorQuda("Test type %d is not supported for the Laplace operator.\n", test_type); }
    solve_type = QUDA_DIRECT_SOLVE;
    solution_type = QUDA_MAT_SOLUTION;
  } else {
    // STAGGERED operator path
    if (solve_type == QUDA_INVALID_SOLVE) {
      if (test_type == 0) {
        solve_type = QUDA_DIRECT_SOLVE;
      } else {
        solve_type = QUDA_DIRECT_PC_SOLVE;
      }
    }
    // If test type is not 3, it is 4 or 0. If 0, the matpc type is irrelevant
    if (test_type == 3)
      matpc_type = QUDA_MATPC_EVEN_EVEN;
    else
      matpc_type = QUDA_MATPC_ODD_ODD;
 
    if (test_type == 0) {
      solution_type = QUDA_MAT_SOLUTION;
    } else {
      solution_type = QUDA_MATPC_SOLUTION;
    }
  }
 
  // Set n_naiks to 2 if eps_naik != 0.0
  if (dslash_type == QUDA_ASQTAD_DSLASH) {
    if (eps_naik != 0.0) {
      if (compute_fatlong) {
        n_naiks = 2;
        printfQuda("Note: epsilon-naik != 0, testing epsilon correction links.\n");
      } else {
        eps_naik = 0.0;
        printfQuda("Not computing fat-long, ignoring epsilon correction.\n");
      }
    } else {
      printfQuda("Note: epsilon-naik = 0, testing original HISQ links.\n");
    }
  }
}