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

PetscErrorCode  TSGLAdaptCreate(MPI_Comm comm,TSGLAdapt *inadapt){  PetscErrorCode ierr;  TSGLAdapt      adapt;  PetscFunctionBegin;  *inadapt = NULL;  ierr     = PetscHeaderCreate(adapt,TSGLADAPT_CLASSID,"TSGLAdapt","General Linear adaptivity","TS",comm,TSGLAdaptDestroy,TSGLAdaptView);CHKERRQ(ierr);  *inadapt = adapt;  PetscFunctionReturn(0);}

PetscErrorCode  PetscViewerFileSetName_ASCII(PetscViewer viewer,const char name[]){  PetscErrorCode    ierr;  size_t            len;  char              fname[PETSC_MAX_PATH_LEN],*gz;  PetscViewer_ASCII *vascii = (PetscViewer_ASCII*)viewer->data;  PetscBool         isstderr,isstdout;  PetscMPIInt       rank;  PetscFunctionBegin;  ierr = PetscViewerFileClose_ASCII(viewer);CHKERRQ(ierr);  if (!name) PetscFunctionReturn(0);  ierr = PetscStrallocpy(name,&vascii->filename);CHKERRQ(ierr);  /* Is this file to be compressed */  vascii->storecompressed = PETSC_FALSE;  ierr = PetscStrstr(vascii->filename,".gz",&gz);CHKERRQ(ierr);  if (gz) {    ierr = PetscStrlen(gz,&len);CHKERRQ(ierr);    if (len == 3) {      *gz = 0;      vascii->storecompressed = PETSC_TRUE;    }  }  ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)viewer),&rank);CHKERRQ(ierr);  if (!rank) {    ierr = PetscStrcmp(name,"stderr",&isstderr);CHKERRQ(ierr);    ierr = PetscStrcmp(name,"stdout",&isstdout);CHKERRQ(ierr);    /* empty filename means stdout */    if (name[0] == 0)  isstdout = PETSC_TRUE;    if (isstderr)      vascii->fd = PETSC_STDERR;    else if (isstdout) vascii->fd = PETSC_STDOUT;    else {      ierr = PetscFixFilename(name,fname);CHKERRQ(ierr);      switch (vascii->mode) {      case FILE_MODE_READ:        vascii->fd = fopen(fname,"r");        break;      case FILE_MODE_WRITE:        vascii->fd = fopen(fname,"w");        break;      case FILE_MODE_APPEND:        vascii->fd = fopen(fname,"a");        break;      case FILE_MODE_UPDATE:        vascii->fd = fopen(fname,"r+");        if (!vascii->fd) vascii->fd = fopen(fname,"w+");        break;      case FILE_MODE_APPEND_UPDATE:        /* I really want a file which is opened at the end for updating,           not a+, which opens at the beginning, but makes writes at the end.        */        vascii->fd = fopen(fname,"r+");        if (!vascii->fd) vascii->fd = fopen(fname,"w+");        else {          ierr     = fseek(vascii->fd, 0, SEEK_END);CHKERRQ(ierr);        }        break;      default:        SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Invalid file mode %d", vascii->mode);      }      if (!vascii->fd) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open PetscViewer file: %s",fname);    }  }#if defined(PETSC_USE_LOG)  PetscLogObjectState((PetscObject)viewer,"File: %s",name);#endif  PetscFunctionReturn(0);}

PetscErrorCode PetscViewerDestroy_ASCII_Subcomm(PetscViewer viewer){  PetscViewer_ASCII *vascii = (PetscViewer_ASCII*)viewer->data;  PetscErrorCode    ierr;  PetscFunctionBegin;  ierr = PetscViewerRestoreSubcomm(vascii->subviewer,PetscObjectComm((PetscObject)viewer),&viewer);CHKERRQ(ierr);  vascii->subviewer = NULL;  PetscFunctionReturn(0);}

EXTERN_C_END#undef __FUNCT__#define __FUNCT__ "PetscViewerRegisterAll"/*@C  PetscViewerRegisterAll - Registers all of the graphics methods in the PetscViewer package.  Not Collective   Level: developer.seealso:  PetscViewerRegisterDestroy()@*/PetscErrorCode  PetscViewerRegisterAll(const char *path){  PetscErrorCode ierr;  PetscFunctionBegin;  ierr = PetscViewerRegisterDynamic(PETSCVIEWERASCII,      path,"PetscViewerCreate_ASCII",      PetscViewerCreate_ASCII);CHKERRQ(ierr);  ierr = PetscViewerRegisterDynamic(PETSCVIEWERBINARY,     path,"PetscViewerCreate_Binary",     PetscViewerCreate_Binary);CHKERRQ(ierr);  ierr = PetscViewerRegisterDynamic(PETSCVIEWERSTRING,     path,"PetscViewerCreate_String",     PetscViewerCreate_String);CHKERRQ(ierr);  ierr = PetscViewerRegisterDynamic(PETSCVIEWERDRAW,       path,"PetscViewerCreate_Draw",       PetscViewerCreate_Draw);CHKERRQ(ierr);#if defined(PETSC_USE_SOCKET_VIEWER)  ierr = PetscViewerRegisterDynamic(PETSCVIEWERSOCKET,     path,"PetscViewerCreate_Socket",     PetscViewerCreate_Socket);CHKERRQ(ierr);#endif#if defined(PETSC_HAVE_MATHEMATICA)  ierr = PetscViewerRegisterDynamic(PETSCVIEWERMATHEMATICA,path,"PetscViewerCreate_Mathematica",PetscViewerCreate_Mathematica);CHKERRQ(ierr);#endif  ierr = PetscViewerRegisterDynamic(PETSCVIEWERVU,         path,"PetscViewerCreate_VU",         PetscViewerCreate_VU);CHKERRQ(ierr);#if defined(PETSC_HAVE_HDF5)  ierr = PetscViewerRegisterDynamic(PETSCVIEWERHDF5,       path,"PetscViewerCreate_HDF5",       PetscViewerCreate_HDF5);CHKERRQ(ierr);#endif#if defined(PETSC_HAVE_MATLAB_ENGINE)  ierr = PetscViewerRegisterDynamic(PETSCVIEWERMATLAB,     path,"PetscViewerCreate_Matlab",     PetscViewerCreate_Matlab);CHKERRQ(ierr);#endif#if defined(PETSC_HAVE_AMS)  ierr = PetscViewerRegisterDynamic(PETSCVIEWERAMS,        path,"PetscViewerCreate_AMS",        PetscViewerCreate_AMS);CHKERRQ(ierr);#endif  ierr = PetscViewerRegisterDynamic(PETSCVIEWERVTK,        path,"PetscViewerCreate_VTK",        PetscViewerCreate_VTK);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*@C    PetscViewerASCIIPrintf - Prints to a file, only from the first    processor in the PetscViewer    Not Collective, but only first processor in set has any effect    Input Parameters:+    viewer - optained with PetscViewerASCIIOpen()-    format - the usual printf() format string    Level: developer    Fortran Note:    The call sequence is PetscViewerASCIIPrintf(PetscViewer, character(*), int ierr) from Fortran.    That is, you can only pass a single character string from Fortran.  Concepts: PetscViewerASCII^printing  Concepts: printing^to file  Concepts: printf.seealso: PetscPrintf(), PetscSynchronizedPrintf(), PetscViewerASCIIOpen(),          PetscViewerASCIIPushTab(), PetscViewerASCIIPopTab(), PetscViewerASCIISynchronizedPrintf(),          PetscViewerCreate(), PetscViewerDestroy(), PetscViewerSetType(), PetscViewerASCIIGetPointer(), PetscViewerASCIISynchronizedAllow()@*/PetscErrorCode  PetscViewerASCIIPrintf(PetscViewer viewer,const char format[],...){  PetscViewer_ASCII *ascii = (PetscViewer_ASCII*)viewer->data;  PetscMPIInt       rank;  PetscInt          tab,intab = ascii->tab;  PetscErrorCode    ierr;  FILE              *fd = ascii->fd;  PetscBool         iascii,issingleton = PETSC_FALSE;  int               err;  PetscFunctionBegin;  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,1);  PetscValidCharPointer(format,2);  ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);  if (!iascii) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not ASCII PetscViewer");  if (ascii->bviewer) {    viewer      = ascii->bviewer;    ascii       = (PetscViewer_ASCII*)viewer->data;    fd          = ascii->fd;    issingleton = PETSC_TRUE;  }  ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)viewer),&rank);CHKERRQ(ierr);  if (!rank) {    va_list Argp;    tab = intab;    while (tab--) {      ierr = PetscFPrintf(PETSC_COMM_SELF,fd,"  ");CHKERRQ(ierr);    }    va_start(Argp,format);    ierr = (*PetscVFPrintf)(fd,format,Argp);CHKERRQ(ierr);    err  = fflush(fd);    if (err) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SYS,"fflush() failed on file");    if (petsc_history) {      va_start(Argp,format);      tab = intab;      while (tab--) {        ierr = PetscFPrintf(PETSC_COMM_SELF,petsc_history,"  ");CHKERRQ(ierr);      }      ierr = (*PetscVFPrintf)(petsc_history,format,Argp);CHKERRQ(ierr);      err  = fflush(petsc_history);      if (err) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SYS,"fflush() failed on file");    }    va_end(Argp);  } else if (issingleton) {    char        *string;    va_list     Argp;    size_t      fullLength;    PrintfQueue next;    ierr = PetscNew(&next);CHKERRQ(ierr);    if (petsc_printfqueue) {      petsc_printfqueue->next = next;      petsc_printfqueue       = next;    } else {      petsc_printfqueuebase = petsc_printfqueue = next;    }    petsc_printfqueuelength++;    next->size = QUEUESTRINGSIZE;    ierr       = PetscCalloc1(next->size, &next->string);CHKERRQ(ierr);    string     = next->string;    tab        = intab;    tab       *= 2;    while (tab--) {      *string++ = ' ';    }    va_start(Argp,format);    ierr = PetscVSNPrintf(string,next->size-2*ascii->tab,format,&fullLength,Argp);CHKERRQ(ierr);    va_end(Argp);  }  PetscFunctionReturn(0);}

int main(int argc, char *argv[]){  PetscErrorCode ierr;  DM             dau,dak,pack;  const PetscInt *lxu;  PetscInt       *lxk,m,sizes;  User           user;  SNES           snes;  Vec            X,F,Xu,Xk,Fu,Fk;  Mat            B;  IS             *isg;  PetscBool      view_draw,pass_dm;  PetscInitialize(&argc,&argv,0,help);  ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-10,1,1,NULL,&dau);CHKERRQ(ierr);  ierr = DMSetOptionsPrefix(dau,"u_");CHKERRQ(ierr);  ierr = DMSetFromOptions(dau);CHKERRQ(ierr);  ierr = DMDAGetOwnershipRanges(dau,&lxu,0,0);CHKERRQ(ierr);  ierr = DMDAGetInfo(dau,0, &m,0,0, &sizes,0,0, 0,0,0,0,0,0);CHKERRQ(ierr);  ierr = PetscMalloc1(sizes,&lxk);CHKERRQ(ierr);  ierr = PetscMemcpy(lxk,lxu,sizes*sizeof(*lxk));CHKERRQ(ierr);  lxk[0]--;  ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,m-1,1,1,lxk,&dak);CHKERRQ(ierr);  ierr = DMSetOptionsPrefix(dak,"k_");CHKERRQ(ierr);  ierr = DMSetFromOptions(dak);CHKERRQ(ierr);  ierr = PetscFree(lxk);CHKERRQ(ierr);  ierr = DMCompositeCreate(PETSC_COMM_WORLD,&pack);CHKERRQ(ierr);  ierr = DMSetOptionsPrefix(pack,"pack_");CHKERRQ(ierr);  ierr = DMCompositeAddDM(pack,dau);CHKERRQ(ierr);  ierr = DMCompositeAddDM(pack,dak);CHKERRQ(ierr);  ierr = DMDASetFieldName(dau,0,"u");CHKERRQ(ierr);  ierr = DMDASetFieldName(dak,0,"k");CHKERRQ(ierr);  ierr = DMSetFromOptions(pack);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(pack,&X);CHKERRQ(ierr);  ierr = VecDuplicate(X,&F);CHKERRQ(ierr);  ierr = PetscNew(&user);CHKERRQ(ierr);  user->pack = pack;  ierr = DMCompositeGetGlobalISs(pack,&isg);CHKERRQ(ierr);  ierr = DMCompositeGetLocalVectors(pack,&user->Uloc,&user->Kloc);CHKERRQ(ierr);  ierr = DMCompositeScatter(pack,X,user->Uloc,user->Kloc);CHKERRQ(ierr);  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Coupled problem options","SNES");CHKERRQ(ierr);  {    user->ptype = 0; view_draw = PETSC_FALSE; pass_dm = PETSC_TRUE;    ierr = PetscOptionsInt("-problem_type","0: solve for u only, 1: solve for k only, 2: solve for both",0,user->ptype,&user->ptype,NULL);CHKERRQ(ierr);    ierr = PetscOptionsBool("-view_draw","Draw the final coupled solution regardless of whether only one physics was solved",0,view_draw,&view_draw,NULL);CHKERRQ(ierr);    ierr = PetscOptionsBool("-pass_dm","Pass the packed DM to SNES to use when determining splits and forward into splits",0,pass_dm,&pass_dm,NULL);CHKERRQ(ierr);  }  ierr = PetscOptionsEnd();CHKERRQ(ierr);  ierr = FormInitial_Coupled(user,X);CHKERRQ(ierr);  ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);  switch (user->ptype) {  case 0:    ierr = DMCompositeGetAccess(pack,X,&Xu,0);CHKERRQ(ierr);    ierr = DMCompositeGetAccess(pack,F,&Fu,0);CHKERRQ(ierr);    ierr = DMCreateMatrix(dau,&B);CHKERRQ(ierr);    ierr = SNESSetFunction(snes,Fu,FormFunction_All,user);CHKERRQ(ierr);    ierr = SNESSetJacobian(snes,B,B,FormJacobian_All,user);CHKERRQ(ierr);    ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);    ierr = SNESSetDM(snes,dau);CHKERRQ(ierr);    ierr = SNESSolve(snes,NULL,Xu);CHKERRQ(ierr);    ierr = DMCompositeRestoreAccess(pack,X,&Xu,0);CHKERRQ(ierr);    ierr = DMCompositeRestoreAccess(pack,F,&Fu,0);CHKERRQ(ierr);    break;  case 1:    ierr = DMCompositeGetAccess(pack,X,0,&Xk);CHKERRQ(ierr);    ierr = DMCompositeGetAccess(pack,F,0,&Fk);CHKERRQ(ierr);    ierr = DMCreateMatrix(dak,&B);CHKERRQ(ierr);    ierr = SNESSetFunction(snes,Fk,FormFunction_All,user);CHKERRQ(ierr);    ierr = SNESSetJacobian(snes,B,B,FormJacobian_All,user);CHKERRQ(ierr);    ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);    ierr = SNESSetDM(snes,dak);CHKERRQ(ierr);    ierr = SNESSolve(snes,NULL,Xk);CHKERRQ(ierr);    ierr = DMCompositeRestoreAccess(pack,X,0,&Xk);CHKERRQ(ierr);    ierr = DMCompositeRestoreAccess(pack,F,0,&Fk);CHKERRQ(ierr);    break;  case 2:    ierr = DMCreateMatrix(pack,&B);CHKERRQ(ierr);    /* This example does not correctly allocate off-diagonal blocks. These options allows new nonzeros (slow). */    ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);    ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);    ierr = SNESSetFunction(snes,F,FormFunction_All,user);CHKERRQ(ierr);    ierr = SNESSetJacobian(snes,B,B,FormJacobian_All,user);CHKERRQ(ierr);    ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);    if (!pass_dm) {             /* Manually provide index sets and names for the splits */      KSP ksp;      PC  pc;      ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);      ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);      ierr = PCFieldSplitSetIS(pc,"u",isg[0]);CHKERRQ(ierr);      ierr = PCFieldSplitSetIS(pc,"k",isg[1]);CHKERRQ(ierr);    } else {//.........这里部分代码省略.........

int main(int argc,char **argv){  TS             ts;            /* time integration context */  Vec            X;             /* solution, residual vectors */  Mat            J;             /* Jacobian matrix */  PetscErrorCode ierr;  PetscScalar    *x;  PetscReal      ftime;  PetscInt       i,steps,nits,lits;  PetscBool      view_final;  Ctx            ctx;  PetscInitialize(&argc,&argv,(char *)0,help);  ctx.n = 3;  ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&ctx.n,PETSC_NULL);CHKERRQ(ierr);  if (ctx.n < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"The dimension specified with -n must be at least 2");  view_final = PETSC_FALSE;  ierr = PetscOptionsGetBool(PETSC_NULL,"-view_final",&view_final,PETSC_NULL);CHKERRQ(ierr);  ctx.monitor_short = PETSC_FALSE;  ierr = PetscOptionsGetBool(PETSC_NULL,"-monitor_short",&ctx.monitor_short,PETSC_NULL);CHKERRQ(ierr);  /*     Create Jacobian matrix data structure and state vector  */  ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr);  ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,ctx.n,ctx.n);CHKERRQ(ierr);  ierr = MatSetFromOptions(J);CHKERRQ(ierr);  ierr = MatSetUp(J);CHKERRQ(ierr);  ierr = MatGetVecs(J,&X,PETSC_NULL);CHKERRQ(ierr);  /* Create time integration context */  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);  ierr = TSSetType(ts,TSPSEUDO);CHKERRQ(ierr);  ierr = TSSetIFunction(ts,PETSC_NULL,FormIFunction,&ctx);CHKERRQ(ierr);  ierr = TSSetIJacobian(ts,J,J,FormIJacobian,&ctx);CHKERRQ(ierr);  ierr = TSSetDuration(ts,1000,1e14);CHKERRQ(ierr);  ierr = TSSetInitialTimeStep(ts,0.0,1e-3);CHKERRQ(ierr);  ierr = TSMonitorSet(ts,MonitorObjective,&ctx,PETSC_NULL);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Customize time integrator; set runtime options   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Evaluate initial guess; then solve nonlinear system   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = VecSet(X,0.0);CHKERRQ(ierr);  ierr = VecGetArray(X,&x);CHKERRQ(ierr);#if 1  x[0] = 5.;  x[1] = -5.;  for (i=2; i<ctx.n; i++) x[i] = 5.;#else  x[0] = 1.0;  x[1] = 15.0;  for (i=2; i<ctx.n; i++) x[i] = 10.0;#endif  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);  ierr = TSSolve(ts,X);CHKERRQ(ierr);  ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);  ierr = TSGetTimeStepNumber(ts,&steps);CHKERRQ(ierr);  ierr = TSGetSNESIterations(ts,&nits);CHKERRQ(ierr);  ierr = TSGetKSPIterations(ts,&lits);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_WORLD,"Time integrator took (%D,%D,%D) iterations to reach final time %G/n",steps,nits,lits,ftime);CHKERRQ(ierr);  if (view_final) {    ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);  }  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Free work space.  All PETSc objects should be destroyed when they     are no longer needed.   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = VecDestroy(&X);CHKERRQ(ierr);  ierr = MatDestroy(&J);CHKERRQ(ierr);  ierr = TSDestroy(&ts);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

PetscErrorCode KSPSolve_STCG(KSP ksp){#ifdef PETSC_USE_COMPLEX  SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP, "STCG is not available for complex systems");#else  KSP_STCG       *cg = (KSP_STCG *)ksp->data;  PetscErrorCode ierr;  MatStructure   pflag;  Mat            Qmat, Mmat;  Vec            r, z, p, d;  PC             pc;  PetscReal      norm_r, norm_d, norm_dp1, norm_p, dMp;  PetscReal      alpha, beta, kappa, rz, rzm1;  PetscReal      rr, r2, step;  PetscInt       max_cg_its;  PetscBool      diagonalscale;  PetscFunctionBegin;  /***************************************************************************/  /* Check the arguments and parameters.                                     */  /***************************************************************************/  ierr = PCGetDiagonalScale(ksp->pc, &diagonalscale);CHKERRQ(ierr);  if (diagonalscale) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);  if (cg->radius < 0.0) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");  /***************************************************************************/  /* Get the workspace vectors and initialize variables                      */  /***************************************************************************/  r2 = cg->radius * cg->radius;  r  = ksp->work[0];  z  = ksp->work[1];  p  = ksp->work[2];  d  = ksp->vec_sol;  pc = ksp->pc;  ierr = PCGetOperators(pc, &Qmat, &Mmat, &pflag);CHKERRQ(ierr);  ierr = VecGetSize(d, &max_cg_its);CHKERRQ(ierr);  max_cg_its = PetscMin(max_cg_its, ksp->max_it);  ksp->its = 0;  /***************************************************************************/  /* Initialize objective function and direction.                            */  /***************************************************************************/  cg->o_fcn = 0.0;  ierr = VecSet(d, 0.0);CHKERRQ(ierr);			/* d = 0             */  cg->norm_d = 0.0;  /***************************************************************************/  /* Begin the conjugate gradient method.  Check the right-hand side for     */  /* numerical problems.  The check for not-a-number and infinite values     */  /* need be performed only once.                                            */  /***************************************************************************/  ierr = VecCopy(ksp->vec_rhs, r);CHKERRQ(ierr);	/* r = -grad         */  ierr = VecDot(r, r, &rr);CHKERRQ(ierr);		/* rr = r^T r        */  if (PetscIsInfOrNanScalar(rr)) {    /*************************************************************************/    /* The right-hand side contains not-a-number or an infinite value.       */    /* The gradient step does not work; return a zero value for the step.    */    /*************************************************************************/    ksp->reason = KSP_DIVERGED_NAN;    ierr = PetscInfo1(ksp, "KSPSolve_STCG: bad right-hand side: rr=%g/n", rr);CHKERRQ(ierr);    PetscFunctionReturn(0);  }  /***************************************************************************/  /* Check the preconditioner for numerical problems and for positive        */  /* definiteness.  The check for not-a-number and infinite values need be   */  /* performed only once.                                                    */  /***************************************************************************/  ierr = KSP_PCApply(ksp, r, z);CHKERRQ(ierr);		/* z = inv(M) r      */  ierr = VecDot(r, z, &rz);CHKERRQ(ierr);		/* rz = r^T inv(M) r */  if (PetscIsInfOrNanScalar(rz)) {    /*************************************************************************/    /* The preconditioner contains not-a-number or an infinite value.        */    /* Return the gradient direction intersected with the trust region.      */    /*************************************************************************/    ksp->reason = KSP_DIVERGED_NAN;    ierr = PetscInfo1(ksp, "KSPSolve_STCG: bad preconditioner: rz=%g/n", rz);CHKERRQ(ierr);    if (cg->radius != 0) {      if (r2 >= rr) {        alpha = 1.0;        cg->norm_d = PetscSqrtReal(rr);      }      else {        alpha = PetscSqrtReal(r2 / rr);        cg->norm_d = cg->radius;//.........这里部分代码省略.........

PetscErrorCode MatSetUpMultiply_MPISBAIJ(Mat mat){  Mat_MPISBAIJ   *sbaij = (Mat_MPISBAIJ*)mat->data;  Mat_SeqBAIJ    *B = (Mat_SeqBAIJ*)(sbaij->B->data);    PetscErrorCode ierr;  PetscInt       Nbs = sbaij->Nbs,i,j,*indices,*aj = B->j,ec = 0,*garray,*sgarray;  PetscInt       bs = mat->rmap->bs,*stmp,mbs=sbaij->mbs, vec_size,nt;  IS             from,to;  Vec            gvec;  PetscMPIInt    rank=sbaij->rank,lsize,size=sbaij->size;   PetscInt       *owners=sbaij->rangebs,*sowners,*ec_owner,k;   PetscScalar    *ptr;  PetscFunctionBegin;  ierr = VecScatterDestroy(&sbaij->sMvctx);CHKERRQ(ierr);    /* For the first stab we make an array as long as the number of columns */  /* mark those columns that are in sbaij->B */  ierr = PetscMalloc(Nbs*sizeof(PetscInt),&indices);CHKERRQ(ierr);  ierr = PetscMemzero(indices,Nbs*sizeof(PetscInt));CHKERRQ(ierr);  for (i=0; i<mbs; i++) {    for (j=0; j<B->ilen[i]; j++) {      if (!indices[aj[B->i[i] + j]]) ec++;       indices[aj[B->i[i] + j] ] = 1;    }  }  /* form arrays of columns we need */  ierr = PetscMalloc(ec*sizeof(PetscInt),&garray);CHKERRQ(ierr);  ierr = PetscMalloc2(2*ec,PetscInt,&sgarray,ec,PetscInt,&ec_owner);CHKERRQ(ierr);    ec = 0;  for (j=0; j<size; j++){    for (i=owners[j]; i<owners[j+1]; i++){      if (indices[i]) {        garray[ec]   = i;        ec_owner[ec] = j;        ec++;      }    }  }  /* make indices now point into garray */  for (i=0; i<ec; i++) indices[garray[i]] = i;  /* compact out the extra columns in B */  for (i=0; i<mbs; i++) {    for (j=0; j<B->ilen[i]; j++) aj[B->i[i] + j] = indices[aj[B->i[i] + j]];  }  B->nbs      = ec;  sbaij->B->cmap->n = sbaij->B->cmap->N = ec*mat->rmap->bs;  ierr = PetscLayoutSetUp((sbaij->B->cmap));CHKERRQ(ierr);  ierr = PetscFree(indices);CHKERRQ(ierr);  /* create local vector that is used to scatter into */  ierr = VecCreateSeq(PETSC_COMM_SELF,ec*bs,&sbaij->lvec);CHKERRQ(ierr);  /* create two temporary index sets for building scatter-gather */  ierr = PetscMalloc(2*ec*sizeof(PetscInt),&stmp);CHKERRQ(ierr);  ierr = ISCreateBlock(PETSC_COMM_SELF,bs,ec,garray,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);     for (i=0; i<ec; i++) { stmp[i] = i; }   ierr = ISCreateBlock(PETSC_COMM_SELF,bs,ec,stmp,PETSC_COPY_VALUES,&to);CHKERRQ(ierr);  /* generate the scatter context      -- Mvctx and lvec are not used by MatMult_MPISBAIJ(), but usefule for some applications */  ierr = VecCreateMPIWithArray(((PetscObject)mat)->comm,1,mat->cmap->n,mat->cmap->N,PETSC_NULL,&gvec);CHKERRQ(ierr);  ierr = VecScatterCreate(gvec,from,sbaij->lvec,to,&sbaij->Mvctx);CHKERRQ(ierr);   ierr = VecDestroy(&gvec);CHKERRQ(ierr);  sbaij->garray = garray;  ierr = PetscLogObjectParent(mat,sbaij->Mvctx);CHKERRQ(ierr);  ierr = PetscLogObjectParent(mat,sbaij->lvec);CHKERRQ(ierr);  ierr = PetscLogObjectParent(mat,from);CHKERRQ(ierr);  ierr = PetscLogObjectParent(mat,to);CHKERRQ(ierr);  ierr = ISDestroy(&from);CHKERRQ(ierr);  ierr = ISDestroy(&to);CHKERRQ(ierr);  /* create parallel vector that is used by SBAIJ matrix to scatter from/into */  lsize = (mbs + ec)*bs;  ierr = VecCreateMPI(((PetscObject)mat)->comm,lsize,PETSC_DETERMINE,&sbaij->slvec0);CHKERRQ(ierr);  ierr = VecDuplicate(sbaij->slvec0,&sbaij->slvec1);CHKERRQ(ierr);  ierr = VecGetSize(sbaij->slvec0,&vec_size);CHKERRQ(ierr);  sowners = sbaij->slvec0->map->range;   /* x index in the IS sfrom */  for (i=0; i<ec; i++) {     j = ec_owner[i];    sgarray[i]  = garray[i] + (sowners[j]/bs - owners[j]);    }   /* b index in the IS sfrom */  k = sowners[rank]/bs + mbs;  for (i=ec,j=0; i< 2*ec; i++,j++) sgarray[i] = k + j;  ierr = ISCreateBlock(PETSC_COMM_SELF,bs,2*ec,sgarray,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);      /* x index in the IS sto */  k = sowners[rank]/bs + mbs;  for (i=0; i<ec; i++) stmp[i] = (k + i);    /* b index in the IS sto *///.........这里部分代码省略.........

int main(int argc,char **argv){  PetscErrorCode ierr;  UserCtx        user;  DM             red,da;  SNES           snes;  DM             packer;  PetscBool      use_monitor = PETSC_FALSE;  PetscInitialize(&argc,&argv,NULL,help);  ierr = PetscOptionsSetFromOptions();CHKERRQ(ierr);  /* Hardwire several options; can be changed at command line */  ierr = PetscOptionsInsertString(common_options);CHKERRQ(ierr);  ierr = PetscOptionsInsertString(matrix_free_options);CHKERRQ(ierr);  ierr = PetscOptionsInsert(&argc,&argv,NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetBool(NULL,"-use_monitor",&use_monitor,PETSC_IGNORE);CHKERRQ(ierr);  /* Create a global vector that includes a single redundant array and two da arrays */  ierr = DMCompositeCreate(PETSC_COMM_WORLD,&packer);CHKERRQ(ierr);  ierr = DMRedundantCreate(PETSC_COMM_WORLD,0,1,&red);CHKERRQ(ierr);  ierr = DMSetOptionsPrefix(red,"red_");CHKERRQ(ierr);  ierr = DMCompositeAddDM(packer,red);CHKERRQ(ierr);  ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-5,2,1,NULL,&da);CHKERRQ(ierr);  ierr = DMSetOptionsPrefix(red,"da_");CHKERRQ(ierr);  ierr = DMCompositeAddDM(packer,(DM)da);CHKERRQ(ierr);  ierr = DMSetApplicationContext(packer,&user);CHKERRQ(ierr);  packer->ops->creatematrix = DMCreateMatrix_MF;  /* create nonlinear multi-level solver */  ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);  ierr = SNESSetDM(snes,packer);CHKERRQ(ierr);  ierr = SNESSetFunction(snes,NULL,ComputeFunction,NULL);CHKERRQ(ierr);  ierr = SNESSetJacobian(snes,NULL, NULL,ComputeJacobian_MF,NULL);CHKERRQ(ierr);  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);  if (use_monitor) {    /* create graphics windows */    ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u_lambda - state variables and Lagrange multipliers",-1,-1,-1,-1,&user.u_lambda_viewer);CHKERRQ(ierr);    ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu_lambda - derivate w.r.t. state variables and Lagrange multipliers",-1,-1,-1,-1,&user.fu_lambda_viewer);CHKERRQ(ierr);    ierr = SNESMonitorSet(snes,Monitor,0,0);CHKERRQ(ierr);  }  ierr = SNESSolve(snes,NULL,NULL);CHKERRQ(ierr);  ierr = SNESDestroy(&snes);CHKERRQ(ierr);  ierr = DMDestroy(&red);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = DMDestroy(&packer);CHKERRQ(ierr);  if (use_monitor) {    ierr = PetscViewerDestroy(&user.u_lambda_viewer);CHKERRQ(ierr);    ierr = PetscViewerDestroy(&user.fu_lambda_viewer);CHKERRQ(ierr);  }  PetscFinalize();  return 0;}

PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy){  UserCtx        *user;  PetscErrorCode ierr;  PetscInt       m,N;  PetscScalar    *w,*dw;  Vec            u_lambda,U,F,Uexact;  DM             packer;  PetscReal      norm;  DM             da;  PetscFunctionBeginUser;  ierr = SNESGetDM(snes,&packer);CHKERRQ(ierr);  ierr = DMGetApplicationContext(packer,&user);CHKERRQ(ierr);  ierr = SNESGetSolution(snes,&U);CHKERRQ(ierr);  ierr = DMCompositeGetAccess(packer,U,&w,&u_lambda);CHKERRQ(ierr);  ierr = VecView(u_lambda,user->u_lambda_viewer);  ierr = DMCompositeRestoreAccess(packer,U,&w,&u_lambda);CHKERRQ(ierr);  ierr = SNESGetFunction(snes,&F,0,0);CHKERRQ(ierr);  ierr = DMCompositeGetAccess(packer,F,&w,&u_lambda);CHKERRQ(ierr);  /* ierr = VecView(u_lambda,user->fu_lambda_viewer); */  ierr = DMCompositeRestoreAccess(packer,U,&w,&u_lambda);CHKERRQ(ierr);  ierr = DMCompositeGetEntries(packer,&m,&da);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  ierr = VecDuplicate(U,&Uexact);CHKERRQ(ierr);  ierr = ExactSolution(packer,Uexact);CHKERRQ(ierr);  ierr = VecAXPY(Uexact,-1.0,U);CHKERRQ(ierr);  ierr = DMCompositeGetAccess(packer,Uexact,&dw,&u_lambda);CHKERRQ(ierr);  ierr = VecStrideNorm(u_lambda,0,NORM_2,&norm);CHKERRQ(ierr);  norm = norm/PetscSqrtReal((PetscReal)N-1.);  if (dw) ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Error at x = 0 %g/n",(double)norm,(double)PetscRealPart(dw[0]));CHKERRQ(ierr);  ierr = VecView(u_lambda,user->fu_lambda_viewer);  ierr = DMCompositeRestoreAccess(packer,Uexact,&dw,&u_lambda);CHKERRQ(ierr);  ierr = VecDestroy(&Uexact);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*      Evaluates FU = Gradiant(L(w,u,lambda))     This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and   DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeGetAccess()).*/PetscErrorCode ComputeFunction(SNES snes,Vec U,Vec FU,void *ctx){  PetscErrorCode ierr;  PetscInt       xs,xm,i,N;  ULambda        *u_lambda,*fu_lambda;  PetscScalar    d,h,*w,*fw;  Vec            vw,vfw,vu_lambda,vfu_lambda;  DM             packer,red,da;  PetscFunctionBeginUser;  ierr = VecGetDM(U, &packer);CHKERRQ(ierr);  ierr = DMCompositeGetEntries(packer,&red,&da);CHKERRQ(ierr);  ierr = DMCompositeGetLocalVectors(packer,&vw,&vu_lambda);CHKERRQ(ierr);  ierr = DMCompositeScatter(packer,U,vw,vu_lambda);CHKERRQ(ierr);  ierr = DMCompositeGetAccess(packer,FU,&vfw,&vfu_lambda);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  ierr = VecGetArray(vw,&w);CHKERRQ(ierr);  ierr = VecGetArray(vfw,&fw);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,vu_lambda,&u_lambda);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,vfu_lambda,&fu_lambda);CHKERRQ(ierr);  d    = N-1.0;  h    = 1.0/d;  /* derivative of L() w.r.t. w */  if (xs == 0) { /* only first processor computes this */    fw[0] = -2.0*d*u_lambda[0].lambda;  }  /* derivative of L() w.r.t. u */  for (i=xs; i<xs+xm; i++) {    if      (i == 0)   fu_lambda[0].lambda   =    h*u_lambda[0].u   + 2.*d*u_lambda[0].lambda   - d*u_lambda[1].lambda;    else if (i == 1)   fu_lambda[1].lambda   = 2.*h*u_lambda[1].u   + 2.*d*u_lambda[1].lambda   - d*u_lambda[2].lambda;    else if (i == N-1) fu_lambda[N-1].lambda =    h*u_lambda[N-1].u + 2.*d*u_lambda[N-1].lambda - d*u_lambda[N-2].lambda;    else if (i == N-2) fu_lambda[N-2].lambda = 2.*h*u_lambda[N-2].u + 2.*d*u_lambda[N-2].lambda - d*u_lambda[N-3].lambda;    else               fu_lambda[i].lambda   = 2.*h*u_lambda[i].u   - d*(u_lambda[i+1].lambda - 2.0*u_lambda[i].lambda + u_lambda[i-1].lambda);  }  /* derivative of L() w.r.t. lambda */  for (i=xs; i<xs+xm; i++) {    if      (i == 0)   fu_lambda[0].u   = 2.0*d*(u_lambda[0].u - w[0]);    else if (i == N-1) fu_lambda[N-1].u = 2.0*d*u_lambda[N-1].u;    else               fu_lambda[i].u   = -(d*(u_lambda[i+1].u - 2.0*u_lambda[i].u + u_lambda[i-1].u) - 2.0*h);  }  ierr = VecRestoreArray(vw,&w);CHKERRQ(ierr);  ierr = VecRestoreArray(vfw,&fw);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,vu_lambda,&u_lambda);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,vfu_lambda,&fu_lambda);CHKERRQ(ierr);  ierr = DMCompositeRestoreLocalVectors(packer,&vw,&vu_lambda);CHKERRQ(ierr);  ierr = DMCompositeRestoreAccess(packer,FU,&vfw,&vfu_lambda);CHKERRQ(ierr);  ierr = PetscLogFlops(13.0*N);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*@  TaoComputeObjectiveAndGradient - Computes the objective function value at a given point  Collective on Tao  Input Parameters:+ tao - the Tao context- X - input vector  Output Parameter:+ f - Objective value at X- g - Gradient vector at X  Notes: TaoComputeObjectiveAndGradient() is typically used within minimization implementations,  so most users would not generally call this routine themselves.  Level: advanced.seealso: TaoComputeGradient(), TaoComputeObjectiveAndGradient(), TaoSetObjectiveRoutine()@*/PetscErrorCode TaoComputeObjectiveAndGradient(Tao tao, Vec X, PetscReal *f, Vec G){    PetscErrorCode ierr;    PetscFunctionBegin;    PetscValidHeaderSpecific(tao,TAO_CLASSID,1);    PetscValidHeaderSpecific(X,VEC_CLASSID,2);    PetscValidHeaderSpecific(G,VEC_CLASSID,4);    PetscCheckSameComm(tao,1,X,2);    PetscCheckSameComm(tao,1,G,4);    if (tao->ops->computeobjectiveandgradient) {        ierr = PetscLogEventBegin(Tao_ObjGradientEval,tao,X,G,NULL);        CHKERRQ(ierr);        PetscStackPush("Tao user objective/gradient evaluation routine");        ierr = (*tao->ops->computeobjectiveandgradient)(tao,X,f,G,tao->user_objgradP);        CHKERRQ(ierr);        PetscStackPop;        if (tao->ops->computegradient == TaoDefaultComputeGradient) {            /* Overwrite gradient with finite difference gradient */            ierr = TaoDefaultComputeGradient(tao,X,G,tao->user_objgradP);            CHKERRQ(ierr);        }        ierr = PetscLogEventEnd(Tao_ObjGradientEval,tao,X,G,NULL);        CHKERRQ(ierr);        tao->nfuncgrads++;    } else if (tao->ops->computeobjective && tao->ops->computegradient) {        ierr = PetscLogEventBegin(Tao_ObjectiveEval,tao,X,NULL,NULL);        CHKERRQ(ierr);        PetscStackPush("Tao user objective evaluation routine");        ierr = (*tao->ops->computeobjective)(tao,X,f,tao->user_objP);        CHKERRQ(ierr);        PetscStackPop;        ierr = PetscLogEventEnd(Tao_ObjectiveEval,tao,X,NULL,NULL);        CHKERRQ(ierr);        tao->nfuncs++;        ierr = PetscLogEventBegin(Tao_GradientEval,tao,X,G,NULL);        CHKERRQ(ierr);        PetscStackPush("Tao user gradient evaluation routine");        ierr = (*tao->ops->computegradient)(tao,X,G,tao->user_gradP);        CHKERRQ(ierr);        PetscStackPop;        ierr = PetscLogEventEnd(Tao_GradientEval,tao,X,G,NULL);        CHKERRQ(ierr);        tao->ngrads++;    } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"TaoSetObjectiveRoutine() or TaoSetGradientRoutine() not set");    ierr = PetscInfo1(tao,"TAO Function evaluation: %14.12e/n",(double)(*f));    CHKERRQ(ierr);    PetscFunctionReturn(0);}

PetscErrorCode MatSetUpMultiply_MPIDense(Mat mat){  Mat_MPIDense   *mdn = (Mat_MPIDense*)mat->data;  PetscErrorCode ierr;  IS             from,to;  Vec            gvec;  PetscFunctionBegin;  /* Create local vector that is used to scatter into */  ierr = VecCreateSeq(PETSC_COMM_SELF,mat->cmap->N,&mdn->lvec);CHKERRQ(ierr);  /* Create temporary index set for building scatter gather */  ierr = ISCreateStride(PetscObjectComm((PetscObject)mat),mat->cmap->N,0,1,&from);CHKERRQ(ierr);  ierr = ISCreateStride(PETSC_COMM_SELF,mat->cmap->N,0,1,&to);CHKERRQ(ierr);  /* Create temporary global vector to generate scatter context */  /* n    = mdn->cowners[mdn->rank+1] - mdn->cowners[mdn->rank]; */  ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat),1,mdn->nvec,mat->cmap->N,NULL,&gvec);CHKERRQ(ierr);  /* Generate the scatter context */  ierr = VecScatterCreate(gvec,from,mdn->lvec,to,&mdn->Mvctx);CHKERRQ(ierr);  ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)mdn->Mvctx);CHKERRQ(ierr);  ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)mdn->lvec);CHKERRQ(ierr);  ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)from);CHKERRQ(ierr);  ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)to);CHKERRQ(ierr);  ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)gvec);CHKERRQ(ierr);  ierr = ISDestroy(&to);CHKERRQ(ierr);  ierr = ISDestroy(&from);CHKERRQ(ierr);  ierr = VecDestroy(&gvec);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode MatSetUpMultiply_MPISBAIJ_2comm(Mat mat){  Mat_MPISBAIJ       *baij = (Mat_MPISBAIJ*)mat->data;  Mat_SeqBAIJ        *B = (Mat_SeqBAIJ*)(baij->B->data);    PetscErrorCode     ierr;  PetscInt           i,j,*aj = B->j,ec = 0,*garray;  PetscInt           bs = mat->rmap->bs,*stmp;  IS                 from,to;  Vec                gvec;#if defined (PETSC_USE_CTABLE)  PetscTable         gid1_lid1;  PetscTablePosition tpos;  PetscInt           gid,lid;#else  PetscInt           Nbs = baij->Nbs,*indices;#endif    PetscFunctionBegin;#if defined (PETSC_USE_CTABLE)  /* use a table - Mark Adams */  PetscTableCreate(B->mbs,baij->Nbs+1,&gid1_lid1);   for (i=0; i<B->mbs; i++) {    for (j=0; j<B->ilen[i]; j++) {      PetscInt data,gid1 = aj[B->i[i]+j] + 1;      ierr = PetscTableFind(gid1_lid1,gid1,&data);CHKERRQ(ierr);      if (!data) {        /* one based table */         ierr = PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES);CHKERRQ(ierr);       }    }  }   /* form array of columns we need */  ierr = PetscMalloc(ec*sizeof(PetscInt),&garray);CHKERRQ(ierr);  ierr = PetscTableGetHeadPosition(gid1_lid1,&tpos);CHKERRQ(ierr);   while (tpos) {      ierr = PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid);CHKERRQ(ierr);     gid--; lid--;    garray[lid] = gid;   }  ierr = PetscSortInt(ec,garray);CHKERRQ(ierr);  ierr = PetscTableRemoveAll(gid1_lid1);CHKERRQ(ierr);  for (i=0; i<ec; i++) {    ierr = PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES);CHKERRQ(ierr);   }  /* compact out the extra columns in B */  for (i=0; i<B->mbs; i++) {    for (j=0; j<B->ilen[i]; j++) {      PetscInt gid1 = aj[B->i[i] + j] + 1;      ierr = PetscTableFind(gid1_lid1,gid1,&lid);CHKERRQ(ierr);      lid --;      aj[B->i[i]+j] = lid;    }  }  B->nbs     = ec;  baij->B->cmap->n = baij->B->cmap->N = ec*mat->rmap->bs;  ierr = PetscLayoutSetUp((baij->B->cmap));CHKERRQ(ierr);  ierr = PetscTableDestroy(&gid1_lid1);CHKERRQ(ierr);#else  /* For the first stab we make an array as long as the number of columns */  /* mark those columns that are in baij->B */  ierr = PetscMalloc(Nbs*sizeof(PetscInt),&indices);CHKERRQ(ierr);  ierr = PetscMemzero(indices,Nbs*sizeof(PetscInt));CHKERRQ(ierr);  for (i=0; i<B->mbs; i++) {    for (j=0; j<B->ilen[i]; j++) {      if (!indices[aj[B->i[i] + j]]) ec++;       indices[aj[B->i[i] + j] ] = 1;    }  }  /* form array of columns we need */  ierr = PetscMalloc(ec*sizeof(PetscInt),&garray);CHKERRQ(ierr);  ec = 0;  for (i=0; i<Nbs; i++) {    if (indices[i]) {      garray[ec++] = i;    }  }  /* make indices now point into garray */  for (i=0; i<ec; i++) {    indices[garray[i]] = i;  }  /* compact out the extra columns in B */  for (i=0; i<B->mbs; i++) {    for (j=0; j<B->ilen[i]; j++) {      aj[B->i[i] + j] = indices[aj[B->i[i] + j]];    }  }  B->nbs       = ec;  baij->B->cmap->n   = ec*mat->rmap->bs;  ierr = PetscFree(indices);CHKERRQ(ierr);#endif      /* create local vector that is used to scatter into */  ierr = VecCreateSeq(PETSC_COMM_SELF,ec*bs,&baij->lvec);CHKERRQ(ierr);  /* create two temporary index sets for building scatter-gather */  ierr = ISCreateBlock(PETSC_COMM_SELF,bs,ec,garray,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);   //.........这里部分代码省略.........

static PetscErrorCode FormJacobian_All(SNES snes,Vec X,Mat J,Mat B,void *ctx){  User           user = (User)ctx;  DM             dau,dak;  DMDALocalInfo  infou,infok;  PetscScalar    *u,*k;  PetscErrorCode ierr;  Vec            Uloc,Kloc;  PetscFunctionBeginUser;  ierr = DMCompositeGetEntries(user->pack,&dau,&dak);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dau,&infou);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dak,&infok);CHKERRQ(ierr);  ierr = DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc);CHKERRQ(ierr);  switch (user->ptype) {  case 0:    ierr = DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd  (dau,X,INSERT_VALUES,Uloc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,user->Kloc,&k);CHKERRQ(ierr);    ierr = FormJacobianLocal_U(user,&infou,u,k,B);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,user->Kloc,&k);CHKERRQ(ierr);    break;  case 1:    ierr = DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd  (dak,X,INSERT_VALUES,Kloc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,user->Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,Kloc,&k);CHKERRQ(ierr);    ierr = FormJacobianLocal_K(user,&infok,u,k,B);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,user->Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,Kloc,&k);CHKERRQ(ierr);    break;  case 2: {    Mat Buu,Buk,Bku,Bkk;    IS  *is;    ierr = DMCompositeScatter(user->pack,X,Uloc,Kloc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,Kloc,&k);CHKERRQ(ierr);    ierr = DMCompositeGetLocalISs(user->pack,&is);CHKERRQ(ierr);    ierr = MatGetLocalSubMatrix(B,is[0],is[0],&Buu);CHKERRQ(ierr);    ierr = MatGetLocalSubMatrix(B,is[0],is[1],&Buk);CHKERRQ(ierr);    ierr = MatGetLocalSubMatrix(B,is[1],is[0],&Bku);CHKERRQ(ierr);    ierr = MatGetLocalSubMatrix(B,is[1],is[1],&Bkk);CHKERRQ(ierr);    ierr = FormJacobianLocal_U(user,&infou,u,k,Buu);CHKERRQ(ierr);    ierr = FormJacobianLocal_UK(user,&infou,&infok,u,k,Buk);CHKERRQ(ierr);    ierr = FormJacobianLocal_KU(user,&infou,&infok,u,k,Bku);CHKERRQ(ierr);    ierr = FormJacobianLocal_K(user,&infok,u,k,Bkk);CHKERRQ(ierr);    ierr = MatRestoreLocalSubMatrix(B,is[0],is[0],&Buu);CHKERRQ(ierr);    ierr = MatRestoreLocalSubMatrix(B,is[0],is[1],&Buk);CHKERRQ(ierr);    ierr = MatRestoreLocalSubMatrix(B,is[1],is[0],&Bku);CHKERRQ(ierr);    ierr = MatRestoreLocalSubMatrix(B,is[1],is[1],&Bkk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,Kloc,&k);CHKERRQ(ierr);    ierr = ISDestroy(&is[0]);CHKERRQ(ierr);    ierr = ISDestroy(&is[1]);CHKERRQ(ierr);    ierr = PetscFree(is);CHKERRQ(ierr);  } break;  }  ierr = DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc);CHKERRQ(ierr);  ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  ierr = MatAssemblyEnd  (B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  if (J != B) {    ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);    ierr = MatAssemblyEnd  (J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

PetscErrorCode MatDisAssemble_MPISBAIJ(Mat A){  Mat_MPISBAIJ   *baij = (Mat_MPISBAIJ*)A->data;  Mat            B = baij->B,Bnew;  Mat_SeqBAIJ    *Bbaij = (Mat_SeqBAIJ*)B->data;  PetscErrorCode ierr;  PetscInt       i,j,mbs=Bbaij->mbs,n = A->cmap->N,col,*garray=baij->garray;  PetscInt       k,bs=A->rmap->bs,bs2=baij->bs2,*rvals,*nz,ec,m=A->rmap->n;  MatScalar      *a = Bbaij->a;  PetscScalar    *atmp;#if defined(PETSC_USE_REAL_MAT_SINGLE)  PetscInt       l;#endif  PetscFunctionBegin;#if defined(PETSC_USE_REAL_MAT_SINGLE)  ierr = PetscMalloc(A->rmap->bs*sizeof(PetscScalar),&atmp);#endif  /* free stuff related to matrix-vec multiply */  ierr = VecGetSize(baij->lvec,&ec);CHKERRQ(ierr); /* needed for PetscLogObjectMemory below */  ierr = VecDestroy(&baij->lvec);CHKERRQ(ierr);   ierr = VecScatterDestroy(&baij->Mvctx);CHKERRQ(ierr);   ierr = VecDestroy(&baij->slvec0);CHKERRQ(ierr);  ierr = VecDestroy(&baij->slvec0b);CHKERRQ(ierr);   ierr = VecDestroy(&baij->slvec1);CHKERRQ(ierr);  ierr = VecDestroy(&baij->slvec1a);CHKERRQ(ierr);  ierr = VecDestroy(&baij->slvec1b);CHKERRQ(ierr);   if (baij->colmap) {#if defined (PETSC_USE_CTABLE)    ierr = PetscTableDestroy(&baij->colmap);CHKERRQ(ierr);#else    ierr = PetscFree(baij->colmap);CHKERRQ(ierr);    ierr = PetscLogObjectMemory(A,-Bbaij->nbs*sizeof(PetscInt));CHKERRQ(ierr);#endif  }  /* make sure that B is assembled so we can access its values */  ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  /* invent new B and copy stuff over */  ierr = PetscMalloc(mbs*sizeof(PetscInt),&nz);CHKERRQ(ierr);  for (i=0; i<mbs; i++) {    nz[i] = Bbaij->i[i+1]-Bbaij->i[i];  }  ierr = MatCreate(PETSC_COMM_SELF,&Bnew);CHKERRQ(ierr);  ierr = MatSetSizes(Bnew,m,n,m,n);CHKERRQ(ierr);  ierr = MatSetType(Bnew,((PetscObject)B)->type_name);CHKERRQ(ierr);  ierr = MatSeqBAIJSetPreallocation(Bnew,B->rmap->bs,0,nz);CHKERRQ(ierr);  ((Mat_SeqSBAIJ*)Bnew->data)->nonew = Bbaij->nonew; /* Inherit insertion error options. */  ierr = PetscFree(nz);CHKERRQ(ierr);    ierr = PetscMalloc(bs*sizeof(PetscInt),&rvals);CHKERRQ(ierr);  for (i=0; i<mbs; i++) {    rvals[0] = bs*i;    for (j=1; j<bs; j++) { rvals[j] = rvals[j-1] + 1; }    for (j=Bbaij->i[i]; j<Bbaij->i[i+1]; j++) {      col = garray[Bbaij->j[j]]*bs;      for (k=0; k<bs; k++) {#if defined(PETSC_USE_REAL_MAT_SINGLE)        for (l=0; l<bs; l++) atmp[l] = a[j*bs2+l];#else        atmp = a+j*bs2 + k*bs;#endif        ierr = MatSetValues_SeqSBAIJ(Bnew,bs,rvals,1,&col,atmp,B->insertmode);CHKERRQ(ierr);        col++;      }    }  }#if defined(PETSC_USE_REAL_MAT_SINGLE)  ierr = PetscFree(atmp);CHKERRQ(ierr);#endif  ierr = PetscFree(baij->garray);CHKERRQ(ierr);  baij->garray = 0;  ierr = PetscFree(rvals);CHKERRQ(ierr);  ierr = PetscLogObjectMemory(A,-ec*sizeof(PetscInt));CHKERRQ(ierr);  ierr = MatDestroy(&B);CHKERRQ(ierr);  ierr = PetscLogObjectParent(A,Bnew);CHKERRQ(ierr);  baij->B = Bnew;  A->was_assembled = PETSC_FALSE;  PetscFunctionReturn(0);}

static PetscErrorCode FormFunction_All(SNES snes,Vec X,Vec F,void *ctx){  User           user = (User)ctx;  DM             dau,dak;  DMDALocalInfo  infou,infok;  PetscScalar    *u,*k;  PetscScalar    *fu,*fk;  PetscErrorCode ierr;  Vec            Uloc,Kloc,Fu,Fk;  PetscFunctionBeginUser;  ierr = DMCompositeGetEntries(user->pack,&dau,&dak);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dau,&infou);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dak,&infok);CHKERRQ(ierr);  ierr = DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc);CHKERRQ(ierr);  switch (user->ptype) {  case 0:    ierr = DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd  (dau,X,INSERT_VALUES,Uloc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,user->Kloc,&k);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,F,&fu);CHKERRQ(ierr);    ierr = FormFunctionLocal_U(user,&infou,u,k,fu);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,F,&fu);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,user->Kloc,&k);CHKERRQ(ierr);    break;  case 1:    ierr = DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd  (dak,X,INSERT_VALUES,Kloc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,user->Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,Kloc,&k);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,F,&fk);CHKERRQ(ierr);    ierr = FormFunctionLocal_K(user,&infok,u,k,fk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,F,&fk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,user->Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,Kloc,&k);CHKERRQ(ierr);    break;  case 2:    ierr = DMCompositeScatter(user->pack,X,Uloc,Kloc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,Kloc,&k);CHKERRQ(ierr);    ierr = DMCompositeGetAccess(user->pack,F,&Fu,&Fk);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,Fu,&fu);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,Fk,&fk);CHKERRQ(ierr);    ierr = FormFunctionLocal_U(user,&infou,u,k,fu);CHKERRQ(ierr);    ierr = FormFunctionLocal_K(user,&infok,u,k,fk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,Fu,&fu);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,Fk,&fk);CHKERRQ(ierr);    ierr = DMCompositeRestoreAccess(user->pack,F,&Fu,&Fk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,Kloc,&k);CHKERRQ(ierr);    break;  }  ierr = DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr){  User           user = (User)ptr;  DM             da;  DMDALocalInfo  info;  PetscInt       i;  Field          *x,*xdot,*f;  PetscReal      hx;  Vec            Xloc;  PetscErrorCode ierr;  PetscFunctionBeginUser;  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);  hx   = 1.0/(PetscReal)(info.mx-1);  /*     Scatter ghost points to local vector,using the 2-step process        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().     By placing code between these two statements, computations can be     done while messages are in transition.  */  ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  /* Get pointers to vector data */  ierr = DMDAVecGetArrayRead(da,Xloc,&x);CHKERRQ(ierr);  ierr = DMDAVecGetArrayRead(da,Xdot,&xdot);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);  /* Compute function over the locally owned part of the grid */  for (i=info.xs; i<info.xs+info.xm; i++) {    if (i == 0) {      f[i].u = hx * (x[i].u - user->uleft);      f[i].v = hx * (x[i].v - user->vleft);    } else if (i == info.mx-1) {      f[i].u = hx * (x[i].u - user->uright);      f[i].v = hx * (x[i].v - user->vright);    } else {      f[i].u = hx * xdot[i].u - user->alpha * (x[i-1].u - 2.*x[i].u + x[i+1].u) / hx;      f[i].v = hx * xdot[i].v - user->alpha * (x[i-1].v - 2.*x[i].v + x[i+1].v) / hx;    }  }  /* Restore vectors */  ierr = DMDAVecRestoreArrayRead(da,Xloc,&x);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayRead(da,Xdot,&xdot);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode MonitorObjective(TS ts,PetscInt step,PetscReal t,Vec X,void *ictx){  Ctx               *ctx = (Ctx*)ictx;  PetscErrorCode    ierr;  const PetscScalar *x;  PetscScalar       f;  PetscReal         dt,gnorm;  PetscInt          i,snesit,linit;  SNES              snes;  Vec               Xdot,F;  PetscFunctionBeginUser;  /* Compute objective functional */  ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);  f = 0;  for (i=0; i<ctx->n-1; i++) {    f += PetscSqr(1. - x[i]) + 100. * PetscSqr(x[i+1] - PetscSqr(x[i]));  }  ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);  /* Compute norm of gradient */  ierr = VecDuplicate(X,&Xdot);CHKERRQ(ierr);  ierr = VecDuplicate(X,&F);CHKERRQ(ierr);  ierr = VecZeroEntries(Xdot);CHKERRQ(ierr);  ierr = FormIFunction(ts,t,X,Xdot,F,ictx);CHKERRQ(ierr);  ierr = VecNorm(F,NORM_2,&gnorm);CHKERRQ(ierr);  ierr = VecDestroy(&Xdot);CHKERRQ(ierr);  ierr = VecDestroy(&F);CHKERRQ(ierr);  ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);  ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);  ierr = SNESGetIterationNumber(snes,&snesit);CHKERRQ(ierr);  ierr = SNESGetLinearSolveIterations(snes,&linit);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_WORLD,                     (ctx->monitor_short                      ? "%3D t=%10.1e  dt=%10.1e  f=%10.1e  df=%10.1e  it=(%2D,%3D)/n"                      : "%3D t=%10.4e  dt=%10.4e  f=%10.4e  df=%10.4e  it=(%2D,%3D)/n"),                     step,(double)t,(double)dt,(double)PetscRealPart(f),(double)gnorm,snesit,linit);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ptr){  User           user = (User)ptr;  PetscErrorCode ierr;  DMDALocalInfo  info;  PetscInt       i;  PetscReal      hx;  DM             da;  Field          *x,*xdot;  PetscFunctionBeginUser;  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);  hx   = 1.0/(PetscReal)(info.mx-1);  /* Get pointers to vector data */  ierr = DMDAVecGetArrayRead(da,X,&x);CHKERRQ(ierr);  ierr = DMDAVecGetArrayRead(da,Xdot,&xdot);CHKERRQ(ierr);  /* Compute function over the locally owned part of the grid */  for (i=info.xs; i<info.xs+info.xm; i++) {    if (i == 0 || i == info.mx-1) {      const PetscInt    row        = i,col = i;      const PetscScalar vals[2][2] = {{hx,0},{0,hx}};      ierr = MatSetValuesBlocked(Jpre,1,&row,1,&col,&vals[0][0],INSERT_VALUES);CHKERRQ(ierr);    } else {      const PetscInt    row           = i,col[] = {i-1,i,i+1};      const PetscScalar dxxL          = -user->alpha/hx,dxx0 = 2.*user->alpha/hx,dxxR = -user->alpha/hx;      const PetscScalar vals[2][3][2] = {{{dxxL,0},{a *hx+dxx0,0},{dxxR,0}},                                         {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}};      ierr = MatSetValuesBlocked(Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES);CHKERRQ(ierr);    }  }  /* Restore vectors */  ierr = DMDAVecRestoreArrayRead(da,X,&x);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayRead(da,Xdot,&xdot);CHKERRQ(ierr);  ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  if (J != Jpre) {    ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);    ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

PetscErrorCode Brusselator(int argc,char **argv,PetscInt cycle){  TS                ts;         /* nonlinear solver */  Vec               X;          /* solution, residual vectors */  Mat               J;          /* Jacobian matrix */  PetscInt          steps,maxsteps,mx;  PetscErrorCode    ierr;  DM                da;  PetscReal         ftime,hx,dt,xmax,xmin;  struct _User      user;       /* user-defined work context */  TSConvergedReason reason;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create distributed array (DMDA) to manage parallel grid and vectors  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-11,2,2,NULL,&da);CHKERRQ(ierr);  /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Extract global vectors from DMDA;   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMCreateGlobalVector(da,&X);CHKERRQ(ierr);  /* Initialize user application context */  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Advection-reaction options","");  {    user.A      = 1;    user.B      = 3;    user.alpha  = 0.1;    user.uleft  = 1;    user.uright = 1;    user.vleft  = 3;    user.vright = 3;    ierr        = PetscOptionsReal("-A","Reaction rate","",user.A,&user.A,NULL);CHKERRQ(ierr);    ierr        = PetscOptionsReal("-B","Reaction rate","",user.B,&user.B,NULL);CHKERRQ(ierr);    ierr        = PetscOptionsReal("-alpha","Diffusion coefficient","",user.alpha,&user.alpha,NULL);CHKERRQ(ierr);    ierr        = PetscOptionsReal("-uleft","Dirichlet boundary condition","",user.uleft,&user.uleft,NULL);CHKERRQ(ierr);    ierr        = PetscOptionsReal("-uright","Dirichlet boundary condition","",user.uright,&user.uright,NULL);CHKERRQ(ierr);    ierr        = PetscOptionsReal("-vleft","Dirichlet boundary condition","",user.vleft,&user.vleft,NULL);CHKERRQ(ierr);    ierr        = PetscOptionsReal("-vright","Dirichlet boundary condition","",user.vright,&user.vright,NULL);CHKERRQ(ierr);  }  ierr = PetscOptionsEnd();CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create timestepping solver context     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);  ierr = TSSetDM(ts,da);CHKERRQ(ierr);  ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr);  ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);CHKERRQ(ierr);  ierr = TSSetIFunction(ts,NULL,FormIFunction,&user);CHKERRQ(ierr);  ierr = DMSetMatType(da,MATAIJ);CHKERRQ(ierr);  ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr);  ierr = TSSetIJacobian(ts,J,J,FormIJacobian,&user);CHKERRQ(ierr);  ftime    = 1.0;  maxsteps = 10000;  ierr     = TSSetDuration(ts,maxsteps,ftime);CHKERRQ(ierr);  ierr     = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Set initial conditions   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = FormInitialSolution(ts,X,&user);CHKERRQ(ierr);  ierr = TSSetSolution(ts,X);CHKERRQ(ierr);  ierr = VecGetSize(X,&mx);CHKERRQ(ierr);  hx = 1.0/(PetscReal)(mx/2-1);  dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */  dt *= PetscPowRealInt(0.2,cycle);     /* Shrink the time step in convergence study. */  ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);  ierr = TSSetTolerances(ts,1e-3*PetscPowRealInt(0.5,cycle),NULL,1e-3*PetscPowRealInt(0.5,cycle),NULL);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Set runtime options   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Solve nonlinear system     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSSolve(ts,X);CHKERRQ(ierr);  ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);  ierr = TSGetTimeStepNumber(ts,&steps);CHKERRQ(ierr);  ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr);  ierr = VecMin(X,NULL,&xmin);CHKERRQ(ierr);  ierr = VecMax(X,NULL,&xmax);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after % 3D steps. Range [%6.4f,%6.4f]/n",TSConvergedReasons[reason],(double)ftime,steps,(double)xmin,(double)xmax);CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Free work space.   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = MatDestroy(&J);CHKERRQ(ierr);  ierr = VecDestroy(&X);CHKERRQ(ierr);  ierr = TSDestroy(&ts);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  return ierr;}

/*@C     PetscViewerFileGetName - Gets the name of the file the PetscViewer uses.    Not Collective  Input Parameter:.  viewer - the PetscViewer; either ASCII or binary  Output Parameter:.  name - the name of the file it is using    Level: advanced.seealso: PetscViewerCreate(), PetscViewerSetType(), PetscViewerASCIIOpen(), PetscViewerBinaryOpen(), PetscViewerFileSetName()@*/PetscErrorCode  PetscViewerFileGetName(PetscViewer viewer,const char **name){  PetscErrorCode ierr;  PetscFunctionBegin;  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,1);  ierr = PetscTryMethod(viewer,"PetscViewerFileGetName_C",(PetscViewer,const char**),(viewer,name));CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*@C    PetscViewerASCIISynchronizedPrintf - Prints synchronized output to the specified file from    several processors.  Output of the first processor is followed by that of the    second, etc.    Not Collective, must call collective PetscViewerFlush() to get the results out    Input Parameters:+   viewer - the ASCII PetscViewer-   format - the usual printf() format string    Level: intermediate    Notes: You must have previously called PetscViewerASCIISynchronizeAllow() to allow this routine to be called.    Fortran Note:      Can only print a single character* string.seealso: PetscSynchronizedPrintf(), PetscSynchronizedFlush(), PetscFPrintf(),          PetscFOpen(), PetscViewerFlush(), PetscViewerASCIIGetPointer(), PetscViewerDestroy(), PetscViewerASCIIOpen(),          PetscViewerASCIIPrintf(), PetscViewerASCIISynchronizedAllow()@*/PetscErrorCode  PetscViewerASCIISynchronizedPrintf(PetscViewer viewer,const char format[],...){  PetscViewer_ASCII *vascii = (PetscViewer_ASCII*)viewer->data;  PetscErrorCode    ierr;  PetscMPIInt       rank,size;  PetscInt          tab = vascii->tab;  MPI_Comm          comm;  FILE              *fp;  PetscBool         iascii;  int               err;  PetscFunctionBegin;  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,1);  PetscValidCharPointer(format,2);  ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);  if (!iascii) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not ASCII PetscViewer");  ierr = MPI_Comm_size(PetscObjectComm((PetscObject)viewer),&size);CHKERRQ(ierr);  if (size > 1 && !vascii->allowsynchronized) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"First call PetscViewerASCIISynchronizedAllow() to allow this call");  if (!viewer->ops->flush) PetscFunctionReturn(0); /* This viewer obtained via PetscViewerGetSubcomm_ASCII(), should not participate. */  ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);  fp   = vascii->fd;  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);  /* First processor prints immediately to fp */  if (!rank) {    va_list Argp;    while (tab--) {      ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"  ");CHKERRQ(ierr);    }    va_start(Argp,format);    ierr = (*PetscVFPrintf)(fp,format,Argp);CHKERRQ(ierr);    err  = fflush(fp);    if (err) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SYS,"fflush() failed on file");    if (petsc_history) {      va_start(Argp,format);      ierr = (*PetscVFPrintf)(petsc_history,format,Argp);CHKERRQ(ierr);      err  = fflush(petsc_history);      if (err) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SYS,"fflush() failed on file");    }    va_end(Argp);  } else { /* other processors add to local queue */    char        *string;    va_list     Argp;    size_t      fullLength;    PrintfQueue next;    ierr = PetscNew(&next);CHKERRQ(ierr);    if (petsc_printfqueue) {      petsc_printfqueue->next = next;      petsc_printfqueue       = next;    } else {      petsc_printfqueuebase = petsc_printfqueue = next;    }    petsc_printfqueuelength++;    next->size = QUEUESTRINGSIZE;    ierr       = PetscMalloc1(next->size, &next->string);CHKERRQ(ierr);    ierr       = PetscMemzero(next->string,next->size);CHKERRQ(ierr);    string     = next->string;    tab       *= 2;    while (tab--) {      *string++ = ' ';    }    va_start(Argp,format);    ierr = PetscVSNPrintf(string,next->size-2*vascii->tab,format,&fullLength,Argp);CHKERRQ(ierr);    va_end(Argp);  }  PetscFunctionReturn(0);}

PetscErrorCode MatGetSubMatrices_MPIDense_Local(Mat C,PetscInt ismax,const IS isrow[],const IS iscol[],MatReuse scall,Mat *submats){  Mat_MPIDense   *c = (Mat_MPIDense*)C->data;  Mat            A  = c->A;  Mat_SeqDense   *a = (Mat_SeqDense*)A->data,*mat;  PetscErrorCode ierr;  PetscMPIInt    rank,size,tag0,tag1,idex,end,i;  PetscInt       N = C->cmap->N,rstart = C->rmap->rstart,count;  const PetscInt **irow,**icol,*irow_i;  PetscInt       *nrow,*ncol,*w1,*w3,*w4,*rtable,start;  PetscInt       **sbuf1,m,j,k,l,ct1,**rbuf1,row,proc;  PetscInt       nrqs,msz,**ptr,*ctr,*pa,*tmp,bsz,nrqr;  PetscInt       is_no,jmax,**rmap,*rmap_i;  PetscInt       ctr_j,*sbuf1_j,*rbuf1_i;  MPI_Request    *s_waits1,*r_waits1,*s_waits2,*r_waits2;  MPI_Status     *r_status1,*r_status2,*s_status1,*s_status2;  MPI_Comm       comm;  PetscScalar    **rbuf2,**sbuf2;  PetscBool      sorted;  PetscFunctionBegin;  ierr = PetscObjectGetComm((PetscObject)C,&comm);CHKERRQ(ierr);  tag0 = ((PetscObject)C)->tag;  size = c->size;  rank = c->rank;  m    = C->rmap->N;  /* Get some new tags to keep the communication clean */  ierr = PetscObjectGetNewTag((PetscObject)C,&tag1);CHKERRQ(ierr);  /* Check if the col indices are sorted */  for (i=0; i<ismax; i++) {    ierr = ISSorted(isrow[i],&sorted);CHKERRQ(ierr);    if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted");    ierr = ISSorted(iscol[i],&sorted);CHKERRQ(ierr);    if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted");  }  ierr = PetscMalloc5(ismax,const PetscInt*,&irow,ismax,const PetscInt*,&icol,ismax,PetscInt,&nrow,ismax,PetscInt,&ncol,m,PetscInt,&rtable);CHKERRQ(ierr);  for (i=0; i<ismax; i++) {    ierr = ISGetIndices(isrow[i],&irow[i]);CHKERRQ(ierr);    ierr = ISGetIndices(iscol[i],&icol[i]);CHKERRQ(ierr);    ierr = ISGetLocalSize(isrow[i],&nrow[i]);CHKERRQ(ierr);    ierr = ISGetLocalSize(iscol[i],&ncol[i]);CHKERRQ(ierr);  }  /* Create hash table for the mapping :row -> proc*/  for (i=0,j=0; i<size; i++) {    jmax = C->rmap->range[i+1];    for (; j<jmax; j++) rtable[j] = i;  }  /* evaluate communication - mesg to who,length of mesg, and buffer space     required. Based on this, buffers are allocated, and data copied into them*/  ierr = PetscMalloc3(2*size,PetscInt,&w1,size,PetscInt,&w3,size,PetscInt,&w4);CHKERRQ(ierr);  ierr = PetscMemzero(w1,size*2*sizeof(PetscInt));CHKERRQ(ierr); /* initialize work vector*/  ierr = PetscMemzero(w3,size*sizeof(PetscInt));CHKERRQ(ierr); /* initialize work vector*/  for (i=0; i<ismax; i++) {    ierr   = PetscMemzero(w4,size*sizeof(PetscInt));CHKERRQ(ierr); /* initialize work vector*/    jmax   = nrow[i];    irow_i = irow[i];    for (j=0; j<jmax; j++) {      row  = irow_i[j];      proc = rtable[row];      w4[proc]++;    }    for (j=0; j<size; j++) {      if (w4[j]) { w1[2*j] += w4[j];  w3[j]++;}    }  }  nrqs       = 0;              /* no of outgoing messages */  msz        = 0;              /* total mesg length (for all procs) */  w1[2*rank] = 0;              /* no mesg sent to self */  w3[rank]   = 0;  for (i=0; i<size; i++) {    if (w1[2*i])  { w1[2*i+1] = 1; nrqs++;} /* there exists a message to proc i */  }  ierr = PetscMalloc((nrqs+1)*sizeof(PetscInt),&pa);CHKERRQ(ierr); /*(proc -array)*/  for (i=0,j=0; i<size; i++) {    if (w1[2*i]) { pa[j] = i; j++; }  }  /* Each message would have a header = 1 + 2*(no of IS) + data */  for (i=0; i<nrqs; i++) {    j        = pa[i];    w1[2*j] += w1[2*j+1] + 2* w3[j];    msz     += w1[2*j];  }  /* Do a global reduction to determine how many messages to expect*/  ierr = PetscMaxSum(comm,w1,&bsz,&nrqr);CHKERRQ(ierr);  /* Allocate memory for recv buffers . Make sure rbuf1[0] exists by adding 1 to the buffer length */  ierr = PetscMalloc((nrqr+1)*sizeof(PetscInt*),&rbuf1);CHKERRQ(ierr);  ierr = PetscMalloc(nrqr*bsz*sizeof(PetscInt),&rbuf1[0]);CHKERRQ(ierr);  for (i=1; i<nrqr; ++i) rbuf1[i] = rbuf1[i-1] + bsz;  /* Post the receives */  ierr = PetscMalloc((nrqr+1)*sizeof(MPI_Request),&r_waits1);CHKERRQ(ierr);  for (i=0; i<nrqr; ++i) {//.........这里部分代码省略.........

int main(int argc,char *argv[]){  PetscErrorCode ierr;  DM             da;  PetscInt       dim = 2,m,n,p,i;  const PetscInt *lx,*ly,*lz;  PetscMPIInt    rank,size;  PetscInitialize(&argc,&argv,0,help);  ierr = PetscOptionsGetInt(0,"-dim",&dim,0);CHKERRQ(ierr);  switch (dim) {  case 2:    ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, -3,-5,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);CHKERRQ(ierr);    break;  case 3:    ierr = DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, -3,-5,-7,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,NULL,&da);CHKERRQ(ierr);    break;  default: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"No support for %D dimensions",dim);  }  ierr = DMDAGetInfo(da, 0, 0,0,0, &m,&n,&p, 0,0, 0,0,0,0);CHKERRQ(ierr);  ierr = DMDAGetOwnershipRanges(da,&lx,&ly,&lz);CHKERRQ(ierr);  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);  for (i=0; i<size; i++) {    ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);    if (i == rank) {      ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF,"[%d] lx ly%s/n",rank,dim>2 ? " lz" : "");CHKERRQ(ierr);      ierr = PetscIntView(m,lx,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);      ierr = PetscIntView(n,ly,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);      if (dim > 2) {ierr = PetscIntView(n,lz,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}    }    ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr);  }  ierr = DMDestroy(&da);CHKERRQ(ierr);  PetscFinalize();  return 0;}

PetscErrorCode  TSGLAdaptChoose(TSGLAdapt adapt,PetscInt n,const PetscInt orders[],const PetscReal errors[],const PetscReal cost[],PetscInt cur,PetscReal h,PetscReal tleft,PetscInt *next_sc,PetscReal *next_h,PetscBool  *finish){  PetscErrorCode ierr;  PetscFunctionBegin;  PetscValidHeaderSpecific(adapt,TSGLADAPT_CLASSID,1);  PetscValidIntPointer(orders,3);  PetscValidPointer(errors,4);  PetscValidPointer(cost,5);  PetscValidIntPointer(next_sc,9);  PetscValidPointer(next_h,10);  PetscValidIntPointer(finish,11);  ierr = (*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish);CHKERRQ(ierr);  PetscFunctionReturn(0);}