BV: modify BVGetSplit() to prefer creating instead of resizing an exi…

…sting BV This partly reverts b44b2f3. In the previous implementation, the subsidiary objects (Vec in BVSVEC, Mat in BVMAT) were not resized and hence there could be side effects, e.g. in block orthogonalization. Reported-by: Pierre Jolivet
firedrakeproject · Jul 7, 2020 · c91a490 · c91a490
1 parent 40b9ec6
commit c91a490
Show file tree

Hide file tree

Showing 6 changed files with 94 additions and 143 deletions.
diff --git a/src/sys/classes/bv/impls/contiguous/contig.c b/src/sys/classes/bv/impls/contiguous/contig.c
@@ -245,45 +245,37 @@ PetscErrorCode BVResize_Contiguous(BV bv,PetscInt m,PetscBool copy)
 {
   PetscErrorCode ierr;
   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)bv->data;
-  PetscInt       j,bs,lsplit;
+  PetscInt       j,bs;
   PetscScalar    *newarray;
   Vec            *newV;
   char           str[50];
-  BV             parent;
 
   PetscFunctionBegin;
-  if (bv->issplit==2) {
-    parent = bv->splitparent;
-    lsplit = parent->lsplit;
-    ctx->V = ((BV_CONTIGUOUS*)parent->data)->V+lsplit;
-    ctx->array = ((BV_CONTIGUOUS*)parent->data)->array+lsplit*bv->n;
-  } else if (!bv->issplit) {
-    ierr = VecGetBlockSize(bv->t,&bs);CHKERRQ(ierr);
-    ierr = PetscMalloc1(m*bv->n,&newarray);CHKERRQ(ierr);
-    ierr = PetscArrayzero(newarray,m*bv->n);CHKERRQ(ierr);
-    ierr = PetscMalloc1(m,&newV);CHKERRQ(ierr);
-    for (j=0;j<m;j++) {
-      if (ctx->mpi) {
-        ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv->t),bs,bv->n,PETSC_DECIDE,newarray+j*bv->n,newV+j);CHKERRQ(ierr);
-      } else {
-        ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv->t),bs,bv->n,newarray+j*bv->n,newV+j);CHKERRQ(ierr);
-      }
-    }
-    ierr = PetscLogObjectParents(bv,m,newV);CHKERRQ(ierr);
-    if (((PetscObject)bv)->name) {
-      for (j=0;j<m;j++) {
-        ierr = PetscSNPrintf(str,sizeof(str),"%s_%d",((PetscObject)bv)->name,(int)j);CHKERRQ(ierr);
-        ierr = PetscObjectSetName((PetscObject)newV[j],str);CHKERRQ(ierr);
-      }
+  ierr = VecGetBlockSize(bv->t,&bs);CHKERRQ(ierr);
+  ierr = PetscMalloc1(m*bv->n,&newarray);CHKERRQ(ierr);
+  ierr = PetscArrayzero(newarray,m*bv->n);CHKERRQ(ierr);
+  ierr = PetscMalloc1(m,&newV);CHKERRQ(ierr);
+  for (j=0;j<m;j++) {
+    if (ctx->mpi) {
+      ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv->t),bs,bv->n,PETSC_DECIDE,newarray+j*bv->n,newV+j);CHKERRQ(ierr);
+    } else {
+      ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv->t),bs,bv->n,newarray+j*bv->n,newV+j);CHKERRQ(ierr);
     }
-    if (copy) {
-      ierr = PetscArraycpy(newarray,ctx->array,PetscMin(m,bv->m)*bv->n);CHKERRQ(ierr);
+  }
+  ierr = PetscLogObjectParents(bv,m,newV);CHKERRQ(ierr);
+  if (((PetscObject)bv)->name) {
+    for (j=0;j<m;j++) {
+      ierr = PetscSNPrintf(str,sizeof(str),"%s_%d",((PetscObject)bv)->name,(int)j);CHKERRQ(ierr);
+      ierr = PetscObjectSetName((PetscObject)newV[j],str);CHKERRQ(ierr);
     }
-    ierr = VecDestroyVecs(bv->m,&ctx->V);CHKERRQ(ierr);
-    ctx->V = newV;
-    ierr = PetscFree(ctx->array);CHKERRQ(ierr);
-    ctx->array = newarray;
   }
+  if (copy) {
+    ierr = PetscArraycpy(newarray,ctx->array,PetscMin(m,bv->m)*bv->n);CHKERRQ(ierr);
+  }
+  ierr = VecDestroyVecs(bv->m,&ctx->V);CHKERRQ(ierr);
+  ctx->V = newV;
+  ierr = PetscFree(ctx->array);CHKERRQ(ierr);
+  ctx->array = newarray;
   PetscFunctionReturn(0);
 }
 

diff --git a/src/sys/classes/bv/impls/mat/bvmat.c b/src/sys/classes/bv/impls/mat/bvmat.c
@@ -285,38 +285,27 @@ PetscErrorCode BVResize_Mat(BV bv,PetscInt m,PetscBool copy)
   PetscErrorCode ierr;
   BV_MAT         *ctx = (BV_MAT*)bv->data;
   PetscScalar    *pA,*pnew;
-  PetscInt       lsplit;
   Mat            A;
   char           str[50];
-  PetscScalar    *array;
-  BV             parent;
 
   PetscFunctionBegin;
-  if (bv->issplit==2) {
-    parent = bv->splitparent;
-    lsplit = parent->lsplit;
-    ierr = MatDenseGetArray(((BV_MAT*)parent->data)->A,&array);CHKERRQ(ierr);
-    ierr = MatDensePlaceArray(ctx->A,array+lsplit*bv->n);CHKERRQ(ierr);
-    ierr = MatDenseRestoreArray(((BV_MAT*)parent->data)->A,&array);CHKERRQ(ierr);
-  } else if (!bv->issplit) {
-    ierr = MatCreateDense(PetscObjectComm((PetscObject)bv->t),bv->n,PETSC_DECIDE,PETSC_DECIDE,m,NULL,&A);CHKERRQ(ierr);
-    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-    ierr = PetscLogObjectParent((PetscObject)bv,(PetscObject)A);CHKERRQ(ierr);
-    if (((PetscObject)bv)->name) {
-      ierr = PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name);CHKERRQ(ierr);
-      ierr = PetscObjectSetName((PetscObject)A,str);CHKERRQ(ierr);
-    }
-    if (copy) {
-      ierr = MatDenseGetArray(ctx->A,&pA);CHKERRQ(ierr);
-      ierr = MatDenseGetArray(A,&pnew);CHKERRQ(ierr);
-      ierr = PetscArraycpy(pnew,pA,PetscMin(m,bv->m)*bv->n);CHKERRQ(ierr);
-      ierr = MatDenseRestoreArray(ctx->A,&pA);CHKERRQ(ierr);
-      ierr = MatDenseRestoreArray(A,&pnew);CHKERRQ(ierr);
-    }
-    ierr = MatDestroy(&ctx->A);CHKERRQ(ierr);
-    ctx->A = A;
+  ierr = MatCreateDense(PetscObjectComm((PetscObject)bv->t),bv->n,PETSC_DECIDE,PETSC_DECIDE,m,NULL,&A);CHKERRQ(ierr);
+  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = PetscLogObjectParent((PetscObject)bv,(PetscObject)A);CHKERRQ(ierr);
+  if (((PetscObject)bv)->name) {
+    ierr = PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name);CHKERRQ(ierr);
+    ierr = PetscObjectSetName((PetscObject)A,str);CHKERRQ(ierr);
+  }
+  if (copy) {
+    ierr = MatDenseGetArray(ctx->A,&pA);CHKERRQ(ierr);
+    ierr = MatDenseGetArray(A,&pnew);CHKERRQ(ierr);
+    ierr = PetscArraycpy(pnew,pA,PetscMin(m,bv->m)*bv->n);CHKERRQ(ierr);
+    ierr = MatDenseRestoreArray(ctx->A,&pA);CHKERRQ(ierr);
+    ierr = MatDenseRestoreArray(A,&pnew);CHKERRQ(ierr);
   }
+  ierr = MatDestroy(&ctx->A);CHKERRQ(ierr);
+  ctx->A = A;
   PetscFunctionReturn(0);
 }
 

diff --git a/src/sys/classes/bv/impls/svec/svec.c b/src/sys/classes/bv/impls/svec/svec.c
@@ -284,41 +284,30 @@ PetscErrorCode BVResize_Svec(BV bv,PetscInt m,PetscBool copy)
   BV_SVEC           *ctx = (BV_SVEC*)bv->data;
   PetscScalar       *pnew;
   const PetscScalar *pv;
-  PetscInt          bs,lsplit;
-  Vec               vnew,vpar;
+  PetscInt          bs;
+  Vec               vnew;
   char              str[50];
-  BV                parent;
 
   PetscFunctionBegin;
-  if (bv->issplit==2) {
-    parent = bv->splitparent;
-    lsplit = parent->lsplit;
-    vpar = ((BV_SVEC*)parent->data)->v;
-    ierr = VecGetArrayRead(vpar,&pv);CHKERRQ(ierr);
-    ierr = VecResetArray(ctx->v);CHKERRQ(ierr);
-    ierr = VecPlaceArray(ctx->v,pv+lsplit*bv->n);CHKERRQ(ierr);
-    ierr = VecRestoreArrayRead(vpar,&pv);CHKERRQ(ierr);
-  } else if (!bv->issplit) {
-    ierr = VecGetBlockSize(bv->t,&bs);CHKERRQ(ierr);
-    ierr = VecCreate(PetscObjectComm((PetscObject)bv->t),&vnew);CHKERRQ(ierr);
-    ierr = VecSetType(vnew,((PetscObject)bv->t)->type_name);CHKERRQ(ierr);
-    ierr = VecSetSizes(vnew,m*bv->n,PETSC_DECIDE);CHKERRQ(ierr);
-    ierr = VecSetBlockSize(vnew,bs);CHKERRQ(ierr);
-    ierr = PetscLogObjectParent((PetscObject)bv,(PetscObject)vnew);CHKERRQ(ierr);
-    if (((PetscObject)bv)->name) {
-      ierr = PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name);CHKERRQ(ierr);
-      ierr = PetscObjectSetName((PetscObject)vnew,str);CHKERRQ(ierr);
-    }
-    if (copy) {
-      ierr = VecGetArrayRead(ctx->v,&pv);CHKERRQ(ierr);
-      ierr = VecGetArray(vnew,&pnew);CHKERRQ(ierr);
-      ierr = PetscArraycpy(pnew,pv,PetscMin(m,bv->m)*bv->n);CHKERRQ(ierr);
-      ierr = VecRestoreArrayRead(ctx->v,&pv);CHKERRQ(ierr);
-      ierr = VecRestoreArray(vnew,&pnew);CHKERRQ(ierr);
-    }
-    ierr = VecDestroy(&ctx->v);CHKERRQ(ierr);
-    ctx->v = vnew;
+  ierr = VecGetBlockSize(bv->t,&bs);CHKERRQ(ierr);
+  ierr = VecCreate(PetscObjectComm((PetscObject)bv->t),&vnew);CHKERRQ(ierr);
+  ierr = VecSetType(vnew,((PetscObject)bv->t)->type_name);CHKERRQ(ierr);
+  ierr = VecSetSizes(vnew,m*bv->n,PETSC_DECIDE);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(vnew,bs);CHKERRQ(ierr);
+  ierr = PetscLogObjectParent((PetscObject)bv,(PetscObject)vnew);CHKERRQ(ierr);
+  if (((PetscObject)bv)->name) {
+    ierr = PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name);CHKERRQ(ierr);
+    ierr = PetscObjectSetName((PetscObject)vnew,str);CHKERRQ(ierr);
   }
+  if (copy) {
+    ierr = VecGetArrayRead(ctx->v,&pv);CHKERRQ(ierr);
+    ierr = VecGetArray(vnew,&pnew);CHKERRQ(ierr);
+    ierr = PetscArraycpy(pnew,pv,PetscMin(m,bv->m)*bv->n);CHKERRQ(ierr);
+    ierr = VecRestoreArrayRead(ctx->v,&pv);CHKERRQ(ierr);
+    ierr = VecRestoreArray(vnew,&pnew);CHKERRQ(ierr);
+  }
+  ierr = VecDestroy(&ctx->v);CHKERRQ(ierr);
+  ctx->v = vnew;
   PetscFunctionReturn(0);
 }
 

diff --git a/src/sys/classes/bv/impls/svec/sveccuda/sveccuda.cu b/src/sys/classes/bv/impls/svec/sveccuda/sveccuda.cu
@@ -617,43 +617,32 @@ PetscErrorCode BVResize_Svec_CUDA(BV bv,PetscInt m,PetscBool copy)
   PetscErrorCode    ierr;
   BV_SVEC           *ctx = (BV_SVEC*)bv->data;
   const PetscScalar *d_pv;
-  PetscScalar       *d_pnew,*d_ptr;
-  PetscInt          bs,lsplit;
-  Vec               vnew,vpar;
+  PetscScalar       *d_pnew;
+  PetscInt          bs;
+  Vec               vnew;
   char              str[50];
   cudaError_t       cerr;
-  BV                parent;
 
   PetscFunctionBegin;
-  if (bv->issplit==2) {
-    parent = bv->splitparent;
-    lsplit = parent->lsplit;
-    vpar = ((BV_SVEC*)parent->data)->v;
-    ierr = VecCUDAResetArray(ctx->v);CHKERRQ(ierr);
-    ierr = VecCUDAGetArray(vpar,&d_ptr);CHKERRQ(ierr);
-    ierr = VecCUDAPlaceArray(ctx->v,d_ptr+lsplit*bv->n);CHKERRQ(ierr);
-    ierr = VecCUDARestoreArray(vpar,&d_ptr);CHKERRQ(ierr);
-  } else if (!bv->issplit) {
-    ierr = VecGetBlockSize(bv->t,&bs);CHKERRQ(ierr);
-    ierr = VecCreate(PetscObjectComm((PetscObject)bv->t),&vnew);CHKERRQ(ierr);
-    ierr = VecSetType(vnew,((PetscObject)bv->t)->type_name);CHKERRQ(ierr);
-    ierr = VecSetSizes(vnew,m*bv->n,PETSC_DECIDE);CHKERRQ(ierr);
-    ierr = VecSetBlockSize(vnew,bs);CHKERRQ(ierr);
-    ierr = PetscLogObjectParent((PetscObject)bv,(PetscObject)vnew);CHKERRQ(ierr);
-    if (((PetscObject)bv)->name) {
-      ierr = PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name);CHKERRQ(ierr);
-      ierr = PetscObjectSetName((PetscObject)vnew,str);CHKERRQ(ierr);
-    }
-    if (copy) {
-      ierr = VecCUDAGetArrayRead(ctx->v,&d_pv);CHKERRQ(ierr);
-      ierr = VecCUDAGetArrayWrite(vnew,&d_pnew);CHKERRQ(ierr);
-      cerr = cudaMemcpy(d_pnew,d_pv,PetscMin(m,bv->m)*bv->n*sizeof(PetscScalar),cudaMemcpyDeviceToDevice);CHKERRCUDA(cerr);
-      ierr = VecCUDARestoreArrayRead(ctx->v,&d_pv);CHKERRQ(ierr);
-      ierr = VecCUDARestoreArrayWrite(vnew,&d_pnew);CHKERRQ(ierr);
-    }
-    ierr = VecDestroy(&ctx->v);CHKERRQ(ierr);
-    ctx->v = vnew;
+  ierr = VecGetBlockSize(bv->t,&bs);CHKERRQ(ierr);
+  ierr = VecCreate(PetscObjectComm((PetscObject)bv->t),&vnew);CHKERRQ(ierr);
+  ierr = VecSetType(vnew,((PetscObject)bv->t)->type_name);CHKERRQ(ierr);
+  ierr = VecSetSizes(vnew,m*bv->n,PETSC_DECIDE);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(vnew,bs);CHKERRQ(ierr);
+  ierr = PetscLogObjectParent((PetscObject)bv,(PetscObject)vnew);CHKERRQ(ierr);
+  if (((PetscObject)bv)->name) {
+    ierr = PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name);CHKERRQ(ierr);
+    ierr = PetscObjectSetName((PetscObject)vnew,str);CHKERRQ(ierr);
+  }
+  if (copy) {
+    ierr = VecCUDAGetArrayRead(ctx->v,&d_pv);CHKERRQ(ierr);
+    ierr = VecCUDAGetArrayWrite(vnew,&d_pnew);CHKERRQ(ierr);
+    cerr = cudaMemcpy(d_pnew,d_pv,PetscMin(m,bv->m)*bv->n*sizeof(PetscScalar),cudaMemcpyDeviceToDevice);CHKERRCUDA(cerr);
+    ierr = VecCUDARestoreArrayRead(ctx->v,&d_pv);CHKERRQ(ierr);
+    ierr = VecCUDARestoreArrayWrite(vnew,&d_pnew);CHKERRQ(ierr);
   }
+  ierr = VecDestroy(&ctx->v);CHKERRQ(ierr);
+  ctx->v = vnew;
   PetscFunctionReturn(0);
 }
 

diff --git a/src/sys/classes/bv/impls/vecs/vecs.c b/src/sys/classes/bv/impls/vecs/vecs.c
@@ -345,32 +345,25 @@ PetscErrorCode BVResize_Vecs(BV bv,PetscInt m,PetscBool copy)
   PetscErrorCode ierr;
   BV_VECS        *ctx = (BV_VECS*)bv->data;
   Vec            *newV;
-  PetscInt       j,lsplit;
+  PetscInt       j;
   char           str[50];
-  BV             parent;
 
   PetscFunctionBegin;
-  if (bv->issplit==2) {
-    parent = bv->splitparent;
-    lsplit = parent->lsplit;
-    ctx->V = ((BV_VECS*)parent->data)->V+lsplit;
-  } else if (!bv->issplit) {
-    ierr = VecDuplicateVecs(bv->t,m,&newV);CHKERRQ(ierr);
-    ierr = PetscLogObjectParents(bv,m,newV);CHKERRQ(ierr);
-    if (((PetscObject)bv)->name) {
-      for (j=0;j<m;j++) {
-        ierr = PetscSNPrintf(str,sizeof(str),"%s_%D",((PetscObject)bv)->name,j);CHKERRQ(ierr);
-        ierr = PetscObjectSetName((PetscObject)newV[j],str);CHKERRQ(ierr);
-      }
+  ierr = VecDuplicateVecs(bv->t,m,&newV);CHKERRQ(ierr);
+  ierr = PetscLogObjectParents(bv,m,newV);CHKERRQ(ierr);
+  if (((PetscObject)bv)->name) {
+    for (j=0;j<m;j++) {
+      ierr = PetscSNPrintf(str,sizeof(str),"%s_%D",((PetscObject)bv)->name,j);CHKERRQ(ierr);
+      ierr = PetscObjectSetName((PetscObject)newV[j],str);CHKERRQ(ierr);
     }
-    if (copy) {
-      for (j=0;j<PetscMin(m,bv->m);j++) {
-        ierr = VecCopy(ctx->V[j],newV[j]);CHKERRQ(ierr);
-      }
+  }
+  if (copy) {
+    for (j=0;j<PetscMin(m,bv->m);j++) {
+      ierr = VecCopy(ctx->V[j],newV[j]);CHKERRQ(ierr);
     }
-    ierr = VecDestroyVecs(bv->m,&ctx->V);CHKERRQ(ierr);
-    ctx->V = newV;
   }
+  ierr = VecDestroyVecs(bv->m,&ctx->V);CHKERRQ(ierr);
+  ctx->V = newV;
   PetscFunctionReturn(0);
 }
 

diff --git a/src/sys/classes/bv/interface/bvbasic.c b/src/sys/classes/bv/interface/bvbasic.c
@@ -1812,15 +1812,14 @@ static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
 
   PetscFunctionBegin;
   ncols = left? bv->nc+bv->l: bv->m-bv->l;
+  if (*split && ncols!=(*split)->m) { ierr = BVDestroy(split);CHKERRQ(ierr); }
   if (!*split) {
     ierr = BVCreate(PetscObjectComm((PetscObject)bv),split);CHKERRQ(ierr);
     ierr = PetscLogObjectParent((PetscObject)bv,(PetscObject)*split);CHKERRQ(ierr);
     (*split)->issplit = left? 1: 2;
     (*split)->splitparent = bv;
     ierr = BVSetSizesFromVec(*split,bv->t,ncols);CHKERRQ(ierr);
     ierr = BVDuplicate_Private(bv,*split);CHKERRQ(ierr);
-  } else {
-    ierr = BVResize(*split,ncols,PETSC_FALSE);CHKERRQ(ierr);
   }
   (*split)->l  = 0;
   (*split)->k  = left? bv->l: bv->k-bv->l;