lambda_high_aux[0] = sqrt(fabs(lambda_high_aux[0]));
lambda_high_aux[0] = fabs(lambda_high_aux[0]);
void RitzCode5D(Handle< const LinearOperator<multi1d<LatticeFermion> > >& MM,		const ChromaWilsonRitz_t& input,		XMLWriter& xml_out){  // Try and get lowest eigenvalue of MM  int N5 = MM->size();  const OrderedSubset& s = MM->subset();    multi1d<Real> lambda(input.ritz_params.Neig+input.ritz_params.Ndummy);  multi1d<Real> check_norm(input.ritz_params.Neig);  multi2d<LatticeFermion> psi(input.ritz_params.Neig			      +input.ritz_params.Ndummy, N5);          for(int i =0; i < input.ritz_params.Neig + input.ritz_params.Ndummy; i++){    for(int n=0; n < N5; n++) {       psi[i][n] = zero;      gaussian(psi[i][n],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);    /*  EigSpecRitzCG( *MM,		 lambda,		 psi,		 input.ritz_params.Neig,		 input.ritz_params.Nrenorm,		 input.ritz_params.MinKSIter,		 input.ritz_params.MaxCG,		 input.ritz_params.RsdR,		 input.ritz_params.RsdA,		 input.ritz_params.RsdZero,		 input.ritz_params.ProjApsiP,		 n_CG_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++) {     multi1d<LatticeFermion> Me(N5);    multi1d<LatticeFermion> lambda_e(N5);    (*MM)(Me, psi[i], PLUS);    for(int n =0; n < N5; n++) {       lambda_e[n][s] = lambda[i]*psi[i][n];    }        multi1d<LatticeFermion> r_norm(N5);    for(int n=0; n < N5; n++) {       r_norm[n][s] = Me[n] - lambda_e[n];    }        check_norm[i] = norm2(r_norm[0],s);    for(int n=1; n < N5; n++) {      check_norm[i] += norm2(r_norm[n],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);    // 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;  Handle<const LinearOperator<multi1d<LatticeFermion> > > MinusMM = new lopscl<multi1d<LatticeFermion>, Real>(MM, Real(-1.0));    multi1d<Real> lambda_high_aux(1);  multi2d<LatticeFermion> lambda_high_vec(1,N5);  for(int n=0; n < N5; n++) {    gaussian(lambda_high_vec[0][n],s);    lambda_high_vec[0][n][s] /= sqrt(norm2(lambda_high_vec[0],s));  }    int n_cg_high;  XMLBufferWriter high_xml;    // 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);}