write(xml_out, "lambda_Msq_hi", lambda_high_aux[0]);
write(xml_out, "lambda_hi", lambda_high_aux[0]);
void RitzCode4DHw(Handle< const LinearOperator<LatticeFermion> >& MM,		  Handle< const LinearOperator<LatticeFermion> >& H,		  const ChromaWilsonRitz_t& input,		  XMLWriter& xml_out){  // Try and get lowest eigenvalue of MM  const OrderedSubset& s = MM->subset();    multi1d<Real> lambda(input.ritz_params.Neig+input.ritz_params.Ndummy);  multi1d<Real> check_norm(input.ritz_params.Neig);  multi1d<LatticeFermion> psi(input.ritz_params.Neig			      +input.ritz_params.Ndummy);          for(int i =0; i < input.ritz_params.Neig + input.ritz_params.Ndummy; i++){    psi[i] = zero;    gaussian(psi[i],s);    lambda[i] = Real(1);  }      int n_CG_count;  Real delta_cycle = Real(1);  XMLBufferWriter eig_spec_xml;  int n_KS_count = 0;  int n_jacob_count = 0;  EigSpecRitzKS(*MM, 		lambda, 		psi, 		input.ritz_params.Neig,		input.ritz_params.Ndummy,                // No of dummies		input.ritz_params.Nrenorm, 		input.ritz_params.MinKSIter, 		input.ritz_params.MaxKSIter,             // Max iters / KS cycle		input.ritz_params.MaxKS,            // Max no of KS cycles		input.ritz_params.GammaFactor,       // Gamma factor		input.ritz_params.MaxCG,		input.ritz_params.RsdR,		input.ritz_params.RsdA,  		input.ritz_params.RsdZero,		input.ritz_params.ProjApsiP,		n_CG_count,		n_KS_count,		n_jacob_count,		eig_spec_xml);    // Dump output  xml_out << eig_spec_xml;  write(xml_out, "lambda_Msq", lambda);     // Check norms  for(int i=0; i < input.ritz_params.Neig; i++) {     LatticeFermion Me;    LatticeFermion lambda_e;    (*MM)(Me, psi[i], PLUS);    lambda_e[s] = lambda[i]*psi[i];            LatticeFermion r_norm;    r_norm[s] = Me - lambda_e;        check_norm[i] = norm2(r_norm,s);    check_norm[i] = sqrt(check_norm[i]);  }  write(xml_out, "check_norm", check_norm);    for(int i=0; i < input.ritz_params.Neig; i++) {    check_norm[i] /= fabs(lambda[i]);  }  write(xml_out, "check_norm_rel", check_norm);    // Fix to ev-s of gamma_5 wilson...  // Try to get one:  multi1d<bool> valid_eig(input.ritz_params.Neig);  int n_valid;  int n_jacob;    fixMMev2Mev(*H, 	      lambda, 	      psi, 	      input.ritz_params.Neig, 	      input.ritz_params.RsdR,	      input.ritz_params.RsdA, 	      input.ritz_params.RsdZero, 	      valid_eig, 	      n_valid, 	      n_jacob);    push(xml_out, "eigFix");  write(xml_out, "lambda_Hw", lambda);  write(xml_out, "n_valid", n_valid);  write(xml_out, "valid_eig", valid_eig);    for(int i=0; i < input.ritz_params.Neig; i++) {     LatticeFermion Me;    (*H)(Me, psi[i], PLUS);        bool zeroP = toBool( fabs(lambda[i]) <  input.ritz_params.RsdZero );    if( zeroP ) {      check_norm[i] = norm2(Me,s);      check_norm[i] = sqrt(check_norm[i]);    }    else {      LatticeFermion lambda_e;      LatticeFermion r_norm;            lambda_e[s] = lambda[i]*psi[i];      r_norm[s] = Me - lambda_e;                  check_norm[i] = norm2(r_norm,s);      check_norm[i] = sqrt(check_norm[i]);    }       QDPIO::cout << "lambda_lo[" << i << "] = " << lambda[i] << "  ";     QDPIO::cout << "check_norm["<<i<<"] = " << check_norm[i] << endl;  }  write(xml_out, "check_norm", check_norm);    for(int i=0; i < input.ritz_params.Neig; i++) {     check_norm[i] /= fabs(lambda[i]);    QDPIO::cout << "check_norm_rel["<< i <<"] = " << check_norm[i] << endl;  }  QDPIO::cout << flush ;  write(xml_out, "check_norm_rel", check_norm);  pop(xml_out);      // Now get the absolute value of the  highest e-value  // Work with H^{dag}H = M^{dag}M  Real hi_RsdR = 1.0e-4;  Real hi_RsdA = 1.0e-4;    multi1d<Real> lambda_high_aux(1);  multi1d<LatticeFermion> lambda_high_vec(1);  gaussian(lambda_high_vec[0],s);  lambda_high_vec[0][s] /= sqrt(norm2(lambda_high_vec[0],s));  int n_cg_high;  XMLBufferWriter high_xml;    Handle<const LinearOperator<LatticeFermion> > MinusMM = new lopscl<LatticeFermion, Real>(MM, Real(-1.0));  // Initial guess -- upper bound on spectrum  lambda_high_aux[0] = Real(8);      push(high_xml, "LambdaHighRitz");        // Minus MM ought to produce a negative e-value  // since MM is herm_pos_def  // ie minus MM is hermitian -ve definite  EigSpecRitzCG( *MinusMM,		 lambda_high_aux,		 lambda_high_vec,		 1,		 input.ritz_params.Nrenorm,		 input.ritz_params.MinKSIter,		 input.ritz_params.MaxCG,		 hi_RsdR,		 hi_RsdA,		 input.ritz_params.RsdZero,		 input.ritz_params.ProjApsiP,		 n_cg_high,		 high_xml);    lambda_high_aux[0] = sqrt(fabs(lambda_high_aux[0]));  QDPIO::cout << "|| lambda_hi || = " << lambda_high_aux[0]  << " hi_Rsd_r = " << hi_RsdR << endl;    pop(high_xml);  xml_out << high_xml;    push(xml_out, "Highest");  write(xml_out, "lambda_Msq_hi", lambda_high_aux[0]);  pop(xml_out);  QDPIO::cout << "Writing low eigenvalues/vectors" << endl;  writeEigen(input, lambda, psi, lambda_high_aux[0], QDPIO_SERIAL);}