Week 5 (2021/11/3)

今日やること

• ポインタのポインタ
• ポインタの配列
• 二次元配列とポインタ
• コマンドライン引数

ポインタのポインタ

今日はポインタの応用編をやっていきます。頑張っていきましょう。まずはポインタのポインタです。

素朴な例

「int型の変数を指すポインタ」を「intへのポインタ」などと呼びました。では、 「『intへのポインタ』を指すポインタ」は考えられるのでしょうか？考えられます。それがポインタのポインタです。 次の素朴な例を考えてみましょう。

int a = 10;
printf("   a: %d\n", a);     // 10

int *p = &a;
printf("   p: %p\n",  p);    // 0xff5c   (以下、アドレスは全て人工的な例です)
printf("  *p: %d\n", *p);    // 10 

int **pp = &p;
printf("  pp: %p\n",   pp);  // 0x51a5
printf(" *pp: %p\n",  *pp);  // 0xff5c
printf("**pp: %d\n", **pp);  // 10

このように、int **ppとアスタリスクを２つ重ねることで、ポインタへのポインタを表現することができます。 あるポインタint *pに対し、ppにはポインタpのアドレス&pを代入します。上の表記の繰り返しですが、a, p, ppの関係は次のようになります。

• a と *p と **pp は同じ値：10. 型は int
• p と *pp は同じ値：0xff5c. 型はint *
• ppの値：0x51a5. 型はint **

関数の引数としてポインタそのものを受け取る場合

さて、素朴なポインタのポインタを使う例を見てみましょう。ここではまず、「main側のポインタを、関数から編集する」という例を見てみましょう。 まず、配列中から一番大きい値を返すfind_largest関数を 考えます。ここでは、ポインタを用いて最大値を返しています。 ここで変数が実際にどのようにメモリ中に表現されるかを見てみましょう。 以下では、通常の変数は黒、配列は赤、ポインタは青、と色分けしています Reese, p99

ここではいくつかのポイントがあります

• mainが扱う領域と、find_largestが扱う領域は違います。そのため、main側の変数をfind_largestから変更するには、ポインタを経由する必要があります。
• int a[]とint *maxは、find_largest側ではどちらもポインタになっています。これは先週習った、「関数引数におけるa[]は*aに等しい」ということですね。
• 関数の呼び出し時に引数の変数３つがコピーされます。このとき、配列全てをコピーしたりするわけではないので注意してください。
• 関数引数（a, n, max）は、あくまで、関数側のメモリ領域中でのローカル変数です

上図はmaxを経由して値を更新するときの動作です。 maxはmain側のansを指しているので、*maxという表記でansの中身を書き換えることができます。

上のfind_largestは、main側のint型変数ansに、arr[3]中の最大値を記述する関数でした。 次にfind_largest2を考えましょう。これは、 main側のint型ポインタp_ansが、arr[3]中の最大要素を指すようにする関数です。ここでは、ポインタのポインタを使います。 下図を見てみましょう。ここではポインタのポインタを緑色で表現しました。

ここでは次のポイントがあります

• 基本的な構造は先ほどのfind_largestと変わりません。
• find_largestでは関数側のint *maxがmain側のint ansを編集しました。つまり、intを編集するためにポインタを使いました。一方、find_largest2では、関数側のint **p_maxがmain側のint *p_ansを編集します。つまり、intポインタを編集するので、intポインタのポインタを使いました。

実際にp_ansを編集するときの挙動が上の図になります。p_ansには、arr要素中で最大のものアドレス（&arr[1]）を代入しています。

このようにして、「main側のポインタ」を編集する際に、関数側では引数として「ポインタのポインタ」を使う場合があります。 ポインタのポインタを扱う例は他にもあり、最も使う例として「ポインタの配列」があります。それを次章で見ていきましょう。

やってみよう（20分）

上記を写経してみましょう

クイズ

• int *max を引数にとる

  void find_largest(int a[], int n, int *max) { ... }
  

  に対し、最大値を直接、値としてreturnする

  int find_largest_(int a[], int n) { ... }
  

  を書いてみましょう。つまり、次のような使い方になります。

  int ans = find_largest(arr, 3);
  

  これは簡単ですね。
• int **p_max を引数にとる

  void find_largest2(int a[], int n, int **p_max) { ... }
  

  に対し、最大要素へのポインタ（アドレス）を直接値としてreturnする

  int *find_largest2_(int a[], int n) { ... }
  

  を書いてみましょう。

答え

• 例えば以下。

  int find_largest_(int a[], int n) {
      int max = a[0];
      for (int i = 0; i < n; ++i) {
          if (max < a[i]) {
              max = a[i];
          }
      }
      return max;
  }
  int main() {
      int arr[] = {3, 6, 2};
      int ans1, ans2;
      find_largest(arr, 3, &ans1);
      ans2 = find_largest_(arr, 3);
      printf("%d %d\n", ans1, ans2);  // 6 6 
  }
  

  ここでポイントは、ある型Tの変数xを変更したいとき、次の二通りがあるということです。
  • void f(T *x) { ... } ポインタとして変更
  • T f() { ... } 直接return
• 例えば以下。

  int *find_largest2_(int a[], int n) {
      int *p_max = &a[0];  // 最初の要素へのポインタ
      for (int i = 0; i < n; ++i) {
          if (*p_max < a[i]) {
              p_max = &a[i];  // i番目の要素へのポインタ
          }
      }
      return p_max;
  }
  int main() {
      int arr[] = {3, 6, 2};
      int *p_ans1, *p_ans2;
      find_largest2(arr, 3, &p_ans1);
      p_ans2 = find_largest2_(arr, 3);
      printf("%d %d\n", *p_ans1, *p_ans2);  // 6 6 
  }
  

  ここでポイントは、ある型Tへのポインタ変数pを変更したいとき、次の二通りがあるということです。
  • void f(T **p) { ... } ポインタとして変更。つまりポインタのポインタを用いる
  • T *f() { ... } 直接return

このように思考を整理しておくと、ソフトウェア２でmallocや線形リストといった、ポインタを駆使する概念を習うときに役にたちます。

ポインタの配列

さて、それではポインタの配列を考えてみましょう。

ポインタ配列の基本

char s1[] = "toyama";
char s2[] = "ishikawa";
char s3[] = "fukui";

char *hokuriku1[3];
hokuriku1[0] = s1;  // &s1[0]といっしょ。s1の先頭アドレスを代入
hokuriku1[1] = s2;
hokuriku1[2] = s3;

ここで、hokuriku1は配列です。そして、その要素は「char型へのポインタ」です。 そして、hokuriku1の各要素には文字列の先頭アドレスを代入にしています。これを図示すると次のようになります。

このように、文字列s1, s2, s3がどこかで確保されます。そして、その先頭を、hokuriku1の各要素が指します。 hokuriku1[0]という表記がs1と同じであることを考えると、hokuriku1[0][0]は't'であることがわかります。 要素のアクセスは次のように書けます。

for (int i = 0; i < 3; ++i) {
    hokuriku1[i][0] = 'A';
    printf("%s\n", hokuriku1[i]);
}
// 結果：
// Aoyama
// Ashikawa
// Aukui

さて、このような「文字列を複数確保 + それらの先頭を指すポインタの配列」は、配列の初期化記法を用いて次のようにも書けます。

char *hokuriku1[] = {"toyama", "ishikawa", "fukui"};

注意として、この記法で変数を宣言した場合は、要素の変更は出来ません。

コラム

なぜ初期化記法を使って作った文字列は編集できないのでしょうか？これはweek4で習った文字列の初期化 がそのまま適用されているからです。char s1[] = "toyama";のように明示的に配列を作った場合はそれらは自由に編集できます。 一方、初期化記法を用いた場合は、「自分が編集できる配列」が確保されているのではなく、プログラムのどこかに確保された「編集不可能なリテラル」を指しているだけなのです。 以下を実行してエラーになることを確認しましょう。

char *hokuriku1[] = {"toyama", "ishikawa", "fukui"};
hokuriku1[0][0] = 'a';  // Segmentation Fault

再度確認しておくと、例えば次のようにすると、最初の"zaku"と最後の"zaku"は、同じリテラルなので、同じアドレスになったりします。

char *ms[] = {"zaku", "gufu", "zaku"};
printf("%p %p %p\n", ms[0], ms[1], ms[2]);
// 0x55a31d9b5004 0x55a31d9b5009 0x55a31d9b5004

上記の文字列の場合は書き込みはできないものの、読み込みはできます。 一方で、文字列でない場合は、そもそも配列を確保できないので注意してください。

int *a[] = {{0, 1}, {2, 3, 4}, {5}};  // この書き方はダメ
printf("%d", a[0][0]);  // Segmentation Fault

なぜなら、{0, 1}のような書き方は、文字列リテラルのようにどこかに領域を確保したりしないからです。この書き方は配列を初期化するときに使えるというだけです。

次のように、明示的に要素の配列を確保すれば、大丈夫です。

int arr1[] = {0, 1};
int arr2[] = {2, 3, 4};
int arr3[] = {5};
int *b[] = {arr1, arr2, arr3};  // これは、通常の「配列の初期化」してるだけです。その要素として、各arrの先頭アドレスが入ります。
printf("%d\n", b[0][0]);   // 0

さて、ポインタ配列の利点は何でしょうか。プログラミングを行う上で、もし全ての要素を配列や多重配列に保持できるのであれば話は楽なのですが、 ときにはそのようなことが難しい場面があります。その際に、配列の入り口だけを用意しておき、必要に応じて配列が指す先を指定する、 ということをしたい場面があります。そのときに、ポインタ配列は最適です。 その例としてスワップがあげられます。ポインタ配列の要素であるポインタは値としてアドレスを持つだけなので、要素を簡単にスワップすることができます。 これは、後述するように、二次元配列では簡単にはできません。

char *tmp = hokuriku[1];
hokuriku1[1] = hokuriku1[0];  // hokuriku[0]: "toyama" になった
hokuriku1[0] = tmp;           // hokuriku[1]: "ishiwaka" になった

これを図示すると次のようになります。ここではわかりやすくするため文字列を横に並べて表示しています。

ここでは、右側の文字列たちには変更が加わっていないことに注意してください。 あくまでポインタ配列の中の値（アドレス）だけが更新されています。ポインタ配列は、特に、ソフト２でmallocを習う際に必須になります。

ポインタ配列を引数にもつ関数

さて、上記のhokuriku1のようなポインタ配列を引数にもつ関数はどうすれば作れるでしょうか？ ここでは、通常の配列 int a[]を引数に持つ関数はvoid f(int a[]) {...}とそのまま書けたことを思い出しましょう。 ポインタ配列はあくまで普通の配列です。その要素がポインタというだけです。なので、通常の配列と同様に、そのまま書けます。

void print_each(char *prefectures[], int n) {
    for (int i = 0; i < n; ++i) {
        printf("%s-ken\n", prefectures[i]);
    }
}

この例では、prefecturesの中身に"-ken"を追加して表示しています。 ここでは、要素がint型の配列を引数にとるときはint a[]と書いたのと同様に、要素がchar *型である配列を引数にとっているのでchar *prefectures[]となっているというだけで、話はシンプルです。 さて、今日の最初に再確認したように、関数の引数に配列が与えられたときは、実はそれはポインタであるという点を思い出しましょう。 すなわち、関数引数中でT a[]の形が現れたときは、それは常にT *aでした。これを適用すると、上記の関数は実は次と同じです。

void print_each(char **prefectures, int n) { ... }

ここではポインタのポインタが出てきました。

ここで知っておくべきことは、ポインタのポインタが出現したとき、それは様々な使い方があるということです。

• ポインタのポインタはポインタ単体を指し、そのポインタは変数単体を指す： void find_largest2(int a[], int n, int **p_max)
• ポインタのポインタはポインタ配列を指し、その要素のポインタは変数の配列（文字列）を指す： void print_each(char **prefectures, int n)

{単体, 配列}×{単体, 配列}になっているので、残りの２パターンを考えることもできます。

• ポインタのポインタはポインタ単体を指し、そのポインタは変数の配列を指す
• ポインタのポインタはポインタ配列を指し、その要素のポインタは変数単体を指す

現実では、find_largest2のように素朴なポインタの例か、あるいはprefecturesのようにポインタ配列が文字列を指す例が多いと思いますが、 複雑なものが来ても驚かないようにしてください。

やってみよう（30分）

上記を写経してみましょう

クイズ

• int *v[]を、「int型配列の先頭アドレスを集めた配列」とします。このvをポインタで指すとすると、どういう型になりますか？また、要素にアクセスするとき、どのようになりますか？
• ポインタ配列char *prefectures[]、整数i、整数jが与えられたとき、prefectures[i]とprefectures[j]をswapする関数を書いてください。
• 上で勉強した、次のパターンの「ポインタのポインタ」の使用例を書いてください。型はcharでお願いします。
  • ポインタのポインタはポインタ単体を指し、そのポインタは変数の配列（すなわち文字列）を指す
  • ポインタのポインタはポインタ配列を指し、その要素のポインタは変数単体（すなわちchar一個）を指す

答え

• 例えば次のようになります。

  int arr1[] = {1, 2, 3};
  int arr2[] = {4, 5};
  int *v[] = {arr1, arr2};  // ポインタ配列 
  int **p;
  p = v;  // ポインタのポインタで指し示す
  printf("v[0][1]: %d\n", p[0][1]);  // v[0][1]: 2
  printf("v[1][0]: %d\n", p[1][0]);  // v[1][0]: 4
  

  これは、型Tに対して、T v[]を考えると、そのポインタはT *pでしたね。このTにint *を代入したと思うとわかります。
• 例えば以下。

  void swap(char *prefectures[], int i, int j) {
      char *tmp = prefectures[i];
      prefectures[i] = prefectures[j];
      prefectures[j] = tmp;
  }
  int main() {
      char *hokuriku1[] = {"toyama", "ishikawa", "fukui"};
      print_each(hokuriku1, 3);  // toyama-ken, ishikawa-ken, fukui-ken
      swap(hokuriku1, 1, 2);     // ishikawaとfukuiをswap
      print_each(hokuriku1, 3);  // toyama-ken, fukui-ken, ishikawa-ken
  }
  

• 例えば以下が「ポインタのポインタはポインタ単体を指し、そのポインタは変数の配列（文字列）を指す」例。
  

  char *animal = "dog";
  char **p;
  p = &animal;
  printf("%s\n", *p);  // dog
  

  例えば以下が「ポインタのポインタはポインタ配列を指し、その要素のポインタは変数単体（char一個）を指す」例です。ちなみにここではalphabetsの各要素は文字列ではないことに注意しましょう（最後にNULL文字がないので）

  char c1 = 'a', c2 = 'b', c3 = 'c';
  char *alphabets[] = {&c1, &c2, &c3};
  char **p;
  p = alphabets;
  for (int i = 0; i < 3; ++i) {
      printf("%c, ", *(p[i]));  // a, b, c
  }
  printf("\n");
  

二次元配列とポインタ

さて、次は二次元配列のことを思い出しましょう。二次元配列は、その使用目的はポインタ配列と似ていますが、中身は違うものです。

二次元配列の基本

以下のような初期化表記で二次元配列を確保できましたね。

char hokuriku2[][9] = {"toyama", "ishikawa", "fukui"};

これは、先ほどの「ポインタ配列」と同じようにアクセスできます。

for (int i = 0; i < 3; ++i) {
    printf("%s\n", hokuriku2[i]);
}
// toyama
// ishikawa
// fukui

二次元配列を可視化すると次のようになります。二次元配列の実際のメモリレイアウトについては既に勉強しましたね。

ここで、ポインタ配列hokuriku1と、二次元配列hokuriku2は、どちらも同じ情報を表現していますが違いがあります。

• 二次元配列は全ての要素が同じ長さである必要があります。すなわち、最長のものに合わせる必要があります。この場合はishikawaと\0の合計9個のcharを、全ての行について確保する必要があります。これは、例えばfukuiの後ろには無駄な領域を確保していることを意味します。
• 一方で、二次元配列はメモリ上では連続して確保されるため、わかりやすいという面があります。例えばhokuriku2の場合は 3 x 9 の二重forループで全体を舐めることが出来ます。

また、二次元配列では、文字列以外も初期化記法で書けることを思い出しましょう。

int arr[][3] = {
    {1, 2, 3},
    {4, 5, 6}
};

この書き方はポインタ配列ではできませんでしたね。 このように、ポインタ配列と二次元配列は長所短所があるので、必要に応じて使い分けましょう。

二次元配列をポイントする

さて、「二次元配列のポインタ」はどうなるのでしょうか？ここは非常に難しいので、わからなければ、今のところは結果だけ知っておいてください。 次のような奇妙な表記になります。Reese, p104

int (*p)[3];  // 「int 3つ分」を1単位とするポインタ
p = arr;      // p と p + 1 の間の値（アドレス）の差は12（つまりint3つ分） 

ここでint (*p)[3];という新しい表記が登場しました。ここでは、いつも通りpというポインタを１つ作っています。 そしてそのpは「int ３つ分」を一単位とするポインタです。カッコを忘れないでください。 カッコを忘れてint *p[3]としてしまうと上で習ったポインタ配列になってしまいます。 pはポインタなので、p = arrという表記でarrの先頭アドレスを代入出来ます。

さて、「int３つ分を単位とする」とはどういうことでしょうか？普通のintポインタ（int *ip;）の場合、「int１つを単位とする」ので、 ip + 1 とipの値（アドレス）の差は4でした。ここでは、pはint３つを単位とするため、p + 1とpの値（アドレス）の差が12になっています。

さて、なぜそんな奇妙なポインタが出てくるのでしょうか？これは二次元配列の構造を考えるとわかります。 2x3の二次元配列int arr[2][3]について、&arr[0]とは、要素３のint配列（の先頭アドレス）ですね。 そして、&arr[1]は、次の「要素３のint配列（の先頭アドレス）」ですね。 以上より、&arr[0]とarr[1]のアドレスの差は、4 * 3 = 12になることがわかります。4バイト * 三つ分です。 すなわち、&arr[0]と&arr[1]の差は12です。

さて、pがarrを指すポインタだとしましょう。ここでポインタと配列の基本ルールを思い出しましょう。

• ポインタpが配列arrを指しているとする
• p + 1、&p[1]、arr + 1、&arr[1]は全て同じものを指す。すなわち、pもarrも同じように扱える

この基本ルールにより、&arr[0]と&arr[1]の差が12である以上、pとp + 1の値（アドレス）の差も12である必要があります。 なので、arrを指すポインタは、「12バイト = int * 3」を単位とするポインタである、ということになります。 この関係を下図にまとめます。

printf("%p %p\n", &arr[0], &arr[1]); // 0x7fff571015f0 0x7fff571015fc  つまり差は12
printf("%p %p\n", p, p + 1);         // 0x7fff571015f0 0x7fff571015fc

ここでの結論は、

• 型Tで幅Wの二次元配列arrを考える。T arr[][W];
• このarrへのポインタpは、T (*p)[W] = arr;と書ける
• このポインタはsizeof(T) * Wを単位とする。
• Wは省略できない

ということです。

二次元配列を引数に持つ関数

さて、二次元配列を関数引数にとるとどうなるでしょうか？その例を以下に示します。

double avg_mat(int height, int width, int mat[][width]) {
    double sum = 0.0;
    for (int h = 0; h < height; ++h) {
        for (int w = 0; w < width; ++w) {
            sum += mat[h][w];
        }
    }
    return sum / (height * width);
}

int main () {
    int arr[][3] = {
        {1, 2, 3},
        {4, 5, 6}
    };
    printf("avg of arr: %f\n", avg_mat(2, 3, arr));  // 3.5
}

avg_matは二次元配列の要素の平均値を出すものです。ここでは以下のポイントがあります

• int mat[][width]のように、一つ目の添字（高さ）は省略できるが、二つ目の添え字（幅）は省略できない
• しかし、幅パラメータは同時にint widthといった形で与えることができる（C99のVariable-Length Arrayの機能, King, p200）

ここでは、結論としては、上記のようにポインタを明示的に表記せずに二次元配列を引数にとれます。 なので実用上はこの場合はそこまで混乱することはないです。

コラム

C99以前では、上記を達成するためには例えば次のようにdefineなどを用いるのが普通であり、非常にめんどくさかったのです。

#define WIDTH 3

double avg_sum(int height, int mat[][WIDTH]) { ... }

さて、関数引数の配列は実はポインタであることを何度も勉強しましたね。関数引数においてint x[]はint *xであるというルールを適用すると、上記の関数は次と同じです。

double avg_mat(int height, int width, int (*mat)[width]) { ... }

ということで、ここでは先ほど勉強した、「sizeof(int) * width」を単位としたポインタが登場しました。 ここでは二次元配列をポインタで受け取っているということなので、まさにその通りになりましたね。

さて、二次元配列をそのまま扱うのではなく、一旦それをほぐし、単純な一次元配列だと解釈して処理することもできます。

double avg_mat2(int *mat, int num_element) {
    double sum = 0.0;
    for (int i = 0; i < num_element; ++i) {
        sum += mat[i];
    }
    return sum / num_element;
}

このavg_mat2は処理の結果はavg_matと同じですが、データの与え方が単純なintポインタになっています。 そして、総合計要素分だけ値をスキャンして、平均を計算します。

printf("%f\n", avg_mat2(&arr[0][0], 2 * 3));  // OK

ここではarr[0][0]は二次元配列の先頭要素であり、型はintですね。そのアドレスである&arr[0][0]は、intへのポインタであるint *で指し示せます。 なので、これを引数に与えることで計算できます。これは定義通りですね。

さて、実はこれは以下のように書くこともできます。

printf("%f\n", avg_mat2(arr[0], 2 * 3));      // OK

なぜ&arr[0][0]ではなくarr[0]でもOKかというと、

• ショートアンサー：
  • arr[0]は、二次元配列であるarrの最初の行を表します。つまり、intの配列を表します。intの配列はintポインタであるint *に渡せるので、これでOKです。
• 別の説明：
  • 任意の配列vについて、&v[3]のような「添え字アクセス + &」は、v + 3のような「変数名 + オフセット」と同じ、というルールを思い出しましょう。
  • これを適用すると、&arr[0][0]はarr[0] + 0と同じということがわかります。この最後の0を消すことで、arr[0]となります。
  • ちなみに細かいですが、arr[0] + 0はポインタですが、arr[0]は配列です。この辺を詳しく知りたい人は次章を見てください。

コラム

さて、&arr, arr, &arr[0], arr[0]は全て、値としては&arr[0][0]と同じ（配列の先頭要素のアドレス）ですが、型が違います。 なので、以下を実行すると計算は実行できて答えも正しいのですが、「ポインタの型が違っているよ」と警告が出ます。 これについて詳しく知りたい人は次章をよく読んでみてください。

printf("%f\n", avg_mat2(&arr, 2 * 3));
printf("%f\n", avg_mat2(arr, 2 * 3));    
printf("%f\n", avg_mat2(&arr[0], 2 * 3));    

やってみよう（20分）

• 上記を写経してみましょう
• double avg_mat(int mat[][width], int height, int width) { ... }のように書くとエラーになります。これはwidthを宣言する前にmat中でwidthを使ってしまっているからです。これを実際に確認してみてください。

クイズ

ポインタ配列のクイズでswapを書きました。この二次元配列版を書いてみましょう。

答え

例えば以下。配列ポインタ版に比べ、多次元配列版では幅の情報（max_len）を追加で入力する必要があります。また、ポインタのコピーといったことも出来ないので、下では愚直に一文字ずつスワップしています。

void swap2(int max_len, char prefectures[][max_len], int i, int j) {
    for (int len = 0; len < max_len; ++len) {
        char tmp = prefectures[i][len];
        prefectures[i][len] = prefectures[j][len];
        prefectures[j][len] = tmp;
    }
}

int main () {
    char hokuriku2[][9] = {"toyama", "ishikawa", "fukui"};
    swap2(9, hokuriku2, 1, 2);  // "toyama", "fukui", "ishikawa"
}

二次元配列とポインタのチートシート

二次元配列の型は複雑なので、以下のようにまとめました。 以下をmain.cppとして保存し、コンパイルすると、手元で実行出来ます。 型情報は例えばデバッガを使えば逐一調べられますが、今回は明示的にプリントすることでまとめました。これらを暗記する必要はないですが、困ったときには調べられるようにしておいてください。

ちなみに、Google Cloud Shell Editor上では、型を間違えて警告（変数の下に破線）が出来るときは、そこにマウスオーバーすることで型の情報を得ることができると思います。

#include <typeinfo>
#include <stdio.h>

// 一次元配列と二次元配列について、型の詳細。
// typeinfoというc++の機能を使えば、以下のように型を表示できる。
// typeid(変数).name() で型が表示される。
// そのままだと見にくいので、c++filtというコマンドで見やすく出来る
//
// How to run:
// (1) コンパイル。gccではなくてg++なので注意。つまりこれはcではなくc++プログラム
//   $ g++ main.cpp
// (2) 実行。パイプでc++filtというプログラムに渡すと、見やすくなる。
//    （ただし、文字列置換を行うので、普通のprintfの出力がおかしくなるかもしれないので注意）
//   $ ./a.out | c++filt -t 


int main(int argc, char *argv[]) {
    // =================== 一次元配列の場合 =================== 
    int a[] = {1, 2, 3};
    int *p;
    p = a;

    // 次のように、配列の型は「int [3]」、int型へのポインタは「int*」となる
    printf("=== 0 ===\n");
    printf("%s\n", typeid(a).name());  // int [3]
    printf("%s\n", typeid(p).name());  // int*

    // 配列とポインタは同じように扱えると習ったが、それでも、
    // int [3] と int*は型が違うじゃないかと思うかもしれない。
    // これはルールなので仕方ないのだが、一つの説明として、
    // a + 0 を考えてみる。 a + 0 の型は実はポインタと同じ int* である。
    // a + 0 も a も似たようなものだと考えると、この事実から、
    // 配列をポインタと同じように扱えることがわかるかもしれない。
    // 配列が関数に渡されるときは、こっちになっていると思えばOK。
    // 逆に、複雑な配列に対してそれをポイントする型がわからないときは、
    // "+0"して型を調べてみるといいかもしれない。
    printf("%s\n", typeid(a + 0).name());  // int*

    // 「a」だけの場合、それは実体としての配列（intを3つ分確保）を指している。
    // sizeofも12になる。 一方、「a + 0」とした場合、
    // それをポインタだとコンパイラは解釈する。sizeofも8になる。
    // この二つはポインタを経由して扱う場合は同じもの。なので、
    // int *p = aと出来る。
    // 結論として、型が「int [3]」のものは、「int*」でポイント出来る。
    printf("%p %p\n", a, a + 0);  // 値は同じ。 0x7ffd864060b4 0x7ffd864060b4
    printf("%lu\n", sizeof(a));      // sizeofは違う：12
    printf("%lu\n", sizeof(a + 0));  // sizeofは違う：8


    // 「a + 1」のように指定すると、「p + 1」と全く同じにふるまう
    printf("=== 1 ===\n");
    printf("%s\n", typeid(a + 1).name());  // int* 
    printf("%s\n", typeid(p + 1).name());  // int*
    printf("%s\n", typeid(*(a + 1)).name());  // int
    printf("%s\n", typeid(*(p + 1)).name());  // int


    // *(a + i) とa[i]は同じ
    printf("=== 2 ===\n");
    printf("%s\n", typeid(&a[1]).name());  // int*  ちなみにこれは &(a[1])のこと
    printf("%s\n", typeid(&p[1]).name());  // int*
    printf("%s\n", typeid(a[1]).name());  // int
    printf("%s\n", typeid(p[1]).name());  // int


    printf("=== 3 ===\n");
    // &aと&pは型も値も違う（week 4を参照）

    // &aは何かを考えることは難しい。結論はint (*)[3] である。
    // すなわち、「int 3つ」を単位とするポインタである。
    // それはなぜだろうか？例えば変数一つの場合を考える。
    // int n;
    // int *p;
    // p = &n;
    // このように、「int 1つ」であるnのアドレス（&n）を代入すべきポインタは、
    // int *である。すなわち、「int 1つ」を単位とするポインタである。
    // この議論を拡張し、int a[3]; とすると、
    // 「int 3つ」であるaのアドレス（&a）を代入すべきポインタは、
    // int (*)[3] すなわち、「int 3つ」を単位とするポインタである、と言える
    printf("%s\n", typeid(&a).name());  // int (*)[3]

    // &pは単純にポインタのさらにアドレスなので、ポインタのポインタ
    printf("%s\n", typeid(&p).name());  // int**

    // &aはaと同じ。&pはポインタ変数そのもののアドレス
    printf("%p %p\n", &a, &p);  // 違う: 0x7ffd864060b4 0x7ffd86406098


    // =================== 二次元配列の場合 =================== 
    int arr[][3] = {{1, 2, 3}, {4, 5, 6}};
    int (*p2)[3];
    p2 = arr; // p2はarrを二次元配列として指すポインタ。

    // 授業で習った通り、arrは二次元配列 int [2][3] で、
    // p2は「int 3つ分を単位とするポインタ」
    // では、int [2][3] を指すポインタが int (*)[3]だと
    // どうすればわかるだろうか？
    printf("=== 4 ===\n");
    printf("%s\n", typeid(arr).name());  // int [2][3]
    printf("%s\n", typeid(p2).name());   // int (*)[3]

    // 一次元配列の例で見た通り、0を足すと、コンパイラは
    // arrをポインタとして解釈する。そのときの型を見ると、
    // それは、arrをポインタとして指すときの型。
    // なので、この場合、p2と型が同じなので、OKということになる。
    printf("%s\n", typeid(arr + 0).name());  // int (*)[3]


    // 添え字を使って各行にアクセスする場合。これは一次元配列そのものになる
    printf("=== 5 ===\n");    
    printf("%s\n", typeid(arr[0]).name());  // int [3]
    printf("%s\n", typeid(p2[0]).name());   // int [3]
    printf("%s\n", typeid(arr[0] + 0).name()); // int*
    printf("%s\n", typeid(p2[0] + 0).name());  // int*


    // arrやp2に&をすると次のようになる。値も異なる。
    printf("=== 6 ===\n");    
    printf("%s\n", typeid(&arr).name());  // int (*)[2][3]
    printf("%s\n", typeid(&p2).name());  // (**)[3]



    printf("=== 7 ===\n");    
    // 以下は全てアドレスの値は同じだが、型は異なる
    // 0x7ffcc93804a0 0x7ffcc93804a0 0x7ffcc93804a0 0x7ffcc93804a0 0x7ffcc93804a0
    printf("%p %p %p %p %p\n", &arr, arr, &arr[0], arr[0], &arr[0][0]);

    // &arrは、二次元配列そのものを一単位として扱う。よって、
    // 「intが2 * 3個」を一つの単位とする。
    // これは一次元配列名に&をしたときと同じ。難しいが、普段は使わない
    printf("%s\n", typeid(&arr).name());       // int (*)[2][3]   

    // arrは二次元配列そのもの。これはそのまんま。
    printf("%s\n", typeid(arr).name());        // int [2][3]

    // 再度の注意として、arrに0を足すと、これをポインタとして解釈してくれる
    printf("%s\n", typeid(arr + 0).name());    // int (*)[3]

    // arr[0]は一行目の一次元配列そのもの
    printf("%s\n", typeid(arr[0]).name());    // int [3]

    // 再度の注意として、arr[0]に0を足すと、これをポインタとして解釈してくれる
    printf("%s\n", typeid(arr[0] + 0).name());    // int*

    // &arr[0]は難しい。考え方として、次の二通りがあると思う
    // (1) &arr[0]（アンド＋添え字アクセス）は arr + 0（変数名＋オフセット）と同じ
    //     というルールを適用すると、上のarr + 0になる、と覚える
    // (2) arr[0]は一次元配列そのものだった。一次元配列そのものに&をつけると、
    //     配列長を単位とするポインタになったことを思い出す
    //     （上のほうの === 3 === のところを参照）
    printf("%s\n", typeid(&arr[0]).name());    // int (*)[3]

    // arr[0][0]はintなので、そのポインタ。
    printf("%s\n", typeid(&arr[0][0]).name()); // int*


    // =================== 三次元配列の場合も同様 =================== 
    int arr2[2][3][4] = {
        {{ 1,  2,  3,  4}, { 5,  6,  7,  8}, { 9, 10, 11, 12}},
        {{13, 14, 15, 16}, {17, 18, 19, 20}, {21, 22, 23, 24}}
    };
    int (*p3)[3][4];
    p3 = arr2; // p3はarr2を三次元配列として指すポインタ。
    printf("=== 8 ===\n");
    printf("%s\n", typeid(arr2).name());      // int [2][3][4]
    printf("%s\n", typeid(p3).name());        // int (*)[3][4]
    printf("%s\n", typeid(arr2 + 0).name());  // int (*)[3][4]

    return 0;
}

コマンドライン引数

さて、次はコマンドライン引数について見ていきましょう。

コマンドライン引数の基本

コマンドライン引数とは、プログラムを実行するときに追加で与える情報です。例えば次を考えましょう。

#include <stdio.h>

int main (int argc, char *argv[]) {
    return 0;
}

ここで、main関数の引数として次の二つを与えています

• int argc: コマンドライン引数の個数
• char *argv[]：コマンドライン引数の内容を表すポインタ配列

これをコンパイルして次のように実行したとします。

$ ./a.out hoge 123

すると、argcおよびargvには次のようにコマンドライン引数情報が入ります。

すなわち、argvの各要素には引数が入り、その個数がargcで与えられます。 例えばargv[1]とすると、それは"hoge"を意味します。 ルールとして、argv[0]にはプログラム名そのものが入ります。 また、argvの最後の要素はNULLになっていると決まっています。

このargcとargvを使うと、プログラムに対して簡単に情報を入力することができます。

コラム

argc, argvの変数名は実は任意です。以下でも大丈夫です

int main (int kosuu, char *hikisuu[]) { ... }

また、これまでに習った通り、関数引数では*argv[]は**argvと同じなので、以下のようにも書けます。

int main (int argc, char **argv) { ... }

ここでargvをポインタ配列にするか、ポインタのポインタにするかは、好みです。

コマンドライン引数を使う例

それでは例を見てみましょう。まず、入力をそのまま出力する、echoコマンドを思い出しましょう。

$ echo hoge 123
hoge 123

これを真似したmyechoを作ってみましょう。

以下をmyecho.cとします。

#include <stdio.h>
int main(int argc, char *argv[]) {
    for (int i = 1; i < argc; ++i) {
        printf("%s ", argv[i]);
    }
    printf("\n");
}

forの開始が1からになっている点に注意してください。これはargv[0]はコマンドそのものになるため省いているからです。 これをmyechoという名前でコンパイルします

$ gcc myecho.c -o myecho

すると、次のように使うことができます。

$ ./myecho hoge 123
hoge 123

ここでは、argv[0]は必ずプログラム名なので、それはプリントしていません。

また、足し算の例も見てみましょう。次をsum.cとして保存します。

#include <stdio.h>
#include <stdlib.h>

int main(int argc, char *argv[]) {
    int nums[argc - 1];
    for (int i = 1; i < argc; ++i) {
        nums[i - 1] = atoi(argv[i]);
    }

    int sum = 0;
    for (int i = 0; i < argc - 1; ++i) {
        sum += nums[i];
    }
    printf("%d\n", sum);
}

ここでは整数を読み取り、合計値を返しています。atoiはstdlib.hで定義されている便利関数で、文字列が整数と解釈出来る場合は整数として読み込んでくれます。 これをコンパイルすると、次のように使えます。

$ ./sum 12 23 34
69

ここで、整数以外を入力するとatoiでおかしくなるので注意してください。そのほかに、少数を読み取ってくれるatofなどもあります。

やってみよう（20分）

上記を写経してみましょう。

クイズ

コマンドライン引数の中に"toyama", "ishikawa", "fukui"のどれかが含まれていたら、"hokuriku!"出力するプログラムを書いてください。仮に引数中にたくさん北陸の県が含まれていても、一回だけ"hokuriku!"と言うようにしてください。この際、文字列の一致を判定するstrcmpは使ってもOKです（string.hをincludeする必要があります）

答え

例えば以下です。

#include <stdio.h>
#include <string.h>

int is_hokuriku(const char *prefecture) {
    char *hokuriku[] = {"toyama", "ishikawa", "fukui"};
    for (int i = 0; i < 3; ++i) {
        if (strcmp(prefecture, hokuriku[i]) == 0) {
            return 1;
        }
    } 
    return 0;
}

int main(int argc, char *argv[]) {
    for (int i = 1; i < argc; ++i) {
        if (is_hokuriku(argv[i])) {
            printf("hokuriku!\n");
            return 0;
        }
    }
    return 0;
}

これを使った例が以下です。

$ a.out aichi shizuoka hokkaido
$ a.out nagano ishikawa
hokuriku!
$ a.out fukui ishikawa
hokuriku!

宿題

それでは今週の宿題です。

• 締切は次の授業の前日の深夜23:59までです（今回の場合、2022/11/9, 23:59）
• 宿題リンクをslackで配付します。このリンクは公開しないでください。

week5_1

LU分解を実装しましょう。 LU分解とは、ある正則行列\(\mathbf{A} \in \mathbb{R}^{D \times D}\)を、以下のように「単位下三角行列\(\mathbf{L}\)」と「上三角行列\(\mathbf{U}\)」の積に分解することです。

\[ \mathbf{A} = \mathbf{L}\mathbf{U} =\begin{bmatrix} 1 & & & \\ l_{21} & 1 & & \\ \vdots & \ddots & \ddots & \\ l_{n1} & \cdots & l_{n, n-1} & 1 \\ \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ & u_{22} & \cdots & u_{2n} \\ & & \ddots & \vdots \\ & & & u_{nn} \\ \end{bmatrix} \]

ここで、下三角行列とは、 対角成分からみて右上が全て0になっている行列です。 単位下三角行列とは、下三角行列でかつ対角成分が1のものです。上三角行列も同様に定義され、左下が0のものです。 LU分解とは、上のように、行列を二つの三角行列に分ける処理を言います。

LU分解は、実はみなさんが大学１年生のときにならった「ガウスの消去法」そのものを行列の形で表記したものになっています。 ということでまずガウスの消去法を復習しましょう。変数\(\mathbf{x}\)に関して、 \(\mathbf{A}\mathbf{x} = \mathbf{b}\)という形で方程式が与えられたとします。\(\mathbf{A}\)と\(\mathbf{b}\)は既知で、\(\mathbf{x}\)は未知です。 ここで、逐次的に\(\mathbf{A}\)および\(\mathbf{b}\)を変形させていき\(\mathbf{A}\)を三角形の形にして、そのうえで\(\mathbf{x}\)を求めていくのがガウスの消去法です。 ここでは繰り返し処理をするので、\(\mathbf{A}\)をまず\(\mathbf{A}^{(1)}\)と呼びましょう。

\[ \mathbf{A}^{(1)}= \begin{bmatrix} 3 & 2 & 1 \\ 2 & 3 & 1 \\ 1 & 1 & 4 \\ \end{bmatrix} \]

ここで、

• 第1行目を2/3倍（\(a_{21} / a_{11}\)倍）したものを第2行目から引きます
• 第1行目を1/3倍（\(a_{31} / a_{11}\)倍）したものを第3行目から引きます

この処理をした後のものを\(\mathbf{A}^{(2)}\)と呼びましょう。

\[ \mathbf{A}^{(2)}=\begin{bmatrix} 3 & 2 & 1 \\ 2 - 3 \cdot 2/3 & 3 - 2 \cdot 2/3 & 1 - 1 \cdot 2/3\\ 1 - 3 \cdot 1/3 & 1 - 2 \cdot 1/3 & 4 - 1 \cdot 1/3\\ \end{bmatrix}=\begin{bmatrix} 3 & 2 & 1 \\ 0 & 5/3 & 1/3 \\ 0 & 1/3 & 11/3 \\ \end{bmatrix} \]

次に、

• 第2行目を1/5倍（\(a_{32} / a_{22}\)倍）したものを第3行目から引きます

この処理をした後のものを\(\mathbf{A}^{(3)}\)と呼びましょう。

\[ \mathbf{A}^{(3)}=\begin{bmatrix} 3 & 2 & 1 \\ 0 & 5/3 & 1/3 \\ 0 & 1/3 - 5/3 \cdot 1/5 & 11/3 - 1/3 \cdot 1/5\\ \end{bmatrix}=\begin{bmatrix} 3 & 2 & 1 \\ 0 & 5/3 & 1/3 \\ 0 & 0 & 54/15\\ \end{bmatrix} \]

というわけで、\(\mathbf{A}\)を三角形にすることが出来ました。

ここで「行を定数倍して引く」という処理はそのような行列を左からかければ達成できます。なので、\(\mathbf{b}\)も同様 に肩に載せる数字で変形回数を表現するとすると、 \(\mathbf{A}^{(1)}\mathbf{x} = \mathbf{b}^{(1)}\) に対し「定数倍して引く」行列を左からかけて \(\mathbf{A}^{(2)}\mathbf{x} = \mathbf{b}^{(2)}\)とできます。 さらにもう一度「定数倍して引く」行列を左からかけて \(\mathbf{A}^{(3)}\mathbf{x} = \mathbf{b}^{(3)}\) という形の式を得ます。\(\mathbf{A}^{(3)}\)は三角形なので、三行目に関して変数は\(x_3\)のみです。 よってこれは簡単に求めることができます。それを順番に代入していくことで\(x_2, x_1\)もわかります。結果として\(\mathbf{x}\)を得ることができます。

さて、LU分解に話を戻します。実はなんともうLU分解が完了しています。

\[ \mathbf{A}^{(1)}= \begin{bmatrix} 3 & 2 & 1 \\ 2 & 3 & 1 \\ 1 & 1 & 4 \\ \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 2/3 & 1 & 0 \\ 1/3 & 1/5 & 1 \\ \end{bmatrix}\begin{bmatrix} 3 & 2 & 1 \\ 0 & 5/3 & 1/3 \\ 0 & 0 & 54/15 \\ \end{bmatrix} \]

上記のように、ガウスの分解の過程で現れた係数をマイナスをつけて並べることで、\(\mathbf{L}\)を構成することができます。

それはなぜでしょうか？まず、\(\mathbf{A}^{(1)}\)を\(\mathbf{A}^{(2)}\)に変形した部分を思い出しましょう。 この変形をじーっと眺めていると、以下のように行列積の形で書き下せることがわかります。

\[ \mathbf{A}^{(1)}= \begin{bmatrix} 3 & 2 & 1 \\ 2 & 3 & 1 \\ 1 & 1 & 4 \\ \end{bmatrix} = \begin{bmatrix} 1 \\ 2/3 \\ 1/3 \\ \end{bmatrix}\begin{bmatrix} 3 & 2 & 1 \\ \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 \\ 0 & 5/3 & 1/3 \\ 0 & 1/3 & 11/3 \\ \end{bmatrix} \]

ここで

\[ \begin{bmatrix} 1 \\ 2/3 \\ 1/3 \\ \end{bmatrix}\begin{bmatrix} 3 & 2 & 1 \\ \end{bmatrix} \]

の部分はベクトルの外積であり、行列を構成します。 上で述べていた

• 第1行目を2/3倍（\(a_{21} / a_{11}\)倍）したものを第2行目から引く
• 第1行目を1/3倍（\(a_{31} / a_{11}\)倍）したものを第3行目から引く

というのは、この外積の形で表現できています。 そして、

\[ \begin{bmatrix} 0 & 0 & 0 \\ 0 & 5/3 & 1/3 \\ 0 & 1/3 & 11/3 \\ \end{bmatrix} \]

の部分は、\(\mathbf{A}^{(2)}\)の右下の2x2の部分です。外積部分が引き算部分を構成するので、\(\mathbf{A}^{(1)}\)からそれを引いたものはまさに\(\mathbf{A}^{(2)}\)の残りの部分そのものということです。 ここは全然難しいことは言っておらず、よく眺めるとそうなっているということです。

あとはこれを再帰的に解くことができます。この2x2の部分について、

\[ \begin{bmatrix} 0 & 0 & 0 \\ 0 & 5/3 & 1/3 \\ 0 & 1/3 & 11/3 \\ \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \\ 1/5 \end{bmatrix} \begin{bmatrix} 0 & 5/3 & 1/3 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 54/15 \\ \end{bmatrix} \]

となっています。これは

• 第2行目を1/5倍（\(a_{32} / a_{22}\)倍）したものを第3行目から引く

に相当します。

これらをまとめると、

\[ \mathbf{A}^{(1)}= \begin{bmatrix} 3 & 2 & 1 \\ 2 & 3 & 1 \\ 1 & 1 & 4 \\ \end{bmatrix} = \begin{bmatrix} 1 \\ 2/3 \\ 1/3 \\ \end{bmatrix}\begin{bmatrix} 3 & 2 & 1 \\ \end{bmatrix} + \begin{bmatrix} 0 \\ 1 \\ 1/5 \end{bmatrix} \begin{bmatrix} 0 & 5/3 & 1/3 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 54/15 \\ \end{bmatrix}　= \begin{bmatrix} 1 & 0 & 0 \\ 2/3 & 1 & 0 \\ 1/3 & 1/5 & 1 \\ \end{bmatrix}\begin{bmatrix} 3 & 2 & 1 \\ 0 & 5/3 & 1/3 \\ 0 & 0 & 54/15 \\ \end{bmatrix} \]

となります（「外積の和」の形は、上記のように束ねて一つの行列積にすることができます）

ひとたび行列\(\mathbf{A}\)を\(\mathbf{A}=\mathbf{L}\mathbf{U}\)と分解できると、次のようにすることで方程式\(\mathbf{A}\mathbf{x}=\mathbf{b}\)を解くことができます。

• LU分解できているので、解きたい方程式は\(\mathbf{L}\mathbf{U}\mathbf{x} = \mathbf{b}\)である。
• ここで、新たに変数\(\mathbf{y}\)を考える。\(\mathbf{L}\mathbf{y} = \mathbf{b}\)を解くことで、\(\mathbf{y}\)を得る。これは\(\mathbf{L}\)が三角形なので順番に簡単に求まる。
• 得られた\(\mathbf{y}\)を元に、\(\mathbf{U}\mathbf{x}=\mathbf{y}\)を解く。これも、\(\mathbf{U}\)が三角形なので、代入していって簡単に解ける。

というわけで、LU分解を用いると、上記のように方程式を簡単に解くことができます。

ガウスの消去法に対しLU分解が優れているのは、「同じ\(\mathbf{A}\)に対し、何度も方程式を解く場合。すなわち、たくさんの「\(\mathbf{x}, \mathbf{b}\)」が与えられる場合」です。 実はガウスの消去法の計算量（計算回数をおおざっぱに見積もったものです。３年生で習います）と、LU分解を実行する計算量は同じです（\(O(D^3)\)という量になります）。 そのうえで、上記のように\(\mathbf{L}\mathbf{U}\mathbf{x} = \mathbf{b}\)を解く計算量は分解の計算量よりずっと小さいです（\(O(D^2)\)という量になります）。

なので、同じ\(\mathbf{A}\)に対し、いろいろな\(\mathbf{x}, \mathbf{b}\)を解く場合、

• ガウスの消去法：毎回\(O(D^3)\)の時間が必要
• LU分解：最初に\(O(D^3)\)で分解を実行する必要があるが、二度目以降は\(O(D^2)\)のコストでOK

となり、便利なわけです。

ちなみに、LU分解に関する書籍や論文や解説記事はたくさんあるのですが、以下の本がコンパクトでよかったです。

• 20世紀のトップ10アルゴリズム、金田行雄・笹井理生 監修、張紹良 編、共立出版、2022、「第四章　行列計算に対する分解アプローチ」

さて、本題の宿題に移りましょう。

今日は、関数を0から書いてみましょう。main関数は変更しないでください。 main.cをコンパイルすると、コンパイルが通りません。

$ gcc main.c
main.c: In function ‘main’:
main.c:98:5: warning: implicit declaration of function ‘compute_L_U’ [-Wimplicit-function-declaration]
   98 |     compute_L_U(A, L, U);
      |     ^~~~~~~~~~~
main.c:106:5: warning: implicit declaration of function ‘matmul’ [-Wimplicit-function-declaration]
  106 |     matmul(L, U, A_reconstructed);
      |     ^~~~~~
/usr/bin/ld: /tmp/ccqDnEsr.o: in function `main':
main.c:(.text+0x368): undefined reference to `compute_L_U'
/usr/bin/ld: main.c:(.text+0x3c0): undefined reference to `matmul'
collect2: error: ld returned 1 exit status

これは、関数compute_L_Uおよびmatmulが書かれていないからです。まずこれらの関数の雛形を書いてコンパイルを通るようにしましょう。 そのうえで、以下のように実行します。この表記により、mat1.txtというテキストファイルを読み込み、そのなかに入っている情報を行列double A[D][D]の中にいれてくれます。 以下、今回の宿題では行列はすべてdoubleの二次元配列で表されるとします。

$ gcc main.c
$ ./a.out mat1.txt
A:
[3.00, 2.00, 1.00, 
 2.00, 3.00, 1.00, 
 1.00, 1.00, 4.00, ]
L:
[0.00, 0.00, 0.00, 
 0.00, 0.00, 0.00, 
 0.00, 0.00, 0.00, ]
U:
[0.00, 0.00, 0.00, 
 0.00, 0.00, 0.00, 
 0.00, 0.00, 0.00, ]
A_reconstructed:
[0.00, 0.00, 0.00, 
 0.00, 0.00, 0.00, 
 0.00, 0.00, 0.00, ]

あとは、compute_L_Uとmatmulを書いてみてください。

• compute_L_U: 行列Aを受け取り、LU分解を実行する。その結果をL, Uに書き込む。今回は上記の説明をもとに愚直に実装してOK
• matmul: 二つの行列をうけとり、行列積を計算し、書き込む。

もしLU分解が成功していれば、最後のA_reconstructedは元のAと同じになるはずです。

自動採点結果は以下です。

• 入力

  $ ./a.out mat1.txt
  

• 出力

  A:
  [3.00, 2.00, 1.00, 
   2.00, 3.00, 1.00, 
   1.00, 1.00, 4.00, ]
  L:
  [1.00, 0.00, 0.00, 
   0.67, 1.00, 0.00, 
   0.33, 0.20, 1.00, ]
  U:
  [3.00, 2.00, 1.00, 
   0.00, 1.67, 0.33, 
   0.00, 0.00, 3.60, ]
  A_reconstructed:
  [3.00, 2.00, 1.00, 
   2.00, 3.00, 1.00, 
   1.00, 1.00, 4.00, ]
  

本採点では、別の行列\(\mathbf{A}\)（mat2.txt）を使います。

発展：本文中の#define D 3を書き換えて大きな行列に対して実行するかもしれないので、\(D\)が3以外でも実行できるようにしておくと良いです。

答えの例

答え

#include <stdio.h>
#include <stdlib.h>

#define D 3

// 関数は自由に追加していいです。


// 行列を表示
void print_mat(const double mat[][D]) {
    printf("[");
    for(int j = 0; j < D; ++j) {
        if (j != 0) {
            printf(" ");
        }
        for (int i = 0; i < D; ++i) {
            printf("%3.2f, ", mat[j][i]);
        }
        if (j != D - 1) {
            printf("\n");
        }
    }
    printf("]\n");
}



void delete_one_col(double A[][D], double L[][D], int n) {
    // n行目を定数倍したものを各行から引いて消していく。ガウス消去法そのもの
    for (int j = n + 1; j < D; ++j) { // n行目を定数倍したものを適用するので、直後のn+1行目から始まる
        double m = A[j][n] / A[n][n]; // 定数倍の値
        for (int i = 0; i < D; ++i) {
            A[j][i] -= A[n][i] * m;   // 引いていく
        }
        L[j][n] = m; // 構成のルールより、定数倍の値をLに埋めていく
    }
    L[n][n] = 1.0;  // 対角は1にしておく

}

void compute_L_U(const double A[][D], double L[][D], double U[][D]) {
    // Aはオリジナルなので変更しない。Aの内容を全てUに移し、Uを変形させていく。
    for(int j = 0; j < D; ++j) {
        for (int i = 0; i < D; ++i) {
            U[j][i] = A[j][i];
        }
    }

    // 列を消す（下のほうにゼロ詰めして階段にする）、を繰り返す
    for (int n = 0; n < D; ++n) {
        delete_one_col(U, L, n);
    }

}


void matmul(const double src1[][D], const double src2[][D], double dst[][D]) {
    for (int j = 0; j < D; ++j) {
        for (int i = 0; i < D; ++i) {
            dst[j][i] = 0.0;
            for (int n = 0; n < D; ++n) {
                dst[j][i] += src1[j][n] * src2[n][i];
            } 
        }
    }
}



int main(int argc, char *argv[]) {
    // ==================== main関数は触らない ==================
    if (argc != 2) {
        printf("Error. Select an argment. For example: $ ./a.out mat1.txt\n");
        return 1;
    }

    // matを読み込む
    FILE *fp;
    fp = fopen(argv[1], "r");
    if (fp == NULL) {
        printf("Cannot open %s\n", argv[1]);
        return 1;
    }
    double A[D][D];
    for (int j = 0; j < D; ++j) {
        for(int i = 0; i < D; ++i) {
            fscanf(fp, "%lf", &A[j][i]);
        }
    }
    fclose(fp);

    // 入力の表示
    printf("A:\n");
    print_mat(A);

    // LU分解を実行して表示
    double L[D][D];
    double U[D][D];
    compute_L_U(A, L, U);
    printf("L:\n");
    print_mat(L);
    printf("U:\n");
    print_mat(U);

    // 計算したL, UをもとにAを再構成して検算
    double A_reconstructed[D][D];
    matmul(L, U, A_reconstructed);
    printf("A_reconstructed:\n");
    print_mat(A_reconstructed);

    return 0;
}