DOLFIN

Extract the diagonal of a PETSc matrix.

Asked by Lizao (Larry) Li

Let A be a Dolfin or PETSc matrix. Is there a way to implement D=Diag(A) to extract its diagonal which works in parallel and does not resort to numpy?

I know this can be done with the Epetra backend. But I am currently having some trouble with trilinos (compiling and parallel performance). So I want to stick with PETSc.

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Johan Hake
Solved:
Last query:
Last reply:
Revision history for this message
Johan Hake (johan-hake) said : 		#1

There are no support for dealing with diagonals in the GenericMatrix
interface. However, I have fiddled around with this a while ago.

  https://code.launchpad.net/~johan-hake/dolfin/dolfin-diagonal

Not sure I implemented exactly what you need but it should be easy to
add something like:

  shared_ptr<GenericVector> dolfin::GenericMatrix::diag() const;

I never pushed it as it turned out I did not need it and I did not
bother implementing the unit tests.

Johan

On 10/09/2012 10:21 PM, Lizao (Larry) Li wrote:
> New question #210777 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/210777
>
> Let A be a Dolfin or PETSc matrix. Is there a way to implement
> D=Diag(A) to extract its diagonal which works in parallel and does
> not resort to numpy?
>
> I know this can be done with the Epetra backend. But I am currently
> having some trouble with trilinos (compiling and parallel
> performance). So I want to stick with PETSc.
>

Revision history for this message
Lizao (Larry) Li (creatorlarryli) said : 		#2

Thank you for the answer. I look through the code. It seems it only sets the diagonal or multiplies the matrix by the diagonal matrix.

Given a matrix A, I need something like A.diagpart() which returns a PETSc matrix contains the diagonal of A. The inverse of the diagonal can be nice to have as a part of some more complicated preconditioners. In PETSc, we have the function MatGetDiagonalBlock (or MatGetDiagonalBlock and MatDiagonalSet).

 I can extract using Trilinos because downcast(A) gives a PyTrilinos matrix. But this does not work using the PETSc backend. I see there is a blue print for petsc4py (https://blueprints.launchpad.net/dolfin/+spec/add-petsc4py-support). It seems that interface is not there yet though.

Revision history for this message
Best Johan Hake (johan-hake) said : 		#3

On 10/10/2012 09:55 PM, Lizao (Larry) Li wrote:
> Question #210777 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/210777
>
> Status: Answered => Open
>
> Lizao (Larry) Li is still having a problem:
> Thank you for the answer. I look through the code. It seems it only sets
> the diagonal or multiplies the matrix by the diagonal matrix.

Yes, but it uses each backend to do the lifting. I also mentioned that
it does not do what you want, rather the opposite as you point out, but
what you want should be easy to implement along the same lines.

> Given a matrix A, I need something like A.diagpart() which returns a
> PETSc matrix contains the diagonal of A. The inverse of the diagonal can
> be nice to have as a part of some more complicated preconditioners. In
> PETSc, we have the function MatGetDiagonalBlock (or MatGetDiagonalBlock
> and MatDiagonalSet).

Yes, and these should be used for the PETSc backend.

> I can extract using Trilinos because downcast(A) gives a PyTrilinos
> matrix. But this does not work using the PETSc backend. I see there is a
> blue print for petsc4py (https://blueprints.launchpad.net/dolfin/+spec
> /add-petsc4py-support). It seems that interface is not there yet though.

No nothing has yet been done for this blueprint.

Johan

Revision history for this message
Lizao (Larry) Li (creatorlarryli) said : 		#4

Thanks Johan Hake, that solved my question.