[petsc-users] DMCreateMatrix from one DMDA to an another

[Dave May]

On 6 December 2015 at 11:19, Timothée Nicolas <timothee.nicolas at gmail.com>
wrote:

> Hi,
>
> The answer to the first question is yes. I was puzzled about the second
> part of the answer, but I realize I was not clear enough. I don't want to
> just transfer information between 2 DMDAs, I want to perform an operation
> on vectors with a rectangular matrix.
>

Okay - I misunderstood.



> To be more explicit, the input vector of the operation is a vector with 3
> components (hence dof=3), and the results of the operator is put in one
> vector equation + 2 scalar equation, which makes it dof = 5.
>
> So, either I cannot fill my matrix (if the number of rows exceeds what is
> allowed by the DMDA), or I can fill my matrix but then the matrix and
> vectors are not compatible, as in this error :
>


Do you need to explicitly assemble this rectangular operator?
Expressing this operator via a matrix-free representation would be
relatively straight forward.

However, since the DMDA's have the same parallel layout you are in good
shape to define the rectangular matrix if you really want the assembled
thing. Your row space would be defined by the DMDA with block size 5 and
the column space would be defined by the DMDA with block size 3. For
example the matrix would be created like in the following way (assuming a
2D DMDA)

PetscInt m,n,M,N;

DMDAGetInfo(dm3,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
nodes = M *N;
DMDAGetCorners(dm3,NULL,NULL,NULL,&m,&n,NULL);
lnodes = m * n;
MatCreate(PETSC_COMM_WORLD,&A);
MatSetSizes(A,lnodes*5,lnodes*3,nodes*5,nodes*3);

The sparsity pattern is defined by the connectivity of the DMDA points AND
the different block sizes.
To define the non-zero structure, you need to traverse the part of the dmda
defined by the natural indices given by si <= i < si+m, sj <= j < sj+n and
check whether the the 5x3 block associated with each i,j is living within
the diagonal or off-diagonal block of the matrix. For a given i,j the
off-diagonal is defined by any i',j' associated your stencil which is not
in the range si <= i' <= si + m, sj <= j' < sj + n. (e.g. it is a ghost
node).

Your arrays counting the non-zero entries on the diagonal (nnz[]) and off
diag (nnz_od[]) can be indexed using i + j * m. This corresponds to the
indexing used by the DMDA. Your non-zero counting function would look
something like the pseudo code below

PetscInt *nnz,*nnz_od;

PetscMalloc(sizeof(PetscInt)*m*n*5,&nnz);
PetscMalloc(sizeof(PetscInt)*m*n*5,&nnz_od);

PetscMemzero(nnz,sizeof(PetscInt)*m*n*5);
PetscMemzero(nnz_od,sizeof(PetscInt)*m*n*5);

for (j=0; j<n; j++) {
for (i=0; i<m; i++) {
  PetscInt idx,d,cnti,cntg;
  PetscBool is_ghost;

  idx = i + j * m;

  cnti = 0;
  cntg = 0;
  for (all_points_in_stencil) { /* User logic to define this loop */

   is_ghost = PETSC_FALSE;
   /*User logic goes here to test if stencil point is a ghost point */

    /* accumulate counts for the stencil */
    if (is_ghost) { cntg++; } /* values living on a neighbour rank */
    else { cnti++; } /* values local to current rank */
  }

  /* assume all dofs in row and col space in block are coupled to each
other */
  for (d=0; d<5; d++) {
    nnz[5*idx+d]       = 3 * cnti;
    nnz_od[5*idx+d] = 3 * cntg;
  }

}
}

Thanks,
  Dave


>
>
> [0]PETSC ERROR: --------------------- Error Message
> --------------------------------------------------------------
> [0]PETSC ERROR: Nonconforming object sizes
> [0]PETSC ERROR: Mat mat,Vec y: global dim 4900 2940
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
> for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015
> [0]PETSC ERROR: ./miips on a arch-linux2-c-debug named Carl-9000 by
> timothee Sun Dec  6 19:17:20 2015
> [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++
> --with-fc=gfortran --download-fblaslapack --download-mpich
> [0]PETSC ERROR: #1 MatMult() line 2215 in
> /home/timothee/Documents/petsc-3.6.1/src/mat/interface/matrix.c
>
> So I assume it means the matrix is assumed to be square. Is there a way to
> change this ?
>
> Best
>
> Timothée
>
>
>
> 2015-12-05 20:41 GMT+09:00 Dave May <dave.mayhem23 at gmail.com>:
>
>>
>>
>> On 5 December 2015 at 11:15, Timothée Nicolas <timothee.nicolas at gmail.com
>> > wrote:
>>
>>> Hi,
>>>
>>> I need to create a rectangular matrix which takes vector structured by a
>>> DMDA and transforms them to vectors structured with an other DMDA. It does
>>> not look like there is a routine to do this since apparently DMCreateMatrix
>>> assumes square matrix. What is the clean way to do this ?
>>>
>>
>> Do the DMDA's have exactly the same parallel layout and only differ by
>> the number of DOFs?
>>
>> If yes, and you only want to move entries between vectors associated with
>> the two DMDA's, you can perform the transfer between vectors locally
>> (rank-wise) in a very simple manner using the local number of points from
>> DMDA and the known block sizes (DOFs) associated with your DMDA and the
>> size of the subspace you ultimately want to construct.
>>
>> Thanks,
>>   Dave
>>
>>
>>> The problem is I have a DMDA with 8 degrees of freedom (dof) and my
>>> matrix sends vectors living in a subspace with 3 dof to the other 5 dof
>>> subspace. I can define a square matrix living in the large, 8 dof subspace,
>>> but this is of course very inefficient in terms of memory and time.
>>>
>>> Thanks
>>>
>>> Best
>>>
>>> Timothée
>>>
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151206/8eab51fa/attachment-0001.html>