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

[Elias Karabelas]

Hi Junchao,

This seems to work.

Thanks

Elias

On 08/04/2020 19:02, Junchao Zhang wrote:
> 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 <mailto: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
>     <mailto: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
>     <mailto: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/20200415/79dc3ebf/attachment-0001.html>