Commits

Mark Adams  committed a0dea87

PCGAMG: fix corruption for multiple solves when MatSetNearNullSpace is used

Member pc_gamg->data_sz was not being updated and thus caused corrupt
memory access. New test runex56_nns demonstrated the failure and now
runs correctly.

Reported-by: Garth N. Wells <gnw20@cam.ac.uk>
Helped-by: Jed Brown <jedbrown@mcs.anl.gov>

  • Participants
  • Parent commits 06e133e

Comments (0)

Files changed (4)

File src/ksp/ksp/examples/tutorials/ex56.c

   PetscErrorCode ierr;
   PetscInt       m,nn,M,Istart,Iend,i,j,k,ii,jj,kk,ic,ne=4,id;
   PetscReal      x,y,z,h,*coords,soft_alpha=1.e-3;
-  PetscBool      two_solves = PETSC_FALSE,test_nonzero_cols = PETSC_FALSE;
+  PetscBool      two_solves=PETSC_FALSE,test_nonzero_cols=PETSC_FALSE,use_nearnullspace=PETSC_FALSE;
   Vec            xx,bb;
   KSP            ksp;
   MPI_Comm       comm;
     ierr = PetscOptionsReal("-alpha","material coefficient inside circle","",soft_alpha,&soft_alpha,NULL);CHKERRQ(ierr);
     ierr = PetscOptionsBool("-two_solves","solve additional variant of the problem","",two_solves,&two_solves,NULL);CHKERRQ(ierr);
     ierr = PetscOptionsBool("-test_nonzero_cols","nonzero test","",test_nonzero_cols,&test_nonzero_cols,NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsBool("-use_mat_nearnullspace","MatNearNullSpace API test","",use_nearnullspace,&use_nearnullspace,NULL);CHKERRQ(ierr);
   }
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
 
 
   /* finish KSP/PC setup */
   ierr = KSPSetOperators(ksp, Amat, Amat, SAME_NONZERO_PATTERN);CHKERRQ(ierr);
-  ierr = PCSetCoordinates(pc, 3, m/3, coords);CHKERRQ(ierr);
+  if (use_nearnullspace) {
+    MatNullSpace matnull;
+    Vec vec_coords;
+    ierr = VecCreateSeqWithArray(MPI_COMM_SELF,3,m,coords,&vec_coords);CHKERRQ(ierr);
+    ierr = MatNullSpaceCreateRigidBody(vec_coords,&matnull);CHKERRQ(ierr);
+    ierr = MatSetNearNullSpace(Amat,matnull);CHKERRQ(ierr);
+    ierr = MatNullSpaceDestroy(&matnull);CHKERRQ(ierr);
+    ierr = VecDestroy(&vec_coords);CHKERRQ(ierr);
+  } else {
+    ierr = PCSetCoordinates(pc, 3, m/3, coords);CHKERRQ(ierr);
+  }
 
   ierr = MaybeLogStagePush(stage[0]);CHKERRQ(ierr);
 

File src/ksp/ksp/examples/tutorials/makefile

          ${DIFF} output/ex56_ml.out ex56.tmp || echo ${PWD} "\nPossible problem with with ex56_2, diffs above \n========================================="; \
         ${RM} -f ex56.tmp
 
+runex56_nns:
+	-@${MPIEXEC} -n 1 ./ex56 -ne 9 -alpha 1.e-3 -ksp_monitor_short -ksp_type cg -ksp_max_it 50 -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 1000 -mg_levels_ksp_type chebyshev -mg_levels_pc_type sor -pc_gamg_reuse_interpolation true -two_solves -use_mat_nearnullspace > ex56.tmp 2>&1;	\
+         ${DIFF} output/ex56_nns.out ex56.tmp || printf "${PWD}\nPossible problem with with ex56_nns, diffs above \n=========================================\n"; \
+         ${RM} -f ex56.tmp
+
 runex58:
 	-@${MPIEXEC} -n 1 ./ex58 -mat_type aij > ex58.tmp 2>&1;         \
 	${DIFF} output/ex58.out ex58.tmp || echo ${PWD} "\nPossible problem with with ex58, diffs above \n========================================="; \
                                  ex43.PETSc runex43 runex43_2 runex43_3 ex43.rm \
                                  ex45.PETSc runex45 runex45_2 ex45.rm \
                                  ex49.PETSc runex49 runex49_2 runex49_3 ex49.rm ex53.PETSc runex53 ex53.rm ex55.PETSc runex55_SA ex55.rm\
+                                 ex56.PETSc runex56_nns runex56 ex56.rm \
                                  ex58.PETSc runex58 runex58_baij runex58_sbaij ex58.rm
 TESTEXAMPLES_C_X	       = ex2.PETSc runex2_5 ex2.rm ex5.PETSc runex5_5 ex5.rm ex8.PETSc ex8.rm ex28.PETSc runex28 ex28.rm
 TESTEXAMPLES_FORTRAN	       = ex1f.PETSc runex1f ex1f.rm ex2f.PETSc runex2f ex2f.rm ex6f.PETSc ex6f.rm \

File src/ksp/ksp/examples/tutorials/output/ex56_nns.out

+  0 KSP Residual norm 761.484 
+  1 KSP Residual norm 101.691 
+  2 KSP Residual norm 83.8991 
+  3 KSP Residual norm 46.6149 
+  4 KSP Residual norm 3.81038 
+  5 KSP Residual norm 2.39808 
+  6 KSP Residual norm 0.683167 
+  7 KSP Residual norm 0.154152 
+  8 KSP Residual norm 0.00730291 
+  0 KSP Residual norm 0.00761484 
+  1 KSP Residual norm 0.00101691 
+  2 KSP Residual norm 0.000838991 
+  3 KSP Residual norm 0.000466149 
+  4 KSP Residual norm 3.81038e-05 
+  5 KSP Residual norm 2.39808e-05 
+  6 KSP Residual norm 6.83167e-06 
+  7 KSP Residual norm 1.54152e-06 
+  8 KSP Residual norm 7.30291e-08 
+  0 KSP Residual norm 7.61484e-08 
+  1 KSP Residual norm 1.01691e-08 
+  2 KSP Residual norm 8.38991e-09 
+  3 KSP Residual norm 4.66149e-09 
+  4 KSP Residual norm 3.810e-10 
+  5 KSP Residual norm 2.398e-10 
+  6 KSP Residual norm 6.832e-11 
+  7 KSP Residual norm 1.542e-11 
+  8 KSP Residual norm < 1.e-11
+[0]main |b-Ax|/|b|=2.271909e-04, |b|=5.391826e+00, emax=9.999167e-01

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

     const Vec *vecs; const PetscScalar *v;
     ierr = MatGetLocalSize(a_A,&mlocal,NULL);CHKERRQ(ierr);
     ierr = MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);CHKERRQ(ierr);
+    pc_gamg->data_sz = (nvec+!!has_const)*mlocal;
     ierr = PetscMalloc((nvec+!!has_const)*mlocal*sizeof(*nullvec),&nullvec);CHKERRQ(ierr);
     if (has_const) for (i=0; i<mlocal; i++) nullvec[i] = 1.0;
     for (i=0; i<nvec; i++) {