自学教程：C++ CHKERRABORT函数代码示例

void ApplyInvDavidsonDiagPrecPETSc(void *x, PRIMME_INT *ldx, void *y,      PRIMME_INT *ldy, int *blockSize, primme_params *primme, int *err) {   int i, j;   double shift, d, minDenominator;   SCALAR *xvec, *yvec;   const int nLocal = primme->nLocal, bs = *blockSize;   const PetscScalar *diag;   Vec vec;   PetscErrorCode ierr;      vec = *(Vec *)primme->preconditioner;   xvec = (SCALAR *)x;   yvec = (SCALAR *)y;   minDenominator = 1e-14*(primme->aNorm >= 0.0L ? primme->aNorm : 1.);   ierr = VecGetArrayRead(vec, &diag); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);   for (i=0; i<bs; i++) {      shift = primme->ShiftsForPreconditioner[i];      for (j=0; j<nLocal; j++) {         d = diag[j] - shift;         d = (fabs(d) > minDenominator) ? d : copysign(minDenominator, d);         yvec[*ldy*i+j] = xvec[*ldx*i+j]/d;      }   }   ierr = VecRestoreArrayRead(vec, &diag); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);   *err = 0;}

 void prod( sparse_matrix_type const& A,            vector_type const& x,            vector_type& b ) const {     int ierr = 0;     petsc_sparse_matrix_type const& _A = dynamic_cast<petsc_sparse_matrix_type const&>( A );     petsc_vector_type const& _x = dynamic_cast<petsc_vector_type const&>( x );     petsc_vector_type const& _b = dynamic_cast<petsc_vector_type const&>( b );     if ( _A.mapCol().worldComm().globalSize() == x.map().worldComm().globalSize() )     {         //std::cout << "BackendPetsc::prod STANDART"<< std::endl;         ierr = MatMult( _A.mat(), _x.vec(), _b.vec() );         CHKERRABORT( _A.comm().globalComm(),ierr );     }     else     {         //std::cout << "BackendPetsc::prod with convert"<< std::endl;         auto x_convert = petscMPI_vector_type(_A.mapColPtr());         x_convert.duplicateFromOtherPartition(x);         x_convert.close();         ierr = MatMult( _A.mat(), x_convert.vec(), _b.vec() );         CHKERRABORT( _A.comm().globalComm(),ierr );     }     b.close(); }

voidPetscOutputter::timestepSetupInternal(){// Only execute if PETSc exists#ifdef LIBMESH_HAVE_PETSC  // Extract the non-linear and linear solvers from PETSc  NonlinearSystem & nl = _problem_ptr->getNonlinearSystem();  PetscNonlinearSolver<Number> * petsc_solver = dynamic_cast<PetscNonlinearSolver<Number> *>(nl.sys().nonlinear_solver.get());  SNES snes = petsc_solver->snes();  KSP ksp;  SNESGetKSP(snes, &ksp);  // Update the pseudo times  _nonlinear_time = _time_old;                   // non-linear time starts with the previous time step  _nonlinear_dt = _dt/_nonlinear_dt_divisor;     // set the pseudo non-linear timestep  _linear_dt = _nonlinear_dt/_linear_dt_divisor; // set the pseudo linear timestep  // Set the PETSc monitor functions  if (_output_nonlinear || (_time >= _nonlinear_start_time - _t_tol && _time <= _nonlinear_end_time + _t_tol) )  {    PetscErrorCode ierr = SNESMonitorSet(snes, petscNonlinearOutput, this, PETSC_NULL);    CHKERRABORT(libMesh::COMM_WORLD,ierr);  }  if (_output_linear || (_time >= _linear_start_time - _t_tol && _time <= _linear_end_time + _t_tol) )  {    PetscErrorCode ierr = KSPMonitorSet(ksp, petscLinearOutput, this, PETSC_NULL);    CHKERRABORT(libMesh::COMM_WORLD,ierr);  }#endif}

void PetscDiffSolver::init (){  START_LOG("init()", "PetscDiffSolver");  Parent::init();  int ierr=0;#if PETSC_VERSION_LESS_THAN(2,1,2)  // At least until Petsc 2.1.1, the SNESCreate had a different  // calling syntax.  The second argument was of type SNESProblemType,  // and could have a value of either SNES_NONLINEAR_EQUATIONS or  // SNES_UNCONSTRAINED_MINIMIZATION.  ierr = SNESCreate(libMesh::COMM_WORLD, SNES_NONLINEAR_EQUATIONS, &_snes);  CHKERRABORT(libMesh::COMM_WORLD,ierr);#else  ierr = SNESCreate(libMesh::COMM_WORLD,&_snes);  CHKERRABORT(libMesh::COMM_WORLD,ierr);#endif#if PETSC_VERSION_LESS_THAN(2,3,3)  ierr = SNESSetMonitor (_snes, __libmesh_petsc_diff_solver_monitor,                         this, PETSC_NULL);#else  // API name change in PETSc 2.3.3  ierr = SNESMonitorSet (_snes, __libmesh_petsc_diff_solver_monitor,                         this, PETSC_NULL);#endif  CHKERRABORT(libMesh::COMM_WORLD,ierr);  ierr = SNESSetFromOptions(_snes);  CHKERRABORT(libMesh::COMM_WORLD,ierr);  STOP_LOG("init()", "PetscDiffSolver");}

 PetscTimeStepper(FEProblem &feproblem) : TimeStepper(feproblem) {   PetscErrorCode ierr;   ierr = TSCreate(libMesh::COMM_WORLD, &this->_ts);   CHKERRABORT(libMesh::COMM_WORLD,ierr);   ierr = TSSetApplicationContext(this->_ts,this);   CHKERRABORT(libMesh::COMM_WORLD,ierr); };

std::pair<unsigned int, unsigned int>SlepcEigenSolver<T>::solve_standard (ShellMatrix<T> &shell_matrix,				     int nev,                  // number of requested eigenpairs				     int ncv,                  // number of basis vectors				     const double tol,         // solver tolerance				     const unsigned int m_its) // maximum number of iterations{  this->init ();  int ierr=0;  // Prepare the matrix.  Mat mat;  ierr = MatCreateShell(libMesh::COMM_WORLD,            shell_matrix.m(), // Specify the number of local rows            shell_matrix.n(), // Specify the number of local columns            PETSC_DETERMINE,            PETSC_DETERMINE,            const_cast<void*>(static_cast<const void*>(&shell_matrix)),            &mat);  /* Note that the const_cast above is only necessary because PETSc     does not accept a const void*.  Inside the member function     _petsc_shell_matrix() below, the pointer is casted back to a     const ShellMatrix<T>*.  */  CHKERRABORT(libMesh::COMM_WORLD,ierr);  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));  CHKERRABORT(libMesh::COMM_WORLD,ierr);  return _solve_standard_helper(mat, nev, ncv, tol, m_its);}

int TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void *ctx){  TS              ts = (TS) ctx;  MPI_Comm        comm = ((PetscObject)ts)->comm;  TS_Sundials     *cvode = (TS_Sundials*)ts->data;  Vec             yy = cvode->w1,yyd = cvode->w2,yydot = cvode->ydot;  PetscScalar     *y_data,*ydot_data;  PetscErrorCode  ierr;  PetscFunctionBegin;  /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/  y_data     = (PetscScalar *) N_VGetArrayPointer(y);  ydot_data  = (PetscScalar *) N_VGetArrayPointer(ydot);  ierr = VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr);  ierr = VecPlaceArray(yyd,ydot_data); CHKERRABORT(comm,ierr);  /* now compute the right hand side function */  if (!ts->userops->ifunction) {    ierr = TSComputeRHSFunction(ts,t,yy,yyd);CHKERRQ(ierr);  } else {                      /* If rhsfunction is also set, this computes both parts and shifts them to the right */    ierr = VecZeroEntries(yydot);CHKERRQ(ierr);    ierr = TSComputeIFunction(ts,t,yy,yydot,yyd,PETSC_FALSE); CHKERRABORT(comm,ierr);    ierr = VecScale(yyd,-1.);CHKERRQ(ierr);  }  ierr = VecResetArray(yy); CHKERRABORT(comm,ierr);  ierr = VecResetArray(yyd); CHKERRABORT(comm,ierr);  PetscFunctionReturn(0);}

int TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void *ctx){  TS             ts = (TS) ctx;  DM             dm;  DMTS           tsdm;  TSIFunction    ifunction;  MPI_Comm       comm;  TS_Sundials    *cvode = (TS_Sundials*)ts->data;  Vec            yy     = cvode->w1,yyd = cvode->w2,yydot = cvode->ydot;  PetscScalar    *y_data,*ydot_data;  PetscErrorCode ierr;  PetscFunctionBegin;  ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr);  /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/  y_data    = (PetscScalar*) N_VGetArrayPointer(y);  ydot_data = (PetscScalar*) N_VGetArrayPointer(ydot);  ierr      = VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr);  ierr      = VecPlaceArray(yyd,ydot_data);CHKERRABORT(comm,ierr);  /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);  ierr = DMGetDMTS(dm,&tsdm);CHKERRQ(ierr);  ierr = DMTSGetIFunction(dm,&ifunction,NULL);CHKERRQ(ierr);  if (!ifunction) {    ierr = TSComputeRHSFunction(ts,t,yy,yyd);CHKERRQ(ierr);  } else {                      /* If rhsfunction is also set, this computes both parts and shifts them to the right */    ierr = VecZeroEntries(yydot);CHKERRQ(ierr);    ierr = TSComputeIFunction(ts,t,yy,yydot,yyd,PETSC_FALSE);CHKERRABORT(comm,ierr);    ierr = VecScale(yyd,-1.);CHKERRQ(ierr);  }  ierr = VecResetArray(yy);CHKERRABORT(comm,ierr);  ierr = VecResetArray(yyd);CHKERRABORT(comm,ierr);  PetscFunctionReturn(0);}

voidSolverLinearPetsc<T>::clear (){    PetscBool pinit;    PetscInitialized( &pinit );    if ( pinit && this->initialized() )    {        this->setInitialized( false );        int ierr=0;        // 2.1.x & earlier style#if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR <= 1)        ierr = SLESDestroy( M_sles );        CHKERRABORT( this->worldComm().globalComm(),ierr );        // 2.2.0 & newer style#else        FEELPP_ASSERT( M_ksp != 0 ).error( "invalid ksp" );        ierr = PETSc::KSPDestroy( M_ksp );        CHKERRABORT( this->worldComm().globalComm(),ierr );#endif        // Mimic PETSc default solver and preconditioner        this->setSolverType(  GMRES );        if ( this->worldComm().globalComm().size() == 1 )            this->setPreconditionerType( LU_PRECOND );        else            this->setPreconditionerType( BLOCK_JACOBI_PRECOND );    }}

 virtual void setStepLimits(Real dtmin,Real dtmax) {   PetscErrorCode ierr;   TSAdapt adapt;   ierr = TSGetAdapt(this->_ts,&adapt);   CHKERRABORT(libMesh::COMM_WORLD,ierr);   ierr = TSAdaptSetStepLimits(adapt,dtmin,dtmax);   CHKERRABORT(libMesh::COMM_WORLD,ierr); }

inline void PetscVector::init (const int n,                               const int n_local,                               const std::vector<int>& ghost,                               const bool fast,                               const ParallelType type) {  int ierr=0;  PetscInt petsc_n=static_cast<int>(n);  PetscInt petsc_n_local=static_cast<int>(n_local);  PetscInt petsc_n_ghost=static_cast<int>(ghost.size());  // If the mesh is not disjoint, every processor will either have  // all the dofs, none of the dofs, or some non-zero dofs at the  // boundary between processors.  //  // However we can't assert this, because someone might want to  // construct a GHOSTED vector which doesn't include neighbor element  // dofs.  Boyce tried to do so in user code, and we're going to want  // to do so in System::project_vector().  //  // libmesh_assert(n_local == 0 || n_local == n || !ghost.empty());  assert(sizeof(PetscInt) == sizeof(int));  // If the mesh is disjoint, the following assertion will fail.  // If the mesh is not disjoint, every processor will either have  // all the dofs, none of the dofs, or some non-zero dofs at the  // boundary between processors.  //assert(n_local == 0 || n_local == n || !ghost.empty());  PetscInt* petsc_ghost = ghost.empty() ? PETSC_NULL :                     const_cast<int*>(reinterpret_cast<const PetscInt*>(&ghost[0]));  // Clear initialized vectors  if (this->initialized())   this->clear();  assert(type == AUTOMATIC || type == GHOSTED);  this->_type = GHOSTED;  /* Make the global-to-local ghost cell map.  */  for (int i=0; i<(int)ghost.size(); i++){    _global_to_local_map[ghost[i]] = i;  }  /* Create vector.  */  ierr = VecCreateGhost (MPI_COMM_WORLD, petsc_n_local, petsc_n,                         petsc_n_ghost, petsc_ghost, &_vec);  CHKERRABORT(MPI_COMM_WORLD,ierr);  ierr = VecSetFromOptions (_vec);  CHKERRABORT(MPI_COMM_WORLD,ierr);  this->_is_initialized = true;  this->_is_closed = true;  if (fast == false)    this->zero ();}

  std::pair<unsigned int, Real>  PetscDMNonlinearSolver<T>::solve (SparseMatrix<T>& jac_in,  // System Jacobian Matrix				    NumericVector<T>& x_in,   // Solution vector				    NumericVector<T>& r_in,   // Residual vector				    const double,             // Stopping tolerance				    const unsigned int)  {    START_LOG("solve()", "PetscNonlinearSolver");    this->init ();    // Make sure the data passed in are really of Petsc types    libmesh_cast_ptr<PetscMatrix<T>*>(&jac_in);    libmesh_cast_ptr<PetscVector<T>*>(&r_in);    // Extract solution vector    PetscVector<T>* x = libmesh_cast_ptr<PetscVector<T>*>(&x_in);    int ierr=0;    int n_iterations =0;    // Should actually be a PetscReal, but I don't know which version of PETSc first introduced PetscReal    Real final_residual_norm=0.;    if (this->user_presolve)      this->user_presolve(this->system());    //Set the preconditioning matrix    if (this->_preconditioner)      this->_preconditioner->set_matrix(jac_in);    ierr = SNESSolve (this->_snes, PETSC_NULL, x->vec());    CHKERRABORT(libMesh::COMM_WORLD,ierr);    ierr = SNESGetIterationNumber(this->_snes,&n_iterations);    CHKERRABORT(libMesh::COMM_WORLD,ierr);    ierr = SNESGetLinearSolveIterations(this->_snes, &this->_n_linear_iterations);    CHKERRABORT(libMesh::COMM_WORLD,ierr);    ierr = SNESGetFunctionNorm(this->_snes,&final_residual_norm);    CHKERRABORT(libMesh::COMM_WORLD,ierr);    // Get and store the reason for convergence    SNESGetConvergedReason(this->_snes, &this->_reason);    //Based on Petsc 2.3.3 documentation all diverged reasons are negative    this->converged = (this->_reason >= 0);    this->clear();    STOP_LOG("solve()", "PetscNonlinearSolver");    // return the # of its. and the final residual norm.    return std::make_pair(n_iterations, final_residual_norm);  }

 virtual TimeStepperStatus step(Real *ftime) {   PetscErrorCode ierr;   TSConvergedReason reason;   ierr = TSStep(_ts);   CHKERRABORT(libMesh::COMM_WORLD,ierr);   ierr = TSGetConvergedReason(_ts,&reason);   CHKERRABORT(libMesh::COMM_WORLD,ierr);   ierr = TSGetTime(_ts,ftime);   CHKERRABORT(libMesh::COMM_WORLD,ierr);   return (TimeStepperStatus)reason; }

void PETScMatvec(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize, primme_params *primme, int *err) {   Mat *matrix;   PetscInt m, n, mLocal, nLocal;   PetscErrorCode ierr;   matrix = (Mat *)primme->matrix;   ierr = MatGetSize(*matrix, &m, &n);  CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);   ierr = MatGetLocalSize(*matrix, &mLocal, &nLocal);  CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);   assert(m == primme->n && n == primme->n && mLocal == primme->nLocal         && nLocal == primme->nLocal);   PETScMatvecGen(x, *ldx, y, *ldy, *blockSize, 0, *matrix, *(MPI_Comm*)primme->commInfo);   *err = 0;}

  void PetscDMRegister()  {    if (PetscDMRegistered)      return;    PetscErrorCode ierr;#if PETSC_RELEASE_LESS_THAN(3,4,0)    ierr = DMRegister(DMLIBMESH, PETSC_NULL, "DMCreate_libMesh", DMCreate_libMesh); CHKERRABORT(libMesh::COMM_WORLD,ierr);#else    ierr = DMRegister(DMLIBMESH, DMCreate_libMesh); CHKERRABORT(libMesh::COMM_WORLD,ierr);#endif    PetscDMRegistered = PETSC_TRUE;  }

    /**     * Call the assemble functions     */    void close ()    {        FEELPP_ASSERT ( this->isInitialized() ).error( "VectorPetsc<> not initialized" );        int ierr=0;        ierr = VecAssemblyBegin( M_vec );        CHKERRABORT( this->comm(),ierr );        ierr = VecAssemblyEnd( M_vec );        CHKERRABORT( this->comm(),ierr );        this->M_is_closed = true;    }

 virtual void setupInternal(NumericVector<Number> &X) {   PetscVector<Number> *pX = cast_ptr<PetscVector<Number> *>(&X);   PetscErrorCode ierr;   ierr = TSSetIFunction(this->_ts,PETSC_NULL,this->_computeIFunction,this);   CHKERRABORT(libMesh::COMM_WORLD,ierr);   PetscMatrix<Number> *mat = cast_ptr<PetscMatrix<Number> *>(this->_fe_problem.getNonlinearSystem().sys().matrix);   Mat pmat = mat->mat();   ierr = TSSetIJacobian(this->_ts,pmat,pmat,this->_computeIJacobian,this);   CHKERRABORT(libMesh::COMM_WORLD,ierr);   ierr = TSSetFromOptions(this->_ts);   CHKERRABORT(libMesh::COMM_WORLD,ierr);   ierr = TSSetSolution(_ts,pX->vec());   CHKERRABORT(libMesh::COMM_WORLD,ierr); }

inline void PetscVector::close () {  this->_restore_array();  int ierr=0;  ierr = VecAssemblyBegin(_vec);  					CHKERRABORT(MPI_COMM_WORLD,ierr);  ierr = VecAssemblyEnd(_vec);  					CHKERRABORT(MPI_COMM_WORLD,ierr);  if (this->type() == GHOSTED) {    ierr = VecGhostUpdateBegin(_vec,INSERT_VALUES,SCATTER_FORWARD);  	CHKERRABORT(MPI_COMM_WORLD,ierr);    ierr = VecGhostUpdateEnd(_vec,INSERT_VALUES,SCATTER_FORWARD);  	CHKERRABORT(MPI_COMM_WORLD,ierr);  }  this->_is_closed = true;}

 void PetscPreconditioner<T>::set_petsc_subpreconditioner_type(const PCType type, PC& pc)#endif{  // For catching PETSc error return codes  int ierr = 0;  // get the communicator from the PETSc object  Parallel::communicator comm;  PetscObjectGetComm((PetscObject)pc, &comm);  Parallel::Communicator communicator(comm);  // All docs say must call KSPSetUp or PCSetUp before calling PCBJacobiGetSubKSP.  // You must call PCSetUp after the preconditioner operators have been set, otherwise you get the:  //  // "Object is in wrong state!"  // "Matrix must be set first."  //  // error messages...  ierr = PCSetUp(pc);  CHKERRABORT(comm,ierr);  // To store array of local KSP contexts on this processor  KSP* subksps;  // the number of blocks on this processor  PetscInt n_local;  // The global number of the first block on this processor.  // This is not used, so we just pass PETSC_NULL instead.  // int first_local;  // Fill array of local KSP contexts  ierr = PCBJacobiGetSubKSP(pc, &n_local, PETSC_NULL, &subksps);  CHKERRABORT(comm,ierr);  // Loop over sub-ksp objects, set ILU preconditioner  for (PetscInt i=0; i<n_local; ++i)    {      // Get pointer to sub KSP object's PC      PC subpc;      ierr = KSPGetPC(subksps[i], &subpc);      CHKERRABORT(comm,ierr);      // Set requested type on the sub PC      ierr = PCSetType(subpc, type);      CHKERRABORT(comm,ierr);    }}

void stdNormalArray(double *eps, int size, PetscRandom ran){  int            i;  double         u1,u2,t;  PetscErrorCode ierr;  for (i=0;i<size;i+=2){    ierr = PetscRandomGetValue(ran,(PetscScalar*)&u1);CHKERRABORT(PETSC_COMM_WORLD,ierr);    ierr = PetscRandomGetValue(ran,(PetscScalar*)&u2);CHKERRABORT(PETSC_COMM_WORLD,ierr);    t = sqrt(-2*log(u1));    eps[i] = t * cos(2*PI*u2);    eps[i+1] = t * sin(2*PI*u2);  }}

voidSolverLinearPetsc<T>::getResidualHistory( std::vector<double>& hist ){    int ierr = 0;    int its  = 0;    // Fill the residual history vector with the residual norms    // Note that GetResidualHistory() does not copy any values, it    // simply sets the pointer p.  Note that for some Krylov subspace    // methods, the number of residuals returned in the history    // vector may be different from what you are expecting.  For    // example, TFQMR returns two residual values per iteration step.    double* p;    ierr = KSPGetResidualHistory( M_ksp, &p, &its );    CHKERRABORT( this->worldComm().globalComm(),ierr );    // Check for early return    if ( its == 0 ) return;    // Create space to store the result    hist.resize( its );    // Copy history into the vector provided by the user.    for ( int i=0; i<its; ++i )    {        hist[i] = *p;        p++;    }}

// =======================================================/// This function initializes the libraries if it is parallelFemusInit::FemusInit(    int & argc,            // integer program input    char** & argv,         // char program input    MPI_Comm comm_world_in // communicator for MPI direct) {// ======================================================#ifdef HAVE_PETSC      int ierr = PetscInitialize (&argc, &argv, NULL, NULL);    CHKERRABORT(PETSC_COMM_WORLD,ierr);#endif#ifdef HAVE_MPI    // redirect libMesh::out to nothing on all    // other processors unless explicitly told    // not to via the --keep-cout command-line argument.    int i;    MPI_Comm_rank(MPI_COMM_WORLD, &i);    if ( i != 0 ) {        std::cout.rdbuf(NULL);    }#endif   std::cout << " FemusInit(): PETSC_COMM_WORLD initialized" << std::endl << std::endl;    return;}