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

[Junchao Zhang]

Hi, Elias,
   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
      D_ij = L_ij * (u[i] - u[j])
  Since D has the same sparsity pattern as L, we may have (not tested),

Mat A,B;
const PetscInt *garray,*cols;
VecScatter vscat;
PetscInt i,j,m,n,ncols;
Vec ur;
PetscScalar *ulocal,*uremote,val,*vals;
IS from,to;

ierr = MatMPIAIJGetSeqAIJ(L,&A,&B,&garray);CHKERRQ(ierr);
ierr = MatGetSize(B,&NULL,&n);CHKERRQ(ierr); /* garray[]'s length = n */
ierr = VecCreateSeq(PETSC_COMM_SELF,n,&ur);CHKERRQ(ierr); /* ur stores
needed off-proc entries of u */
ierr = ISCreateStride(PETSC_COMM_SELF,n,0,1,&to);
ierr =
ISCreateGeneral(PETSC_COMM_SELF,n,garray,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
ierr = VecScatterCreate(u,from,ur,to,&vscat);CHKERRQ(ierr); /* vscat is D's
Mvctx, which however is not exposed to users */
ierr =
VecScatterBegin(vscat,u,ur,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr =
VecScatterEnd(vscat,u,ur,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);

ierr = VecGetArrayRead(u,&ulocal);CHKERRQ(ierr);
ierr = VecGetArrayRead(ur,&uremote);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(D,&rstart,NULL);CHKERRQ(ierr);
ierr = MatGetOwnershipRangeColumn(D,&cstart,NULL);
ierr = MatDuplicate(L,&D,MAT_DO_NOT_COPY_VALUES);CHKERRQ(ierr);
ierr =
MatSetOption(D,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
ierr = MatAssemblyBegin(D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatGetSize(A,&m,NULL);CHKERRQ(ierr);
for (i=0; i<m; i++) {
ierr = MatGetRow(A,i,ncols,cols,vals);
for (j=0; j<ncols; j++) {
val = vals[j]*(ulocal[i] - ulocal[cols[j]]);
ierr =
MatSetValue(D,rstart+i,cstart+cols[j],val,INSERT_VALUES);CHKERRQ(ierr);
}
ierr = MatGetRow(B,i,ncols,cols,vals);
for (j=0; j<ncols; j++) {
val = vals[j]*(ulocal[i] - uremote[cols[j]]);
ierr =
MatSetValue(D,rstart+i,garray[cols[j]],val,INSERT_VALUES);CHKERRQ(ierr);
}
}
ierr = MatAssemblyEnd(D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = VecRestoreArrayRead(u,&ulocal);CHKERRQ(ierr);
ierr = VecRestoreArrayRead(ur,&uremote);CHKERRQ(ierr);

--Junchao Zhang

On Wed, Apr 8, 2020 at 2:17 AM Elias Karabelas <karabelaselias at gmail.com>
wrote:

> 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);
> >>> }
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200408/06e3b7dd/attachment.html>