CODEWITHSUNDEEP
matrixfortran
LM_O
LM_O

Reputation: 131

Build a block tri-diagonal matrix

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 ?

Upvotes: 2

Views: 1865

Answers (1)

High Performance Mark
High Performance Mark

Reputation: 78354

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.

Upvotes: 4

Related Questions