Use of static condensation to solve a linear system

When using the spectral element method (SEM) to solve partial differential equations (PDE), it is common for the basis function to have modes with support on boundaries (of elements) and modes with support only inside an element. This structure can be exploited when solving the linear system.

Here, as an exercise for implementing the solver using this technique, we will use it to solve a symmetric linear system that is decomposed as

$$ A\cdot x = \left[ \begin{matrix} A_{bb} & A_{bi} \\ A_{ib} & A_{ii} \\ \end{matrix}\right] \cdot \left\{ \begin{matrix} x_b \\ x_i \\ \end{matrix}\right\} = \left\{\begin{matrix} f_b \\ f_i \\\end{matrix}\right\} $$

where $A_{bb}$, $A_{bi}$, $A_{ib}$ and $A_{ii}$ are submatrices.

The second row can be solved for $x_i$:

$$ x_i = A_{ii}^{-1} f_i - A_{ii}^{-1} A_{ib} x_b $$

Substituting this on the first row:

$$ \left( A_{bb} - A_{bi}A_{ii}^{-1} A_{ib} \right) x_b = f_b - A_{bi}A_{ii}^{-1} f_i $$

This can be rewritten as

$$ A'_{bb} \cdot x_b = f'_b $$

where $A'_{bb} = A_{bb} - A_{bi}A_{ii}^{-1} A_{ib}$ and $f'_b = f_b - A_{bi}A_{ii}^{-1} f_i$.

With these transormations, if the linear system is symmetric and positive and definite, to solve a linear system, we need to use the following steps:

Compute $A_{ii}^{-1}$ (in reality do a Cholesky decomposition).
Compute $A_{bi} A_{ii}^{-1}$
Compute the Choleksy factorization of $A'_{bb} = A_{bb} - A_{bi}A_{ii}^{-1} \cdot A_{ib}$

To solve the linear system given a RHS (f):

Compute $f'_b = f_b - A_{bi}A_{ii}^{-1} f_i$.
Using the Cholesky factorization of $A'_{bb}$ compute $x_b$.
Now, using the Choleksy factorization of $A_{ii}$ compute $x_i = A_{ii}^{-1} f_i - A_{ii}^{-1} A_{ib} x_b$

A function to generate random symmetric matrices



In [ ]:

    
function random_matrix(n, symm=true, diag=4.0)
    A = randn(n,n)
    for i = 1:n
        A[i,i] += diag
    end
    
    if symm
        A = (A + A') / 2
    end
    
    return A
end

Creating the matrix



In [ ]:

    
N = 12
nb = 6
ni = N - nb



In [ ]:

    
A = [
 3.73295   0.351419  -0.534708   0.301579   0.558047;
  0.351419  6.48719    0.31895    0.198893   0.400289;
 -0.534708  0.31895    4.81739    0.309232  -0.103873;
  0.301579  0.198893   0.309232   5.411     -0.37619 ;
    0.558047  0.400289  -0.103873  -0.37619    3.44621 ]


A = random_matrix(N);

Obtaining the sub-matrices



In [ ]:

    
Abb = A[1:nb, 1:nb];



In [ ]:

    
Abi = A[1:nb, (nb+1):N];



In [ ]:

    
Aib = A[(nb+1):N, 1:nb];



In [ ]:

    
Aii = A[(nb+1):N, (nb+1):N];

Loading the appropriate modules



In [ ]:

    
using Base.LinAlg.BLAS.gemm!
using Base.LinAlg.BLAS.gemv!
using Base.LinAlg.LAPACK.potrf!
using Base.LinAlg.LAPACK.potrs!

Cholesky decompostion of $A_{ii}$



In [ ]:

    
potrf!('L', Aii);

Computing the matrix $A_{bi}A_{ii}^{-1}$. This is kindy of tricky. Remembering that $(A\cdot B)^T = A^T\cdot B^T$, and that the system is symmetric, if we compute $M = A_{ii}^{-1}\cdot A_{ib}$, then $M^T = \left(A_{ii}^{-1}\cdot A_{ib}\right)^T = A_{bi}\cdot A_{ii}^{-1}$



In [ ]:

    
M = copy(Aib);
potrs!('L', Aii, M);

Now we need to compute $A'_{bb} = A_{bb} - A_{bi}A_{ii}^{-1} A_{ib} = A_{bb} - M^T\cdot A_{ib}$. BLAS makes it simple:



In [ ]:

    
gemm!('T', 'N', -1.0, M, Aib, 1.0, Abb);

Cholesky decomposition of $A'_{bb}$:



In [ ]:

    
potrf!('L', Abb);

Solving for a RHS

This is just a test. Let's just create any RHS:



In [ ]:

    
f = [1.0:N;]



In [ ]:

    
fbi = copy(f);



In [ ]:

    
fb = view(fbi, 1:nb)
fi = view(fbi, (nb+1):N);

Solving for boundary modes. First we need to correct the RHS: $f'_b = f_b - A_{bi}A_{ii}^{-1} f_i$.



In [ ]:

    
gemv!('T', -1.0, M, fi, 1.0, fb);

Solve the boundary linear system:



In [ ]:

    
potrs!('L', Abb, fb);

Variable fb contains $x_b$. Now solve the equation for $x_i$:



In [ ]:

    
potrs!('L', Aii, fi)
gemv!('N', -1.0, M, fb, 1.0, fi);



In [ ]:

    
x = copy(fbi)

Typical solution:



In [ ]:

    
x0 = A\f



In [ ]:

    
x0 - x

Non symmetric matrices

In this case,

$$ A_{bi}^T \ne A_{ib} $$

Generating a nonrandom matrix B



In [ ]:

    
B = random_matrix(N, false)



In [ ]:

    
Bbb = B[1:nb, 1:nb]
Bbi = B[1:nb, (nb+1):N]
Bib = B[(nb+1):N, 1:nb]
Bii = B[(nb+1):N, (nb+1):N];



In [ ]:

    
using Base.LinAlg.LAPACK.getrf!

Bii, ipiv, info = getrf!(Bii);

We need to calculate the matrix $A_{bi}A_{ii}^{-1}AA_{ib}$. In this case the matrix is not symmetric and therefore we need to calculate the matrix

$$ M = \left(A_{bi}A_{ii}^{-1}\right)^T = \left(A_{ii}^{-1}\right)^T A_{bi}^T $$



In [ ]:

    
using Base.LinAlg.LAPACK.getrs!

M2 = Bbi'

getrs!('T', Bii, ipiv, M2)

When calculating the RHS, we need to compute $M_b = A_{ii}^{-1}A_{ib} \ne M^T$. No problem just pre-calculate it and store for further use.



In [ ]:

    
M2b = copy(Bib)
getrs!('N', Bii, ipiv, M2b);



In [ ]:

    
gemm!('T', 'N', -1.0, M2, Bib, 1.0, Bbb);



In [ ]:



In [ ]:

    
Bbb, ipiv2, info = getrf!(Bbb)

Solving for a RHS



In [ ]:

    
g = float([1:N;])
gbi = copy(g)



In [ ]:

    
gb = view(gbi, 1:nb)
gi = view(gbi, (nb+1):N)



In [ ]:

    
gemv!('T', -1.0, M2, gi, 1.0, gb);



In [ ]:

    
getrs!('N', Bbb, ipiv2, gb);
gb



In [ ]:

    
getrs!('N', Bii, ipiv, gi)
gemv!('N', -1.0, M2b, gb, 1.0, gi);



In [ ]:

    
y = copy(gbi)
y0 = B\g

hcat(y0, y, y-y0)



In [ ]:



In [ ]: