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

[Elias Karabelas]

Dear Jed,

Thanks I will try to keep you updated. So I managed to get my stuff 
running with one process but it fails with more than one. The matrix 
assembly as you described it works fine but I forgot that in the FCT 
algorithm one also needs to assemble a matrix D that looks kinda like

D_ij = L_ij * (u[i] - u[j])

where u is a vector with solution data from a previous time step. So the 
problem is that I have to now access off-process values of the Vector u, 
and doing this for each row of D. My basic first idea was to utilize the 
VecScatter Example from the Petsc Manual to get the off-process values 
per row of L to build this matrix D. But somehow this does not really 
work. Any ideas on how to scatter the Vector-data so that I can get the 
correct Vector values ?

At the moment I'm using something along this lines (in a loop where I 
get the local array of u with VecGetArray --> u_arr)

MatGetRow(L, row, &ncolsL, &colsL, &valsL);

double * Dvals = new double[ncolsM]();
int * idx_to = new int[ncolsL]();
for(int k=0; k < ncolsM; k++)

     idx_to[k] = k;
MPI_Barrier(PETSC_COMM_WORLD);
//Scatter the right stuff
Vec x;
VecScatter scatter;
IS from, to;
double *vals;
VecCreateSeq(PETSC_COMM_SELF, ncolsL, &x);
ISCreateGeneral(PETSC_COMM_SELF,ncolsL,colsL,PETSC_COPY_VALUES,&from);
ISCreateGeneral(PETSC_COMM_SELF,ncolsL,idx_to,PETSC_COPY_VALUES,&to);
VecScatterCreate(u,from,x,to,&scatter);
VecScatterBegin(scatter,u,x,INSERT_VALUES,SCATTER_FORWARD);
VecScatterEnd(scatter,u,x,INSERT_VALUES,SCATTER_FORWARD);
VecGetArray(x,&vals);
for(int c=0; c < ncolsL; c++)
{
     Dvals[c] = valsL[c] * (vals[c] -  u_arr[row])

}

MatSetValues(D, 1, &row, ncolsL, colsL, Dvals, ADD_VALUES);

and so on...

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