<div dir="ltr"><div><div><div>Wow, Great answer, thanks, I should probably be able to do it like this, I'll try my best !<br><br></div>I actually do want the assembled matrix now. I have coded the matrix-free version already and it works fine, but I want to use multigrid and I figured that in order to be be most efficient at the coarse level it is best to use -ksp_type preonly -pc_type lu, which requires a true matrix (the matrix I am talking about is the product between a matrix 3dof --> 5dof with a matrix 5dof --> 3dof, which ends up square). But I want to retain a memory light implementation so I use matrix-free formulation at all the levels but the last, which is also light by definition, so I should not suffer much in terms of memory.<br><br></div>Thanks<br><br></div>Timothée<br><br><br></div><div class="gmail_extra"><br><div class="gmail_quote">2015-12-06 22:29 GMT+09:00 Dave May <span dir="ltr"><<a href="mailto:dave.mayhem23@gmail.com" target="_blank">dave.mayhem23@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote"><span class="">On 6 December 2015 at 11:19, Timothée Nicolas <span dir="ltr"><<a href="mailto:timothee.nicolas@gmail.com" target="_blank">timothee.nicolas@gmail.com</a>></span> wrote:<br></span><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div><div><div><div>Hi,<br><br></div><span class="">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. </span></div></div></div></div></div></blockquote><div><br></div><div>Okay - I misunderstood.<br></div><span class=""><div><br> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div><div><div>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.<br><br></div>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 :<br></div></div></div></div></blockquote><div><br></div><br></span><div>Do you need to explicitly assemble this rectangular operator?<br></div><div>Expressing this operator via a matrix-free representation would be relatively straight forward.<br><br><div>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)<br></div><div><br></div><div>PetscInt m,n,M,N;<br></div><br>DMDAGetInfo(dm3,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);</div><div><div>nodes = M *N;<br>DMDAGetCorners(dm3,NULL,NULL,NULL,&m,&n,NULL);<br></div></div><div>lnodes = m * n;<br>MatCreate(PETSC_COMM_WORLD,&A);<br></div>MatSetSizes(A,lnodes*5,lnodes*3,nodes*5,nodes*3);<br><div><br>The sparsity pattern is defined by the connectivity of the DMDA points AND the different block sizes. <br></div><div>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).<br><br></div><div>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<br></div><div><br></div><div>PetscInt *nnz,*nnz_od;<br></div><div><br></div><div>PetscMalloc(sizeof(PetscInt)*m*n*5,&nnz);<br></div><div><div>PetscMalloc(sizeof(PetscInt)*m*n*5,&nnz_od);<br></div><br>PetscMemzero(nnz,sizeof(PetscInt)*m*n*5);</div><div>PetscMemzero(nnz_od,sizeof(PetscInt)*m*n*5);<div><br></div></div><div>for (j=0; j<n; j++) {<br></div><div>for (i=0; i<m; i++) {<br></div><div> PetscInt idx,d,cnti,cntg;<br></div><div> PetscBool is_ghost;<br></div><div><br></div><div> idx = i + j * m;<br></div><div><br></div><div><div> cnti = 0;<br></div> cntg = 0;<br></div><div> for (all_points_in_stencil) { /* User logic to define this loop */<br></div><div><br></div><div> is_ghost = PETSC_FALSE;<br></div><div> /*User logic goes here to test if stencil point is a ghost point */<br><br></div><div> /* accumulate counts for the stencil */<br></div><div> if (is_ghost) { cntg++; } /* values living on a neighbour rank */<br></div><div> else { cnti++; } /* values local to current rank */<br></div><div> }<br></div><div><br></div><div> /* assume all dofs in row and col space in block are coupled to each other */<br></div><div> for (d=0; d<5; d++) {<br></div><div> nnz[5*idx+d] = 3 * cnti;</div><div><div> nnz_od[5*idx+d] = 3 * cntg;</div> }<br><br>}<br>}<br></div><br><div>Thanks,<br></div><div> Dave<br></div><span class=""><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div><div><br><br>[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>[0]PETSC ERROR: Nonconforming object sizes<br>[0]PETSC ERROR: Mat mat,Vec y: global dim 4900 2940<br>[0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br>[0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015 <br>[0]PETSC ERROR: ./miips on a arch-linux2-c-debug named Carl-9000 by timothee Sun Dec 6 19:17:20 2015<br>[0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich<br>[0]PETSC ERROR: #1 MatMult() line 2215 in /home/timothee/Documents/petsc-3.6.1/src/mat/interface/matrix.c<br><br></div>So I assume it means the matrix is assumed to be square. Is there a way to change this ?<br><br></div>Best<span><font color="#888888"><br><br></font></span></div><span><font color="#888888">Timothée<br><div><div><div><br><br></div></div></div></font></span></div><div><div><div class="gmail_extra"><br><div class="gmail_quote">2015-12-05 20:41 GMT+09:00 Dave May <span dir="ltr"><<a href="mailto:dave.mayhem23@gmail.com" target="_blank">dave.mayhem23@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote"><span>On 5 December 2015 at 11:15, Timothée Nicolas <span dir="ltr"><<a href="mailto:timothee.nicolas@gmail.com" target="_blank">timothee.nicolas@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div><div><div>Hi,<br><br></div>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 ?<br></div></div></div></div></blockquote><div><br></div></span><div>Do the DMDA's have exactly the same parallel layout and only differ by the number of DOFs?<br></div><div><br>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. </div><div><br></div><div>Thanks,<br></div><div> Dave <br></div><span><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div><div><br></div><div>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.<br></div><div><br></div>Thanks<br><br></div>Best<span><font color="#888888"><br><br></font></span></div><span><font color="#888888">Timothée<br></font></span></div>
</blockquote></span></div><br></div></div>
</blockquote></div><br></div>
</div></div></blockquote></span></div><br></div></div>
</blockquote></div><br></div>