构建块三对角矩阵

Question

我正在尝试在 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 中构建块对角矩阵？

Answer 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，但这需要对第一次和最后一次迭代进行特殊处理，并且可能看起来要复杂得多。至于性能，如果对您很重要，请运行一些测试。

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