ソフトウェアII 第4回（2023/12/21）

本日のメニュー

  • 関数へのポインタ
  • 分割コンパイル
  • Makefile
  • 連続関数の最適化

準備

必要なファイルをダウンロードしましょう。以下はwget での例ですが、curlの場合はcurl -O に置き換えてください。

# wget
wget https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/optimize.tar.gz
wget https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/qsort.tar.gz
wget https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/scalar_calc.tar.gz
wget https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/vector_calc.tar.gz
wget https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/func1.c
wget https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/data.csv

関数へのポインタ

関数のアドレス

今日は関数へのポインタを見ていきます。一般に、プログラムを実行するとその構成要素が全てメモリ上に展開されます。つまり変数に加えて、関数の実行コードもメモリ上のどこかに置かれるので、関数が置かれているアドレスというものを考えることができます。以下のコードはプログラム中の色々な要素のアドレスを全て表示しています。

// 関数を含めて、いろいろなアドレスを見てみる
#include <stdio.h>
#include <stdlib.h>

int global = 0;

int add_one(int x){

    // 静的領域に変数を確保する
    static int count = 0;
    
    // 関数内のstatic変数	
    printf("count_in_add_one(%d): %p\n", count,&count);
    // add_one関数の引数x
    printf("x_in_add_one(%d): %p\n",count,&x);
    
    count++; // 何回呼び出されたかを記録
    return x + 1;
}

void call(int x, int (*fp)(int)){
    // call関数の引数x
    printf("x_in_call: %p\n",&x);
    // call関数の引数fp
    printf("fp_in_call: %p\n", &fp);
    
    fp(x);
}

int main(int argc, char **argv){
    int x = 10;
    int *a = (int*) malloc(sizeof(int));
    
    // main関数のアドレス
    printf("function_main: %p\n", &main);
    
    // main内の変数x
    printf("x_in_main: %p\n", &x);
    
    // argc
    printf("argc_in_main: %p\n", &argc);
    
    // argv
    printf("argv_in_main: %p\n", &argv);
    
    // グローバル変数global
    printf("global: %p\n", &global);
    
    // malloc されたa
    printf("allocated_a: %p\n", a);
    // ポインタそのもののアドレス
    
    printf("pointer_a_in_main: %p\n", &a);
    
    add_one(x);
    
    // add_one関数のアドレス
    printf("function_add_one: %p\n", &add_one);
    
    // もう一度呼ぶ
    add_one(x);
    
    // call 経由でadd_one を呼ぶ
    call(x,add_one);

    // call関数のアドレス
    printf("function_call: %p\n", &call);
    
    free(a);
    return 0;
}
Let's try
  • 上記を写経して実行してみましょう（print_address.c）。

一般的なUNIX環境において、実行する際に以下のようにsort コマンドに結果をパイプすることで、アドレス順に並べ替えができます。

# コンパイル (-O0 [オーゼロ]は最適化がかからないようにするオプション)
gcc -O0 -o print_address print_address.c
# 実行とソート
./print_address | sort -k 2

実行毎に結果は変化しますが、プログラム中の要素がそれぞれが置かれているメモリ領域に応じてある程度固まって配置されていることが読み取れます。

関数へのポインタ

ポインタはアドレスを格納する変数なので、関数に対するポインタを考えることができます。演算子の優先順位から、関数へのポインタは以下のような形式で宣言されます。

// fp は返り値がintで引数がint一つの関数へのポインタ
int (*fp)(int);
// fp は返り値がdoubleで、引数がdouble二つの関数へのポインタ
double (*fp)(double, double);

関数ポインタを宣言する際は、アドレスを受け取りたい関数の返り値の型と、引数の型を指定します。これにより同じ入力、出力を持つ関数をこのポインタでひとまとめに扱うことができます。

関数ポインタを実際に用いて関数を実行してみましょう。

#include <stdio.h>

int add_one(int x){
    return x+1;
}

int main(){
    int (*fp)(int);
    // fp に add_one関数のアドレスを代入
    fp = &add_one;
    
    int x = 1;
    // fp で指されている関数を呼び出す
    int y = (*fp)(x);
    
    printf("%d\n", y);
    
    return 0;
}

なお、関数のアドレス表示および関数ポインタは「暗黙的な型変換」が適応されるため、かなりいろいろな書き方が許容されます。例えばアドレスを表示や代入の際に以下のように書いても正しく動きます。

// 表示の例
printf("function_add_one: %p\n", add_one);
// 代入の例
fp = add_one;

これは引数なしで関数名を扱うと暗黙的に関数ポインタ型=アドレスに変換されるためです。つまり& の有無が結果に影響をあたえません。

また実行する際も以下のような記述でも動きます。

// 上記のコードで int (*fp)(x) に int add_one(int) のアドレスを代入している前提
int y = fp(x);
// 暗黙的な型変換が繰り返された結果、以下もOK（使えとはいっていない...）
int y = (*******fp)(x);

関数ポインタの使い方の一つに、関数に引数として関数のアドレスを渡して、同じ型の動作を関数ポインタに渡すアドレスを変更することで切り替えるような使い方があります。このように渡される関数をコールバック関数と呼びます。

#include <stdio.h>

int add(int x, int y) {
    return x + y;
}

int mul(int x, int y) {
    return x * y;
}

void apply(int v[], int n, int y, int (*fp)(int,int)) {
    for (int i = 0 ; i < n ; i++){
        v[i] = (*fp)(v[i],y);
    }
}

int main() {
    int v[10];
    const int n = sizeof(v)/sizeof(int);
    for (int i = 0 ; i < n ; i++){
        v[i] = i;
    }
    
    apply(v, n, 2, add);
    apply(v, n, 3, mul);
    
    for (int i = 0 ; i < n ; i++){
        printf("v[%d] = %d\n", i, v[i]);
    }
}
Let's try

上記のプログラムを写経して実行してみましょう。

クイックソート

関数ポインタによるコールバック関数をうまく使った標準ライブラリ関数の一つにqsort があります。qsort はデータの型に依存せずに、配列を決められた基準で並べ替える機能を提供します。

#include <stdlib.h>

// 先頭アドレスがbase、個々の要素がsize、個数membの配列を、比較関数compar の基準で並べ替える
void qsort(void *base, size_t memb, size_t size, int (*compar)(const void *, const void*));

関数ポインタで、要素を比較した時に、同じならば0、大小関係を正負の値で返す関数を指定します。

サンプルコード（整数配列）

標準ライブラリ関数qsort を使ったサンプルコードqsort.cをダウンロードして、コンパイル・実行してみましょう。

# wgetの場合
wget https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/qsort.c
# curlの場合
curl -O https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/qsort.c
# コンパイル
gcc -o qsort qsort.c
# 実行
./qsort 100

サンプルコード（構造体配列）

qsort はデータ型に依存しないので、例えば構造体の特定のメンバを基準にした構造体配列の並べ替えなども可能です。構造体配列に対するqsort の例をダウンロードしてみましょう。

# wgetの場合
wget https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/qsort_struct.c
# curlの場合
curl -O https://www.gavo.t.u-tokyo.ac.jp/~dsk_saito/lecture/software2/resource/qsort_struct.c
Let's try

以下に取り組んでみましょう。

  • 整数配列のソートで、整数値を降順にソートするように修正してください。
  • 構造体配列のソートで、comp_name関数を完成させ、Student を名前の辞書順でソートするように修正してください。

Let's try 一例（ソート）

[講義後追記] ソートの一例です。

整数配列のソート 基本的には昇順の時と返り値の符号を入れ替えるだけです。簡潔版のみ示します。 
int comp_int(const void *x0, const void *x1) {
    return *(int*)x1 - *(int*)x0;
}
構造体配列のソート 入力された構造体の文字列に着目し比較する。`string.h` の比較関数である`strcmp`がこの用途に合致する。 
int comp_name(const void *x0, const void *x1) {
    const Student *p0 = (Student*)x0;
    const Student *p1 = (Student*)x1;
	
    return strcmp(p0->name, p1->name);
}

分割コンパイル

これまではプログラムを作成する際、一つのファイルに全ての要素を詰め込んで記述していました。中規模以上のプログラムを作成する場合、複数のソースファイルに分けて関数等をコーディングし、コンパイル時にこれらをマージする方式でコンパイルすることがあります。これを分割コンパイルといいます。ソースファイルを分割することで、例えば汎用的な機能を別のプログラムでも使いまわしたり、またバグがあった場合にも修正が容易になるなど、保守性・再利用性の面で大きな利点があります。

コンパイルの実際

これまでのプログラムのコンパイルの流れを再確認します。例えばqsort.c をコンパイルする場合、以下のようにコンパイルしていました。

# qsort.c をコンパイル
gcc -o qsort qsort.c

実はこのコンパイルは以下に示す2つの手順を一度に行っていたことになります。

# オブジェクトファイルの作成: 個々の関数の実体ができる
gcc -c qsort.c
# リンク: オブジェクトファイルを結合して実行バイナリを作成する
gcc -o qsort qsort.o

ソースコードを分割する

まずは以下のようなコードを考えます（scalar_calc/main.c）。

# tar で解凍し、ディレクトリ移動
tar xvzf scalar_calc.tar.gz
cd scalar_calc

// 四則演算のプログラム
#include <stdio.h>

double add(double a,double b);
double sub(double a,double b);
double mul(double a, double b);
double div(double a, double b);

int main() {
    double a = 1.2;
    double b = 1.4;

    printf("a+b: %f\n",add(a,b));
    printf("a-b: %f\n",sub(a,b));
    printf("a*b: %f\n",mul(a,b));
    printf("a/b: %f\n",div(a,b));
    return 0;
}

double add(double a, double b){
    return a + b;
}

double sub(double a, double b){
    return a - b;
}

double mul(double a, double b){
    return a * b;
}

double div(double a, double b){
    return a / b;
}

このコードについて、まずは実装部分をcalc.cという別ファイルにうつします。

// main.c
#include <stdio.h>

double add(double a, double b);
double sub(double a, double b);
double mul(double a, double b);
double div(double a, double b);

int main(){
    double a = 1.2;
    double b = 1.4;

    printf("a+b: %f\n",add(a,b));
    printf("a-b: %f\n",sub(a,b));
    printf("a*b: %f\n",mul(a,b));
    printf("a/b: %f\n",div(a,b));
    return 0;
}

// calc.c

double add(double a, double b){
    return a + b;
}

double sub(double a, double b){
    return a - b;
}

double mul(double a, double b){
    return a * b;
}

double div(double a, double b){
    return a / b;
}

現状はmain.c にプロトタイプ宣言とmain関数が、calc.c に関数の実際の実装が記述されている状態です。この状態で以下のようにコンパイルするとどうなるでしょうか。

# エラーがでるはず (undefined reference / undefined symbols)
gcc main.c

またcalc.c もコンパイルしてみます。

# mainがないことについてエラーがでるはず
gcc calc.c

main.cのコンパイルについては、プロトタイプ宣言はあるけど関数の実体が定義されていないということは、calc.c ではmain の実体がないことがエラーの原因です。

これを解消するために、まずは他の関数の実体が定義されていない状態でも個々のファイルをコンパイルできるようにします。

# -c オプションを追加する
gcc -c main.c
gcc -c calc.c

このようにするとmain.o calc.o というファイルが生成されます。

最後にこれらのオブジェクトファイルを結びつけることで最終的なバイナリを作ります。これをリンクと呼びます。

# -o オプションがないので a.out が生成される
gcc main.o calc.o
# 実行する
./a.out

さらにプロトタイプ宣言を一つのファイルにまとめます。これが皆さんがよく使うヘッダファイルの基本形です。

#include <stdio.h>
#include "calc.h"

int main(){
    double a = 1.2;
    double b = 1.4;

    printf("a+b: %f\n",add(a,b));
    printf("a-b: %f\n",sub(a,b));
    printf("a*b: %f\n",mul(a,b));
    printf("a/b: %f\n",div(a,b));
    return 0;
}


// calc.h 
// 実際は後述のインクルードガードがあった方がよいが、一旦省略
double add(double a, double b);
double sub(double a, double b);
double mul(double a, double b);
double div(double a, double b);
Let's try

qsort.c をヘッダファイルintlib.hとソースファイルmain.cおよびintlib.c に適切に分割し、分割コンパイルを実行してみましょう

  • main.c には main関数と#include文のみ
  • intlib.h には 関数のプロトタイプ宣言 (load_int, comp_int)
  • intlib.c に実際の実装をかく
  • intlib.h で宣言される関数が使われる場合はincludeする

Let's try 一例（分割）

[講義後追記] 分割されたファイルの一例です。

main.c 
// main.c 
#include <stdio.h> 
#include <stdlib.h>
#include "intlib.h"

int main(int argc, char **argv){
    if(argc != 2){
        fprintf(stderr, "usage: %s <number of elements>\n",argv[0]);
        exit(EXIT_FAILURE);
    }
    
    int n = load_int(argv[1]);
    
    int *v = (int*) malloc(sizeof(int) * n);
    for(int i = 0 ; i < n ; i++){
        v[i] = rand()%1024;// 0 - 1023 で値を入れる
    }
    
    // sort前表示
    printf("[");
    for(int i = 0 ; i < n ; i++){
        printf("%d ",v[i]);
    }
    printf("]\n");
    
    // qsort実行
    qsort(v, n, sizeof(int), comp_int);
    
    // sort後表示
    for (int i = 0; i < n; i++) {
        printf("v[%d] = %d\n", i, v[i]);
    }
    
    free(v);
    return 0;
}
intlib.c 
// intlib.c
#include <stdio.h>
#include <stdlib.h>
#include <errno.h> 
#include <string.h>
#include "intlib.h" // これはこのファイルのコンパイル自体には不要だがライブラリとペアということで設置する

int comp_int(const void *x0, const void *x1)
{
    return *(int*)x0 - *(int*)x1;
}

int load_int(const char *str)
{
    errno = 0;
    char *e;
    long n = strtol(str,&e,10);
    
    // エラーチェック
    if (errno == ERANGE){
	fprintf(stderr, "%s: %s\n", str, strerror(errno));
	exit(1);
    }
    if( *e != '\0'){
	fprintf(stderr, "%s: an irregular character '%c' is detected.\n",str,*e);
	exit(EXIT_FAILURE);
    }
    return (int)n;
}
intlib.h 
// intlib.h
#ifndef QSORT_QSORT_H
#define QSORT_QSORT_H

int comp_int(const void *x0, const void *x1);

int load_int(const char *str);

#endif

インクルードガード

さきほどはヘッダファイルに特に情報を加えませんでした。先程の例の場合、このヘッダファイルが複数回呼ばれても「宣言」しか存在しないため影響しません（コンパイルが遅くなる可能性はあります）。

// 
#include <stdio.h>
#include "calc.h"
#include "calc.h"
#include "calc.h"

int main(){
    double a = 1.2;
    double b = 1.4;

    printf("a+b: %f\n",add(a,b));
    printf("a-b: %f\n",sub(a,b));
    printf("a*b: %f\n",mul(a,b));
    printf("a/b: %f\n",div(a,b));
    return 0;
}

上記のコードのうち、プリプロセッサが展開した状態のコードから"double add(" という部分を抜き出してみましょう。

# コンパイルせずにプリプロセッサだけ実行するには
# -E オプションをつける
gcc -E main.c | grep "double add("

プロトタイプ宣言が3回出てきているのがわかります。

しかし定義ファイルが含まれる場合には後述のように事情が異なるため、一般にはインクルードガードというテクニックで、ヘッダファイルが複数回呼ばれるのを抑制します。

// よくあるインクルードガードの構文
// calc.h 
#ifndef _CALC_H_
#define _CALC_H_
double add(double a, double b);
double sub(double a, double b);
double mul(double a, double b);
double div(double a, double b);

#endif

上記の例ではプリプロセッセマクロの_CALC_H_ が定義されていなければ、これを定義して#endifまでのブロックを読み込むという意味になります。一度このブロックに入ると必ずマクロが定義されるので、もう一度calc.h がincludeされても読み込まれません。

# コンパイルせずにプリプロセッサだけ実行するには
# -E オプションをつける
gcc -E main.c | grep "double add("

なお上記は非常に昔から書かれている方法ですが、最近では以下のような書き方も多くのコンパイラで使えます。

// calc.h
#pragma once
double add(double a, double b);
double sub(double a, double b);
double mul(double a, double b);
double div(double a, double b);

定義（特に構造体定義）を含む場合の分割コンパイル

インクルードガードが必要になる例として行列とベクトルの演算を考えてみましょう。例えばベクトルの掛け算だけであればベクトルの定義だけで十分ですが、行列とベクトルの掛け算などの場合、両方の定義が必要となります。簡単のためここでは2次元の演算のみを考えます（構造体ポインタを使わないようにハードコーディングされています）。

今

  • matrix.h : 行列の構造体定義とプロトタイプ宣言（行列とベクトルの積を含む）
  • vector.h : ベクトルの構造体定義とプロトタイプ宣言
  • matrix.c : 行列関係の関数実装
  • vector.c : ベクトル関係の関数実装

のような分割を考えます。

# tar で解凍し、ディレクトリ移動
tar xvzf vector_calc.tar.gz
cd vector_calc

// main.c
#include "matrix.h"
#include "vector.h"
		    
int main()
{
    Vector x = (Vector){ .v = { 1.0, 2.0 } };
    Vector y = (Vector){ .v = { 3.0, 4.0 } };
    Matrix A = (Matrix){ .v = { 1.0, -1.0, -1.0, 1.0} };
    
    // c = 2.0 x + y
    Vector c = axpy(2.0, x, y);
    // d = 2.0 Ax + y
    Vector d = gemv(2.0, A, x, 1.0, y);

    print_vector(c);
    print_vector(d);
    return 0;
}

// vector.c
#include <stdio.h>
#include "vector.h"

Vector axpy(double alpha, Vector a, Vector b)
{
    return (Vector){ .v = { alpha * a.v[0] + b.v[0], alpha * a.v[1] + b.v[1]}};
}

void print_vector(Vector a)
{
    printf("[%f %f]\n",a.v[0],a.v[1]);
}

// matrix.c
#include <stdio.h>
#include "matrix.h"
#include "vector.h"

Vector gemv(double alpha, Matrix A, Vector x, double beta, Vector y)
{
    return (Vector){ .v = {alpha * (A.v[0]*x.v[0]+A.v[1]*x.v[1]) + beta * y.v[0], alpha * (A.v[2]*x.v[0]+A.v[3]*x.v[1]) + beta * y.v[1]}};
}

void print_matrix(Matrix m)
{
    printf("[%f %f ;%f %f]\n",m.v[0],m.v[1],m.v[2],m.v[3]);    
}

// vector.h
typedef struct {
    double v[2];
} Vector;

// ax + y の形（BLASレベル1と同じ名前のつけかた）
Vector axpy(double alpha, Vector a, Vector b);
void print_vector(Vector a);

// matrix.h
#include "vector.h"

typedef struct{
    double v[4]; // 2 x 2
} Matrix;

// プロトタイプ宣言でVector型を知っている必要がある
// alpha * Ax + beta * y（BLASレベル2と同じ名前のつけかた）
Vector gemv(double alpha, Matrix A, Vector x, double beta, Vector y);
void print_matrix(Matrix m);

この状態でmain.c をコンパイルしてみましょう。

gcc -c main.c

コンパイラがvector.h が複数回呼ばれている旨エラーを出すと思います。これはmatrix.h でVector型の宣言が必要なので、vector.h をinclude しているために、結果的にVector型が複数回定義されてしまっています。

ヘッダファイルを呼び出す側としては、matrix.h を呼べばvector.h は呼ばれるので呼ぶ必要がないという考え方ではコードを書かないため、複数回呼ばれても大丈夫なように衝突回避する必要があります。

Let's try

matrix.h および vector.h にインクルードガードを施して全てのコードをコンパイルしてください。

さらにモジュールを隠蔽する: 宣言と定義の分離

これまで構造体の宣言と定義は一緒に行うことが多かったと思います。この場合、ヘッダファイルをインクルードすることで各プログラムは構造体の定義を知ることができ、そのメンバへのアクセス等が可能になっていました。しかしモジュールを細かく分割していった場合、究極的には対象となる構造体がどのようなメンバを持つかはわからないほうが望ましい です。そのようにすると構造体のメンバーをその構造体を使う別の関数からは見えなくなります。どうやって使うのかという話になりますが、メンバー変数にアクセスするための関数のプロトタイプをヘッダファイルに記載してあげます。

今回の例の場合、

  • matrix.h, vector.h ではtypedef により構造体の宣言のみ行う。
  • matrix.c, vector.c にて構造体の定義を行い、それぞれのファイル内ではメンバを直接参照可能。他の関数にはアクセスできる関数を提供する。
  • Matrix, Vector は基本的にポインタだけでやりとりする。

としてモジュールを隠蔽することができます。

Let's try

高難度ですが、適宜関数を追加した上で、構造体の定義部分をヘッダファイルから隠蔽してみてください。

Let's try 一例（宣言と定義の分離）

[講義後追記] 構造体の定義をmatrix.c, vector.c 内で行います。matrix.c 内でVectorの値を取り出したい場合は、定義の情報がないので、vector.h でアクセスのための関数を提供します。

main.c 
#include "matrix.h"
#include "vector.h"

int main(){
    // 複合リテラルでdouble型配列を用意して引数とする。
    // 配列が引数の場合は、関数内部ではポインタ扱いで処理される。
    Vector *x = init_vector((double []){1.0, 2.0});
    Vector *y = init_vector((double []){3.0, 4.0});
    Matrix *A = init_matrix((double []){1.0, -1.0, -1.0, 1.0});

    // c = 2.0 x + y
    Vector *c = axpy(2.0, x, y); 
    
    // d = 2.0 Ax + y
    Vector *d = gemv(2.0, A, x, 1.0, y);

    print_vector(c);
    print_vector(d);
    
    // 以下のfree関数群の結果、ポインタにはNULLが代入されるので、init関数でmallocされたときのアドレス情報は残らない
    free_vector(&x);
    free_vector(&y);

    free_vector(&c);
    free_vector(&d); 
    free_matrix(&A); 

    return 0;
}
vector.h 
#pragma once
#include <stddef.h> // size_t
// 宣言のみ行う
typedef struct vector Vector;

// ax + y の形（BLASレベル1と同じ名前のつけかた）
// ポインタでのやり取りに変更
Vector* axpy(double alpha, const Vector *a, const Vector *b);
void print_vector(const Vector *a);

// Vector型を動的に確保して初期化する関数
Vector *init_vector(const double v[]);
// Vector 用の free / Vector型ポインタへのポインタとすることで、nullリセット可能にする
void free_vector(Vector **v);
// 要素の取り出し
double get_vec_element(const Vector *v, size_t i);
vector.c 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "vector.h"

// ここで構造体の定義をかく。このメンバーに直接触れるのはこのファイルの中の関数だけ
struct vector {
    double v[2];
};

Vector* axpy(double alpha, const Vector *a, const Vector *b)
{
    Vector *v = malloc(sizeof(Vector));

    // このvector.c 内に限り、アロー演算子が有効（構造体の定義を知っているため）
    *v = (Vector){ .v = { alpha * a->v[0] + b->v[0], alpha * a->v[1] + b->v[1]} };
    
    // 以下は上記の置き換えパターンのうち、かなり裏技的な書き方。Vector型は配列のみのメンバなので、doubleのポインタにキャストして
    // バナナをつけると、doubleとして値を取り出せる。
    // memcpy(v,(double []){alpha *  *(double*)a + *(double*)b, alpha * *((double*)a + 1) + *((double*)b + 1)}, sizeof(double)*2);

    return v;
}

void print_vector(const Vector *a)
{
    printf("[%f %f]\n",a->v[0],a->v[1]);
}

Vector *init_vector(const double v[]){
    Vector *ret = malloc(sizeof(Vector));
    memcpy(ret->v, v, sizeof(double)*2);
    return ret;
}

void free_vector(Vector **v){
    free(*v);
    *v = NULL;
}

double get_vec_element(const Vector *v, size_t i){
    // 不正なインデックスアクセスを抑制
    if ( i != 0 && i != 1) {
        fprintf(stderr, "%s: index should be 0 or 1.",__func__);
        exit(1);
    }
    return v->v[i];
}
matrix.h 
#pragma once

#include <stddef.h> // size_t
#include "vector.h"
typedef struct matrix Matrix;

// プロトタイプ宣言ではVector型があることをを知っていればよい
// alpha * Ax + beta * y
Vector* gemv(double alpha, const Matrix *A, const Vector *x, double beta, const Vector *y);
void print_matrix(const Matrix *m);

// Matrix型を動的確保し、初期化する関数
Matrix *init_matrix(const double m[]);
// Matrix用の free / 引数をMatrix型ポインタへのポインタとすることで、nullリセット可能にする
void free_matrix(Matrix **m);
double get_mat_element(const Matrix *m, size_t i, size_t j);
matrix.c 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "matrix.h"
#include "vector.h"

struct matrix{
    double v[4];
};

Vector *gemv(double alpha, const Matrix *A, const Vector *x, double beta, const Vector *y)
{
    double xx, xy, yx, yy; // x = [xx, xy], y = [yx, yy];
    xx = get_vec_element(x, 0);
    xy = get_vec_element(x, 1);
    yx = get_vec_element(y, 0);
    yy = get_vec_element(y, 1);

    double axx, axy, ayx, ayy;
    axx = get_mat_element(A, 0, 0);
    axy = get_mat_element(A, 0, 1);
    ayx = get_mat_element(A, 1, 0);
    ayy = get_mat_element(A, 1, 1);

    double a = alpha * (axx * xx + axy * xy) + beta * yx; 
    double b = alpha * (ayx * xx + ayy * xy) + beta * yy;

    return init_vector((double []){a,b});
}

void print_matrix(const Matrix *m)
{
    printf("[%f %f ;%f %f]\n",m->v[0],m->v[1],m->v[2],m->v[3]);    
}

Matrix *init_matrix(const double m[]){
    Matrix *ret = (Matrix *)malloc(sizeof(Matrix));
    memcpy(ret->v, m, sizeof(double) * 4);
    return ret;
}

void free_matrix(Matrix **m){
    free(*m);
    *m = NULL;
}

double get_mat_element(const Matrix *m, size_t i, size_t j){
    if ( i != 0 && i != 1){
        fprintf(stderr, "%s: index out of range %zu. It should be [0 1].\n",__func__,i);
        exit(1);
    }
     if ( j != 0 && j != 1){
        fprintf(stderr, "%s: index out of range %zu. It should be [0 1].\n",__func__,j);
        exit(1);
    }
    return m->v[2 * i + j];
}

ライブラリ

これまで最初の方の課題や、数学関数を含むコンパイルなどで例えば以下のようなコンパイルをしていたかと思います。

# ライフゲーム用のプログラムのコンパイル
gcc -o lifegame -Wall life.c -L. -lgol
# 数学関数を使う場合（特にLinuxやWSL2）
gcc -o hoge hoge.c -lm

実はライブラリとは基本的複数のオブジェクトファイルを一つに固めたものです。関数の仕様はヘッダファイルに記述しておくことで、ソースコードの実装を意識せずに関数を使用することができます。ライブラリには基本的に以下の2種類があります

  • 静的ライブラリ: プログラム構築時に組み込まれる (lib*.a)

    • ライブラリ部分の再コンパイルは不要
    • 実行バイナリは大きくなりがち

  • 共有ライブラリ: プログラム実行時にメモリに展開 (lib*.{so,dylib})

    • ライブラリが更新されても本体プログラムの再コンパイルが不要
    • 実行バイナリは比較的小さい
    • ライブラリが置かれた情報を実行時に適切に知る必要がある

前者の静的ライブラリは全ての関数オブジェクトがライブラリに含まれており、コンパイル時にそれらをリンクします。一方で後者の共有ライブラリではコンパイル時には、ファイルレベルでどこにライブラリがあるかをコンパイル時に教えておき、実行時にそれを探します。今後みなさんが、何らかのオープンソースのプログラムをコンパイルした際に、最終的に実行時にerror while loading shared libraries のようなエラーと遭遇した場合は、後者のライブラリの置かれているとされるところにアクセスできていないことが原因です。

静的ライブラリの作成と使用

今日は静的ライブラリの作成を行ってみましょう。ライブラリの作成にはar というコマンドを用います。

ar rsv libinteger.a intlib.o

ここで

  • r: 既存ファイルの更新 or 新規作成
  • s: 複数のオブジェクトファイルをランダムアクセス可能にする
  • v: 詳細を表示

という意味です。このようにできたライブラリをコンパイル時に使用する場合は2通りの方法があります。ひとつはオブジェクトファイルと同様に直接指定する方法です。

gcc -o qsort main.c libinteger.a

もう一つはgcc において、-l というオプションを用いることでリンクするライブラリを指定する方法です。これはlib????.a のリンクを -l???? の形式で指定します。

gcc -o qsort main.c -L. -linteger

ソースコードのディレクトリを整頓する

分割コンパイルを行うことで、プログラム中の役割に応じて異なるファイルが生成されることになります。これを単一のディレクトリで扱うと非常にごちゃごちゃします。オープンソースのプロジェクトなどでもそうですが、ヘッダファイル、ソースコード、ライブラリ、実行コードを階層的に配置しておくことで、保守性を高めることができます。

  • src: ソースコードを格納する
  • include: ヘッダファイルを格納（-I./include で参照）
  • lib: ライブラリファイルを格納（-L./lib で参照）
  • bin: 実行バイナリを保存するディレクトリ

分割コンパイル時、gcc -c はデフォルトでは実行したディレクトリにオブジェクトファイルを出力します。-o は -c と組み合わせても用いることが可能で、この場合はオブジェクトファイル名を指定できます。そのため例えば、src以下にオブジェクトを生成したい場合は

gcc -c -I./include -o src/intlib.o src/intlib.c

のようになります。あくまでオブジェクトの生成なので、-c を忘れないのが重要です。

最終的にはライブラリをlib に作成した上で、以下のコードでコンパイルします。

gcc -I./include -o bin/qsort src/main.c -L./lib -linteger
Let's try

qsort.h, main.c, intlib.c をqsort.tar.gz を展開したディレクトリに適切に配置し、ライブラリlibinteger.a および 実行ファイルqsort を生成してみてください。

Let's try 一例

[講義後追記] 以下にLet's tryでやるべき手順を示します。 まず、qsort.h 、main.c およびintlib.c を以下のように配置します。

.
├── bin
├── include
│   └── qsort.h
├── lib
└── src
    ├── main.c
    └── intlib.c

次に、intlib.c からオブジェクトファイルを生成し、続いてこれをライブラリにします。

# オブジェクトファイルの生成
gcc -I./include -c -o src/intlib.o src/intlib.c 
# ar でライブラリにする
ar rsv lib/libinteger.a src/intlib.o

最後にこのライブラリとmain.c からバイナリをbin以下に生成します。

# コンパイル
gcc -o bin/qsort -I./include src/main.c -L./lib -linteger 
# 実行
bin/qsort 10

Makefile

Makefileとは

これまでは単一のプログラムを一回コンパイルするのみなので、逐次コンパイルコマンドを叩けばOKでした。しかしプログラムの規模が大きくなり、複数のソースファイルでコードが構成されるようになると、手順が複雑になりマニュアルでコンパイルを叩いていくのは現実的でなくなります。それを解消するのがMakefileの使用です。Makefileはプログラムをコンパイルするのに必要な情報を記述したファイルで、どのファイルがどのファイルで使われるかなどの依存関係を適切に記述することで、適切な順序でコンパイルを行ってくれます。Makefileが存在するディレクトリでmakeコマンドを実行すると、用意されたMakefileを参照してコンパイルを実行してくれます。

target: prerequisites
	recipe

ここで

  • target: 生成したいファイル
  • prerequisites: target を生成するのに必要なファイル
  • recipe: 実行するコマンド

を表します。注意すべき点として個々のrecipeの前にはタブが必要なことです。

ルールの実際

以下はルールの例です。

main.o: main.c
	gcc -c main.c

上記はmain.oを作成するためには、main.c が必要で、作成のためのコマンドはgcc –c main.cという意味を示しています。

実際に実行する場合は以下のようなコマンドになります。

make target

target は実際のターゲットの名前によって変わります。これによってルールが実行されtargetが作成されます。なお、ルールの実行の際にはprerequisitesのタイムスタンプがチェックされ、更新されている場合にのみ recipe が実行されます。すなわちプログラムの再コンパイルの必要がない場合にはコンパイルされず、ファイルが修正された場合にのみルールが実行されます。

Makefileの例

vector_calcにおけるMakefileの例を示します。

calc: main.o vector.o matrix.o
	gcc -o calc main.o vector.o matrix.o

main.o: main.c
	gcc -c main.c

vector.o: vector.c
	gcc -c vector.c

matrix.o: matrix.c
	gcc -c matrix.c

このMakefileが存在するディレクトリ で

make

とするとmain.o, vector.o, matrix.o, calc が順々に作成されます。なおtargetを省略したmakeは最初のターゲットがしてされたものとみなされます。

Phony target（偽ターゲット）

実はmakeで指定するtargetは必ずしもファイルでなくともよいです。以下のようなMakefileを作成すると、make clean と実行することでa.outやオブジェクトファイル等をまとめて消せます。

clean:
	rm *.o *~ a.out

例えば全部コンパイルしなおしたい場合等には

make clean
make

といった実行を行います。なお、この際にもしclean というファイルが存在すると意図した挙動をしないため、ファイルでないコマンド的なターゲットであることを明示します。

.PHONY: clean
clean:
	rm *.o *~ a.out

なお*~ はEmacsなどのエディタが生成するバックアップファイル（ファイル名に~がつく）を表しています。

変数の定義

Makefileでは変数が定義可能です。例えばコンパイラとしてgcc以外を使いたい等、環境に応じて変数への代入を変更することで、異なる挙動を実現できます。

CC = gcc

calc: main.o vector.o matrix.o
	$(CC) main.o vector.o matrix.o

main.o: main.c
	$(CC) -c main.c
vector.o: vector.c
	$(CC) -c vector.c
matrix.o: main.c
	$(CC) -c matrix.c

さらに各recipesで使える変数として以下があります。これを用いることで、ファイル名のハードコーディングを回避できます。

  • $@ : ターゲット
  • $< : prerequisitesの最初
  • $^: 全てのprerequisites

暗黙ルール

実はMakefileはいくつかの言語の基本的なコンパイルパターンを暗黙のルールを持っていて、ターゲットに指定されていないファイルがprerequisites に含まれる場合にも自動的に暗黙ルールによるコンパイルを行ってくれます。例えばfoo.o は foo.c から

$(CC) $(CPPFLAGS) $(CFLAGS) -c

で作成されるというルールは暗黙ルールに含まれています。そのため、先ほどのMakefileはさらに短縮されます。

CC = gcc

calc: main.o vector.o matrix.o
	$(CC) main.o vector.o matrix.o

最終的なMakefile

最終的にcalc をターゲットにする場合は、結局以下のMakefileだけで分割コンパイルが可能です。

CC = gcc
TARGET = calc
OBJS = main.o vector.o matrix.o

$(TARGET): $(OBJS)
	$(CC) -o $@ $^
Let's try

以下をやってみましょう。

  • vector_calc をmakeでコンパイルできるようにMakefileを書いてみましょう。
  • ディレクトリを構成し、整頓したqsortにおいて、Makefileを書いてみましょう。

連続関数の最適化

今日は関数ポインタを使いつつ、一般的な連続関数の最適化を行ってみましょう。ここでは入力がベクトル、出力がスカラーのような関数を考え、この関数の最小値およびそれを実現する入力を計算することを考えます。これは現在の機械学習において、「目的関数の最大化/最小化」と呼ばれる最も基本的なプロセスです。入力ベクトルを変化させたときの目的関数の変動（勾配: gradient）がわかっている場合には、勾配法と呼ばれる方法を用いることができます。

1変数関数の場合

たとえば\( f(x) = (x-3)^2 + 1 \) を最小化する\(x\) を求めたい場合、\(f'(x) = 2(x-3)\) を考えます。この場合は人間は\(f'(x) = 0 \) を簡単に求めることができます（と思います）が、これが求められない場合は

  • 適当に \(x\) を決める
  • \(f'(x) < 0\) なら右に移動
  • \(f'(x) > 0 \) なら左に移動

として繰り返します。

2変数関数の場合

例えば$$f(x,y) = (x-3)^4 + (y-2)^2 $$ のような関数を最適化したい場合勾配として

$$\nabla f = \left(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}\right)=[4(x-3)^3, 2(y-2)]$$

を考えます。この勾配の反対方向に\( (x,y) \) を進めることで最適化を実現します。

最急降下法 (gradient descent)

以上を一般的にまとめると、最急降下法のプロセスは以下のようになります。

  1. 適当に初期位置\( x \)を決める
  2. \(\nabla x \)を計算。ノルムが十分小さければ終了
  3. \(x \leftarrow x - \alpha \nabla x\)として更新　　　　　　　
  4. 2に戻り繰り返す

サンプルプログラム

サンプルプログラムを実行してみましょう。

# 展開
tar xvzf optimize.tar.gz
cd optimize
# コンパイル&実行
make
./optimizer 0.05

Makefile

今回のサンプルにはMakefileが同梱されています。

CC = gcc
CFLAGS = -Wall
LDLIBS = -lm
OBJS = main.o func.o optimize.o
TARGET = optimizer

$(TARGET): $(OBJS)
	$(CC) -o $@ $^ $(LDLIBS)

.PHONY: tmpclean clean
tmpclean:
	rm -f *~
clean: tmpclean
	rm -f $(OBJS) $(TARGET)
Let's try

以下に取り組みましょう。

  • optimize.c を修正し、最適化の終了条件を厳しくしてみましょう

    • make したときに optimize.c しかコンパイルされないことを確認

  • 最適化の各ステップで関数の値 f(x) も表示するように改良してみましょう

    • 関数 optimize() に関数 f_value() のポインタを渡すようにする

Let's try 一例

[講義後追記] 終了条件については optimize.c にある breakを呼び出しているif文を修正すればよいです。一方関数の値も表示する場合は以下のようになります。

// optimize.h
#pragma once

int optimize(const double alpha, const int dim, double x[], 
	     void (*calc_grad)(const double [], double []),
	     double (*calc_value)(const double[]));

// optimize.c
// ヘッダと calc_norm()は省略

int optimize(const double alpha, const int dim, double x[], 
             void (*calc_grad)(const double [], double []),
	     double (*calc_value)(const double[]))
{
    // 勾配ベクトルを記録する領域を確保
    double *g = malloc(dim * sizeof(double));
    
    int iter = 0;
    while (++iter < 10000) {
	
	// 引数で渡された関数を使って勾配ベクトルを計算
	(*calc_grad)(x, g);
	
	// 勾配ベクトルのノルムを評価
	const double norm = calc_norm(dim, g);
	printf("%3d norm = %7.4f", iter, norm);
	for (int i = 0; i < dim; i++) {
	    printf(", x[%d] = %7.4f", i, x[i]);
	}
	printf(", f(x) = %7.4f\n",(*calc_value)(x));
	
	if (norm < 0.01) break;
	
	// 最急降下法による更新
	for (int i = 0; i < dim; i++) {
	    x[i] -= alpha * g[i];
	}
    }
    
    free(g);
    
    return iter;
}

// main.c
// 以下のみ修正
const int iter = optimize(alpha, dim, x, f_gradient,f_value);

応用

最小二乗基準による線形回帰を勾配法を用いて行ってみましょう。

最小二乗法について

最小二乗法による線形回帰はデータセット$${(x_1,y_1),\ldots,(x_n,y_n)}$$ が与えられたときに、二乗誤差の和$$E=\sum_{i=1}^N (y_i - ax_i - b)^2$$ が最小になるように$$a,b$$ を設定する問題です。

大まかな手順

  1. a, b を適当な値で初期化
  2. 与えられた(x, y) と a, bを用いてパラメータa, b に対する勾配を求める
  3. a, b の値を上記勾配に基づいて最急降下法で更新
  4. 平均二乗誤差が十分小さければ終了
  5. そうでない場合は2 - 4 を繰り返す

実際、今回は1次関数の線形回帰のため、実際には決定的に定まりますが、上記の手順は昨今流行のニューラルネットの最適化でも基本的な考えとなります。

課題（締切: 1/10）

  1. 前回、配布された線形リスト版のペイントソフト（こちら）について、複数のヘッダファイル、ソースファイルに分割し、Makefile でコンパイルできるように修正せよ。
    • コンパイルディレクトリを整頓し、make コマンドのみでコンパイルできるようにすること
    • 可能であれば構造体宣言と定義の分離、ポインタによるメンバの隠蔽、ライブラリ化を行って、配列版と線形リスト版のいずれの実装でも同じmain関数でコンパイル可能にせよ。
      • このときMakefile に異なるコンパイルパターンのターゲットを書くとよい。
    • ディレクトリ名はpaint、main関数を含むファイル名はpaintapp.c とする。適宜分割したソースファイルを用いる。
  2. 最小二乗法により標高と気温の関係を推定せよ
    • 得られたパラメータから推測される富士山頂の７月の平均気温を示すこと
    • func1.cの構造体定義 を参考にdata.csv のデータを読みとって利用する
    • あわせて与えられたデータを標高に基づいてソートし、表示せよ
    • プログラムおよびコンパイルするディレクトリを同梱すること（方法は後述）
    • ディレクトリ名はmtfuji、 main関数を含むファイル名は mtfuji.cとする。適宜分割したソースファイルを用いる。
  3. ［発展課題］より高度な回帰プログラムを実装し評価せよ
    • モデルの改良の例
      • 重回帰分析、非線形回帰、カーネル回帰、リッジ回帰、ニューラルネットワーク etc.
    • 最適化の手法を改良しても良い
    • 可能であれば自身で分析対象を探し変更するとよい
      • 回帰に際してはパラメータを推定するデータと、実際に評価するデータを適切に分割すること。交差検証またはクロスバリデーションというキーワードも合わせて調べてみると良い。
    • プログラム、データおよびコンパイルするディレクトリを添付すること
    • ディレクトリ名はadv、main関数を含むファイル名はadv_regression.c とする。適宜分割したソースファイルを用いる。

ディレクトリ構成およびMakefile

  • 基本課題のみの場合はSOFT-12-21-NNNNNNNN 以下にpaintおよびmtfuji というディレクトリ、発展課題も提出の場合は paint、mtfuji および adv という3つのディレクトリを作成する
  • 今回の課題ではコンパイル用のMakefileを作成し同梱せよ
  • それぞれのディレクトリに作成したMakefileを設置する
  • make clean も実装すること
  • 提出前にmake clean しておく（= バイナリは同梱しない）
  • 追加のデータファイルがある場合はそれも同梱せよ

提出方法

  • ITC-LMS にて提出する
  • 全てのプログラム/ファイルをまとめ、zip またはtar.gzで圧縮し提出
  • 必ずREADME.md を同梱し、やった内容を簡潔に記述する
  • ファイル名はSOFT-12-21-NNNNNNNN.zip または SOFT-12-21-NNNNNNNN.tar.gz とする
    • NNNNNNNN部分は学籍番号（ハイフンは除く）
    • JについてはJ??????? のようにする
  • 課題について
  • 基本課題は毎回提出
  • 発展課題は成績計算に全6回中上位3回分を採用する