<div dir="ltr">Hi, Elias, <br> VecScatterToAll is implemented with MPI_Allgatherv. If not large scale, I guess it won't be a problem. I assume you want to assemble a MATMPIAIJ D with<br> D_ij = L_ij * (u[i] - u[j])<br> Since D has the same sparsity pattern as L, we may have (not tested),<div><br></div><div><div style="color:rgb(0,0,0);font-family:Menlo,Monaco,"Courier New",monospace;font-size:14px;line-height:21px;white-space:pre"><div> Mat A,B;</div><div> <span style="color:rgb(0,0,255)">const</span> PetscInt *garray,*cols;</div><div> VecScatter vscat;</div><div> PetscInt i,j,m,n,ncols;</div><div> Vec ur;</div><div> PetscScalar *ulocal,*uremote,val,*vals;</div><div> IS from,to;</div><br><div> ierr = MatMPIAIJGetSeqAIJ(L,&A,&B,&garray);CHKERRQ(ierr);</div><div> ierr = MatGetSize(B,&<span style="color:rgb(0,0,255)">NULL</span>,&n);CHKERRQ(ierr); <span style="color:rgb(0,128,0)">/* garray[]'s length = n */</span></div><div> ierr = VecCreateSeq(PETSC_COMM_SELF,n,&ur);CHKERRQ(ierr); <span style="color:rgb(0,128,0)">/* ur stores needed off-proc entries of u */</span></div><div> ierr = ISCreateStride(PETSC_COMM_SELF,n,<span style="color:rgb(9,134,88)">0</span>,<span style="color:rgb(9,134,88)">1</span>,&to);</div><div> ierr = ISCreateGeneral(PETSC_COMM_SELF,n,garray,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);</div><div> ierr = VecScatterCreate(u,from,ur,to,&vscat);CHKERRQ(ierr); <span style="color:rgb(0,128,0)">/* vscat is D's Mvctx, which however is not exposed to users */</span></div><div> ierr = VecScatterBegin(vscat,u,ur,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);</div><div> ierr = VecScatterEnd(vscat,u,ur,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);</div><br><div> ierr = VecGetArrayRead(u,&ulocal);CHKERRQ(ierr);</div><div> ierr = VecGetArrayRead(ur,&uremote);CHKERRQ(ierr);</div><div> ierr = MatGetOwnershipRange(D,&rstart,<span style="color:rgb(0,0,255)">NULL</span>);CHKERRQ(ierr);</div><div> ierr = MatGetOwnershipRangeColumn(D,&cstart,<span style="color:rgb(0,0,255)">NULL</span>);</div><div> ierr = MatDuplicate(L,&D,MAT_DO_NOT_COPY_VALUES);CHKERRQ(ierr);</div><div> ierr = MatSetOption(D,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);</div><div> ierr = MatAssemblyBegin(D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);</div><div> ierr = MatGetSize(A,&m,<span style="color:rgb(0,0,255)">NULL</span>);CHKERRQ(ierr);</div><div> <span style="color:rgb(0,0,255)">for</span> (i=<span style="color:rgb(9,134,88)">0</span>; i<m; i++) {</div><div> ierr = MatGetRow(A,i,ncols,cols,vals);</div><div> <span style="color:rgb(0,0,255)">for</span> (j=<span style="color:rgb(9,134,88)">0</span>; j<ncols; j++) {</div><div> val = vals[j]*(ulocal[i] - ulocal[cols[j]]);</div><div> ierr = MatSetValue(D,rstart+i,cstart+cols[j],val,INSERT_VALUES);CHKERRQ(ierr);</div><div> }</div><div> ierr = MatGetRow(B,i,ncols,cols,vals);</div><div> <span style="color:rgb(0,0,255)">for</span> (j=<span style="color:rgb(9,134,88)">0</span>; j<ncols; j++) {</div><div> val = vals[j]*(ulocal[i] - uremote[cols[j]]);</div><div> ierr = MatSetValue(D,rstart+i,garray[cols[j]],val,INSERT_VALUES);CHKERRQ(ierr);</div><div> }</div><div> }</div><div> ierr = MatAssemblyEnd(D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);</div><div> ierr = VecRestoreArrayRead(u,&ulocal);CHKERRQ(ierr);</div><div> ierr = VecRestoreArrayRead(ur,&uremote);CHKERRQ(ierr);</div></div><div><br></div><div><div><div><div>--Junchao Zhang<br></div></div></div></div></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Apr 8, 2020 at 2:17 AM Elias Karabelas <<a href="mailto:karabelaselias@gmail.com" target="_blank">karabelaselias@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Dear Jed,<br>
<br>
I'm done implementing my FCT-Solver and it works fairly well on a small <br>
amount of MPI-Procs. Additionally to your little snippet I have used a <br>
VecScatterToAll.<br>
<br>
Reason is that the flux correction f_i takes the form Sum_{j} alpha_ij <br>
r_ij where r_ij is defined on the sparsity pattern of my FEM Matrix and <br>
alpha_ij is based on two vectors Rp and Rm. So basically I need <br>
off-process values of these vectors to construct alpha which I made with <br>
a VecScatterToAll. However I guess this will slow down my overall <br>
program quite significant. Any Ideas?<br>
<br>
Best regards<br>
<br>
Elias<br>
<br>
On 23/03/2020 15:53, Jed Brown wrote:<br>
> Thanks; please don't drop the list.<br>
><br>
> I'd be curious whether this operation is common enough that we should<br>
> add it to PETSc. My hesitance has been that people may want many<br>
> different variants when working with systems of equations, for example.<br>
><br>
> Elias Karabelas <<a href="mailto:karabelaselias@gmail.com" target="_blank">karabelaselias@gmail.com</a>> writes:<br>
><br>
>> Dear Jed,<br>
>><br>
>> Yes the Matrix A comes from assembling a FEM-convection-diffusion<br>
>> operator over a tetrahedral mesh. So my matrix graph should be<br>
>> symmetric. Thanks for the snippet<br>
>><br>
>> On 23/03/2020 15:42, Jed Brown wrote:<br>
>>> Elias Karabelas <<a href="mailto:karabelaselias@gmail.com" target="_blank">karabelaselias@gmail.com</a>> writes:<br>
>>><br>
>>>> Dear Users,<br>
>>>><br>
>>>> I want to implement a FCT (flux corrected transport) scheme with PETSc.<br>
>>>> To this end I have amongst other things create a Matrix whose entries<br>
>>>> are given by<br>
>>>><br>
>>>> L_ij = -max(0, A_ij, A_ji) for i neq j<br>
>>>><br>
>>>> L_ii = Sum_{j=0,..n, j neq i} L_ij<br>
>>>><br>
>>>> where Mat A is an (non-symmetric) Input Matrix created beforehand.<br>
>>>><br>
>>>> I was wondering how to do this. My first search brought me to<br>
>>>> <a href="https://www.mcs.anl.gov/petsc/petsc-current/src/mat/examples/tutorials/ex16.c.html" rel="noreferrer" target="_blank">https://www.mcs.anl.gov/petsc/petsc-current/src/mat/examples/tutorials/ex16.c.html</a><br>
>>>><br>
>>>><br>
>>>> but this just goes over the rows of one matrix to set new values and now<br>
>>>> I would need to run over the rows and columns of the matrix. My Idea was<br>
>>>> to just create a transpose of A and do the same but then the row-layout<br>
>>>> will be different and I can't use the same for loop for A and AT and<br>
>>>> thus also won't be able to calculate the max's above.<br>
>>> Does your matrix have symmetric nonzero structure? (It's typical for<br>
>>> finite element methods.)<br>
>>><br>
>>> If so, all the indices will match up so I think you can do something like:<br>
>>><br>
>>> for (row=rowstart; row<rowend; row++) {<br>
>>> PetscScalar Lvals[MAX_LEN];<br>
>>> PetscInt diag;<br>
>>> MatGetRow(A, row, &ncols, &cols, &vals);<br>
>>> MatGetRow(At, row, &ncolst, &colst, &valst);<br>
>>> assert(ncols == ncolst); // symmetric structure<br>
>>> PetscScalar sum = 0;<br>
>>> for (c=0; c<ncols; c++) {<br>
>>> assert(cols[c] == colst[c]); // symmetric structure<br>
>>> if (cols[c] == row) diag = c;<br>
>>> else sum -= (Lvals[c] = -max(0, vals[c], valst[c]));<br>
>>> }<br>
>>> Lvals[diag] = sum;<br>
>>> MatSetValues(L, 1, &row, ncols, cols, Lvals, INSERT_VALUES);<br>
>>> MatRestoreRow(A, row, &ncols, &cols, &vals);<br>
>>> MatRestoreRow(At, row, &ncolst, &colst, &valst);<br>
>>> }<br>
</blockquote></div>