【C＆Fortran演習で学ぶ数値計算(2)】行列の作成

本日も数値計算の勉強をしていきます。
読んでいる参考書はこちら。

C & Fortran 演習で学ぶ数値計算

連立方程式の解法

連立方程式の解法にはいくつかの手法があります。
大きく分けて2種類になります。

https://speakerdeck.com/kamakiri1225/shu-zhi-ji-suan

・直接法：厳密解を求める
・反復法：収束させて近似解を求める

まずは直接法での解法から見ていくことになりますが、直接法のプログラム作成は次回に回します。
直接には「ガウス・ザイデル法」「ガウス消去法」「LU分解」の解説が書かれており、プログラムコードはC言語とFortranで書かれています。

行列の作成

連立1次方程式を解くにあたって、10×10の正方行列と1次元ベクトルを作成する必要があります。

コードはC言語とFortranで書かれているため参考書を見ながら書き写しました。

C言語

#include <stdio.h>

#define N 10

int main(){
    double A[N][N], b[N];
    int i, j, ii;
    
  /* 行列作成 */
  for (j=0; j<N; j++) {
    ii = 0;
    for (i=j; i<N; i++) {
      A[j][i] = (N-j) - ii;
      A[i][j] = A[j][i];
      ii++;
    }
  }
    
    printf("Matrix A:\n");
    for (j = 0; j < N; j++){
        for ( i = 0; i < N; i++){
            printf("%5.2f",A[j][i]);
        }
    printf("\n");
    }

    printf("\n");
    
    // ベクトル作成
    for ( i = 0; i < N; i++){
        for ( j = 0; j < N; j++){
            b[i] += A[i][j];
        }
    }

    printf("Vector b:\n");
    for ( i = 0; i < N; i++){
        printf(" %5.2f",b[i]);
    }
    printf("\n");
    return 0;
}

コンパイルして実行すると以下のようになります。

Fortran

program main
    implicit none

    integer, parameter:: N = 10
    real(kind=8), dimension(N,N)::A
    real(kind=8), dimension(N)::b
    integer:: i, j , ii

    !行列の作成
    do j=1, N
        ii = 0
        do i=j,N
            A(j,i) = (N-(j-1))-ii
            A(i,j) = A(j,i)
            ii = ii +1
        end do
    end do

    !行列の確認
    write(*,*)"Matrix:A"
    do j=1, N
        do i=1,N
            write(*, fmt='(1x, f5.2)', advance='no') A(j,i)
        end do
        write(*,*)""
    end do   

    !行列の作成
    do i=1, N
        b(i)=0.d0
        do j=1,N
            b(i) = b(i) + A(i,j)
        end do
    end do

    !ベクトルの確認
    write(*,*)"Vector:b"
    do i=1,N
        write(*, fmt='(1x, f5.2)', advance='no') b(i)
    end do
    write(*,*)""
end program main

コンパイル後に実行すると以下のようになります。

C言語とFortranの違い

学生時代にはじめてFortranを触って以来ずっとFortranに慣れてきたので、そういえば配列要素を作成したときにFortranのindexは1から開始するということを忘れていました。
また、C言語とFortranでメモリの並び順も違うんですよね。

C言語の場合、2次元配列array[4][4]を作成するとメモリが連続するのは、array[0][0],array[0][1],array[0][2],array[0][3],array[1][0]・・・・
だから連続するメモリアクセスを行うにあたって、C言語は以下のようにfor文を回せばよいのですが・・・・

    for (i = 0; i < N; i++){
        for ( j = 0; j < N; j++){
            printf("%5.2f",A[i][j]);
        }
    }

Fortranの場合、2次元配列array(4,4)を作成するとメモリが連続するのは、array(1,1),array(2,1),array(3,1),array(4,1),array(1,2)・・・・
だから連続するメモリアクセスを行うにあたって、Fortranは以下のようにdo文を回さなくてはならない・・・・

    do j=1, N
        do i=1,N
            write(*, fmt='(1x, f5.2)', advance='no') A(i,j)
        end do
    end do   

不連続なメモリアクセスは実行速度を落とすのでFortranを使う場合は注意が必要ですね。

本日は以上です(^^♪

Fortranのおすすめ参考書

Fortranの文法の基礎を学びたい方にはこちらの参考書がおすすめです。
後半は簡単な偏微分方程式を差分法で解くためのプログラムが乗っているため、基礎からちょっとした応用まで学ぶことができます。
もちろん高度な数値計算（例えば、構造解析や流体解析）をしたい場合は他の参考書を読んで勉強する必要があります。

数値計算のためのFortran90/95プログラミング入門(第2版)

Fortranのコード集として「Fortranハンドブック」がお勧めです。
様々な差分法や数値積分方法などが乗っているのでとても勉強になります。

Fortran ハンドブック

C言語のおすすめ参考書

実はC言語は最近勉強し始めたばかりであまり詳しくなく、個人的にわかりやすいと感じた参考書をひとつ紹介します。

スッキリわかるC言語入門 第2版 (スッキリわかる入門シリーズ)

Twitter➡@t_kun_kamakiri
Instagram➡kamakiri1225
youtube➡https://www.youtube.com/channel/UCbG6_Q9ZRqqVT6YZOpcjDlQ
ブログ➡宇宙に入ったカマキリ（物理ブログ）
ココナラ➡物理の質問サポートサービス
コミュニティ➡製造業ブロガー