# HG changeset patch
# User Vijay Mahadevan <vijaysm@mcs.anl.gov>
# Date 1362105232 21600
# Branch vecnest_fixes
# Node ID f419331c4d51fca26b1b0be9b774df747e01db9a
# Parent  d343163d88abee75ed0b961db42e7f9509cb8fe5
Fixing a problem with the nested vector when VecDuplicate was called in parallel. The component IS were not being updated. The routine now handles all cases correctly.

diff -r d343163d88ab -r f419331c4d51 src/vec/vec/impls/nest/vecnest.c
--- a/src/vec/vec/impls/nest/vecnest.c	Thu Feb 28 13:50:39 2013 -0600
+++ b/src/vec/vec/impls/nest/vecnest.c	Thu Feb 28 20:33:52 2013 -0600
@@ -959,24 +959,20 @@ PetscErrorCode  VecNestGetSubVecs(Vec X,
   PetscFunctionReturn(0);
 }
 
+
 #undef __FUNCT__
 #define __FUNCT__ "VecNestSetSubVec_Private"
 static PetscErrorCode  VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
 {
   Vec_Nest       *bx = (Vec_Nest*)X->data;
-  PetscInt       i,offset=0,n=0,bs;
+  PetscInt       i,bs,offset=0,n=0,N=0;
   IS             is;
   PetscErrorCode ierr;
-  PetscBool      issame = PETSC_FALSE;
-  PetscInt       N=0;
+  PetscBool      issame_v=PETSC_FALSE, issame_is=PETSC_FALSE;
 
+  PetscFunctionBegin;
   /* check if idxm < bx->nb */
   if (idxm >= bx->nb) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",idxm,bx->nb);
-
-  PetscFunctionBegin;
-  ierr = VecDestroy(&bx->v[idxm]);CHKERRQ(ierr);       /* destroy the existing vector */
-  ierr = VecDuplicate(x,&bx->v[idxm]);CHKERRQ(ierr);   /* duplicate the layout of given vector */
-  ierr = VecCopy(x,bx->v[idxm]);CHKERRQ(ierr);         /* copy the contents of the given vector */
 
   /* check if we need to update the IS for the block */
   offset = X->map->rstart;
@@ -986,7 +982,8 @@ static PetscErrorCode  VecNestSetSubVec_
     offset += n;
   }
 
-  /* get the local size and block size */
+  n = bs = 0;
+  /* get the local size and block size of x */
   ierr = VecGetLocalSize(x,&n);CHKERRQ(ierr);
   ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr);
 
@@ -994,36 +991,57 @@ static PetscErrorCode  VecNestSetSubVec_
   ierr = ISCreateStride(PetscObjectComm((PetscObject)x),n,offset,1,&is);CHKERRQ(ierr);
   ierr = ISSetBlockSize(is,bs);CHKERRQ(ierr);
 
-  /* check if they are equal */
-  ierr = ISEqual(is,bx->is[idxm],&issame);CHKERRQ(ierr);
+  /* check if the vectors are equal */
+  ierr = VecEqual(x,bx->v[idxm],&issame_v);CHKERRQ(ierr);
 
-  if (!issame) {
-    /* The IS of given vector has a different layout compared to the existing block vector.
-     Destroy the existing reference and update the IS. */
-    ierr = ISDestroy(&bx->is[idxm]);CHKERRQ(ierr);
-    ierr = ISDuplicate(is,&bx->is[idxm]);CHKERRQ(ierr);
-    ierr = ISCopy(is,bx->is[idxm]);CHKERRQ(ierr);
+  /* if the vector layout is different, modify the object reference appropriately */
+  if (!issame_v) {
+    ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);
+    ierr = PetscObjectDereference((PetscObject)bx->v[i]);CHKERRQ(ierr);
+    bx->v[i] = x;
+  }
 
-    offset += n;
-    /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
-    for (i=idxm+1; i<bx->nb; i++) {
+  /* check if the index sets need to be updated */
+  ierr = ISEqual(is,bx->is[idxm],&issame_is);CHKERRQ(ierr);
+  if (!issame_is) {
+    N = n = 0;
+    /* get the local and global size of X */
+    ierr = VecSize_Nest_Recursive(X,PETSC_TRUE,&N);CHKERRQ(ierr);
+    ierr = VecSize_Nest_Recursive(X,PETSC_FALSE,&n);CHKERRQ(ierr);
+    
+    /* update the attributes of the layout */
+    ierr = PetscLayoutSetSize(X->map,N);CHKERRQ(ierr);
+    ierr = PetscLayoutSetLocalSize(X->map,n);CHKERRQ(ierr);
+    ierr = PetscLayoutSetBlockSize(X->map,1);CHKERRQ(ierr);
+    
+    /* free the allocated memory for range so that the layout can be recomputed in setup */
+    ierr = PetscFree(X->map->range);
+    ierr = PetscLayoutSetUp(X->map);CHKERRQ(ierr);
+
+    /* Since the current IS[idxm] changed, we need to update the offsets for all IS.
+       It is not sufficient to update only subsequent IS since the parallel layout could have
+       been modified by the new block vector reference. */
+    offset = X->map->rstart;
+    for (i=0; i<bx->nb; i++) {
+      n = bs = 0;
+      
       /* get the local size and block size */
       ierr = VecGetLocalSize(bx->v[i],&n);CHKERRQ(ierr);
       ierr = VecGetBlockSize(bx->v[i],&bs);CHKERRQ(ierr);
 
-      /* destroy the old and create the new IS */
-      ierr = ISDestroy(&bx->is[i]);CHKERRQ(ierr);
-      ierr = ISCreateStride(((PetscObject)bx->v[i])->comm,n,offset,1,&bx->is[i]);CHKERRQ(ierr);
-      ierr = ISSetBlockSize(bx->is[i],bs);CHKERRQ(ierr);
-
+      if (i == idxm) {
+        /* The IS of given vector has a different layout compared to the existing block vector.
+         Destroy the existing reference and update the IS. */
+        ierr = ISDestroy(&bx->is[idxm]);CHKERRQ(ierr);
+        ierr = ISCreateStride(((PetscObject)x)->comm,n,offset,1,&bx->is[idxm]);CHKERRQ(ierr);
+        ierr = ISSetBlockSize(bx->is[idxm],bs);CHKERRQ(ierr);
+      }
+      else {
+        /* The parallel layout has not changed. Only update the stride and offset. */
+        ierr = ISStrideSetStride(bx->is[i],n,offset,bs);CHKERRQ(ierr);
+      }
       offset += n;
     }
-
-    n=0;
-    ierr = VecSize_Nest_Recursive(X,PETSC_TRUE,&N);CHKERRQ(ierr);
-    ierr = VecSize_Nest_Recursive(X,PETSC_FALSE,&n);CHKERRQ(ierr);
-    ierr = PetscLayoutSetSize(X->map,N);CHKERRQ(ierr);
-    ierr = PetscLayoutSetLocalSize(X->map,n);CHKERRQ(ierr);
   }
 
   ierr = ISDestroy(&is);CHKERRQ(ierr);