[petsc-users] Using Pmat different from Amat

[Barry Smith]

On Feb 21, 2014, at 2:41 PM, Chung-Kan Huang <ckhuangf at gmail.com> wrote:

> Hello,
> 
> In my application I like to use 
> PetscErrorCode  KSPSetOperators(KSP ksp,Mat Amat,Mat Pmat,MatStructure flag)
> and have Pmat different from Amat
> 
> if Amat = L + D + U
> then Pmat  = Amat  - L* - U* + rowsum(L* + U*)
> where L* and U* is a portion of L and U and rowsum(L* + U*) will add into D
> 
> I know I can explicitly construct Pmat step-by-step but wonder what will be the most easy way to do this in PETSC?

   Kan,

   This is something that is trivial in Matlab but not trivial to do efficiently in PETSc. 

   How do you determine what “portion” of L and U are used? Does it depend on numerical values? Is if fixed? 

   Is Amat an AIJ matrix or a BAIJ matrix?

   Do you want to do it in parallel? 

   I would start by doing it sequentially for a SeqAIJ matrix and #include "src/mat/impls/aij/seq/aij.h” into the source code then you would write a routine that looped over the Amat rows accessing values directly in the CSR storage format 

 PetscInt          *i;               /* pointer to beginning of each row */               \
  PetscInt          *j;               /* column values: j + i[k] - 1 is start of row k */  \

  datatype          *a;               /* nonzero elements */     

  you would then determine how many non zeros you wanted per row, call MatCreateSeqAIJ() with the preallocation from the non zeros per row and then loop over the rows of Amat again calling MatSetValues() to move the values you want moved into the Pmat matrix and added up for the diagonal.

  It is not terribly more difficult for MPIAIJ except that the matrix entries in Amat are stored in two parts so accessing those values values requires a bit more complicated code and I would not recommend doing it until you have the sequential version working perfectly.

  Barry



      


> 
> Thanks,
> 
> 
> Kan
> 
>