C++ VecGetLocalSize函数代码示例

/** Create a cache for Dirichlet part of closure vector, and scatter from global closure to Dirichlet cache.

@arg[in] gvec Global vector
@arg[out] dcache New vector to hold the Dirichlet values
@arg[out] dscat Scatter from global closure to \a dcache

@note This could be local but it doesn't cost anything to make it global.
**/
dErr VecDohpCreateDirichletCache(Vec gvec,Vec *dcache,VecScatter *dscat)
{
  MPI_Comm comm;
  dErr     err;
  dBool    isdohp;
  IS       from;
  Vec      gc;
  dInt     n,nc,crstart;

  dFunctionBegin;
  dValidHeader(gvec,VEC_CLASSID,1);
  dValidPointer(dcache,2);
  dValidPointer(dscat,3);
  err = PetscTypeCompare((PetscObject)gvec,VECDOHP,&isdohp);dCHK(err);
  if (!isdohp) dERROR(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vec type %s",((PetscObject)gvec)->type_name);
  err = PetscObjectGetComm((PetscObject)gvec,&comm);dCHK(err);
  err = VecGetLocalSize(gvec,&n);dCHK(err);
  err = VecDohpGetClosure(gvec,&gc);dCHK(err);
  err = VecGetLocalSize(gc,&nc);dCHK(err);
  err = VecGetOwnershipRange(gc,&crstart,NULL);dCHK(err);
  err = VecCreateMPI(comm,nc-n,PETSC_DECIDE,dcache);dCHK(err);
  err = ISCreateStride(comm,nc-n,crstart+n,1,&from);dCHK(err);
  err = VecScatterCreate(gc,from,*dcache,NULL,dscat);dCHK(err);
  err = VecDohpRestoreClosure(gvec,&gc);dCHK(err);
  err = ISDestroy(&from);dCHK(err);
  /* \todo deal with rotations */
  dFunctionReturn(0);
}

static PetscErrorCode VecRestoreArray_Nest(Vec X,PetscScalar **x)
{
  Vec_Nest       *bx = (Vec_Nest*)X->data;
  PetscInt       i,m,rstart,rend;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = VecGetOwnershipRange(X,&rstart,&rend);CHKERRQ(ierr);
  ierr = VecGetLocalSize(X,&m);CHKERRQ(ierr);
  for (i=0; i<bx->nb; i++) {
    Vec            subvec = bx->v[i];
    IS             isy    = bx->is[i];
    PetscInt       j,sm;
    const PetscInt *ixy;
    PetscScalar    *y;
    ierr = VecGetLocalSize(subvec,&sm);CHKERRQ(ierr);
    ierr = VecGetArray(subvec,&y);CHKERRQ(ierr);
    ierr = ISGetIndices(isy,&ixy);CHKERRQ(ierr);
    for (j=0; j<sm; j++) {
      PetscInt ix = ixy[j];
      if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
      y[j] = (*x)[ix-rstart];
    }
    ierr = ISRestoreIndices(isy,&ixy);CHKERRQ(ierr);
    ierr = VecRestoreArray(subvec,&y);CHKERRQ(ierr);
  }
  ierr = PetscFree(*x);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode  MatDiagonalScaleLocal_MPIAIJ(Mat A,Vec scale)
{
  Mat_MPIAIJ        *a = (Mat_MPIAIJ*) A->data; /*access private part of matrix */
  PetscErrorCode    ierr;
  PetscInt          n,i;
  PetscScalar       *d,*o;
  const PetscScalar *s;

  PetscFunctionBegin;
  if (!auglyrmapd) {
    ierr = MatMPIAIJDiagonalScaleLocalSetUp(A,scale);CHKERRQ(ierr);
  }

  ierr = VecGetArrayRead(scale,&s);CHKERRQ(ierr);

  ierr = VecGetLocalSize(auglydd,&n);CHKERRQ(ierr);
  ierr = VecGetArray(auglydd,&d);CHKERRQ(ierr);
  for (i=0; i<n; i++) {
    d[i] = s[auglyrmapd[i]]; /* copy "diagonal" (true local) portion of scale into dd vector */
  }
  ierr = VecRestoreArray(auglydd,&d);CHKERRQ(ierr);
  /* column scale "diagonal" portion of local matrix */
  ierr = MatDiagonalScale(a->A,NULL,auglydd);CHKERRQ(ierr);

  ierr = VecGetLocalSize(auglyoo,&n);CHKERRQ(ierr);
  ierr = VecGetArray(auglyoo,&o);CHKERRQ(ierr);
  for (i=0; i<n; i++) {
    o[i] = s[auglyrmapo[i]]; /* copy "off-diagonal" portion of scale into oo vector */
  }
  ierr = VecRestoreArrayRead(scale,&s);CHKERRQ(ierr);
  ierr = VecRestoreArray(auglyoo,&o);CHKERRQ(ierr);
  /* column scale "off-diagonal" portion of local matrix */
  ierr = MatDiagonalScale(a->B,NULL,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);
}

void MultigridSolver_UpdateWorkVectors( MultigridSolver* self ) {
	MultigridSolver_Level*	level;
	//unsigned		rowSize, colSize;
	PetscInt		rowSize, colSize, vecSize;
	unsigned		l_i;

	assert( self && Stg_CheckType( self, MultigridSolver ) );

	for( l_i = 0; l_i < self->nLevels - 1; l_i++ ) {
		level = self->levels + l_i;

		//Matrix_GetLocalSize( MatrixSolver_GetMatrix( level->downSolver ), &rowSize, &colSize );
		MatGetLocalSize( level->downSolver->matrix, &rowSize, &colSize );

		VecGetLocalSize( level->workSol, &vecSize );

		//if( !level->workSol || Vector_GetLocalSize( level->workSol ) != rowSize ) {
		if( !level->workSol || vecSize != rowSize ) {
			//if( level->workSol )
			//	Stg_Class_RemoveRef( level->workSol );
			//Vector_Duplicate( self->curSolution, (void**)&level->workSol );
			//Vector_SetLocalSize( level->workSol, rowSize );
			if( level->workSol != PETSC_NULL )
				Stg_VecDestroy(&level->workSol );
			VecCreate( self->mgData->comm, &level->workSol );
			VecSetSizes( level->workSol, rowSize, PETSC_DECIDE );
			VecSetFromOptions( level->workSol );
#if( PETSC_VERSION_MAJOR <= 2 && PETSC_VERSION_MINOR >= 3 && PETSC_VERSION_SUBMINOR >= 3 )
			VecSetOption( level->workSol, VEC_IGNORE_NEGATIVE_INDICES );
#elif( PETSC_VERSION_MAJOR >= 3 )
			VecSetOption( level->workSol, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE );
#endif
		}

		VecGetLocalSize( level->workRHS, &vecSize );

		if( !level->workRHS || /*Vector_GetLocalSize( level->workRHS )*/vecSize != rowSize ) {
			//if( level->workRHS )
			//	Stg_Class_RemoveRef( level->workRHS );
			//Vector_Duplicate( self->curSolution, (void**)&level->workRHS );
			//Vector_SetLocalSize( level->workRHS, rowSize );
			if( level->workRHS != PETSC_NULL )
				Stg_VecDestroy(&level->workRHS );
			VecCreate( self->mgData->comm, &level->workRHS );
			VecSetSizes( level->workRHS, rowSize, PETSC_DECIDE );
			VecSetFromOptions( level->workRHS );
#if( PETSC_VERSION_MAJOR <= 2 && PETSC_VERSION_MINOR >= 3 && PETSC_VERSION_SUBMINOR >= 3 )
			VecSetOption( level->workRHS, VEC_IGNORE_NEGATIVE_INDICES );
#elif( PETSC_VERSION_MAJOR >= 3 )
			VecSetOption( level->workRHS, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE );
#endif
		}
	}
}

PetscErrorCode MatMult_myMat(Mat A,Vec x,Vec y)
{
  PetscFunctionBeginUser;
  PetscErrorCode ierr;
  model *m;
  ierr = MatShellGetContext(A, &m); CHKERRQ(ierr);

  PetscInt len_x, len_y;
  ierr = VecGetLocalSize(x, &len_x);  CHKERRQ(ierr);
  ierr = VecGetLocalSize(y, &len_y); CHKERRQ(ierr);
  //std::cout << "len_x=" << len_x << "len_y=" << len_y << std::endl;

  const PetscScalar *px;
  PetscScalar       *py, *pgx;

  VecScatter ctx;
  Vec vout;
  VecScatterCreateToAll(x,&ctx,&vout);
  // scatter as many times as you need
  VecScatterBegin(ctx,x,vout,INSERT_VALUES,SCATTER_FORWARD);
  VecScatterEnd(ctx,x,vout,INSERT_VALUES,SCATTER_FORWARD);

  ierr = VecGetArray(vout, &pgx); CHKERRQ(ierr);
  ierr = VecGetArrayRead(x, &px); CHKERRQ(ierr);
  ierr = VecGetArray(y, &py); CHKERRQ(ierr);

  PetscInt Istart, Iend;
  ierr = MatGetOwnershipRange(A, &Istart, &Iend); CHKERRQ(ierr);
  PetscInt A_m, n;
  ierr = MatGetSize(A, &A_m, &n); CHKERRQ(ierr);
  PetscScalar* pgy = new PetscScalar[n];
  rokko::heisenberg_hamiltonian::multiply(m->L, m->lattice, pgx, pgy);

  for(int j = 0; j < len_y; ++j) {
    py[j] = 0.;
  }

  int k=0;
  for(int i = Istart; i < Iend; ++i) {
    py[k] = pgy[i];
    ++k;
  }

  // destroy scatter context and local vector when no longer needed
  VecScatterDestroy(&ctx);
  //VecView(vout, PETSC_VIEWER_STDOUT_WORLD);
  VecDestroy(&vout);

  ierr = VecRestoreArrayRead(x,&px); CHKERRQ(ierr);
  ierr = VecRestoreArray(y,&py); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

/*@
   PFApplyVec - Applies the mathematical function to a vector

   Collective on PF

   Input Parameters:
+  pf - the function context
-  x - input vector (or PETSC_NULL for the vector (0,1, .... N-1)

   Output Parameter:
.  y - output vector

   Level: beginner

.keywords: PF, apply

.seealso: PFApply(), PFCreate(), PFDestroy(), PFSetType(), PFSet()
@*/
PetscErrorCode  PFApplyVec(PF pf,Vec x,Vec y)
{
  PetscErrorCode ierr;
  PetscInt       i,rstart,rend,n,p;
  PetscBool      nox = PETSC_FALSE;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pf,PF_CLASSID,1);
  PetscValidHeaderSpecific(y,VEC_CLASSID,3);
  if (x) {
    PetscValidHeaderSpecific(x,VEC_CLASSID,2);
    if (x == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"x and y must be different vectors");
  } else {
    PetscScalar *xx;
    PetscInt    lsize;

    ierr = VecGetLocalSize(y,&lsize);CHKERRQ(ierr);
    lsize = pf->dimin*lsize/pf->dimout;
    ierr = VecCreateMPI(((PetscObject)y)->comm,lsize,PETSC_DETERMINE,&x);CHKERRQ(ierr);
    nox  = PETSC_TRUE;
    ierr = VecGetOwnershipRange(x,&rstart,&rend);CHKERRQ(ierr);
    ierr = VecGetArray(x,&xx);CHKERRQ(ierr);
    for (i=rstart; i<rend; i++) {
      xx[i-rstart] = (PetscScalar)i;
    }
    ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);
  }

  ierr = VecGetLocalSize(x,&n);CHKERRQ(ierr);
  ierr = VecGetLocalSize(y,&p);CHKERRQ(ierr);
  if ((pf->dimin*(n/pf->dimin)) != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local input vector length %D not divisible by dimin %D of function",n,pf->dimin);
  if ((pf->dimout*(p/pf->dimout)) != p) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local output vector length %D not divisible by dimout %D of function",p,pf->dimout);
  if ((n/pf->dimin) != (p/pf->dimout)) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local vector lengths %D %D are wrong for dimin and dimout %D %D of function",n,p,pf->dimin,pf->dimout);

  if (pf->ops->applyvec) {
    ierr = (*pf->ops->applyvec)(pf->data,x,y);CHKERRQ(ierr);
  } else {
    PetscScalar *xx,*yy;

    ierr = VecGetLocalSize(x,&n);CHKERRQ(ierr);
    n    = n/pf->dimin;
    ierr = VecGetArray(x,&xx);CHKERRQ(ierr);
    ierr = VecGetArray(y,&yy);CHKERRQ(ierr);
    if (!pf->ops->apply) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No function has been provided for this PF");
    ierr = (*pf->ops->apply)(pf->data,n,xx,yy);CHKERRQ(ierr);
    ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);
    ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr);
  }
  if (nox) {
    ierr = VecDestroy(&x);CHKERRQ(ierr);
  } 
  PetscFunctionReturn(0);
}

 PetscErrorCode   restrictMatVecType2(Mat R, Vec f, Vec c) {
   TransferOpData *data;
   PetscFunctionBegin;
   iC(MatShellGetContext( R, (void **)&data));
   MPI_Comm comm = data->comm;
   Vec tmp = data->tmp;		
   PetscInt tmpSz;
   PetscInt fSz;
   iC(VecGetLocalSize(tmp,&tmpSz));
   iC(VecGetLocalSize(f,&fSz));
   scatterValues(f, tmp, fSz, tmpSz, data->sendSzR,
       data->sendOffR, data->recvSzR, data->recvOffR, comm);
   restrictMatVecType1(R, tmp, c);
   PetscFunctionReturn(0);
 }

    void TestInterleavedVecScatter()
    {
        // Vectors will be twice PROBLEM_SIZE, since this is to be used in bidomain code.
        const unsigned PROBLEM_SIZE = 10;
        DistributedVectorFactory factory(PROBLEM_SIZE);

        // Source vector = [-1 1 -2 2 ... -PROBLEM_SIZE/2+1 PROBLEM_SIZE/2+1]
        Vec interleaved_vec=factory.CreateVec(2);
        DistributedVector interleaved_dist_vec = factory.CreateDistributedVector(interleaved_vec);
        DistributedVector::Stripe first_variable(interleaved_dist_vec, 0);
        DistributedVector::Stripe second_variable(interleaved_dist_vec, 1);
        for (DistributedVector::Iterator index = interleaved_dist_vec.Begin();
             index!= interleaved_dist_vec.End();
             ++index)
        {
            first_variable[index]  = -1.0 * (index.Global+1);
            second_variable[index] =        (index.Global+1);
        }

        // Destination vectors. It is important they have a compatible parallel layout with the source vector, hence the 2nd param of CreateVec()
        PetscInt local_num_rows;
        VecGetLocalSize(interleaved_vec, &local_num_rows);
        Vec first_variable_vec = PetscTools::CreateVec(PROBLEM_SIZE, local_num_rows/2);
        Vec second_variable_vec = PetscTools::CreateVec(PROBLEM_SIZE, local_num_rows/2);

        // Setup and perform scatter operation
        VecScatter first_variable_context;
        VecScatter second_variable_context;
        PetscVecTools::SetupInterleavedVectorScatterGather(interleaved_vec, first_variable_context, second_variable_context);

#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 6) //PETSc 3.6 or later
        // When used in the preconditioner code within a KSP solve, the main (bidomain) vector will be locked.
        // Lock the vector so that modern PETSc (3.6) won't want to change it
        VecLockPush(interleaved_vec);
#endif

        PetscVecTools::DoInterleavedVecScatter(interleaved_vec, first_variable_context, first_variable_vec, second_variable_context, second_variable_vec);

#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 6) //PETSc 3.6 or later
        // Unlock the vector (for symmetry)
        VecLockPop(interleaved_vec);
#endif

        // Check destination vectors are [-1 -2 -3 ...] and [1 2 3 ...] respectively.
        DistributedVector dist_1st_var_vec = factory.CreateDistributedVector(first_variable_vec);
        DistributedVector dist_2nd_var_vec = factory.CreateDistributedVector(second_variable_vec);
        for (DistributedVector::Iterator index = dist_1st_var_vec.Begin();
             index!= dist_1st_var_vec.End();
             ++index)
        {
            TS_ASSERT_EQUALS(dist_1st_var_vec[index], -1.0 * (index.Global+1));
            TS_ASSERT_EQUALS(dist_2nd_var_vec[index],         index.Global+1);
        }

        PetscTools::Destroy(interleaved_vec);
        PetscTools::Destroy(first_variable_vec);
        PetscTools::Destroy(second_variable_vec);
        VecScatterDestroy(PETSC_DESTROY_PARAM(first_variable_context));
        VecScatterDestroy(PETSC_DESTROY_PARAM(second_variable_context));
    }

PetscErrorCode VecReorder(Vec r, Vec order, Vec or){
    PetscErrorCode err;
    PetscInt n_local, low, high;
    PetscInt *to_idx, *from_idx;
    PetscScalar *o_array;
    IS to, from;
    VecScatter scatter;
    int i;

    err = VecGetLocalSize(r, &n_local); CHKERRQ(err);
    err = PetscMalloc1(n_local, &to_idx); CHKERRQ(err);
    err = PetscMalloc1(n_local, &from_idx); CHKERRQ(err);
    err = VecGetOwnershipRange(r, &low, &high); CHKERRQ(err);
    err = VecGetArray(order, &o_array);
    for(i = 0; i < n_local; i++){
        to_idx[i] = (PetscInt) o_array[low + i];
        from_idx[i] = (PetscInt) low + i;
    }
    err = VecRestoreArray(order, &o_array);

    err = ISCreateGeneral(PETSC_COMM_SELF, n_local, from_idx, PETSC_OWN_POINTER, &to); CHKERRQ(err);
    err = ISCreateGeneral(PETSC_COMM_SELF, n_local, from_idx, PETSC_OWN_POINTER, &from); CHKERRQ(err);

    err = VecScatterCreate(r, from, or, to, &scatter); CHKERRQ(err);
    err = VecScatterBegin(scatter, r, or, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(err);
    err = VecScatterEnd(scatter, r, or, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(err);

    err = PetscFree(to_idx); CHKERRQ(err);
    err = PetscFree(from_idx); CHKERRQ(err);

    return err;
}

int LargeVecCreate(Vec *x, PetscInt nvec, Vec xvec[]){

 PetscInt       size,rank;
 PetscErrorCode ierr;
 PetscInt       lsize,mylsize;
 PetscInt       i;
 PetscScalar    *ptx;
 PetscInt       pstart[nvec], pend[nvec];


 MPI_Comm_size(PETSC_COMM_WORLD,&size);
 MPI_Comm_rank(PETSC_COMM_WORLD,&rank);


 LargeVecGetOwnershipRange(x,nvec,pstart,pend);

 ierr =  VecGetArray(*x,&ptx);CHKERRQ(ierr);
 ierr =  VecGetLocalSize(*x,&lsize);CHKERRQ(ierr);

 for (i=0;i<nvec;i++){
   if(rank>=pstart[i]&&rank<pend[i]){mylsize = lsize;}else{mylsize=0;}
   ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,mylsize,PETSC_DETERMINE,ptx,xvec+i);CHKERRQ(ierr);
 }
 

 return 0;
}

dErr dFSGetBoundingBox(dFS fs,dReal bbox[3][2])
{
    dErr err;
    Vec X;
    const dScalar *x;
    dInt n;

    dFunctionBegin;
    for (dInt i=0; i<3; i++) {
        bbox[i][0] = PETSC_MAX_REAL;
        bbox[i][1] = PETSC_MIN_REAL;
    }
    err = dFSGetGeometryVectorExpanded(fs,&X);
    dCHK(err);
    err = VecGetLocalSize(X,&n);
    dCHK(err);
    err = VecGetArrayRead(X,&x);
    dCHK(err);
    for (dInt i=0; i<n; i++) {
        dInt j = i%3;
        bbox[j][0] = dMin(bbox[j][0],x[i]);
        bbox[j][1] = dMax(bbox[j][1],x[i]);
    }
    err = VecRestoreArrayRead(X,&x);
    dCHK(err);
    dFunctionReturn(0);
}

PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
{
  PetscErrorCode    ierr;
  SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
  PetscInt          *indices;
  PetscInt          i,n,rstart,rend;
  SNESLineSearch    linesearch;

  PetscFunctionBegin;
  ierr = SNESSetUp_VI(snes);CHKERRQ(ierr);

  /* Set up previous active index set for the first snes solve
   vi->IS_inact_prev = 0,1,2,....N */

  ierr = VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);CHKERRQ(ierr);
  ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr);
  ierr = PetscMalloc(n*sizeof(PetscInt),&indices);CHKERRQ(ierr);
  for (i=0;i < n; i++) indices[i] = rstart + i;
  ierr = ISCreateGeneral(((PetscObject)snes)->comm,n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);CHKERRQ(ierr);

  /* set the line search functions */
  if (!snes->linesearch) {
    ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr);
    ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
{
  Mat_Shell *shell = (Mat_Shell*)A->data;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (shell->dshift) {          /* get arrays because there is no VecPointwiseMultAdd() */
    PetscInt          i,m;
    const PetscScalar *x,*d;
    PetscScalar       *y;
    ierr = VecGetLocalSize(X,&m);CHKERRQ(ierr);
    ierr = VecGetArrayRead(shell->dshift,&d);CHKERRQ(ierr);
    ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
    ierr = VecGetArray(Y,&y);CHKERRQ(ierr);
    for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
    ierr = VecRestoreArrayRead(shell->dshift,&d);CHKERRQ(ierr);
    ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
    ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr);
  } else if (PetscAbsScalar(shell->vshift) != 0) {
    ierr = VecAXPBY(Y,shell->vshift,shell->vscale,X);CHKERRQ(ierr);
  } else if (shell->vscale != 1.0) {
    ierr = VecScale(Y,shell->vscale);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

PetscErrorCode Initial(Vec global,void *ctx)
{
  Data           *data = (Data*)ctx;
  PetscInt       m,row,col;
  PetscReal      x,y,dx,dy;
  PetscScalar    *localptr;
  PetscInt       i,mybase,myend,locsize;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  /* make the local  copies of parameters */
  m  = data->m;
  dx = data->dx;
  dy = data->dy;

  /* determine starting point of each processor */
  ierr = VecGetOwnershipRange(global,&mybase,&myend);CHKERRQ(ierr);
  ierr = VecGetLocalSize(global,&locsize);CHKERRQ(ierr);

  /* Initialize the array */
  ierr = VecGetArray(global,&localptr);CHKERRQ(ierr);

  for (i=0; i<locsize; i++) {
    row         = 1+(mybase+i)-((mybase+i)/m)*m;
    col         = (mybase+i)/m+1;
    x           = dx*row;
    y           = dy*col;
    localptr[i] = f_ini(x,y);
  }

  ierr = VecRestoreArray(global,&localptr);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}