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

[Elias Karabelas]

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.

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);
>>> }