C++ PetscMemzero函数代码示例

PETSC_EXTERN PetscErrorCode PetscViewerCreate_Draw(PetscViewer viewer)
{
  PetscInt         i;
  PetscErrorCode   ierr;
  PetscViewer_Draw *vdraw;

  PetscFunctionBegin;
  ierr         = PetscNewLog(viewer,PetscViewer_Draw,&vdraw);CHKERRQ(ierr);
  viewer->data = (void*)vdraw;

  viewer->ops->flush            = PetscViewerFlush_Draw;
  viewer->ops->destroy          = PetscViewerDestroy_Draw;
  viewer->ops->setfromoptions   = PetscViewerSetFromOptions_Draw;
  viewer->ops->getsingleton     = PetscViewerGetSingleton_Draw;
  viewer->ops->restoresingleton = PetscViewerRestoreSingleton_Draw;

  /* these are created on the fly if requested */
  vdraw->draw_max  = 5;
  vdraw->draw_base = 0;
  vdraw->w         = PETSC_DECIDE;
  vdraw->h         = PETSC_DECIDE;

  ierr = PetscMalloc3(vdraw->draw_max,PetscDraw,&vdraw->draw,vdraw->draw_max,PetscDrawLG,&vdraw->drawlg,vdraw->draw_max,PetscDrawAxis,&vdraw->drawaxis);CHKERRQ(ierr);
  ierr = PetscMemzero(vdraw->draw,vdraw->draw_max*sizeof(PetscDraw));CHKERRQ(ierr);
  ierr = PetscMemzero(vdraw->drawlg,vdraw->draw_max*sizeof(PetscDrawLG));CHKERRQ(ierr);
  ierr = PetscMemzero(vdraw->drawaxis,vdraw->draw_max*sizeof(PetscDrawAxis));CHKERRQ(ierr);
  for (i=0; i<vdraw->draw_max; i++) {
    vdraw->draw[i]     = 0;
    vdraw->drawlg[i]   = 0;
    vdraw->drawaxis[i] = 0;
  }
  vdraw->singleton_made = PETSC_FALSE;
  PetscFunctionReturn(0);
}

PetscErrorCode OpForcing(Op op,DM dm,Vec F) {
  PetscErrorCode ierr;
  Vec X,Floc;
  DM dmx;
  const PetscScalar *x;
  PetscScalar *f;
  const PetscReal *B,*D,*w3;
  PetscReal L[3];
  PetscInt nelem,ne = op->ne,P,Q,P3,Q3;

  PetscFunctionBegin;
  ierr = PetscLogEventBegin(OP_Forcing,dm,F,0,0);CHKERRQ(ierr);
  ierr = DMFEGetTensorEval(dm,&P,&Q,&B,&D,NULL,NULL,&w3);CHKERRQ(ierr);
  P3 = P*P*P;
  Q3 = Q*Q*Q;

  ierr = DMFEGetUniformCoordinates(dm,L);CHKERRQ(ierr);
  ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);
  ierr = DMGetCoordinateDM(dm,&dmx);CHKERRQ(ierr);
  ierr = DMGetCoordinatesLocal(dm,&X);CHKERRQ(ierr);
  ierr = DMFEGetNumElements(dm,&nelem);CHKERRQ(ierr);
  ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
  ierr = VecZeroEntries(Floc);CHKERRQ(ierr);
  ierr = VecGetArray(Floc,&f);CHKERRQ(ierr);

  for (PetscInt e=0; e<nelem; e+=ne) {
    PetscScalar fe[op->dof*P3*ne]_align,fq[op->dof][Q3][ne]_align,xe[3*P3*ne]_align,xq[3][Q3][ne]_align,dx[3][3][Q3][ne]_align,wdxdet[Q3][ne]_align;

    ierr = DMFEExtractElements(dmx,x,e,ne,xe);CHKERRQ(ierr);
    ierr = PetscMemzero(xq,sizeof xq);CHKERRQ(ierr);
    ierr = TensorContract(op->Tensor3,B,B,B,TENSOR_EVAL,xe,xq[0][0]);CHKERRQ(ierr);
    ierr = PetscMemzero(dx,sizeof dx);CHKERRQ(ierr);
    ierr = TensorContract(op->Tensor3,D,B,B,TENSOR_EVAL,xe,dx[0][0][0]);CHKERRQ(ierr);
    ierr = TensorContract(op->Tensor3,B,D,B,TENSOR_EVAL,xe,dx[1][0][0]);CHKERRQ(ierr);
    ierr = TensorContract(op->Tensor3,B,B,D,TENSOR_EVAL,xe,dx[2][0][0]);CHKERRQ(ierr);
    ierr = PointwiseJacobianInvert(ne,Q*Q*Q,w3,dx,wdxdet);CHKERRQ(ierr);

    for (PetscInt i=0; i<Q3; i++) {
      for (PetscInt l=0; l<ne; l++) {
        PetscReal xx[] = {xq[0][i][l],xq[1][i][l],xq[2][i][l]};
        PetscScalar fql[op->dof];
        ierr = (op->PointwiseForcing)(op,xx,L,fql);CHKERRQ(ierr);
        for (PetscInt d=0; d<op->dof; d++) fq[d][i][l] = wdxdet[i][l] * fql[d];
      }
    }
    ierr = PetscMemzero(fe,sizeof fe);CHKERRQ(ierr);
    ierr = TensorContract(op->TensorDOF,B,B,B,TENSOR_TRANSPOSE,fq[0][0],fe);CHKERRQ(ierr);
    ierr = DMFESetElements(dm,f,e,ne,ADD_VALUES,DOMAIN_INTERIOR,fe);CHKERRQ(ierr);
  }
  ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
  ierr = VecRestoreArray(Floc,&f);CHKERRQ(ierr);
  ierr = VecZeroEntries(F);CHKERRQ(ierr);
  ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
  ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
  ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);
  ierr = PetscLogEventEnd(OP_Forcing,dm,F,0,0);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

int main(int argc, char *argv[])
{
  PetscErrorCode ierr;
  PetscInt       opts[6] = {0};
  PetscBool      hascl   = PETSC_FALSE,hasstr = PETSC_FALSE;

  ierr = PetscInitialize(&argc,&argv,0,help);CHKERRQ(ierr);
  ierr = PetscOptionsSetValue(NULL,"-zero","0");CHKERRQ(ierr);
  ierr = PetscOptionsPrefixPush(NULL,"a_");CHKERRQ(ierr);
  ierr = PetscOptionsSetValue(NULL,"-one","1");CHKERRQ(ierr);
  ierr = PetscOptionsPrefixPush(NULL,"bb_");CHKERRQ(ierr);
  ierr = PetscOptionsSetValue(NULL,"-two","2");CHKERRQ(ierr);
  ierr = PetscOptionsPrefixPop(NULL);CHKERRQ(ierr);
  ierr = PetscOptionsSetValue(NULL,"-three","3");CHKERRQ(ierr);
  ierr = PetscOptionsPrefixPush(NULL,"cc_");CHKERRQ(ierr);
  ierr = PetscOptionsPrefixPush(NULL,"ddd_");CHKERRQ(ierr);
  ierr = PetscOptionsSetValue(NULL,"-four","4");CHKERRQ(ierr);
  ierr = PetscOptionsPrefixPop(NULL);CHKERRQ(ierr);
  ierr = PetscOptionsPrefixPop(NULL);CHKERRQ(ierr);
  ierr = PetscOptionsPrefixPop(NULL);CHKERRQ(ierr);
  ierr = PetscOptionsSetValue(NULL,"-five","5");CHKERRQ(ierr);

  ierr = PetscOptionsGetInt(NULL,0,"-zero",&opts[0],0);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,0,"-a_one",&opts[1],0);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,0,"-a_bb_two",&opts[2],0);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,0,"-a_three",&opts[3],0);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,0,"-a_cc_ddd_four",&opts[4],0);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,0,"-five",&opts[5],0);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"opts = {%D %D %D %D %D %D}\n",opts[0],opts[1],opts[2],opts[3],opts[4],opts[5]);CHKERRQ(ierr);

  ierr = PetscOptionsGetBool(NULL,0,"-cl",&hascl,0);CHKERRQ(ierr);
  if (hascl) {
    ierr = PetscMemzero(opts,sizeof(opts));CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,0,"-cl_zero",&opts[0],0);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,0,"-cl_a_one",&opts[1],0);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,0,"-cl_a_bb_two",&opts[2],0);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,0,"-cl_a_three",&opts[3],0);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,0,"-cl_a_cc_ddd_four",&opts[4],0);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,0,"-cl_five",&opts[5],0);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"cl_opts = {%D %D %D %D %D %D}\n",opts[0],opts[1],opts[2],opts[3],opts[4],opts[5]);CHKERRQ(ierr);
  }

  ierr = PetscOptionsGetBool(NULL,0,"-str",&hasstr,0);CHKERRQ(ierr);
  if (hasstr) {
    ierr = PetscOptionsInsertString(NULL,"-prefix_push str_ -zero 100 -prefix_push a_ -one 101 -prefix_push bb_ -two 102 -prefix_pop -three 103 -prefix_push cc_ -prefix_push ddd_ -four 104 -prefix_pop -prefix_pop -prefix_pop -five 105 -prefix_pop");CHKERRQ(ierr);
    ierr = PetscMemzero(opts,sizeof(opts));CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,0,"-str_zero",&opts[0],0);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,0,"-str_a_one",&opts[1],0);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,0,"-str_a_bb_two",&opts[2],0);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,0,"-str_a_three",&opts[3],0);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,0,"-str_a_cc_ddd_four",&opts[4],0);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(NULL,0,"-str_five",&opts[5],0);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"str_opts = {%D %D %D %D %D %D}\n",opts[0],opts[1],opts[2],opts[3],opts[4],opts[5]);CHKERRQ(ierr);
  }

  ierr = PetscFinalize();CHKERRQ(ierr);
  return 0;
}

PetscErrorCode MatMPIAIJDiagonalScaleLocalSetUp(Mat inA,Vec scale)
{
  Mat_MPIAIJ     *ina = (Mat_MPIAIJ*) inA->data; /*access private part of matrix */
  PetscErrorCode ierr;
  PetscInt       i,n,nt,cstart,cend,no,*garray = ina->garray,*lindices;
  PetscInt       *r_rmapd,*r_rmapo;
  
  PetscFunctionBegin;
  ierr = MatGetOwnershipRange(inA,&cstart,&cend);CHKERRQ(ierr);
  ierr = MatGetSize(ina->A,PETSC_NULL,&n);CHKERRQ(ierr);
  ierr = PetscMalloc((inA->rmap->mapping->n+1)*sizeof(PetscInt),&r_rmapd);CHKERRQ(ierr);
  ierr = PetscMemzero(r_rmapd,inA->rmap->mapping->n*sizeof(PetscInt));CHKERRQ(ierr);
  nt   = 0;
  for (i=0; i<inA->rmap->mapping->n; i++) {
    if (inA->rmap->mapping->indices[i] >= cstart && inA->rmap->mapping->indices[i] < cend) {
      nt++;
      r_rmapd[i] = inA->rmap->mapping->indices[i] + 1;
    }
  }
  if (nt != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Hmm nt %D n %D",nt,n);
  ierr = PetscMalloc((n+1)*sizeof(PetscInt),&auglyrmapd);CHKERRQ(ierr);
  for (i=0; i<inA->rmap->mapping->n; i++) {
    if (r_rmapd[i]){
      auglyrmapd[(r_rmapd[i]-1)-cstart] = i;
    }
  }
  ierr = PetscFree(r_rmapd);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&auglydd);CHKERRQ(ierr);

  ierr = PetscMalloc((inA->cmap->N+1)*sizeof(PetscInt),&lindices);CHKERRQ(ierr);
  ierr = PetscMemzero(lindices,inA->cmap->N*sizeof(PetscInt));CHKERRQ(ierr);
  for (i=0; i<ina->B->cmap->n; i++) {
    lindices[garray[i]] = i+1;
  }
  no   = inA->rmap->mapping->n - nt;
  ierr = PetscMalloc((inA->rmap->mapping->n+1)*sizeof(PetscInt),&r_rmapo);CHKERRQ(ierr);
  ierr = PetscMemzero(r_rmapo,inA->rmap->mapping->n*sizeof(PetscInt));CHKERRQ(ierr);
  nt   = 0;
  for (i=0; i<inA->rmap->mapping->n; i++) {
    if (lindices[inA->rmap->mapping->indices[i]]) {
      nt++;
      r_rmapo[i] = lindices[inA->rmap->mapping->indices[i]];
    }
  }
  if (nt > no) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Hmm nt %D no %D",nt,n);
  ierr = PetscFree(lindices);CHKERRQ(ierr);
  ierr = PetscMalloc((nt+1)*sizeof(PetscInt),&auglyrmapo);CHKERRQ(ierr);
  for (i=0; i<inA->rmap->mapping->n; i++) {
    if (r_rmapo[i]){
      auglyrmapo[(r_rmapo[i]-1)] = i;
    }
  }
  ierr = PetscFree(r_rmapo);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,nt,&auglyoo);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

PetscErrorCode  DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
{
  PetscErrorCode         ierr;
  PetscInt               m,n,M,N,nDM,i;
  struct DMCompositeLink *nextc;
  struct DMCompositeLink *nextf;
  Vec                    gcoarse,gfine,*vecs;
  DM_Composite           *comcoarse = (DM_Composite*)coarse->data;
  DM_Composite           *comfine = (DM_Composite*)fine->data;
  Mat                    *mats;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(coarse,DM_CLASSID,1);
  PetscValidHeaderSpecific(fine,DM_CLASSID,2);
  ierr = DMSetUp(coarse);CHKERRQ(ierr);
  ierr = DMSetUp(fine);CHKERRQ(ierr);
  /* use global vectors only for determining matrix layout */
  ierr = DMGetGlobalVector(coarse,&gcoarse);CHKERRQ(ierr);
  ierr = DMGetGlobalVector(fine,&gfine);CHKERRQ(ierr);
  ierr = VecGetLocalSize(gcoarse,&n);CHKERRQ(ierr);
  ierr = VecGetLocalSize(gfine,&m);CHKERRQ(ierr);
  ierr = VecGetSize(gcoarse,&N);CHKERRQ(ierr);
  ierr = VecGetSize(gfine,&M);CHKERRQ(ierr);
  ierr = DMRestoreGlobalVector(coarse,&gcoarse);CHKERRQ(ierr);
  ierr = DMRestoreGlobalVector(fine,&gfine);CHKERRQ(ierr);

  nDM = comfine->nDM;
  if (nDM != comcoarse->nDM) SETERRQ2(((PetscObject)fine)->comm,PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
  ierr = PetscMalloc(nDM*nDM*sizeof(Mat),&mats);CHKERRQ(ierr);
  ierr = PetscMemzero(mats,nDM*nDM*sizeof(Mat));CHKERRQ(ierr);
  if (v) {
    ierr = PetscMalloc(nDM*sizeof(Vec),&vecs);CHKERRQ(ierr);
    ierr = PetscMemzero(vecs,nDM*sizeof(Vec));CHKERRQ(ierr);
  }

  /* loop over packed objects, handling one at at time */
  for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
    if (!v) {
      ierr = DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],PETSC_NULL);CHKERRQ(ierr);
    } else {
      ierr = DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);CHKERRQ(ierr);
    }
  }
  ierr = MatCreateNest(((PetscObject)fine)->comm,nDM,PETSC_NULL,nDM,PETSC_NULL,mats,A);CHKERRQ(ierr);
  if (v) {
    ierr = VecCreateNest(((PetscObject)fine)->comm,nDM,PETSC_NULL,vecs,v);CHKERRQ(ierr);
  }
  for (i=0; i<nDM*nDM; i++) {ierr = MatDestroy(&mats[i]);CHKERRQ(ierr);}
  ierr = PetscFree(mats);CHKERRQ(ierr);
  if (v) {
    for (i=0; i<nDM; i++) {ierr = VecDestroy(&vecs[i]);CHKERRQ(ierr);}
    ierr = PetscFree(vecs);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

/*@C
   ISColoringGetIS - Extracts index sets from the coloring context

   Collective on ISColoring

   Input Parameter:
.  iscoloring - the coloring context

   Output Parameters:
+  nn - number of index sets in the coloring context
-  is - array of index sets

   Level: advanced

.seealso: ISColoringRestoreIS(), ISColoringView()
@*/
PetscErrorCode  ISColoringGetIS(ISColoring iscoloring,PetscInt *nn,IS *isis[])
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidPointer(iscoloring,1);

  if (nn) *nn = iscoloring->n;
  if (isis) {
    if (!iscoloring->is) {
      PetscInt        *mcolors,**ii,nc = iscoloring->n,i,base, n = iscoloring->N;
      ISColoringValue *colors = iscoloring->colors;
      IS              *is;

#if defined(PETSC_USE_DEBUG)
      for (i=0; i<n; i++) {
        if (((PetscInt)colors[i]) >= nc) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coloring is our of range index %d value %d number colors %d",(int)i,(int)colors[i],(int)nc);
      }
#endif

      /* generate the lists of nodes for each color */
      ierr = PetscMalloc(nc*sizeof(PetscInt),&mcolors);CHKERRQ(ierr);
      ierr = PetscMemzero(mcolors,nc*sizeof(PetscInt));CHKERRQ(ierr);
      for (i=0; i<n; i++) mcolors[colors[i]]++;

      ierr = PetscMalloc(nc*sizeof(PetscInt*),&ii);CHKERRQ(ierr);
      ierr = PetscMalloc(n*sizeof(PetscInt),&ii[0]);CHKERRQ(ierr);
      for (i=1; i<nc; i++) ii[i] = ii[i-1] + mcolors[i-1];
      ierr = PetscMemzero(mcolors,nc*sizeof(PetscInt));CHKERRQ(ierr);

      if (iscoloring->ctype == IS_COLORING_GLOBAL) {
        ierr = MPI_Scan(&iscoloring->N,&base,1,MPIU_INT,MPI_SUM,iscoloring->comm);CHKERRQ(ierr);
        base -= iscoloring->N;
        for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
      } else if (iscoloring->ctype == IS_COLORING_GHOSTED) {
        for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i;   /* local idx */
      } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this ISColoringType type");

      ierr = PetscMalloc(nc*sizeof(IS),&is);CHKERRQ(ierr);
      for (i=0; i<nc; i++) {
        ierr = ISCreateGeneral(iscoloring->comm,mcolors[i],ii[i],PETSC_COPY_VALUES,is+i);CHKERRQ(ierr);
      }

      iscoloring->is = is;
      ierr = PetscFree(ii[0]);CHKERRQ(ierr);
      ierr = PetscFree(ii);CHKERRQ(ierr);
      ierr = PetscFree(mcolors);CHKERRQ(ierr);
    }
    *isis = iscoloring->is;
  }
  PetscFunctionReturn(0);
}

PetscErrorCode    KSPSetUp_AGMRES(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       hes;
  PetscInt       nloc;
  KSP_AGMRES     *agmres = (KSP_AGMRES*)ksp->data;
  PetscInt       neig    = agmres->neig;
  PetscInt       max_k   = agmres->max_k;
  PetscInt       N       = MAXKSPSIZE;
  PetscInt       lwork   = PetscMax(8 * N + 16, 4 * neig * (N - neig));

  PetscFunctionBegin;
  if (ksp->pc_side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"no symmetric preconditioning for KSPAGMRES");
  max_k = agmres->max_k;
  N     = MAXKSPSIZE;
  /* Preallocate space during the call to KSPSetup_GMRES for the Krylov basis */
  agmres->q_preallocate = PETSC_TRUE; /* No allocation on the fly */
  /* Preallocate space to compute later the eigenvalues in GMRES */
  ksp->calc_sings = PETSC_TRUE;
  agmres->max_k   = N; /* Set the augmented size to be allocated in KSPSetup_GMRES */
  ierr            = KSPSetUp_DGMRES(ksp);CHKERRQ(ierr);
  agmres->max_k   = max_k;
  hes             = (N + 1) * (N + 1);

  /* Data for the Newton basis GMRES */
  ierr = PetscCalloc4(max_k,&agmres->Rshift,max_k,&agmres->Ishift,hes,&agmres->Rloc,((N+1)*4),&agmres->wbufptr);CHKERRQ(ierr);
  ierr = PetscMalloc7((N+1),&agmres->Scale,(N+1),&agmres->sgn,(N+1),&agmres->tloc,(N+1),&agmres->temp,(N+1),&agmres->tau,lwork,&agmres->work,(N+1),&agmres->nrs);CHKERRQ(ierr);
  ierr = PetscMemzero(agmres->Scale, (N+1)*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr = PetscMemzero(agmres->sgn, (N+1)*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr = PetscMemzero(agmres->tloc, (N+1)*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr = PetscMemzero(agmres->temp, (N+1)*sizeof(PetscScalar));CHKERRQ(ierr);

  /* Allocate space for the vectors in the orthogonalized basis*/
  ierr = VecGetLocalSize(agmres->vecs[0], &nloc);CHKERRQ(ierr);
  ierr = PetscMalloc1(nloc*(N+1), &agmres->Qloc);CHKERRQ(ierr);

  /* Init the ring of processors for the roddec orthogonalization */
  ierr = KSPAGMRESRoddecInitNeighboor(ksp);CHKERRQ(ierr);

  if (agmres->neig < 1) PetscFunctionReturn(0);

  /* Allocate space for the deflation */
  ierr = PetscMalloc1(N, &agmres->select);CHKERRQ(ierr);
  ierr = VecDuplicateVecs(VEC_V(0), N, &agmres->TmpU);CHKERRQ(ierr);
  ierr = PetscMalloc2(N*N, &agmres->MatEigL, N*N, &agmres->MatEigR);CHKERRQ(ierr);
  /*  ierr = PetscMalloc6(N*N, &agmres->Q, N*N, &agmres->Z, N, &agmres->wr, N, &agmres->wi, N, &agmres->beta, N, &agmres->modul);CHKERRQ(ierr); */
  ierr = PetscMalloc3(N*N, &agmres->Q, N*N, &agmres->Z, N, &agmres->beta);CHKERRQ(ierr);
  ierr = PetscMalloc2((N+1),&agmres->perm,(2*neig*N),&agmres->iwork);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode VecCreate_Seq(Vec V)
{
  Vec_Seq        *s;
  PetscScalar    *array;
  PetscErrorCode ierr;
  PetscInt       n = PetscMax(V->map->n,V->map->N);
  PetscMPIInt    size;

  PetscFunctionBegin;
  ierr = MPI_Comm_size(PetscObjectComm((PetscObject)V),&size);CHKERRQ(ierr);
  if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot create VECSEQ on more than one process");
#if !defined(PETSC_USE_MIXED_PRECISION)
  ierr = PetscMalloc1(n,&array);CHKERRQ(ierr);
  ierr = PetscLogObjectMemory((PetscObject)V, n*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr = VecCreate_Seq_Private(V,array);CHKERRQ(ierr);

  s                  = (Vec_Seq*)V->data;
  s->array_allocated = array;

  ierr = VecSet(V,0.0);CHKERRQ(ierr);
#else
  switch (((PetscObject)V)->precision) {
  case PETSC_PRECISION_SINGLE: {
    float *aarray;

    ierr = PetscMalloc1(n,&aarray);CHKERRQ(ierr);
    ierr = PetscLogObjectMemory((PetscObject)V, n*sizeof(float));CHKERRQ(ierr);
    ierr = PetscMemzero(aarray,n*sizeof(float));CHKERRQ(ierr);
    ierr = VecCreate_Seq_Private(V,aarray);CHKERRQ(ierr);

    s                  = (Vec_Seq*)V->data;
    s->array_allocated = (PetscScalar*)aarray;
  } break;
  case PETSC_PRECISION_DOUBLE: {
    double *aarray;

    ierr = PetscMalloc1(n,&aarray);CHKERRQ(ierr);
    ierr = PetscLogObjectMemory((PetscObject)V, n*sizeof(double));CHKERRQ(ierr);
    ierr = PetscMemzero(aarray,n*sizeof(double));CHKERRQ(ierr);
    ierr = VecCreate_Seq_Private(V,aarray);CHKERRQ(ierr);

    s                  = (Vec_Seq*)V->data;
    s->array_allocated = (PetscScalar*)aarray;
  } break;
  default: SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"No support for mixed precision %d",(int)(((PetscObject)V)->precision));
  }
#endif
  PetscFunctionReturn(0);
}

PetscErrorCode EPSCreate_XD(EPS eps)
{
  PetscErrorCode ierr;
  EPS_DAVIDSON   *data;

  PetscFunctionBegin;
  eps->ops->solve                = EPSSolve_XD;
  eps->ops->setup                = EPSSetUp_XD;
  eps->ops->reset                = EPSReset_XD;
  eps->ops->backtransform        = EPSBackTransform_Default;
  eps->ops->computevectors       = EPSComputeVectors_XD;
  eps->ops->view                 = EPSView_XD;

  ierr = PetscNewLog(eps,&data);CHKERRQ(ierr);
  eps->data = (void*)data;
  ierr = PetscMemzero(&data->ddb,sizeof(dvdDashboard));CHKERRQ(ierr);

  /* Set default values */
  ierr = EPSXDSetKrylovStart_XD(eps,PETSC_FALSE);CHKERRQ(ierr);
  ierr = EPSXDSetBlockSize_XD(eps,1);CHKERRQ(ierr);
  ierr = EPSXDSetRestart_XD(eps,6,0);CHKERRQ(ierr);
  ierr = EPSXDSetInitialSize_XD(eps,5);CHKERRQ(ierr);
  ierr = EPSJDSetFix_JD(eps,0.01);CHKERRQ(ierr);
  ierr = EPSXDSetBOrth_XD(eps,PETSC_TRUE);CHKERRQ(ierr);
  ierr = EPSJDSetConstCorrectionTol_JD(eps,PETSC_TRUE);CHKERRQ(ierr);
  ierr = EPSXDSetWindowSizes_XD(eps,0,0);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
   PCSetType - Builds PC for a particular preconditioner type

   Collective on PC

   Input Parameter:
+  pc - the preconditioner context.
-  type - a known method

   Options Database Key:
.  -pc_type <type> - Sets PC type

   Use -help for a list of available methods (for instance,
   jacobi or bjacobi)

  Notes:
  See "petsc/include/petscpc.h" for available methods (for instance,
  PCJACOBI, PCILU, or PCBJACOBI).

  Normally, it is best to use the KSPSetFromOptions() command and
  then set the PC type from the options database rather than by using
  this routine.  Using the options database provides the user with
  maximum flexibility in evaluating the many different preconditioners.
  The PCSetType() routine is provided for those situations where it
  is necessary to set the preconditioner independently of the command
  line or options database.  This might be the case, for example, when
  the choice of preconditioner changes during the execution of the
  program, and the user's application is taking responsibility for
  choosing the appropriate preconditioner.  In other words, this
  routine is not for beginners.

  Level: intermediate

  Developer Note: PCRegister() is used to add preconditioner types to PCList from which they
  are accessed by PCSetType().

.keywords: PC, set, method, type

.seealso: KSPSetType(), PCType, PCRegister(), PCCreate(), KSPGetPC()

@*/
PetscErrorCode  PCSetType(PC pc,PCType type)
{
  PetscErrorCode ierr,(*r)(PC);
  PetscBool      match;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc,PC_CLASSID,1);
  PetscValidCharPointer(type,2);

  ierr = PetscObjectTypeCompare((PetscObject)pc,type,&match);CHKERRQ(ierr);
  if (match) PetscFunctionReturn(0);

  ierr =  PetscFunctionListFind(PCList,type,&r);CHKERRQ(ierr);
  if (!r) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested PC type %s",type);
  /* Destroy the previous private PC context */
  if (pc->ops->destroy) {
    ierr             =  (*pc->ops->destroy)(pc);CHKERRQ(ierr);
    pc->ops->destroy = NULL;
    pc->data         = 0;
  }
  ierr = PetscFunctionListDestroy(&((PetscObject)pc)->qlist);CHKERRQ(ierr);
  /* Reinitialize function pointers in PCOps structure */
  ierr = PetscMemzero(pc->ops,sizeof(struct _PCOps));CHKERRQ(ierr);
  /* XXX Is this OK?? */
  pc->modifysubmatrices  = 0;
  pc->modifysubmatricesP = 0;
  /* Call the PCCreate_XXX routine for this particular preconditioner */
  pc->setupcalled = 0;

  ierr = PetscObjectChangeTypeName((PetscObject)pc,type);CHKERRQ(ierr);
  ierr = (*r)(pc);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
   PetscObjectCopyFortranFunctionPointers - Copy function pointers to another object

   Logically Collective on PetscObject

   Input Parameter:
+  src - source object
-  dest - destination object

   Level: developer

   Note:
   Both objects must have the same class.
@*/
PetscErrorCode PetscObjectCopyFortranFunctionPointers(PetscObject src,PetscObject dest)
{
  PetscErrorCode ierr;
  PetscInt       cbtype,numcb[PETSC_FORTRAN_CALLBACK_MAXTYPE];

  PetscFunctionBegin;
  PetscValidHeader(src,1);
  PetscValidHeader(dest,2);
  if (src->classid != dest->classid) SETERRQ(src->comm,PETSC_ERR_ARG_INCOMP,"Objects must be of the same class");

  ierr = PetscFree(dest->fortran_func_pointers);CHKERRQ(ierr);
  ierr = PetscMalloc(src->num_fortran_func_pointers*sizeof(void(*)(void)),&dest->fortran_func_pointers);CHKERRQ(ierr);
  ierr = PetscMemcpy(dest->fortran_func_pointers,src->fortran_func_pointers,src->num_fortran_func_pointers*sizeof(void(*)(void)));CHKERRQ(ierr);

  dest->num_fortran_func_pointers = src->num_fortran_func_pointers;

  ierr = PetscFortranCallbackGetSizes(src->classid,&numcb[PETSC_FORTRAN_CALLBACK_CLASS],&numcb[PETSC_FORTRAN_CALLBACK_SUBTYPE]);CHKERRQ(ierr);
  for (cbtype=PETSC_FORTRAN_CALLBACK_CLASS; cbtype<PETSC_FORTRAN_CALLBACK_MAXTYPE; cbtype++) {
    ierr = PetscFree(dest->fortrancallback[cbtype]);CHKERRQ(ierr);
    ierr = PetscMalloc(numcb[cbtype]*sizeof(PetscFortranCallback),&dest->fortrancallback[cbtype]);CHKERRQ(ierr);
    ierr = PetscMemzero(dest->fortrancallback[cbtype],numcb[cbtype]*sizeof(PetscFortranCallback));CHKERRQ(ierr);
    ierr = PetscMemcpy(dest->fortrancallback[cbtype],src->fortrancallback[cbtype],src->num_fortrancallback[cbtype]*sizeof(PetscFortranCallback));CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,MatStructure *str,void *ctx)
{
  PetscErrorCode ierr;
  PetscInt       i,mx,xm,xs;
  PetscScalar    v[7],Hx;
  MatStencil     row,col[7];
  PetscScalar    lambda;
  DM             da;

  PetscFunctionBeginUser;
  ierr   = KSPGetDM(ksp,&da);CHKERRQ(ierr);
  ierr   = PetscMemzero(col,7*sizeof(MatStencil));CHKERRQ(ierr);
  ierr   = DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
  Hx     = 2.0*PETSC_PI / (PetscReal)(mx);
  ierr   = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
  lambda = 2.0*Hx;
  for (i=xs; i<xs+xm; i++) {
    row.i = i; row.j = 0; row.k = 0; row.c = 0;
    v[0]  = Hx;     col[0].i = i;   col[0].c = 0;
    v[1]  = lambda; col[1].i = i-1;   col[1].c = 1;
    v[2]  = -lambda;col[2].i = i+1; col[2].c = 1;
    ierr  = MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr);

    row.i = i; row.j = 0; row.k = 0; row.c = 1;
    v[0]  = lambda; col[0].i = i-1;   col[0].c = 0;
    v[1]  = Hx;     col[1].i = i;   col[1].c = 1;
    v[2]  = -lambda;col[2].i = i+1; col[2].c = 0;
    ierr  = MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatView(jac,PETSC_VIEWER_BINARY_(PETSC_COMM_SELF));CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
   PetscViewerSetType - Builds PetscViewer for a particular implementation.

   Collective on PetscViewer

   Input Parameter:
+  viewer      - the PetscViewer context
-  type        - for example, PETSCVIEWERASCII

   Options Database Command:
.  -draw_type  <type> - Sets the type; use -help for a list
    of available methods (for instance, ascii)

   Level: advanced

   Notes:
   See "include/petscviewer.h" for available methods (for instance,
   PETSCVIEWERSOCKET)

.seealso: PetscViewerCreate(), PetscViewerGetType(), PetscViewerType, PetscViewerSetFormat()
@*/
PetscErrorCode  PetscViewerSetType(PetscViewer viewer,PetscViewerType type)
{
  PetscErrorCode ierr,(*r)(PetscViewer);
  PetscBool      match;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,1);
  PetscValidCharPointer(type,2);
  ierr = PetscObjectTypeCompare((PetscObject)viewer,type,&match);CHKERRQ(ierr);
  if (match) PetscFunctionReturn(0);

  /* cleanup any old type that may be there */
  if (viewer->data) {
    ierr         = (*viewer->ops->destroy)(viewer);CHKERRQ(ierr);

    viewer->ops->destroy = NULL;
    viewer->data         = 0;
  }
  ierr = PetscMemzero(viewer->ops,sizeof(struct _PetscViewerOps));CHKERRQ(ierr);

  ierr =  PetscFunctionListFind(PetscViewerList,type,&r);CHKERRQ(ierr);
  if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown PetscViewer type given: %s",type);

  ierr = PetscObjectChangeTypeName((PetscObject)viewer,type);CHKERRQ(ierr);
  ierr = (*r)(viewer);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode MatCreate_Preallocator(Mat A)
{
  Mat_Preallocator *p;
  PetscErrorCode    ierr;

  PetscFunctionBegin;
  ierr = PetscNewLog(A, &p);CHKERRQ(ierr);
  A->data = (void *) p;

  p->ht  = NULL;
  p->dnz = NULL;
  p->onz = NULL;

  /* matrix ops */
  ierr = PetscMemzero(A->ops, sizeof(struct _MatOps));CHKERRQ(ierr);
  A->ops->destroy                 = MatDestroy_Preallocator;
  A->ops->setup                   = MatSetUp_Preallocator;
  A->ops->setvalues               = MatSetValues_Preallocator;
  A->ops->assemblybegin           = MatAssemblyBegin_Preallocator;
  A->ops->assemblyend             = MatAssemblyEnd_Preallocator;
  A->ops->view                    = MatView_Preallocator;
  A->ops->setoption               = MatSetOption_Preallocator;

  /* special MATPREALLOCATOR functions */
  ierr = PetscObjectComposeFunction((PetscObject) A, "MatPreallocatorPreallocate_C", MatPreallocatorPreallocate_Preallocator);CHKERRQ(ierr);
  ierr = PetscObjectChangeTypeName((PetscObject) A, MATPREALLOCATOR);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode MatGetDiagonalHermitian_Normal(Mat N,Vec v)
{
  Mat_Normal        *Na = (Mat_Normal*)N->data;
  Mat               A   = Na->A;
  PetscErrorCode    ierr;
  PetscInt          i,j,rstart,rend,nnz;
  const PetscInt    *cols;
  PetscScalar       *diag,*work,*values;
  const PetscScalar *mvalues;

  PetscFunctionBegin;
  ierr = PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);CHKERRQ(ierr);
  ierr = PetscMemzero(work,A->cmap->N*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
  for (i=rstart; i<rend; i++) {
    ierr = MatGetRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr);
    for (j=0; j<nnz; j++) {
      work[cols[j]] += mvalues[j]*PetscConj(mvalues[j]);
    }
    ierr = MatRestoreRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr);
  }
  ierr   = MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));CHKERRQ(ierr);
  rstart = N->cmap->rstart;
  rend   = N->cmap->rend;
  ierr   = VecGetArray(v,&values);CHKERRQ(ierr);
  ierr   = PetscMemcpy(values,diag+rstart,(rend-rstart)*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr   = VecRestoreArray(v,&values);CHKERRQ(ierr);
  ierr   = PetscFree2(diag,work);CHKERRQ(ierr);
  ierr   = VecScale(v,Na->scale);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}