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

[Jed Brown]

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