On Feb 21, 2014, at 2:41 PM, Chung-Kan Huang <[email protected]> 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
>
>
