[petsc-users] Construct Matrix based on row and column values

[Jed Brown]

Elias Karabelas <karabelaselias at gmail.com> writes:

> Dear Jed,
>
> I'm done implementing my FCT-Solver and it works fairly well on a small 
> amount of MPI-Procs. Additionally to your little snippet I have used a 
> VecScatterToAll.

Can you share the code including that VecScatterToAll?

There more local ways to get that information.

> Reason is that the flux correction f_i takes the form Sum_{j} alpha_ij 
> r_ij where r_ij is defined on the sparsity pattern of my FEM Matrix and 
> alpha_ij is based on two vectors Rp and Rm. So basically I need 
> off-process values of these vectors to construct alpha which I made with 
> a VecScatterToAll. However I guess this will slow down my overall 
> program quite significant. Any Ideas?
>
> Best regards
>
> Elias
>
> On 23/03/2020 15:53, Jed Brown wrote:
>> Thanks; please don't drop the list.
>>
>> I'd be curious whether this operation is common enough that we should
>> add it to PETSc.  My hesitance has been that people may want many
>> different variants when working with systems of equations, for example.
>>
>> Elias Karabelas <karabelaselias at gmail.com> writes:
>>
>>> Dear Jed,
>>>
>>> Yes the Matrix A comes from assembling a FEM-convection-diffusion
>>> operator over a tetrahedral mesh. So my matrix graph should be
>>> symmetric. Thanks for the snippet
>>>
>>> On 23/03/2020 15:42, Jed Brown wrote:
>>>> Elias Karabelas <karabelaselias at gmail.com> writes:
>>>>
>>>>> Dear Users,
>>>>>
>>>>> I want to implement a FCT (flux corrected transport) scheme with PETSc.
>>>>> To this end I have amongst other things create a Matrix whose entries
>>>>> are given by
>>>>>
>>>>> L_ij = -max(0, A_ij, A_ji) for i neq j
>>>>>
>>>>> L_ii = Sum_{j=0,..n, j neq i} L_ij
>>>>>
>>>>> where Mat A is an (non-symmetric) Input Matrix created beforehand.
>>>>>
>>>>> I was wondering how to do this. My first search brought me to
>>>>> https://www.mcs.anl.gov/petsc/petsc-current/src/mat/examples/tutorials/ex16.c.html
>>>>>
>>>>>
>>>>> but this just goes over the rows of one matrix to set new values and now
>>>>> I would need to run over the rows and columns of the matrix. My Idea was
>>>>> to just create a transpose of A and do the same but then the row-layout
>>>>> will be different and I can't use the same for loop for A and AT and
>>>>> thus also won't be able to calculate the max's above.
>>>> Does your matrix have symmetric nonzero structure?  (It's typical for
>>>> finite element methods.)
>>>>
>>>> If so, all the indices will match up so I think you can do something like:
>>>>
>>>> for (row=rowstart; row<rowend; row++) {
>>>>     PetscScalar Lvals[MAX_LEN];
>>>>     PetscInt diag;
>>>>     MatGetRow(A, row, &ncols, &cols, &vals);
>>>>     MatGetRow(At, row, &ncolst, &colst, &valst);
>>>>     assert(ncols == ncolst); // symmetric structure
>>>>     PetscScalar sum = 0;
>>>>     for (c=0; c<ncols; c++) {
>>>>       assert(cols[c] == colst[c]); // symmetric structure
>>>>       if (cols[c] == row) diag = c;
>>>>       else sum -= (Lvals[c] = -max(0, vals[c], valst[c]));
>>>>     }
>>>>     Lvals[diag] = sum;
>>>>     MatSetValues(L, 1, &row, ncols, cols, Lvals, INSERT_VALUES);
>>>>     MatRestoreRow(A, row, &ncols, &cols, &vals);
>>>>     MatRestoreRow(At, row, &ncolst, &colst, &valst);
>>>> }