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
\begin{equation} 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\} \end{equation}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:
To solve the linear system given a RHS (f):
In [862]:
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
Out[862]:
In [863]:
N = 12
nb = 6
ni = N - nb
Out[863]:
In [864]:
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);
In [865]:
Abb = A[1:nb, 1:nb];
In [866]:
Abi = A[1:nb, (nb+1):N];
In [867]:
Aib = A[(nb+1):N, 1:nb];
In [868]:
Aii = A[(nb+1):N, (nb+1):N];
In [869]:
using ArrayViews
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 [870]:
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 [871]:
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 [872]:
gemm!('T', 'N', -1.0, M, Aib, 1.0, Abb);
Cholesky decomposition of $A'_{bb}$:
In [873]:
potrf!('L', Abb);
In [874]:
f = float([1:N]);
In [875]:
fbi = copy(f);
In [876]:
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 [877]:
gemv!('T', -1.0, M, fi, 1.0, fb);
Solve the boundary linear system:
In [878]:
potrs!('L', Abb, fb);
Variable fb contains $x_b$. Now solve the equation for $x_i$:
In [879]:
potrs!('L', Aii, fi)
gemv!('N', -1.0, M, fb, 1.0, fi);
In [880]:
x = copy(fbi)
Out[880]:
Typical solution:
In [881]:
x0 = A\f
Out[881]:
In [882]:
x0 - x
Out[882]:
In [883]:
B = random_matrix(N, false)
Out[883]:
In [884]:
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 [886]:
using Base.LinAlg.LAPACK.getrf!
Bii, ipiv, info = getrf!(Bii);
Out[886]:
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 [887]:
using Base.LinAlg.LAPACK.getrs!
M2 = Bbi'
getrs!('T', Bii, ipiv, M2)
Out[887]:
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 [888]:
M2b = copy(Bib)
getrs!('N', Bii, ipiv, M2b);
In [889]:
gemm!('T', 'N', -1.0, M2, Bib, 1.0, Bbb);
In [890]:
M2' * Bib - MM
Bbb - Bbb2
Out[890]:
In [891]:
Bbb, ipiv2, info = getrf!(Bbb)
Out[891]:
In [892]:
g = float([1:N;])
gbi = copy(g)
Out[892]:
In [893]:
gb = view(gbi, 1:nb)
gi = view(gbi, (nb+1):N)
Out[893]:
In [894]:
gemv!('T', -1.0, M2, gi, 1.0, gb);
In [895]:
getrs!('N', Bbb, ipiv2, gb);
gb
Out[895]:
In [896]:
getrs!('N', Bii, ipiv, gi)
gemv!('N', -1.0, M2b, gb, 1.0, gi);
In [897]:
y = copy(gbi)
y0 = B\g
hcat(y0, y, y-y0)
Out[897]:
In [ ]: