C++ VecAssemblyEnd函数代码示例

void PETSc::Solve(void)
{
        
  //start_clock("Before Assemble matrix and vector");
	ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
  	ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
  	
  	ierr = VecAssemblyBegin(x);
  	ierr = VecAssemblyEnd(x);
  	
  	ierr = VecAssemblyBegin(b);
  	ierr = VecAssemblyEnd(b);
	//stop_clock("After Assembly matrix and vector");


        KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
        KSPSetType(ksp,KSPBCGSL);
	KSPBCGSLSetEll(ksp,2);

	//KSPGetPC(ksp, &pc);
	//PCSetType(pc, PCJACOBI);

        KSPSetFromOptions(ksp);
        KSPSetUp(ksp);

	//start_clock("Before KSPSolve");
        KSPSolve(ksp,b,x);
	//stop_clock("After KSPSolve");
}

PetscErrorCode test_vec_ops( void )
{
  Vec            X, a,b;
  Vec            c,d,e,f;
  PetscScalar    val;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscPrintf( PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME);

  ierr = VecCreate( PETSC_COMM_WORLD, &X );CHKERRQ(ierr);
  ierr = VecSetSizes( X, 2, 2 );CHKERRQ(ierr);
  ierr = VecSetType( X, VECNEST );CHKERRQ(ierr);

  ierr = VecCreate( PETSC_COMM_WORLD, &a );CHKERRQ(ierr);
  ierr = VecSetSizes( a, 2, 2 );CHKERRQ(ierr);
  ierr = VecSetType( a, VECNEST );CHKERRQ(ierr);

  ierr = VecCreate( PETSC_COMM_WORLD, &b );CHKERRQ(ierr);
  ierr = VecSetSizes( b, 2, 2 );CHKERRQ(ierr);
  ierr = VecSetType( b, VECNEST );CHKERRQ(ierr);

  /* assemble X */
  ierr = VecNestSetSubVec( X, 0, a );CHKERRQ(ierr);
  ierr = VecNestSetSubVec( X, 1, b );CHKERRQ(ierr);
  ierr = VecAssemblyBegin(X);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(X);CHKERRQ(ierr);

  ierr = VecCreate( PETSC_COMM_WORLD, &c );CHKERRQ(ierr);
  ierr = VecSetSizes( c, 3, 3 );CHKERRQ(ierr);
  ierr = VecSetType( c, VECSEQ );CHKERRQ(ierr);
  ierr = VecDuplicate( c, &d );CHKERRQ(ierr);
  ierr = VecDuplicate( c, &e );CHKERRQ(ierr);
  ierr = VecDuplicate( c, &f );CHKERRQ(ierr);

  ierr = VecSet( c, 1.0 );CHKERRQ(ierr);
  ierr = VecSet( d, 2.0 );CHKERRQ(ierr);
  ierr = VecSet( e, 3.0 );CHKERRQ(ierr);
  ierr = VecSet( f, 4.0 );CHKERRQ(ierr);

  /* assemble a */
  ierr = VecNestSetSubVec( a, 0, c );CHKERRQ(ierr);
  ierr = VecNestSetSubVec( a, 1, d );CHKERRQ(ierr);
  ierr = VecAssemblyBegin(a);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(a);CHKERRQ(ierr);

  /* assemble b */
  ierr = VecNestSetSubVec( b, 0, e );CHKERRQ(ierr);
  ierr = VecNestSetSubVec( b, 1, f );CHKERRQ(ierr);
  ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

  //PetscPrintf( PETSC_COMM_WORLD, "X \n");
  //VecView( X, PETSC_VIEWER_STDOUT_WORLD );

  ierr = VecDot( X,X, &val );CHKERRQ(ierr);
  PetscPrintf( PETSC_COMM_WORLD, "X.X = %f \n", val );

  PetscFunctionReturn(0);
}

PetscErrorCode KSPBuildPressure_CB_Nullspace_BSSCR(KSP ksp)
{
    KSP_BSSCR        *bsscr = (KSP_BSSCR *)ksp->data;
    FeEquationNumber *eq_num = bsscr->solver->st_sle->pSolnVec->feVariable->eqNum;
    FeMesh           *feMesh = bsscr->solver->st_sle->pSolnVec->feVariable->feMesh; /* is the pressure mesh */
    unsigned          ijk[3];
    Vec               t, v;
    int numLocalNodes, globalNodeNumber, i, j, eq;
    MatStokesBlockScaling BA = bsscr->BA;
    PetscErrorCode ierr;
    Mat            Amat,Pmat, G;
    MatStructure   pflag;

    PetscFunctionBegin;
    /* get G matrix from Amat matrix operator on ksp */
    ierr=Stg_PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
    MatNestGetSubMat( Amat, 0,1, &G );/* G should always exist */
    /* now create Vecs t and v to match size of G: i.e. pressure */ /* NOTE: not using "h" vector from ksp->vec_rhs because this part of the block vector doesn't always exist */
    MatGetVecs( G, &t, PETSC_NULL );/* t and v are destroyed in KSPDestroy_BSSCR */
    MatGetVecs( G, &v, PETSC_NULL );/* t and v such that can do G*t */

    numLocalNodes = Mesh_GetLocalSize( feMesh,  MT_VERTEX); /* number of nodes on current proc not counting any shadow nodes */
    for(j=0;j<numLocalNodes;j++){
	i = globalNodeNumber = Mesh_DomainToGlobal( feMesh, MT_VERTEX, j);
	RegularMeshUtils_Element_1DTo3D(feMesh, i, ijk);
	eq = eq_num->destinationArray[j][0];/* get global equation number -- 2nd arg is always 0 because pressure has only one dof */
	if(eq != -1){
	    if( (ijk[0]+ijk[1]+ijk[2])%2 ==0 ){
		VecSetValue(t,eq,1.0,INSERT_VALUES);
	    }
	    else{
		VecSetValue(v,eq,1.0,INSERT_VALUES);	    }}}

    VecAssemblyBegin( t );
    VecAssemblyEnd( t );
    VecAssemblyBegin( v );
    VecAssemblyEnd( v );

    /* Scaling the null vectors here because it easier at the moment *//* maybe should do this in the original scaling function */
    if( BA->scaling_exists == PETSC_TRUE ){
	Vec R2;
	/* Get the scalings out the block mat data */
	VecNestGetSubVec( BA->Rz, 1, &R2 );
	VecPointwiseDivide( t, t, R2); /* x <- x * 1/R2 */
	VecPointwiseDivide( v, v, R2);
    }

    bsscr_writeVec( t, "t", "Writing t vector");
    bsscr_writeVec( v, "v", "Writing v vector");

    bsscr->t=t;
    bsscr->v=v;

    PetscFunctionReturn(0);
}

int main(int argc,char **args)
{
  const PetscScalar xvals[] = {11,13},yvals[] = {17,19};
  const PetscInt    inds[]  = {0,1};
  PetscScalar       avals[] = {2,3,5,7};
  Mat               S1,S2;
  Vec               X,Y;
  User              user;
  PetscErrorCode    ierr;

  ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;

  ierr = PetscNew(&user);CHKERRQ(ierr);
  ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD,2,2,2,NULL,&user->A);CHKERRQ(ierr);
  ierr = MatSetUp(user->A);CHKERRQ(ierr);
  ierr = MatSetValues(user->A,2,inds,2,inds,avals,INSERT_VALUES);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(user->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(user->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_WORLD,2,&X);CHKERRQ(ierr);
  ierr = VecSetValues(X,2,inds,xvals,INSERT_VALUES);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(X);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(X);CHKERRQ(ierr);
  ierr = VecDuplicate(X,&Y);CHKERRQ(ierr);
  ierr = VecSetValues(Y,2,inds,yvals,INSERT_VALUES);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(Y);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(Y);CHKERRQ(ierr);

  ierr = MatCreateShell(PETSC_COMM_WORLD,2,2,2,2,user,&S1);CHKERRQ(ierr);
  ierr = MatSetUp(S1);CHKERRQ(ierr);
  ierr = MatShellSetOperation(S1,MATOP_MULT,(void (*)(void))MatMult_User);CHKERRQ(ierr);
  ierr = MatShellSetOperation(S1,MATOP_COPY,(void (*)(void))MatCopy_User);CHKERRQ(ierr);
  ierr = MatShellSetOperation(S1,MATOP_DESTROY,(void (*)(void))MatDestroy_User);CHKERRQ(ierr);
  ierr = MatCreateShell(PETSC_COMM_WORLD,2,2,2,2,NULL,&S2);CHKERRQ(ierr);
  ierr = MatSetUp(S2);CHKERRQ(ierr);
  ierr = MatShellSetOperation(S2,MATOP_MULT,(void (*)(void))MatMult_User);CHKERRQ(ierr);
  ierr = MatShellSetOperation(S2,MATOP_COPY,(void (*)(void))MatCopy_User);CHKERRQ(ierr);
  ierr = MatShellSetOperation(S2,MATOP_DESTROY,(void (*)(void))MatDestroy_User);CHKERRQ(ierr);

  ierr = MatScale(S1,31);CHKERRQ(ierr);
  ierr = MatShift(S1,37);CHKERRQ(ierr);
  ierr = MatDiagonalScale(S1,X,Y);CHKERRQ(ierr);
  ierr = MatCopy(S1,S2,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = MatMult(S1,X,Y);CHKERRQ(ierr);
  ierr = VecView(Y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = MatMult(S2,X,Y);CHKERRQ(ierr);
  ierr = VecView(Y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  ierr = MatDestroy(&S1);CHKERRQ(ierr);
  ierr = MatDestroy(&S2);CHKERRQ(ierr);
  ierr = VecDestroy(&X);CHKERRQ(ierr);
  ierr = VecDestroy(&Y);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  PetscMPIInt    size;
  PetscInt       n = 9,bs = 3,indices[2],i;
  PetscScalar    values[6];
  Vec            x;

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

  if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"Must be run with one processor");

  /* create vector */
  ierr = VecCreate(PETSC_COMM_SELF,&x);CHKERRQ(ierr);
  ierr = VecSetSizes(x,n,n);CHKERRQ(ierr);
  ierr = VecSetBlockSize(x,bs);CHKERRQ(ierr);
  ierr = VecSetType(x,VECSEQ);CHKERRQ(ierr);

  for (i=0; i<6; i++) values[i] = 4.0*i;
  indices[0] = 0;
  indices[1] = 2;

  ierr = VecSetValuesBlocked(x,2,indices,values,INSERT_VALUES);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(x);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(x);CHKERRQ(ierr);

  /*
      Resulting vector should be 0 4 8  0 0 0 12 16 20
  */
  ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  /* test insertion with negative indices */
  ierr = VecSetOption(x,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr);
  for (i=0; i<6; i++) values[i] = -4.0*i;
  indices[0] = -1;
  indices[1] = 2;

  ierr = VecSetValuesBlocked(x,2,indices,values,ADD_VALUES);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(x);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(x);CHKERRQ(ierr);

  /*
      Resulting vector should be 0 4 8  0 0 0 0 0 0
  */
  ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  ierr = VecDestroy(&x);CHKERRQ(ierr);

  ierr = PetscFinalize();
  return ierr;
}

void PkStruct::finalize() {

  VecAssemblyBegin(pkvec); VecAssemblyEnd(pkvec);
  VecAssemblyBegin(nmodes); VecAssemblyEnd(nmodes);
  VecAssemblyBegin(kvec); VecAssemblyEnd(kvec);

  //Normalize and return 
  VecShift(nmodes, 1.e-20);
  VecPointwiseDivide(pkvec, nmodes);
  VecPointwiseDivide(kvec, nmodes);

  VecGetArray(kvec, &_kvec);
  VecGetArray(pkvec, &_pkvec);
  VecGetArray(nmodes, &_nmodes);
}

/*@
    KSPComputeExplicitOperator - Computes the explicit preconditioned operator.

    Collective on KSP

    Input Parameter:
.   ksp - the Krylov subspace context

    Output Parameter:
.   mat - the explict preconditioned operator

    Notes:
    This computation is done by applying the operators to columns of the
    identity matrix.

    Currently, this routine uses a dense matrix format when 1 processor
    is used and a sparse format otherwise.  This routine is costly in general,
    and is recommended for use only with relatively small systems.

    Level: advanced

.keywords: KSP, compute, explicit, operator

.seealso: KSPComputeEigenvaluesExplicitly(), PCComputeExplicitOperator()
@*/
PetscErrorCode  KSPComputeExplicitOperator(KSP ksp,Mat *mat)
{
  Vec            in,out;
  PetscErrorCode ierr;
  PetscMPIInt    size;
  PetscInt       i,M,m,*rows,start,end;
  Mat            A;
  MPI_Comm       comm;
  PetscScalar    *array,one = 1.0;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
  PetscValidPointer(mat,2);
  comm = ((PetscObject)ksp)->comm;

  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);

  ierr = VecDuplicate(ksp->vec_sol,&in);CHKERRQ(ierr);
  ierr = VecDuplicate(ksp->vec_sol,&out);CHKERRQ(ierr);
  ierr = VecGetSize(in,&M);CHKERRQ(ierr);
  ierr = VecGetLocalSize(in,&m);CHKERRQ(ierr);
  ierr = VecGetOwnershipRange(in,&start,&end);CHKERRQ(ierr);
  ierr = PetscMalloc(m*sizeof(PetscInt),&rows);CHKERRQ(ierr);
  for (i=0; i<m; i++) {rows[i] = start + i;}

  ierr = MatCreate(comm,mat);CHKERRQ(ierr);
  ierr = MatSetSizes(*mat,m,m,M,M);CHKERRQ(ierr);
  if (size == 1) {
    ierr = MatSetType(*mat,MATSEQDENSE);CHKERRQ(ierr);
    ierr = MatSeqDenseSetPreallocation(*mat,PETSC_NULL);CHKERRQ(ierr);
  } else {
    ierr = MatSetType(*mat,MATMPIAIJ);CHKERRQ(ierr);
    ierr = MatMPIAIJSetPreallocation(*mat,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr);
  }
  ierr = MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
  if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
  ierr = PCGetOperators(ksp->pc,&A,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);

  for (i=0; i<M; i++) {

    ierr = VecSet(in,0.0);CHKERRQ(ierr);
    ierr = VecSetValues(in,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(in);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(in);CHKERRQ(ierr);

    ierr = KSP_MatMult(ksp,A,in,out);CHKERRQ(ierr);
    ierr = KSP_PCApply(ksp,out,in);CHKERRQ(ierr);

    ierr = VecGetArray(in,&array);CHKERRQ(ierr);
    ierr = MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecRestoreArray(in,&array);CHKERRQ(ierr);

  }
  ierr = PetscFree(rows);CHKERRQ(ierr);
  ierr = VecDestroy(&in);CHKERRQ(ierr);
  ierr = VecDestroy(&out);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

Ensure(FFT, transforms_constant_into_delta_function)
{
    Vec v;

    DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,4,2,1,NULL,&da);

    DMCreateGlobalVector(da, &v);
    VecZeroEntries(v);
    VecSetValue(v, 0, 0.7, INSERT_VALUES);
    VecAssemblyBegin(v);
    VecAssemblyEnd(v);

    scFftCreate(da, &fft);

    Vec x, y, z;
    scFftCreateVecsFFTW(fft, &x, &y, &z);

    scFftTransform(fft, v, 0, y);
    const PetscScalar *arr;
    VecGetArrayRead(y, &arr);
    //assert_that_double(fabs(arr[0] - 0.7), is_less_than_double(1.0e-6));
    assert_that(fabs(arr[1] - 0.0) < 1.0e-6, is_true);
    assert_that(fabs(arr[1] - 0.7) < 1.0e-6, is_true);
    assert_that(fabs(arr[2] - 0.0) < 1.0e-6, is_true);
    assert_that(fabs(arr[3] - 0.7) < 1.0e-6, is_true);
    assert_that(fabs(arr[4] - 0.0) < 1.0e-6, is_true);
    VecRestoreArrayRead(y, &arr);

    VecDestroy(&x);
    VecDestroy(&y);
    VecDestroy(&z);
    VecDestroy(&v);
    scFftDestroy(&fft);
    DMDestroy(&da);
}

PetscErrorCode BSSPotR1Vec(BSS self, PF pot, Vec V) {
  int order = self->order;
  int nb = self->num_basis;
  int ne = self->num_ele;
  int nq = order*ne;

  PetscErrorCode ierr;
  ierr = BSSCheck(self); CHKERRQ(ierr);

  PetscScalar *vs; PetscMalloc1(nq, &vs);  
  ierr = PFApply(pot, nq, self->Rrs, vs); CHKERRQ(ierr);
  for(int i = 0; i < nb; i++) {
    int k0, k1; Non0QuadIndex(i, i, order, nq, &k0, &k1);
    PetscScalar v = 0.0;
    for(int k = k0; k < k1; k++) {
      //v += self->vals[k+i*nq] * self->ws[k] * vs[k];
      v += self->vals[k+i*nq] * self->ws[k] * self->qrs[k] * vs[k];
    }
    ierr = VecSetValue(V, i, v, INSERT_VALUES); CHKERRQ(ierr);
  }

  PetscFree(vs);
  VecAssemblyBegin(V);
  VecAssemblyEnd(V);

  return 0;
}

/*---------------------------------------------------------
        Function to fill the rhs vector with
        By + g values ****
---------------------------------------------------------*/
PetscErrorCode rhs(AppCtx *obj,PetscScalar *y, PetscInt nz, PetscScalar *z, PetscReal t)
{
  PetscInt       i,j,js,je,jj;
  PetscScalar    val,g[num_z],btri[num_z][3],add_term;
  PetscErrorCode ierr;

  for (i=0; i < nz-2; i++) {
    for (j=0; j <= 2; j++) btri[i][j]=0.0;
    g[i] = 0.0;
  }

  /*  call femBg to set the tri-diagonal b matrix and vector g  */
  femBg(btri,g,nz,z,t);

  /*  setting the entries of the right hand side vector  */
  for (i=0; i < nz-2; i++) {
    val = 0.0;
    js  = 0;
    if (i == 0) js = 1;
    je = 2;
    if (i == nz-2) je = 1;

    for (jj=js; jj <= je; jj++) {
      j    = i+jj-1;
      val += (btri[i][jj])*(y[j]);
    }
    add_term = val + g[i];
    ierr = VecSetValue(obj->ksp_rhs,(PetscInt)i,(PetscScalar)add_term,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(obj->ksp_rhs);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(obj->ksp_rhs);CHKERRQ(ierr);

  /*  return to main driver function  */
  return 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);
}

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  PetscInt       i,blocks[2],nlocal;
  PetscMPIInt    size,rank;
  PetscScalar    value;
  Vec            x,y;
  IS             is1,is2;
  VecScatter     ctx = 0;
  PetscViewer    subviewer;

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

  if (size != 2) SETERRQ(PETSC_COMM_SELF,1,"Must run with 2 processors");

  /* create two vectors */
  if (!rank) nlocal = 8;
  else nlocal = 4;
  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr = VecSetSizes(x,nlocal,12);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr);
  ierr = VecCreate(PETSC_COMM_SELF,&y);CHKERRQ(ierr);
  ierr = VecSetSizes(y,8,PETSC_DECIDE);CHKERRQ(ierr);
  ierr = VecSetFromOptions(y);CHKERRQ(ierr);

  /* create two index sets */
  if (!rank) {
    blocks[0] = 0; blocks[1] = 2;
  } else {
    blocks[0] = 1; blocks[1] = 2;
  }
  ierr = ISCreateBlock(PETSC_COMM_SELF,4,2,blocks,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,8,0,1,&is2);CHKERRQ(ierr);

  for (i=0; i<12; i++) {
    value = i;
    ierr  = VecSetValues(x,1,&i,&value,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(x);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(x);CHKERRQ(ierr);

  ierr = VecScatterCreateWithData(x,is1,y,is2,&ctx);CHKERRQ(ierr);
  ierr = VecScatterBegin(ctx,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd(ctx,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterDestroy(&ctx);CHKERRQ(ierr);

  ierr = PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&subviewer);CHKERRQ(ierr);
  ierr = VecView(y,subviewer);CHKERRQ(ierr);
  ierr = PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&subviewer);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 ierr;
}

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

/*
    FormInitialGuess - Computes an initial approximation to the solution.

    Input Parameters:
.   user - user-defined application context
.   X    - vector

    Output Parameters:
    X    - vector
*/
PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
{
  PetscErrorCode ierr;
  PetscInt       i, j, k, mx = user->mx, my = user->my;
  PetscInt       xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye;
  PetscReal      hx = 1.0/(mx+1), hy = 1.0/(my+1), temp, val;

  PetscFunctionBegin;
  /* Get local mesh boundaries */
  ierr = DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
  ierr = DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);CHKERRQ(ierr);

  /* Compute initial guess over locally owned part of mesh */
  xe = xs+xm;
  ye = ys+ym;
  for (j=ys; j<ye; j++) {  /*  for (j=0; j<my; j++) */
    temp = PetscMin(j+1,my-j)*hy;
    for (i=xs; i<xe; i++) {  /*  for (i=0; i<mx; i++) */
      k   = (j-gys)*gxm + i-gxs;
      val = PetscMin((PetscMin(i+1,mx-i))*hx,temp);
      ierr = VecSetValuesLocal(X,1,&k,&val,ADD_VALUES);CHKERRQ(ierr);
    }
  }
  ierr = VecAssemblyBegin(X);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(X);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}