ブロック三角対角行列の作成

Question

Fortranでブロック三重対角行列を作成しようとしています。今度はA_matrixのメイン対角線に配置された行列だけを扱うこのコードを用意しました。これはすべてのステップの1つの新しい行列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)は、その主対角線に別の3×3行列（impflux）を受け取ります大行列です。 iの各ステップでは、新しいimpflux行列をA_matrixのメイン対角線に配置する必要があることに注意してください。

私はこの問題のより簡単な解決策は考えられません。人々が通常Fortranで対角行列をブロックする方法

Answer 1

ブロック三重対角行列を作成する方法の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 

ここでエラーをチェックすることは一切ありません。私はいくつかのテストを行いました。

1つの明白な選択肢は、同じ繰り返しでamx、bmx、cmxを挿入し、一度だけmat_aの行をループになりますが、これは、最初と最後の反復のための特別な処理を必要とし、おそらく多くのより複雑になります。パフォーマンスに関しては、いくつかのテストを実行することが重要な場合。

これは密度の高いマトリックスを生成することにも注意してください。行列が非常に大きくなると、対角要素だけを格納するアプローチがより役に立ちます。それは、派生した型と操作に向けて私たちを連れて行きます。それはまったく別の質問です。