今週のプログラムその3

今週は二次元配列を使って行列の計算をする。いきなりプログラムが長くなったと思うかもしれないが、短い部品の組み合わせなので、心配しなくてもオッケー

ポイント

  1. 行列は二次元の配列で表現する。二次元配列の書きかたはCとFORTRANで若干違い、2行2列の実数配列matrixの宣言は
    
     float matrix[2][2];  (Cの場合)
          real matrix(2,2)  (FORTRANの場合)
    
    とする。
  2. 行列の次元などのように、扱うデータによって大きさが変わり、ただし計算の間は定数とみなせるような値はFORTRANならパラメタ文で定義すると便利である。

    parameter(nn=3)
    という文がプログラムの初めに書いてあれば、プログラム中のnnはすべてコンパイラによって3という値と置き換えられ、実行時には定数として扱われる。

    同様の効果をCで得るには#define文を使う。

    #define NN 3
    という文がmainの前にあれば、プログラム中のNNという文字はすべて3と置き換えられる。 ただし、こちらは単に文字を置き換えるだけであって、"3"を数として扱っている わけではない。エディタの置換機能を使って文字を置き換えるのと同じこと をやっているだけなので、FORTRANのパラメタ文の働きとは本来違うものである。

  3. program 8で使う入力データは、ふたつの行列の値を並べたもの
    3*3行列の入力データの例(数値には特に意味はありません)。 上側三行分が行列a、下側三行分が行列bに代入される値
    0.1  2.0  1.34
    2.2  0    5.5
    0.01 3.14 2.42
    0.2  2.2  4.34
    0.01 0.02 0.05
    1    3.14 2
    
なお、プログラム自体は

まずはC版

program 7


#include<stdio.h>
/* 行列とベクトルの掛け算 */
main()
{
  float a[2][2],x[2],y[2];       実数配列を宣言
  int i,j;

/*  行列に値をいれる(値自身は意味ないです) */
  a[0][0]=0.5;
  a[0][1]=1.0;
  a[1][0]=0.8;
  a[1][1]=2.0;

/*  ベクトルに値をいれる */
  x[0]=1.0;
  x[1]=2.0;

/*  ベクトルyを初期化(値を0に) */
  for( i=0; i<=1; ++i){          配列のすべての要素を0にする
     y[i]=0;
  }    

/*  a*xによりyを求める */
  for( i=0; i<1; ++i){            二次元配列なので要素を指定するために
                                    i,jのふたつの変数を使う。
                                    ループも二重になる。
    for( j=0; j<=1; ++j){
       y[i] += a[i][j]*x[j];      行列とベクトルの積の計算はいつもこれ。
                                  a[i][j]とx[j]の積をy[i]に足している。
    }
  }

/* 出力  */
  printf("matrix a:\n");
  for( i=0; i<=1; ++i){
    for( j=0; j<=1; ++j){       二次元配列を出力したいので二重ループ
      printf("  %f",a[i][j]);
    }
    printf("\n");                 一行出力したら改行
  }
  printf("vector x:\n");
  for( i=0; i<=1; ++i){         こちらは一次元配列なので一重ループ
    printf("  %f\n",x[i]);
  }
  printf("vector y=a*x:\n");
  for( i=0; i<=1; ++i){ 
    printf("  %f\n",y[i]);
  }
}

program 8


#include<stdio.h>
#define NSIZ 3             プログラム中のNSIZという文字をすべて3に置き換える
/* defineを使った任意サイズの行列同士の掛け算 */
main()
{
  float a[NSIZ][NSIZ],b[NSIZ][NSIZ],c[NSIZ][NSIZ];    配列の宣言
  int i,j,k;

/*  キーボード(またはファイル)から行列に値を読み込む */
  for( i=0; i<=NSIZ-1; ++i){
    for( j=0; j<=NSIZ-1; ++j){        二次元配列なので二重ループ
      scanf("%f",&a[i][j]);
    }
  }
  for( i=0; i<=NSIZ-1; ++i){
    for( j=0; j<=NSIZ-1; ++j){       これも二重ループ
      scanf("%f",&b[i][j]);
   }
  }

/*  行列cを初期化(値を0に) */
  for( i=0; i<=NSIZ-1; ++i){
    for( j=0; j<=NSIZ-1; ++j){      行列の要素すべてに0を代入する
      c[i][j]=0;
    }
  }

/* a*bによりcを求める */
  for( i=0; i<=NSIZ-1; ++i){
    for( j=0; j<=NSIZ-1; ++j){
      for( k=0; k<=NSIZ-1; ++k){
        c[i][j] += a[i][k]*b[k][j];    行列の積の計算はいつもこれ。
                                    a[i][k]とb[k][j]の積をc[i][j]に足している。
                                    三重ループになっていることに注意
      }
    }
  }

/* 出力 */
  printf("matrix a:\n");
  for( i=0; i<=NSIZ-1; ++i){
    for( j=0; j<=NSIZ-1; ++j){      二次元配列なので二重ループ
      printf("  %f",a[i][j]);
    }
    printf("\n");
  }
  printf("matrix b:\n");
  for( i=0; i<=NSIZ-1; ++i){
    for( j=0; j<=NSIZ-1; ++j){
      printf("  %f",b[i][j]);
    }
    printf("\n");
  }
  printf("matrix c=a*b:\n");
  for( i=0; i<=NSIZ-1; ++i){
    for( j=0; j<=NSIZ-1; ++j){
      printf("  %f",c[i][j]);
    }
    printf("\n");
  }
}

FORTRAN版

program 7


C 行列とベクトルの掛け算
      program main
      real a(2,2),x(2),y(2)             実数配列を宣言

C  行列に値をいれる(値自身は意味ないです)
      a(1,1)=0.5
      a(1,2)=1.0
      a(2,1)=0.8
      a(2,2)=2.0

C  ベクトルに値をいれる
      x(1)=1.0
      x(2)=2.0

C  ベクトルyを初期化(値を0に)
      do i=1,2
        y(i)=0           配列の全要素に0を代入する
      enddo

C a*xによりyを求める
      do i=1,2
       do j=1,2                  二次元配列なので要素を指定するために
                                 i,jのふたつの変数を使う。ループも二重になる。
         y(i)=y(i)+a(i,j)*x(j)        行列とベクトルの積の計算はいつもこれ。
                                 a(i,j)とx(j)の積をy(i)に足している。
       enddo
      enddo        

C 出力 
      write(*,*)"matrix a:"
      do i=1,2
        write(*,*)(a(i,j),j=1,2)    この書きかたで一行にj=1と2の値を並べてくれる。
                    二次元配列を出力したいので二重ループになるのだが、
                   二重ループの内側ループが(a(i,j),j=1,2)ですんで
                  しまっている。これはwrite文などいくつかの場合にだけ使える
                   方法。
      enddo
      write(*,*)"vector x:"
      write(*,*)(x(i),i=1,2)    こちらは一次元配列なので一重ループ。
                               write文中であるため、このように簡単に書ける。
      write(*,*)"vector y=a*x:"
      write(*,*)(y(i),i=1,2)
      end

program 8


C パラメタ文を使った任意サイズの行列同士の掛け算
      program main
      parameter(nn=3)                 プログラム中のnnという変数は
                                       すべて3という定数に置き換える
      real a(nn,nn),b(nn,nn),c(nn,nn)

C  ファイルから行列に値を読み込む
      do i=1,nn
        read(*,*)(a(i,j),j=1,nn)      二次元配列なので二重ループ。
                               write文と同様、データを読むほうのread文でも、
                               この書き方でj=1からnnまでを読める
      enddo
      do i=1,nn
        read(*,*)(b(i,j),j=1,nn)    これも二重ループ
      enddo

C  行列cを初期化(値を0に)
      do i=1,nn                     行列の要素すべてに0を代入する二重ループ
        do j=1,nn
          c(i,j)=0
        enddo
      enddo

C a*bによりcを求める
      do i=1,nn
       do j=1,nn
         do k=1,nn
           c(i,j)=c(i,j)+a(i,k)*b(k,j)   行列の積の計算はいつもこれ。
                                 a(i,j)とb(k,j)の積をc(i,j)に足している。
                                   三重ループになっていることに注意

         enddo
       enddo
      enddo        

C 出力 
      write(*,*)"matrix a:"
      do i=1,nn
        write(*,*)(a(i,j),j=1,nn)    二次元配列なので二重ループ。
                                     内側ループは例によって簡単な書き方
      enddo
      write(*,*)"matrix b:"
      do i=1,nn
        write(*,*)(b(i,j),j=1,nn)
      enddo
      write(*,*)"matrix c=a*b:"
      do i=1,nn
        write(*,*)(c(i,j),j=1,nn)
      enddo
      end