<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Hi Junchao,</p>
<p>This seems to work.</p>
<p>Thanks</p>
<p>Elias<br>
</p>
<div class="moz-cite-prefix">On 08/04/2020 19:02, Junchao Zhang
wrote:<br>
</div>
<blockquote type="cite"
cite="mid:CA+MQGp82YhZoSy4hc2CtBHWrdnqRqVcgz0rNM02igT5Q3G42wQ@mail.gmail.com">
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<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>
<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>
<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" moz-do-not-send="true">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"
moz-do-not-send="true">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"
moz-do-not-send="true">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" moz-do-not-send="true">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>
</blockquote>
</body>
</html>