[petsc-users] DaSetGetMatrix

[domenico.borzacchiello at univ-st-etienne.fr]

Hi,

I'm running my code (3D Stokes Solver with MAC arrangement) pretty fine so
far with a FieldSplit/Schur Preconditioning.

I had to write my own DAGetMatrix routine cause it was using too much
memory (5 times the required size) for the matrices. The code works  with
1 2 3 5 etc procs (I presume with any number of procs for which are only
possible 1D cartesian topologies of communicators i.e. any prime numbers) 
then if I run with 4 procs for example it stops when assembling the MPIAIJ
matrix with the following error:

[1]PETSC ERROR: Nonconforming object sizes!
[1]PETSC ERROR: Local scatter sizes don't match!

What could be causing the error?

Thank you,
Domenico

here's the getmatrix function I'm using

#undef __FUNCT__
#define __FUNCT__ "DAGetMatrix_User_2"
PetscErrorCode DAGetMatrix_User_2(DA da,const MatType mtype,Mat *J)
{
  PetscErrorCode ierr;
  Mat A;
  PetscInt xm,ym,zm,dim,dof,starts[3],dims[3];
  const MatType  Atype;
  void          
(*aij)(void)=PETSC_NULL,(*baij)(void)=PETSC_NULL,(*sbaij)(void)=PETSC_NULL;
  ISLocalToGlobalMapping ltog,ltogb;

  PetscFunctionBegin;
  ierr = DAGetInfo(da,&dim, 0,0,0, 0,0,0,&dof,0,0,0);CHKERRQ(ierr);
  if (dim != 3) SETERRQ(PETSC_ERR_ARG_WRONG,"Expected DA to be 3D");

  ierr = DAGetCorners(da,0,0,0,&zm,&ym,&xm);CHKERRQ(ierr);
  ierr = DAGetISLocalToGlobalMapping(da,&ltog);CHKERRQ(ierr);
  ierr = DAGetISLocalToGlobalMappingBlck(da,&ltogb);CHKERRQ(ierr);
  ierr = MatCreate(((PetscObject)da)->comm,&A);CHKERRQ(ierr);
  ierr =
MatSetSizes(A,dof*xm*ym*zm,dof*xm*ym*zm,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
  ierr = MatSetType(A,mtype);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatSeqAIJSetPreallocation(A,17,PETSC_NULL);CHKERRQ(ierr);
  ierr =
MatMPIAIJSetPreallocation(A,17,PETSC_NULL,12,PETSC_NULL);CHKERRQ(ierr);
  ierr = MatSeqBAIJSetPreallocation(A,dof,7,PETSC_NULL);CHKERRQ(ierr);
  ierr =
MatMPIBAIJSetPreallocation(A,dof,7,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr);
  ierr = MatSeqSBAIJSetPreallocation(A,dof,4,PETSC_NULL);CHKERRQ(ierr);
  ierr =
MatMPISBAIJSetPreallocation(A,dof,4,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr);
  ierr = MatSetDA(A,da);
  ierr = MatSetFromOptions(A);
  ierr = MatGetType(A,&Atype);
  ierr = MatSetBlockSize(A,dof);CHKERRQ(ierr);
  ierr = MatSetLocalToGlobalMapping(A,ltog);CHKERRQ(ierr);
  ierr = MatSetLocalToGlobalMappingBlock(A,ltogb);CHKERRQ(ierr);
  ierr =
DAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);CHKERRQ(ierr);
  ierr = MatSetStencil(A,dim,dims,starts,dof);CHKERRQ(ierr);
  *J = A;
  PetscFunctionReturn(0);
}