[petsc-dev] Branches with additional features

[Jed Brown]

Pierre Jolivet <Pierre.Jolivet at enseeiht.fr> writes:

> On Mon, 29 May 2017 04:38:39 -0600, Jed Brown wrote:
>> Pierre Jolivet <pierre.jolivet at enseeiht.fr> writes:
>>
>>>> On 29 May 2017, at 06:49, Jed Brown <jed at jedbrown.org> wrote:
>>>>
>>>> You can't get a callback from this "external Krylov solver"?
>>>
>>> You mean, a callback to let the external solver use the memory 
>>> allocated by PETSc for the Mats? Unfortunately, I need the underlying 
>>> memory to be contiguous between all Mats, and I don't think it is 
>>> possible with the current PETSc operators for matrix creation.
>>
>> You can allocate the full block, then pass portions of the
>> (column-aligned) array in MatCreateDense().
>
> So I was doing:
> MatCreate();
> MatSetSizes();
> MatSetType();
> MatMPIDenseSetPreallocation();
> MatAssemblyBegin();
> MatAssemblyEnd();
> MatMatMult();
>
> and what you are suggesting is:
> MatCreateDense(); // data != nullptr as provided by the external solver 
> so that there is no need for MatAssemblyBegin/MatAssemblyEnd 
> http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/dense/mpi/mpidense.c.html#line1473 
> ?
> MatMatMult();

PetscErrorCode  MatCreateDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscScalar *data,Mat *A)
{
  PetscErrorCode ierr;
  PetscMPIInt    size;

  PetscFunctionBegin;
  ierr = MatCreate(comm,A);CHKERRQ(ierr);
  ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr);
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
  if (size > 1) {
    ierr = MatSetType(*A,MATMPIDENSE);CHKERRQ(ierr);
    ierr = MatMPIDenseSetPreallocation(*A,data);CHKERRQ(ierr);
    if (data) {  /* user provided data array, so no need to assemble */
      ierr = MatSetUpMultiply_MPIDense(*A);CHKERRQ(ierr);
      (*A)->assembled = PETSC_TRUE;
    }
  } else {
    ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr);
    ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

> That makes perfect sense (assuming I'm not wrong about the 
> MatAssemblyBegin/MatAssemblyEnd). MatSetUpMultiply_MPIDense would still 
> be called at each iteration but I doubt this is too costly.

It creates a VecScatter so it isn't nothing (in terms of parallel
semantics), but I'd like to see profiling data before chasing this
around.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 832 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20170529/77302172/attachment.sig>