[petsc-users] Assembling primal Schur matrix in FETI-DP method

[Barry Smith]

On Nov 21, 2011, at 9:11 AM, Jed Brown wrote:

> On Mon, Nov 21, 2011 at 09:00, Thomas Witkowski <Thomas.Witkowski at tu-dresden.de> wrote:
> I have to multiple two matrices, one is implicitly defined by a KSP, the other one is explicitly assembled. So when using KSPSolve(), it must be called for each column of the second matrix? This seems not to be very efficient as I thought the matrices are stored row wise? Is there a way to get the columns of a SEQAIJ matrix?
> 
> In this context, I think it is easier to store the matrix as a collection of column vectors, but you can also use MATDENSE. This thing isn't sparse (and it wouldn't do you any good anyway; there are few algorithms that can solve sparse right hand sides and there isn't much point because the solutions would still be dense).


   From the FAQ I pointed you to earlier

How can I compute the Schur complement, Kbb - Kab * inverse(Kbb) * Kba in PETSc?

It is very expensive to compute the Schur complement of a matrix and very rarely needed in practice. We highly recommend avoiding  algorithms that need it. The Schur complement of a matrix (dense or sparse) is essentially always dense, so begin by
	• forming a dense matrix Kba, 
	• also create another dense matrix T of the same size. 
	• Then factor the matrix Kaa  with MatLUFactor() or MatCholeskyFactor(), call the result A.
	• Then call MatMatSolve(A,Kba,T). 
	• Then call MatMatMult(Kab,T,MAT_INITIAL_MATRIX,1.0,&S).
	• Now call MatAXPY(S,-1.0,Kbb,MAT_SUBSET_NONZERO).
	• Followed by MatScale(S,-1.0);

Note there is never a reason to use KSP to do the solve because for the many solves needed a direct solver will always win over using an iterative solver and since the result is dense it doesn't make sense to do this computation with huge matrices.

   Barry

If you want to use an external LU solver (it will not be faster than PETSc's so why bother). You would use MatGetFactor() then MatLUFactorSymbolic() followed by MatLUFactorNumeric().