构建块三对角矩阵

Posted 2023-02-22

【中文标题】构建块三对角矩阵

【英文标题】：Build a block tri-diagonal matrix

【发布时间】：2016-08-31 17:52:29

【问题描述】：

我正在尝试在 Fortran 中。现在我有了这段代码，它只处理放置在A_matrix 的主对角线中的矩阵，i 中的每一步都有一个新矩阵。

do i = gs+1, total_mesh_points 

    start_line = (3*i)-2
    start_colu = (3*i)-2
    final_line = (3*i)
    final_colu = (3*i)

    do ii = 1, 3
    do jj = 1, 3
        A_matrix(start_line:final_line,start_colu:final_colu) = &
            impflux(ii,jj)
    end do
    end do

end do

这里我的A_matrix(i,j) 是一个大矩阵，它将在其主对角线上接收另一个三乘三矩阵 (impflux)。请注意，对于i 中的每一步，我都会有一个新的impflux 矩阵，该矩阵需要位于A_matrix 的主对角线中。

对于这个问题，我想不出更简单的解决方案。人们通常如何在 Fortran 中构建块对角矩阵？

【问题讨论】：

该代码没有多大意义。在以ii 和jj 为索引的循环内，它反复更新a_matrix 的相同元素，首先是impflux(1,1)，然后是impflux(1,2)，等。如果您向我们展示您正在尝试创建的内容的小示例，我可能会更好地理解。

感谢您的回复。我正在处理偏微分方程的隐式方法。在这些方法中，我们必须为网格中的每个点创建三个矩阵 (3x3)。代码的第一个循环通过网格中的所有点。其余代码将负责正确定位 A_matrix 内的矩阵 impflux（为网格中的每个点生成）。最终的结果将是 A_matrix 的主对角线内的一堆 impflux 矩阵。

理想的情况是，如果我将三个矩阵向量 a(i,j,nm)、b(i,j,nm) 和 c(i,j,nm) 传递给一个函数，然后它返回给我一个非常大的矩阵 (C)，这三个矩阵 (a,b,c) 形成了三个主对角线。

您没有尝试创建三对角矩阵。您正在尝试创建一个block tridiagonal matrix。它们根本不是一回事

【参考方案1】：

这是构建块三对角矩阵的一种方法。我不确定在一些知名库之外是否存在通常的方式。这是一个程序，我会留给你把它变成一个函数。

PROGRAM test

  USE iso_fortran_env

  IMPLICIT NONE

  INTEGER :: k    ! submatrix size
  INTEGER :: n    ! number of submatrices along main diagonal
  INTEGER :: ix   ! loop index

  ! the submatrices, a (lower diagonal) b (main diagonal) c (upper diagonal)
  REAL(real64), DIMENSION(:,:,:), ALLOCATABLE :: amx, bmx, cmx

  ! the block tridiagonal matrix
  REAL(real64), DIMENSION(:,:), ALLOCATABLE :: mat_a

  k = 3  ! set these values as you wish
  n = 4

  ALLOCATE(amx(n-1,k,k), bmx(n,k,k), cmx(n-1,k,k))
  ALLOCATE(mat_a(n*k,n*k))

  mat_a = 0.0

  ! populate these as you wish
  amx = 1.0
  bmx = 2.0
  cmx = 3.0

  ! first the lower diagonal
  DO ix = 1,k*(n-1),k
     mat_a(ix+k:ix+2*k-1,ix:ix+k-1) = amx(CEILING(REAL(ix)/REAL(k)),:,:)
  END DO

  ! now the main diagonal
  DO ix = 1,k*n,k
     mat_a(ix:ix+k-1,ix:ix+k-1) = bmx(CEILING(REAL(ix)/REAL(k)),:,:)
  END DO

  ! finally the upper diagonal
  DO ix = 1,k*(n-1),k
     mat_a(ix:ix+k-1,ix+k:ix+2*k-1) = cmx(CEILING(REAL(ix)/REAL(k)),:,:)
  END DO

END PROGRAM test

请注意，这里根本没有错误检查，我只做了一些测试。

一个明显的替代方法是只循环mat_a 的行一次，在同一迭代中插入amx、bmx、cmx，但这需要对第一次和最后一次迭代进行特殊处理，并且可能看起来要复杂得多。至于性能，如果对您很重要，请运行一些测试。

还要注意，这会产生一个密集矩阵。如果您的矩阵变得非常大，那么仅存储对角线元素的方法可能会更有用。这将我们引向派生类型和对它们的操作，这完全是另一个问题。

【讨论】：

感谢 High Performance Mark 的回复。我会处理你的程序！