[firedrake] square root of matrix

[Patrick Farrell]

On 14/02/17 08:21, Colin Cotter wrote:
> Dear all,
>   For weird reasons I'd like to solve Ux = y where U is the Cholesky square root of a given symmetric positive definite matrix A. Is there a way to abuse the
> PETSc solver interface into doing this?
>

Hi Colin,

I have partial information about your question.

I need to access the Cholesky factorisation data to verify
the second-order sufficiency conditions for constrained minimisation
problems. (I'm using dolfin, but it shouldn't matter very much
for your purposes.) My petsc4py code looks like

         pc = PETSc.PC().create(comm)
         pc.setOperators(A.mat())
         pc.setType("cholesky")
         pc.setFactorSolverPackage("mumps")
         pc.setUp()

to calculate the Cholesky factorisation, then

         F = pc.getFactorMatrix()
         (neg, zero, pos) = F.getInertia()

to compute the inertia of the matrix. (I also need
some MUMPS options, but I presume you don't care about
these, since you're not asking it to compute the inertia.)

Maybe you can get the sparse matrix representation out of
F and solve a linear system with it, who knows.

Cheerio,

Patrick