C++ VecGetOwnershipRange函数代码示例

void Hdf5DataReader::GetVariableOverNodes(Vec data,
                                          const std::string& rVariableName,
                                          unsigned timestep)
{
    if (!mIsDataComplete)
    {
        EXCEPTION("You can only get a vector for complete data");
    }
    if (!mIsUnlimitedDimensionSet && timestep!=0)
    {
        EXCEPTION("The dataset '" << mDatasetName << "' does not contain time dependent data");
    }

    std::map<std::string, unsigned>::iterator col_iter = mVariableToColumnIndex.find(rVariableName);
    if (col_iter == mVariableToColumnIndex.end())
    {
        EXCEPTION("The dataset '" << mDatasetName << "' does not contain data for variable " << rVariableName);
    }
    unsigned column_index = (*col_iter).second;

    // Check for valid timestep
    if (timestep >= mNumberTimesteps)
    {
        EXCEPTION("The dataset '" << mDatasetName << "' does not contain data for timestep number " << timestep);
    }

    int lo, hi, size;
    VecGetSize(data, &size);
    if ((unsigned)size != mDatasetDims[1])
    {
        EXCEPTION("Could not read data because Vec is the wrong size");
    }
    // Get range owned by each processor
    VecGetOwnershipRange(data, &lo, &hi);

    if (hi > lo) // i.e. we own some...
    {
        // Define a dataset in memory for this process
        hsize_t v_size[1] = {(unsigned)(hi-lo)};
        hid_t memspace = H5Screate_simple(1, v_size, NULL);

        // Select hyperslab in the file.
        hsize_t offset[3] = {timestep, (unsigned)(lo), column_index};
        hsize_t count[3]  = {1, (unsigned)(hi-lo), 1};
        hid_t hyperslab_space = H5Dget_space(mVariablesDatasetId);
        H5Sselect_hyperslab(hyperslab_space, H5S_SELECT_SET, offset, NULL, count, NULL);

        double* p_petsc_vector;
        VecGetArray(data, &p_petsc_vector);

        herr_t err = H5Dread(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, H5P_DEFAULT, p_petsc_vector);
        UNUSED_OPT(err);
        assert(err==0);

        VecRestoreArray(data, &p_petsc_vector);

        H5Sclose(hyperslab_space);
        H5Sclose(memspace);
    }
}

/** 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);
}

/*
   ExactSolution - Computes the exact solution at a given time.

   Input Parameters:
   t - current time
   solution - vector in which exact solution will be computed
   appctx - user-defined application context

   Output Parameter:
   solution - vector with the newly computed exact solution
*/
PetscErrorCode ExactSolution(PetscReal t,Vec solution,AppCtx *appctx)
{
  PetscScalar    *s_localptr,h = appctx->h,x;
  PetscInt       i,mybase,myend;
  PetscErrorCode ierr;

  /*
     Determine starting and ending points of each processor's
     range of grid values
  */
  ierr = VecGetOwnershipRange(solution,&mybase,&myend);CHKERRQ(ierr);

  /*
     Get a pointer to vector data.
  */
  ierr = VecGetArray(solution,&s_localptr);CHKERRQ(ierr);

  /*
     Simply write the solution directly into the array locations.
     Alternatively, we could use VecSetValues() or VecSetValuesLocal().
  */
  for (i=mybase; i<myend; i++) {
    x = h*(PetscReal)i;
    s_localptr[i-mybase] = (t + 1.0)*(1.0 + x*x);
  }

  /*
     Restore vector
  */
  ierr = VecRestoreArray(solution,&s_localptr);CHKERRQ(ierr);
  return 0;
}

inline int PetscVector::last_local_index () const {
  assert (this->initialized());
  int ierr=0, petsc_first=0, petsc_last=0;
  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
  CHKERRABORT(MPI_COMM_WORLD,ierr);
  return static_cast<int>(petsc_last);
}

SparseGp *
SparseGp_create (AcfrKern *kern, int rank, gsl_matrix *trainObs, 
                 gsl_vector *trainLabels, gsl_matrix *testObs, 
                 gsl_vector *testLabels)
{
    PetscErrorCode ierr;        /* not used here - not returning int */
    (void) ierr;

    SparseGp *gp = calloc (1, sizeof (*gp));

    gp->kern = kern;

    gp->trainObs = trainObs;
    gp->trainLabels = trainLabels;
    gp->testObs = testObs;
    gp->testLabels = testLabels;

    PetscInt N = gp->trainLabels->size;
    ierr = VecCreate (PETSC_COMM_WORLD, &(gp->_trainLabels)); /* CHKERRQ (ierr); */
    ierr = VecSetSizes (gp->_trainLabels, PETSC_DECIDE, N); /* CHKERRQ (ierr); */
    ierr = VecSetFromOptions (gp->_trainLabels); /* CHKERRQ (ierr); */

    ierr = VecGetOwnershipRange (gp->_trainLabels, &(gp->rstart), &(gp->rend)); /* CHKERRQ (ierr); */
    ierr = VecGetLocalSize (gp->_trainLabels, &(gp->nlocal)); /* CHKERRQ (ierr); */
    petsc_util_fillVec (gp->trainLabels, &gp->_trainLabels, gp->rstart, gp->rend);

    gp->rank = rank;

    return gp;
}

void AbstractContinuumMechanicsSolver<DIM>::AddIdentityBlockForDummyPressureVariables(ApplyDirichletBcsType type)
{
    assert(mCompressibilityType==INCOMPRESSIBLE);

    int lo, hi;
    VecGetOwnershipRange(mResidualVector, &lo, &hi);

    for (unsigned i=0; i<mrQuadMesh.GetNumNodes(); i++)
    {
        if (mrQuadMesh.GetNode(i)->IsInternal())
        {
            unsigned row = (DIM+1)*i + DIM; // DIM+1 is the problem dimension
            if (lo <= (int)row && (int)row < hi)
            {
                if (type!=LINEAR_PROBLEM)
                {
                    PetscVecTools::SetElement(mResidualVector, row, mCurrentSolution[row]-0.0);
                }
                if (type!=NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY) // ie doing a whole linear system
                {
                    double rhs_vector_val = type==LINEAR_PROBLEM ? 0.0 : mCurrentSolution[row]-0.0;
                    PetscVecTools::SetElement(mLinearSystemRhsVector, row, rhs_vector_val);
                    // This assumes the row is already zero, which is should be..
                    PetscMatTools::SetElement(mSystemLhsMatrix, row, row, 1.0);
                    PetscMatTools::SetElement(mPreconditionMatrix, row, row, 1.0);
                }
            }
        }
    }
}

inline PetscErrorCode ScatIEMatMult(Mat mat,Vec xx,Vec yy)
{
  PetscInt       n,s,t;
  PetscScalar    tmp,v;
  PetscInt       xstart,xend;

  PetscFunctionBegin;
  VecGetOwnershipRange(xx,&xstart,&xend);
  VecGetSize(yy,&n);
  VecZeroEntries(yy);

  for(s = 0; s <n ; s++)
  {
      tmp = 0;
      for(t = xstart; t < xend; t++)
      {
          VecGetValues(xx,1,&t,&v);
          tmp += ScatIE.CoefMatFast(s,t)*v;
      }
      VecSetValues(yy,1,&s,&tmp,ADD_VALUES);
  }
  VecAssemblyBegin(yy);
  VecAssemblyEnd(yy);

  PetscFunctionReturn(0);
}

/* Histogram a PETSC vector 
 *
 * x is the vector
 * nbins -- number of bins
 * xmin, xmax -- histogram xmin, xmax -- assume uniform bins
 * hh -- output vector -- assumed to be defined.
 */
void VecHist(const Vec& x, int nbins, double xmin, double xmax, vector<double>& hh) {
  gsl_histogram *h1; 
  double *_x, x1;
  PetscInt lo, hi;
  vector<double> tmp(nbins);

  // Set up the histogram struct
  h1 = gsl_histogram_alloc(nbins);
  gsl_histogram_set_ranges_uniform(h1, xmin, xmax);

  // Get the array
  VecGetOwnershipRange(x, &lo, &hi);
  hi -= lo;
  VecGetArray(x, &_x);
  for (PetscInt ii=0; ii < hi; ++ii) {
    x1 = _x[ii];
    if (x1 < xmin) x1 = xmin;
    if (x1 >= xmax) x1 = xmax - 1.e-10;
    gsl_histogram_increment(h1, x1);
  }
  VecRestoreArray(x, &_x);
  
  // Fill the temporary output vector
  for (int ii =0; ii<nbins; ++ii) 
    tmp[ii] = gsl_histogram_get(h1, ii);

  // MPI Allreduce
  MPI_Allreduce(&tmp[0], &hh[0], nbins, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); 

  // Clean up
  gsl_histogram_free(h1);
}

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  PetscInt       n = 5,N,low,high,iglobal,i;
  PetscMPIInt    size,rank;
  PetscScalar    value,zero = 0.0;
  Vec            x,y;
  IS             is1,is2;
  VecScatter     ctx;

  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr); 
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  /* create two vectors */
  N = size*n;
  ierr = VecCreate(PETSC_COMM_WORLD,&y);CHKERRQ(ierr);
  ierr = VecSetSizes(y,PETSC_DECIDE,N);CHKERRQ(ierr);
  ierr = VecSetFromOptions(y);CHKERRQ(ierr);
  if (!rank) {
    ierr = VecCreateSeq(PETSC_COMM_SELF,N,&x);CHKERRQ(ierr);
  } else {
    ierr = VecCreateSeq(PETSC_COMM_SELF,0,&x);CHKERRQ(ierr);
  }

  /* create two index sets */
  if (!rank) {
    ierr = ISCreateStride(PETSC_COMM_SELF,N,0,1,&is1);CHKERRQ(ierr);
    ierr = ISCreateStride(PETSC_COMM_SELF,N,0,1,&is2);CHKERRQ(ierr);
  } else {
    ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is1);CHKERRQ(ierr);
    ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is2);CHKERRQ(ierr);
  }

  ierr = VecSet(x,zero);CHKERRQ(ierr);
  ierr = VecGetOwnershipRange(y,&low,&high);CHKERRQ(ierr);
  for (i=0; i<n; i++) {
    iglobal = i + low; value = (PetscScalar) (i + 10*rank);
    ierr = VecSetValues(y,1,&iglobal,&value,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(y);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(y);CHKERRQ(ierr);
  ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  ierr = VecScatterCreate(y,is2,x,is1,&ctx);CHKERRQ(ierr);
  ierr = VecScatterBegin(ctx,y,x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd(ctx,y,x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterDestroy(&ctx);CHKERRQ(ierr);
  
  if (!rank) 
    {printf("----\n"); ierr = VecView(x,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}

  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&y);CHKERRQ(ierr);
  ierr = ISDestroy(&is1);CHKERRQ(ierr);
  ierr = ISDestroy(&is2);CHKERRQ(ierr);

  ierr = PetscFinalize();
  return 0;
}

PetscErrorCode VecMin_MPI(Vec xin,PetscInt *idx,PetscReal *z)
{
  PetscErrorCode ierr;
  PetscReal      work;

  PetscFunctionBegin;
  /* Find the local Min */
  ierr = VecMin_Seq(xin,idx,&work);CHKERRQ(ierr);

  /* Find the global Min */
  if (!idx) {
    ierr = MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)xin));CHKERRQ(ierr);
  } else {
    PetscReal work2[2],z2[2];
    PetscInt  rstart;

    ierr = VecGetOwnershipRange(xin,&rstart,NULL);CHKERRQ(ierr);
    work2[0] = work;
    work2[1] = *idx + rstart;
    ierr = MPIU_Allreduce(work2,z2,2,MPIU_REAL,MPIU_MININDEX_OP,PetscObjectComm((PetscObject)xin));CHKERRQ(ierr);
    *z   = z2[0];
    *idx = (PetscInt)z2[1];
  }
  PetscFunctionReturn(0);
}

PetscErrorCode GetRadiusVec(Vec vecRad, int Nr, int Nz, double hr, int m)
{
   PetscErrorCode ierr;
   int i, j, ir, ic, ns, ne;
   double value=1.0;
   ierr = VecGetOwnershipRange(vecRad,&ns,&ne); CHKERRQ(ierr);
   for(i=ns;i<ne; i++)
     {
       j = i;
       ir = (j /= Nz) % Nr;
       ic = (j /= Nr) % 3;
 
       if (ic==0) value = (ir + 0.5)*hr;
       if (ic==1) value = (ir + 0.0)*hr;
       if (ic==2) value = (ir + 0.0)*hr + (ir==0)*(m==0)*hr/8; 
        
       VecSetValue(vecRad, i, value, INSERT_VALUES);

       }

   ierr = VecAssemblyBegin(vecRad); CHKERRQ(ierr);
   ierr = VecAssemblyEnd(vecRad); CHKERRQ(ierr);

   PetscFunctionReturn(0);
}

PetscErrorCode GetUnitVec(Vec ej, int pol, int N)
{
   PetscErrorCode ierr;
   int i, j, ns, ne, ic;
   int Nc=3;
   int Nxyz=N/6;

   ierr = VecGetOwnershipRange(ej,&ns,&ne); CHKERRQ(ierr);

   for(i=ns; i<ne; i++)
     {
       j=i;
       ic = (j /= Nxyz) % Nc;
       
       if (ic==pol)
	 VecSetValue(ej,i,1.0,INSERT_VALUES);
       else
	 VecSetValue(ej,i,0.0,INSERT_VALUES);
     }

   ierr = VecAssemblyBegin(ej); CHKERRQ(ierr);
   ierr = VecAssemblyEnd(ej); CHKERRQ(ierr);
   
   PetscFunctionReturn(0);
  
}

/*@
   DMDANaturalAllToGlobalCreate - Creates a scatter context that maps from a copy
     of the entire vector on each processor to its local part in the global vector.

   Collective on DMDA

   Input Parameter:
.  da - the distributed array context

   Output Parameter:
.  scatter - the scatter context

   Level: advanced

.keywords: distributed array, global to local, begin, coarse problem

.seealso: DMDAGlobalToNaturalEnd(), DMLocalToGlobalBegin(), DMDACreate2d(),
          DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMDACreateNaturalVector()
@*/
PetscErrorCode  DMDANaturalAllToGlobalCreate(DM da,VecScatter *scatter)
{
  PetscErrorCode ierr;
  DM_DA          *dd = (DM_DA*)da->data;
  PetscInt       M,m = dd->Nlocal,start;
  IS             from,to;
  Vec            tmplocal,global;
  AO             ao;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(da,DM_CLASSID,1);
  PetscValidPointer(scatter,2);
  ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr);

  /* create the scatter context */
  ierr = MPI_Allreduce(&m,&M,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)da));CHKERRQ(ierr);
  ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)da),dd->w,m,PETSC_DETERMINE,0,&global);CHKERRQ(ierr);
  ierr = VecGetOwnershipRange(global,&start,NULL);CHKERRQ(ierr);
  ierr = ISCreateStride(PetscObjectComm((PetscObject)da),m,start,1,&from);CHKERRQ(ierr);
  ierr = AOPetscToApplicationIS(ao,from);CHKERRQ(ierr);
  ierr = ISCreateStride(PetscObjectComm((PetscObject)da),m,start,1,&to);CHKERRQ(ierr);
  ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dd->w,M,0,&tmplocal);CHKERRQ(ierr);
  ierr = VecScatterCreate(tmplocal,from,global,to,scatter);CHKERRQ(ierr);
  ierr = VecDestroy(&tmplocal);CHKERRQ(ierr);
  ierr = VecDestroy(&global);CHKERRQ(ierr);
  ierr = ISDestroy(&from);CHKERRQ(ierr);
  ierr = ISDestroy(&to);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

void SetJacobian(Geometry geo, Mat J, Vec v, int jc, int jr, int jh){
// use jc = -1 for last 2 columns, and jc = -2 for blocks (all jc)

	AssembleVec(v);

	int ns, ne, i;
	VecGetOwnershipRange(v, &ns, &ne);
	const double *vals;
	VecGetArrayRead(v, &vals);

	for(i=ns; i<ne; i++){
		if( Last2(geo, i) || vals[i-ns] == 0.0) continue;

		int col, offset = jh*(Nxyzcr(geo)+2),
		ij = i%(Nxyzcr(geo)+2);
	
		Point p;
		CreatePoint_i(&p, ij, &geo->gN);
	
		if(jc == -1) //columns
			col = Nxyzcr(geo)+jr;
		else if(jc == -2) //blocks
			col = jr*Nxyzc(geo) + xyzc(&p);
		else // tensors
			col = jr*Nxyzc(geo) + jc*Nxyz(geo) + xyz(&p);
	
		MatSetValue(J, i, col+offset, vals[i-ns], ADD_VALUES);
	
	}
	VecRestoreArrayRead(v, &vals);
}

void CollectVec(Geometry geo, Vec vN, Vec vM){
	VecSet(vM, 0.0);

	Vecfun f;
	CreateVecfun(&f, geo->vf);

	const double *vals;
	VecGetArrayRead(vN, &vals);
	int ns, ne;
	VecGetOwnershipRange(vN, &ns, &ne);
	if( ne > Nxyzcr(geo)-2) ne = Nxyzcr(geo)-2;

	int i;
	for(i=ns; i<ne; i++){
		if( valr(&f, i) == 0.0) continue;
		// skip if gain profile zero here

		Point p;
		CreatePoint_i(&p, i, &geo->gN);
		if( projectmedium(&p, &geo->gM, geo->LowerPML) )
			VecSetValue(vM, xyz(&p), vals[i-ns], ADD_VALUES);
	}

	VecRestoreArrayRead(vN, &vals);
	DestroyVecfun(&f);
	AssembleVec(vM);
}