Build a block tri-diagonal matrix

Question

I am trying to build a block tridiagonal matrix in Fortran. Now I have this piece of code that would deal with just the matrices that are placed in the main diagonal of the A_matrix, one new matrix for every step in 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

Here my A_matrix(i,j) is a big matrix that will receive another three by three matrix (impflux) in its main diagonal. Note that for each step in i I will have a new impflux matrix that needs to be positioned in the main diagonal of the A_matrix.

I can't think in a more simple solution for this problem. How people usually build block diagonal matrices in Fortran ?

That code doesn't make a lot of sense. Inside the loops indexed over ii and jj it repeatedly updates the same elements of a_matrix, first with impflux(1,1), then with impflux(1,2), etc. I'd probably understand better if you showed us a small example of what you are trying to create. — High Performance Mark
Thank you for the reply. I'm dealing with implicit methods for partial differential equations. In those methods we have to create three matrix (3x3) for each point in the mesh. The first loop of the code pass through all points in the mesh. The rest of the code would take care of the correct positioning of the matrices impflux (generated for each point in the mesh) inside the A_matrix. The final result would be a bunch of impflux matrices inside the main diagonal of A_matrix. — LM_O
The ideal world would be if I passed three vectors of matrices a(i,j,nm), b(i,j,nm) and c(i,j,nm) to a function and it returned to me a very big matrix (C) with these three matrix (a,b,c) forming the three main diagonals. — LM_O
You are not trying to create a tridiagonal matrix. You are trying to create a block tridiagonal matrix. They are not at all the same thing — talonmies

High Performance Mark High Performance Mark · Accepted Answer · 2016-09-02T20:19:32

Here's one way to build a block tridiagonal matrix. I'm not sure that there is, outside some well-known libraries, a usual way. This is a program, I'll leave it up to you to turn it into a function.

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

Be warned, there's no error checking here at all and I've only made a few tests.

One obvious alternative would be to loop over the rows of mat_a only once, inserting amx, bmx, cmx at the same iteration, but this would require special handling for the first and last iterations and probably look a lot more complicated. As for performance, if it matters to you run some tests.

Note also that this produces a dense matrix. If your matrix gets very large then an approach which stores only the diagonal elements might be more useful. That takes us towards derived types and operations on them, and that's a whole other question.

Build a block tri-diagonal matrix

1 Answers