Commits

Mark Adams committed 3ae0bb6

protected for zero auxilary data size, fixed agg's destructor to unregister its setCoordinates method, changed name of a var for clarity

Comments (0)

Files changed (2)

src/ksp/pc/impls/gamg/agg.c

 
   PetscFunctionBegin;
   ierr = PetscFree(pc_gamg->subctx);CHKERRQ(ierr);
+  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 

src/ksp/pc/impls/gamg/gamg.c

   Mat             Cmat,Pold=*a_P_inout;
   MPI_Comm        comm;
   PetscMPIInt     rank,size,new_size,nactive=*a_nactive_proc;
-  PetscInt        ncrs_eq,ncrs_prim,f_bs;
+  PetscInt        ncrs_eq,ncrs,f_bs;
 
   PetscFunctionBegin;
   ierr = PetscObjectGetComm((PetscObject)Amat_fine,&comm);CHKERRQ(ierr);
   /* RAP */
   ierr = MatPtAP(Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat);CHKERRQ(ierr);
 
-  /* set 'ncrs_prim' (nodes), 'ncrs_eq' (equations)*/
-  ncrs_prim = pc_gamg->data_sz/pc_gamg->data_cell_cols/pc_gamg->data_cell_rows;
-  if (pc_gamg->data_sz % (pc_gamg->data_cell_cols*pc_gamg->data_cell_rows)) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"pc_gamg->data_sz %D not divisible by (pc_gamg->data_cell_cols %D *pc_gamg->data_cell_rows %D)",pc_gamg->data_sz,pc_gamg->data_cell_cols,pc_gamg->data_cell_rows);
+  /* set 'ncrs' (nodes), 'ncrs_eq' (equations)*/
   ierr = MatGetLocalSize(Cmat, &ncrs_eq, NULL);CHKERRQ(ierr);
+  if (pc_gamg->data_cell_rows>0) {
+    ncrs = pc_gamg->data_sz/pc_gamg->data_cell_cols/pc_gamg->data_cell_rows;
+  }
+  else {
+    PetscInt  bs;
+    ierr = MatGetBlockSize(Cmat, &bs);CHKERRQ(ierr);
+    ncrs = ncrs_eq/bs;
+  }
 
   /* get number of PEs to make active 'new_size', reduce, can be any integer 1-P */
   {
 
   if (!repart && new_size==nactive) *a_Amat_crs = Cmat; /* output - no repartitioning or reduction - could bail here */
   else {
-    PetscInt       *counts,*newproc_idx,ii,jj,kk,strideNew,*tidx,ncrs_prim_new,ncrs_eq_new,nloc_old;
+    PetscInt       *counts,*newproc_idx,ii,jj,kk,strideNew,*tidx,ncrs_new,ncrs_eq_new,nloc_old;
     IS             is_eq_newproc,is_eq_num,is_eq_num_prim,new_eq_indices;
 
     nloc_old = ncrs_eq/cr_bs;
         PetscInt          *d_nnz, *o_nnz, M, N;
         static PetscInt   llev = 0;
 
-        ierr = PetscMalloc(ncrs_prim*sizeof(PetscInt), &d_nnz);CHKERRQ(ierr);
-        ierr = PetscMalloc(ncrs_prim*sizeof(PetscInt), &o_nnz);CHKERRQ(ierr);
+        ierr = PetscMalloc(ncrs*sizeof(PetscInt), &d_nnz);CHKERRQ(ierr);
+        ierr = PetscMalloc(ncrs*sizeof(PetscInt), &o_nnz);CHKERRQ(ierr);
         ierr = MatGetOwnershipRange(Cmat, &Istart_crs, &Iend_crs);CHKERRQ(ierr);
         ierr = MatGetSize(Cmat, &M, &N);CHKERRQ(ierr);
         for (Ii = Istart_crs, jj = 0; Ii < Iend_crs; Ii += cr_bs, jj++) {
           d_nnz[jj] = ncols/cr_bs;
           o_nnz[jj] = ncols/cr_bs;
           ierr      = MatRestoreRow(Cmat,Ii,&ncols,0,0);CHKERRQ(ierr);
-          if (d_nnz[jj] > ncrs_prim) d_nnz[jj] = ncrs_prim;
-          if (o_nnz[jj] > (M/cr_bs-ncrs_prim)) o_nnz[jj] = M/cr_bs-ncrs_prim;
+          if (d_nnz[jj] > ncrs) d_nnz[jj] = ncrs;
+          if (o_nnz[jj] > (M/cr_bs-ncrs)) o_nnz[jj] = M/cr_bs-ncrs;
         }
 
         ierr = MatCreate(comm, &tMat);CHKERRQ(ierr);
-        ierr = MatSetSizes(tMat, ncrs_prim, ncrs_prim,PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr);
+        ierr = MatSetSizes(tMat, ncrs, ncrs,PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr);
         ierr = MatSetType(tMat,MATAIJ);CHKERRQ(ierr);
         ierr = MatSeqAIJSetPreallocation(tMat,0,d_nnz);CHKERRQ(ierr);
         ierr = MatMPIAIJSetPreallocation(tMat,0,d_nnz,0,o_nnz);CHKERRQ(ierr);
     ierr        = ISPartitioningCount(is_eq_newproc, size, counts);CHKERRQ(ierr);
     ncrs_eq_new = counts[rank];
     ierr        = ISDestroy(&is_eq_newproc);CHKERRQ(ierr);
-    ncrs_prim_new = ncrs_eq_new/cr_bs; /* eqs */
+    ncrs_new = ncrs_eq_new/cr_bs; /* eqs */
 
     ierr = PetscFree(counts);CHKERRQ(ierr);
 #if defined PETSC_GAMG_USE_LOG
     /* move data (for primal equations only) */
     /* Create a vector to contain the newly ordered element information */
     ierr = VecCreate(comm, &dest_crd);CHKERRQ(ierr);
-    ierr = VecSetSizes(dest_crd, node_data_sz*ncrs_prim_new, PETSC_DECIDE);CHKERRQ(ierr);
+    ierr = VecSetSizes(dest_crd, node_data_sz*ncrs_new, PETSC_DECIDE);CHKERRQ(ierr);
     ierr = VecSetType(dest_crd,VECSTANDARD);CHKERRQ(ierr); /* this is needed! */
     /*
      There are 'ndata_rows*ndata_cols' data items per node, (one can think of the vectors of having
      a block size of ...).  Note, ISs are expanded into equation space by 'cr_bs'.
      */
-    ierr = PetscMalloc((ncrs_prim*node_data_sz)*sizeof(PetscInt), &tidx);CHKERRQ(ierr);
+    ierr = PetscMalloc((ncrs*node_data_sz)*sizeof(PetscInt), &tidx);CHKERRQ(ierr);
     ierr = ISGetIndices(is_eq_num_prim, &idx);CHKERRQ(ierr);
-    for (ii=0,jj=0; ii<ncrs_prim; ii++) {
+    for (ii=0,jj=0; ii<ncrs; ii++) {
       PetscInt id = idx[ii*cr_bs]/cr_bs; /* get node back */
       for (kk=0; kk<node_data_sz; kk++, jj++) tidx[jj] = id*node_data_sz + kk;
     }
     ierr = ISRestoreIndices(is_eq_num_prim, &idx);CHKERRQ(ierr);
-    ierr = ISCreateGeneral(comm, node_data_sz*ncrs_prim, tidx, PETSC_COPY_VALUES, &isscat);CHKERRQ(ierr);
+    ierr = ISCreateGeneral(comm, node_data_sz*ncrs, tidx, PETSC_COPY_VALUES, &isscat);CHKERRQ(ierr);
     ierr = PetscFree(tidx);CHKERRQ(ierr);
     /*
      Create a vector to contain the original vertex information for each element
      */
-    ierr = VecCreateSeq(PETSC_COMM_SELF, node_data_sz*ncrs_prim, &src_crd);CHKERRQ(ierr);
+    ierr = VecCreateSeq(PETSC_COMM_SELF, node_data_sz*ncrs, &src_crd);CHKERRQ(ierr);
     for (jj=0; jj<ndata_cols; jj++) {
-      const PetscInt stride0=ncrs_prim*pc_gamg->data_cell_rows;
-      for (ii=0; ii<ncrs_prim; ii++) {
+      const PetscInt stride0=ncrs*pc_gamg->data_cell_rows;
+      for (ii=0; ii<ncrs; ii++) {
         for (kk=0; kk<ndata_rows; kk++) {
           PetscInt    ix = ii*ndata_rows + kk + jj*stride0, jx = ii*node_data_sz + kk*ndata_cols + jj;
           PetscScalar tt = (PetscScalar)pc_gamg->data[ix];
       Put the element vertex data into a new allocation of the gdata->ele
     */
     ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr);
-    ierr = PetscMalloc(node_data_sz*ncrs_prim_new*sizeof(PetscReal), &pc_gamg->data);CHKERRQ(ierr);
+    ierr = PetscMalloc(node_data_sz*ncrs_new*sizeof(PetscReal), &pc_gamg->data);CHKERRQ(ierr);
 
-    pc_gamg->data_sz = node_data_sz*ncrs_prim_new;
-    strideNew        = ncrs_prim_new*ndata_rows;
+    pc_gamg->data_sz = node_data_sz*ncrs_new;
+    strideNew        = ncrs_new*ndata_rows;
 
     ierr = VecGetArray(dest_crd, &array);CHKERRQ(ierr);
     for (jj=0; jj<ndata_cols; jj++) {
-      for (ii=0; ii<ncrs_prim_new; ii++) {
+      for (ii=0; ii<ncrs_new; ii++) {
         for (kk=0; kk<ndata_rows; kk++) {
           PetscInt ix = ii*ndata_rows + kk + jj*strideNew, jx = ii*node_data_sz + kk*ndata_cols + jj;
           pc_gamg->data[ix] = PetscRealPart(array[jx]);
     }
 
     /* stop if one node -- could pull back for singular problems */
-    if (M/pc_gamg->data_cell_cols < 2) {
+    if ( (pc_gamg->data_cell_cols && M/pc_gamg->data_cell_cols < 2) || (!pc_gamg->data_cell_cols && M < 2)) {
       level++;
       break;
     }
   ierr = PetscFunctionListFind(GAMGList,type,&r);CHKERRQ(ierr);
   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown GAMG type %s given",type);
   if (pc_gamg->ops->destroy) {
+    /* there was something here - kill it */
     ierr = (*pc_gamg->ops->destroy)(pc);CHKERRQ(ierr);
     ierr = PetscMemzero(pc_gamg->ops,sizeof(struct _PCGAMGOps));CHKERRQ(ierr);
+    /* cleaning up common data in pc_gamg - this should disapear someday */
+    pc_gamg->data_cell_cols = 0;
+    pc_gamg->data_cell_rows = 0;
+    pc_gamg->orig_data_cell_cols = 0;
+    pc_gamg->orig_data_cell_rows = 0;
+    if (pc_gamg->data_sz) {
+      ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr);
+      pc_gamg->data_sz = 0;
+    }
+    else if (pc_gamg->data) {
+      ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr); /* can this happen ? */
+    }
   }
   ierr = PetscFree(pc_gamg->gamg_type_name);CHKERRQ(ierr);
   ierr = PetscStrallocpy(type,&pc_gamg->gamg_type_name);CHKERRQ(ierr);