[petsc-users] Augmented Linear System

[Karabelas, Elias (elias.karabelas@uni-graz.at)]

Dear Community,

I have a linear system of the form

|K B| du    f1

        =

|C D| dp    f2

where K is a big m x m sparsematrix that comes from some FE 
discretization, B is a coupling matrix (of the form m x 4), C is of the 
form (4 x m) and D is 4x4.

I save B and C as 4 Vecs and D as a 4x4 double array. D might be 
singular so at the moment I use the following schur-complement approach 
to solve this system

1) Calculate the vecs v1 = KSP(K,PrecK) * f1, invB = [ KSP(K, PrecK) * 
B[0], KSP(K, PrecK) * B[1], KSP(K, PrecK) * B[2], KSP(K, PrecK) * B[3] ]

2) Build the schurcomplement S=[C ^ K^-1 B - D] via VecDots (C ^ K^-1 B 
[i, j] = VecDot(C[i], invB[j])

3) invert S (this seems to be mostly non-singular) to get dp

4) calculate du with dp

So counting this, I need 5x to call KSP which can be expensive and I 
thought of somehow doing the Schur-Complement the other way around, 
however due to the (possible) singularity of D this seems like a bad 
idea (even using a pseudoinverse)

Two things puzzle me here still

A) Can this be done more efficiently?

B) In case my above matrix is the Jacobian in a newton method, how do I 
make sure with any form of Schur Complement approach that I hit the 
correct residual reduction?

Thanks

Elias

-- 
Dr. Elias Karabelas
Universitätsassistent | PostDoc

Institut für Mathematik & Wissenschaftliches Rechnen | Institute of Mathematics & Scientific Computing
Universität Graz | University of Graz

Heinrichstraße 36, 8010 Graz
Tel.:   +43/(0)316/380-8546
E-Mail: elias.karabelas at uni-graz.at
Web: 	https://ccl.medunigraz.at
  
Bitte denken Sie an die Umwelt, bevor Sie dieses E-Mail drucken. Danke!
Please consider the environment before printing this e-mail. Thank you!