C++ MultiFab::plus方法代码示例

//.........这里部分代码省略.........
             const Real rel_error = error / norm_to_test_against;
             Spacer(std::cout, level);
             if (using_bnorm)
             {
                 std::cout << "MultiGrid: Iteration   "
                           << nit
                           << " resid/bnorm = "
                           << rel_error << '\n';
             } else {
                 std::cout << "MultiGrid: Iteration   "
                           << nit
                           << " resid/resid0 = "
                           << rel_error << '\n';
             }
         }
     }
  }

  Real run_time = (ParallelDescriptor::second() - strt_time);

  if ( verbose > 0 )
  {
      if ( ParallelDescriptor::IOProcessor() )
      {
          const Real rel_error = error / norm_to_test_against;
          Spacer(std::cout, level);
          if (using_bnorm)
          {
              std::cout << "MultiGrid: Iteration   "
                        << nit-1
                        << " resid/bnorm = "
                        << rel_error << '\n';
          } else {
              std::cout << "MultiGrid: Iteration   "
                        << nit-1
                        << " resid/resid0 = "
                        << rel_error << '\n';
             }
      }

      if ( verbose > 1 )
      {
          Real tmp[2] = { run_time, cg_time };

          ParallelDescriptor::ReduceRealMax(tmp,2,ParallelDescriptor::IOProcessorNumber());

          if ( ParallelDescriptor::IOProcessor() )
              std::cout << ", Solve time: " << tmp[0] << ", CG time: " << tmp[1];
      }

      if ( ParallelDescriptor::IOProcessor() ) std::cout << '\n';
  }

  if ( ParallelDescriptor::IOProcessor() && (verbose > 0) )
  {
      if ( do_fixed_number_of_iters == 1)
      {
          std::cout << "   Did fixed number of iterations: " << maxiter << std::endl;
      } 
      else if ( error < eps_rel*norm_to_test_against )
      {
          std::cout << "   Converged res < eps_rel*max(bnorm,res_norm)\n";
      } 
      else if ( (use_Anorm_for_convergence == 1) && (error < eps_rel*norm_Lp*norm_cor) )
      {
          std::cout << "   Converged res < eps_rel*Anorm*sol\n";
      } 
      else if ( error < eps_abs )
      {
          std::cout << "   Converged res < eps_abs\n";
      }
  }

  //
  // Omit ghost update since maybe not initialized in calling routine.
  // Add to boundary values stored in initialsolution.
  //
  _sol.copy(*cor[level]);
  _sol.plus(*initialsolution,0,_sol.nComp(),0);

  if ( use_Anorm_for_convergence == 1 ) 
  {
     if ( do_fixed_number_of_iters == 1                ||
          error <= eps_rel*(norm_Lp*norm_cor+norm_to_test_against) ||
          error <= eps_abs )
       returnVal = 1;
  } 
  else 
  {
     if ( do_fixed_number_of_iters == 1 ||
          error <= eps_rel*(norm_to_test_against)   ||
          error <= eps_abs )
       returnVal = 1;
  } 

  //
  // Otherwise, failed to solve satisfactorily
  //
  return returnVal;
}

//.........这里部分代码省略.........
                      << rnorm/(rnorm0) << '\n';
        }

#ifdef CG_USE_OLD_CONVERGENCE_CRITERIA
        if ( rnorm < eps_rel*rnorm0 || rnorm < eps_abs ) break;
#else
        sol_norm = norm_inf(sol);
        if ( rnorm < eps_rel*(Lp_norm*sol_norm + rnorm0 ) || rnorm < eps_abs ) break;
#endif
        if ( use_mg_precond )
        {
            sh.setVal(0);
            mg_precond->solve(sh, s, eps_rel, eps_abs, temp_bc_mode);
        }
        else if ( use_jacobi_precond )
        {
            sh.setVal(0);
            Lp.jacobi_smooth(sh, s, lev, temp_bc_mode);
        }
        else
        {
            MultiFab::Copy(sh,s,0,0,1,0);
        }
        Lp.apply(t, sh, lev, temp_bc_mode);
        //
        // This is a little funky.  I want to elide one of the reductions
        // in the following two dotxy()s.  We do that by calculating the "local"
        // values and then reducing the two local values at the same time.
        //
        Real vals[2] = { dotxy(t,t,true), dotxy(t,s,true) };

        ParallelDescriptor::ReduceRealSum(vals,2,color());

        if ( vals[0] )
	{
            omega = vals[1]/vals[0];
	}
        else
	{
            ret = 3; break;
	}
        sxay(sol, sol,  omega, sh);
        sxay(r,     s, -omega,  t);

        rnorm = norm_inf(r);

        if ( verbose > 2 && ParallelDescriptor::IOProcessor(color()) )
        {
            Spacer(std::cout, lev);
            std::cout << "CGSolver_BiCGStab: Iteration "
                      << std::setw(11) << nit
                      << " rel. err. "
                      << rnorm/(rnorm0) << '\n';
        }

#ifdef CG_USE_OLD_CONVERGENCE_CRITERIA
        if ( rnorm < eps_rel*rnorm0 || rnorm < eps_abs ) break;
#else
        sol_norm = norm_inf(sol);
        if ( rnorm < eps_rel*(Lp_norm*sol_norm + rnorm0 ) || rnorm < eps_abs ) break;
#endif
        if ( omega == 0 )
	{
            ret = 4; break;
	}
        rho_1 = rho;
    }

    if ( verbose > 0 && ParallelDescriptor::IOProcessor(color()) )
    {
        Spacer(std::cout, lev);
        std::cout << "CGSolver_BiCGStab: Final: Iteration "
                  << std::setw(4) << nit
                  << " rel. err. "
                  << rnorm/(rnorm0) << '\n';
    }

#ifdef CG_USE_OLD_CONVERGENCE_CRITERIA
    if ( ret == 0 && rnorm > eps_rel*rnorm0 && rnorm > eps_abs)
#else
    if ( ret == 0 && rnorm > eps_rel*(Lp_norm*sol_norm + rnorm0 ) && rnorm > eps_abs )
#endif
    {
        if ( ParallelDescriptor::IOProcessor(color()) )
            BoxLib::Warning("CGSolver_BiCGStab:: failed to converge!");
        ret = 8;
    }

    if ( ( ret == 0 || ret == 8 ) && (rnorm < rnorm0) )
    {
        sol.plus(sorig, 0, 1, 0);
    } 
    else 
    {
        sol.setVal(0);
        sol.plus(sorig, 0, 1, 0);
    }

    return ret;
}

//.........这里部分代码省略.........
        }
        else
        {
            MultiFab::Copy(z,r,0,0,1,0);
        }

        Real rho = dotxy(z,r);

        if (nit == 1)
        {
            MultiFab::Copy(p,z,0,0,1,0);
        }
        else
        {
            Real beta = rho/rho_1;
            sxay(p, z, beta, p);
        }
        Lp.apply(q, p, lev, temp_bc_mode);

        Real alpha;
        if ( Real pw = dotxy(p,q) )
	{
            alpha = rho/pw;
	}
        else
	{
            ret = 1; break;
	}
        
        if ( verbose > 2 && ParallelDescriptor::IOProcessor(color()) )
        {
            Spacer(std::cout, lev);
            std::cout << "CGSolver_cg:"
                      << " nit " << nit
                      << " rho " << rho
                      << " alpha " << alpha << '\n';
        }
        sxay(sol, sol, alpha, p);
        sxay(  r,   r,-alpha, q);
        rnorm = norm_inf(r);
        sol_norm = norm_inf(sol);

        if ( verbose > 2 && ParallelDescriptor::IOProcessor(color()) )
        {
            Spacer(std::cout, lev);
            std::cout << "       CG:       Iteration"
                      << std::setw(4) << nit
                      << " rel. err. "
                      << rnorm/(rnorm0) << '\n';
        }

#ifdef CG_USE_OLD_CONVERGENCE_CRITERIA
        if ( rnorm < eps_rel*rnorm0 || rnorm < eps_abs ) break;
#else
        if ( rnorm < eps_rel*(Lp_norm*sol_norm + rnorm0) || rnorm < eps_abs ) break;
#endif
        if ( rnorm > def_unstable_criterion*minrnorm )
	{
            ret = 2; break;
	}
        else if ( rnorm < minrnorm )
	{
            minrnorm = rnorm;
	}

        rho_1 = rho;
    }
    
    if ( verbose > 0 && ParallelDescriptor::IOProcessor(color()) )
    {
        Spacer(std::cout, lev);
        std::cout << "       CG: Final Iteration"
                  << std::setw(4) << nit
                  << " rel. err. "
                  << rnorm/(rnorm0) << '\n';
    }

#ifdef CG_USE_OLD_CONVERGENCE_CRITERIA
    if ( ret == 0 &&  rnorm > eps_rel*rnorm0 && rnorm > eps_abs )
#else
    if ( ret == 0 && rnorm > eps_rel*(Lp_norm*sol_norm + rnorm0) && rnorm > eps_abs )
#endif
    {
        if ( ParallelDescriptor::IOProcessor(color()) )
            BoxLib::Warning("CGSolver_cg: failed to converge!");
        ret = 8;
    }

    if ( ( ret == 0 || ret == 8 ) && (rnorm < rnorm0) )
    {
        sol.plus(sorig, 0, 1, 0);
    } 
    else 
    {
        sol.setVal(0);
        sol.plus(sorig, 0, 1, 0);
    }

    return ret;
}

//.........这里部分代码省略.........
  const Real error0    = tmp[1];
  int        returnVal = 0;
  Real       error     = error0;

  if ( ParallelDescriptor::IOProcessor() && (verbose > 0) )
  {
      Spacer(std::cout, level);
      std::cout << "MCMultiGrid: Initial rhs                = " << norm_rhs << '\n';
      std::cout << "MCMultiGrid: Initial error (error0)     = " << error0 << '\n';
  }
  
  if ( ParallelDescriptor::IOProcessor() && eps_rel < 1.0e-16 && eps_rel > 0 )
  {
      std::cout << "MCMultiGrid: Tolerance "
                << eps_rel
                << " < 1e-16 is probably set too low" << '\n';
  }
  //
  // Initialize correction to zero at this level (auto-filled at levels below)
  //
  (*cor[level]).setVal(0.0);
  //
  // Note: if eps_rel, eps_abs < 0 then that test is effectively bypassed.
  //
  int        nit         = 1;
  const Real new_error_0 = norm_rhs;
  //const Real norm_Lp     = Lp.norm(0, level);


  for ( ;
        error > eps_abs &&
          error > eps_rel*norm_rhs &&
          nit <= maxiter;
        ++nit)
  {
    relax(*cor[level], *rhs[level], level, eps_rel, eps_abs, bc_mode);

    error = errorEstimate(level,bc_mode);
	
    if ( ParallelDescriptor::IOProcessor() && verbose > 1 )
    {
      const Real rel_error = (error0 != 0) ? error/new_error_0 : 0;
      Spacer(std::cout, level);
      std::cout << "MCMultiGrid: Iteration   "
                << nit
                << " error/error0 = "
                << rel_error << '\n';
    }
  }

  Real run_time = (ParallelDescriptor::second() - strt_time);
  if ( verbose > 0 )
  {
      if ( ParallelDescriptor::IOProcessor() )
      {
          const Real rel_error = (error0 != 0) ? error/error0 : 0;
          Spacer(std::cout, level);
          std::cout << "MCMultiGrid: Final Iter. "
                    << nit-1
                    << " error/error0 = "
                    << rel_error;
      }

      if ( verbose > 1 )
      {
        
        ParallelDescriptor::ReduceRealMax(run_time);

        if ( ParallelDescriptor::IOProcessor() )
          std::cout << ", Solve time: " << run_time << '\n';
      }
  }

  if ( ParallelDescriptor::IOProcessor() && (verbose > 0) )
  {
    if ( error < eps_rel*norm_rhs )
    {
      std::cout << "   Converged res < eps_rel*bnorm\n";
    } 
    else if ( error < eps_abs )
    {
      std::cout << "   Converged res < eps_abs\n";
    }
  }
  //
  // Omit ghost update since maybe not initialized in calling routine.
  // Add to boundary values stored in initialsolution.
  //
  _sol.copy(*cor[level]);
  _sol.plus(*initialsolution,0,_sol.nComp(),0);

  if ( error <= eps_rel*(norm_rhs) ||
       error <= eps_abs )
    returnVal = 1;

  //
  // Otherwise, failed to solve satisfactorily
  //
  return returnVal;
}