Thanks! For some reason I was thinking that -pc_mg_galerkin was a default, which it is not as can be seen from the -help output. Now I have hardcoded it, and it works.
I'll try to apply that rows/cols logic to the distributed interpolation matrices for the parallel case, after I have verified that the problem is solved correctly in the sequential case. Right now I only know that it runs without errors. Best regards Torquil Sørensen On 12/11/13 20:27, Barry Smith wrote: > y = A *x > > The number of local rows in A is the number of local rows in y. The > number of local columns in A is the number of local rows in x. > > To keep it from crashing you need to either > > 1) provide the coarse grid operator or > > 2) run with -pc_mg_galerkin to have the coarse grid operator computed from > the interpolation you provided and the fine grid operator you provided. > > Barry > > On Nov 12, 2013, at 1:16 PM, Torquil Macdonald Sørensen <[email protected]> > wrote: > >> Hi! >> >> I'm having some problems getting PCMG to work when manually populating the >> interpolation matrix, i.e. using MatSetValue to set its values, instead of >> the use of DMCreateInterpolation, as I have seen in an example. I have to >> use MatSetValue to populate my PETSc matrices, because I'm obtaining them in >> a different format from a FEM library. >> >> I have reduced my program to a minimal test case (without involving the FEM >> library), where the system matrix is simply an identity matrix, with two >> multigrid levels of 8 and 4 degrees of freedom, and a very simple >> interpolation matrix that can be seen by looking at the code I'm attaching. >> I'm not including KSPSolve because the problem occurs prior to that, on line >> 56, which corresponds to where I'm calling KSPSetUp(). >> >> I'm actually not certain how to create the matrix that will be passed to >> PCMGSetInterpolation(). At the moment I'm only running sequentially, but how >> should this non-square matrix be distributed across the MPI nodes when >> running in parallel, i.e. what values should I pass to MatSetSizes() for >> that matrix? >> >> In any case, the error message I'm getting (running it with no command line >> options) is: >> >> [0]PETSC ERROR: --------------------- Error Message >> ------------------------------------ >> [0]PETSC ERROR: Object is in wrong state! >> [0]PETSC ERROR: Mat object's type is not set: Argument # 1! >> [0]PETSC ERROR: >> ------------------------------------------------------------------------ >> [0]PETSC ERROR: Petsc Release Version 3.4.3, Oct, 15, 2013 >> [0]PETSC ERROR: See docs/changes/index.html for recent updates. >> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting. >> [0]PETSC ERROR: See docs/index.html for manual pages. >> [0]PETSC ERROR: >> ------------------------------------------------------------------------ >> [0]PETSC ERROR: ./multigrid2.elf on a arch-linux2-cxx-debug named >> enlil.uio.no by tmac Tue Nov 12 19:57:59 2013 >> [0]PETSC ERROR: Libraries linked from >> /mn/anatu/cma-u3/tmac/usr/stow/petsc_complex_dbg/lib >> [0]PETSC ERROR: Configure run at Tue Nov 12 16:27:32 2013 >> [0]PETSC ERROR: Configure options >> CPPFLAGS=-I/mn/anatu/cma-u3/tmac/usr/include >> LDFLAGS=-L/mn/anatu/cma-u3/tmac/usr/lib -L/mn/anatu/cma-u3/tmac/usr/lib64 >> --prefix=/mn/anatu/cma-u3/tmac/usr/stow/petsc_complex_dbg >> --with-clanguage=C++ --with-scalar-type=complex --with-shared-libraries=1 >> --with-debugging=1 --with-superlu=1 >> --with-superlu-lib=/mn/anatu/cma-u3/tmac/usr/lib/libsuperlu.so >> --with-superlu-include=/mn/anatu/cma-u3/tmac/usr/include/superlu/ >> [0]PETSC ERROR: >> ------------------------------------------------------------------------ >> [0]PETSC ERROR: MatGetVecs() line 8131 in >> /tmp/petsc-3.4.3/src/mat/interface/matrix.c >> [0]PETSC ERROR: KSPGetVecs() line 929 in >> /tmp/petsc-3.4.3/src/ksp/ksp/interface/iterativ.c >> [0]PETSC ERROR: PCSetUp_MG() line 691 in >> /tmp/petsc-3.4.3/src/ksp/pc/impls/mg/mg.c >> [0]PETSC ERROR: PCSetUp() line 890 in >> /tmp/petsc-3.4.3/src/ksp/pc/interface/precon.c >> [0]PETSC ERROR: KSPSetUp() line 278 in >> /tmp/petsc-3.4.3/src/ksp/ksp/interface/itfunc.c >> [0]PETSC ERROR: main() line 56 in >> "unknowndirectory/"/mn/anatu/cma-u3/tmac/programming/fem/getfem/source/multigrid2.cpp >> -------------------------------------------------------------------------- >> MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD >> with errorcode 73. >> >> NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes. >> You may or may not see output from other processes, depending on >> exactly when Open MPI kills them. >> -------------------------------------------------------------------------- >> >> The full code is: >> >> #include "petscksp.h" >> >> int main(int argc, char ** argv) >> { >> PetscErrorCode ierr = PetscInitialize(&argc, &argv, 0, 0); CHKERRQ(ierr); >> >> Vec x; >> Mat A; >> PetscInt rstart, rend, nlocal; >> >> PetscInt const n = 8; >> >> ierr = VecCreate(PETSC_COMM_WORLD, &x); CHKERRQ(ierr); >> ierr = VecSetSizes(x, PETSC_DECIDE, n); CHKERRQ(ierr); >> ierr = VecSetFromOptions(x); CHKERRQ(ierr); >> ierr = VecGetOwnershipRange(x, &rstart, &rend); CHKERRQ(ierr); >> ierr = VecGetLocalSize(x, &nlocal); CHKERRQ(ierr); >> >> ierr = MatCreate(PETSC_COMM_WORLD, &A); CHKERRQ(ierr); >> ierr = MatSetSizes(A, nlocal, nlocal, n, n); CHKERRQ(ierr); >> ierr = MatSetFromOptions(A); CHKERRQ(ierr); >> ierr = MatSetUp(A); CHKERRQ(ierr); >> for(PetscInt i = rstart; i < rend; ++i) { >> ierr = MatSetValue(A, i, i, 1.0, INSERT_VALUES); CHKERRQ(ierr); >> } >> ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); >> ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); >> >> KSP ksp; >> ierr = KSPCreate(PETSC_COMM_WORLD, &ksp); CHKERRQ(ierr); >> ierr = KSPSetOperators(ksp, A, A, SAME_NONZERO_PATTERN); CHKERRQ(ierr); >> ierr = KSPSetType(ksp, KSPFGMRES); CHKERRQ(ierr); >> >> PC pc; >> ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr); >> ierr = PCSetType(pc, PCMG); CHKERRQ(ierr); >> ierr = PCMGSetLevels(pc, 2, 0); CHKERRQ(ierr); >> >> Mat interp; >> ierr = MatCreate(PETSC_COMM_WORLD, &interp); CHKERRQ(ierr); >> ierr = MatSetSizes(interp, PETSC_DECIDE, PETSC_DECIDE, n, n/2); >> CHKERRQ(ierr); >> ierr = MatSetFromOptions(interp); CHKERRQ(ierr); >> ierr = MatSetUp(interp); CHKERRQ(ierr); >> ierr = MatGetOwnershipRange(interp, &rstart, &rend); CHKERRQ(ierr); >> for(PetscInt i = rstart; i < rend; ++i) { >> ierr = MatSetValue(interp, i, i/2, 1.0, INSERT_VALUES); >> CHKERRQ(ierr); >> } >> ierr = MatAssemblyBegin(interp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); >> ierr = MatAssemblyEnd(interp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); >> >> ierr = PCMGSetInterpolation(pc, 1, interp); CHKERRQ(ierr); >> ierr = MatDestroy(&interp); CHKERRQ(ierr); >> >> ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr); >> ierr = KSPSetUp(ksp); CHKERRQ(ierr); >> >> ierr = VecDestroy(&x); CHKERRQ(ierr); >> ierr = MatDestroy(&A); CHKERRQ(ierr); >> ierr = KSPDestroy(&ksp); CHKERRQ(ierr); >> ierr = PetscFinalize(); CHKERRQ(ierr); >> >> return 0; >> }
