我正在尝试构建并求解一个如下所示的矩阵系统:
\begin{equation}
\begin{bmatrix}
B & D & 0 & \cdots & & 0 \\
A & B & D & \ddots & & \vdots \\
0 & A & B & D & \ddots & \\
\vdots & \ddots & & \ddots & \\
\vdots & \cdots & \cdots & A & B & D \\
0 & \cdots & \cdots & 0 & A & B
\end{bmatrix}
\
\begin{bmatrix}
\Delta c
\end{bmatrix}
\
=
\
\begin{bmatrix}
G
\end{bmatrix}
\end{equation}
该系统来自用于求解方程组 (PDE) 的一维有限体积方法。
大矩阵 ( BBM
) 由数组构造而成A, B, D
。以下代码可以正常运行,但我有兴趣使该代码的计算效率更高。
我很好奇堆栈溢出社区是否有任何建议。特别是有没有一种方法可以消除for
函数内的循环?
using BlockBandedMatrices
function calculate_BBM(; A, B, D)
# initialize the Block Banded Matrix (BBM)
int_array = map(Int, N*ones(NJ));
BBM = BlockBandedMatrix(Zeros{Float64}(N*NJ, N*NJ), int_array, int_array, (1,1))
# Precompute indices to avoid calculating them in each iteration
ind_begin = [(1 + N * (j - 1)) for j in 1:NJ]
ind_end = [N * j for j in 1:NJ]
for j in 2:NJ - 1
# Use precomputed indices
ib = ind_begin[j]
ie = ind_end[j]
BBM[ib:ie, (ib - N):(ie - N)] .= A[:, :, j - 1]
BBM[ib:ie, ib:ie] .= B[:, :, j]
BBM[ib:ie, (ib + N):(ie + N)] .= D[:, :, j]
end
j = 1
ib = ind_begin[j]
ie = ind_end[j]
BBM[ib:ie, ib:ie] .= B[:, :, j]
BBM[ib:ie, (ib + N):(ie + N)] .= D[:, :, j]
j = NJ
ib = ind_begin[j]
ie = ind_end[j]
BBM[ib:ie, (ib - N):(ie - N)] .= A[:, :, j - 1]
BBM[ib:ie, ib:ie] .= B[:, :, j]
return BBM
end
N = 3
NJ = 102
A = randn(N, N, NJ-1)
B = randn(N, N, NJ)
D = randn(N, N, NJ-1)
# A = -ones(N, N, NJ-1)
# B = ones(N, N, NJ) * 2
# D = ones(N, N, NJ-1)
BBM = calculate_BBM(A=A, B=B, D=D);