質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

90.40%

  • C

    4131questions

    C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

2次元配列によるシュトラッセンアルゴリズムの実装について

解決済

回答 3

投稿

  • 評価
  • クリップ 1
  • VIEW 1,168

Takabatake_Y

score 7

こちらの質問をご覧頂きありがとうございます。

前提

行列積計算の速度向上を目的に、シュトラッセンアルゴリズムの実装を行っています。
現在のプログラムでは、シュトラッセンの関数(str)に入るたびに2重のポインタと逐次mallocを使うことで二次元配列を実現しています。

質問

mallocのメモリ確保・開放分の時間を削減するために、このプログラムを、2重ポインタから固定サイズの2次元配列で書き換え、実現することは可能でしょうか。

ソースコードの方に現在の2重ポインタを使用したプログラムを添付します。
ご回答の程よろしくお願いします。

ソースコード

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <sys/time.h>
#include <sys/resource.h>

#define BASE 2

void input(double **, double **, int n);
void output(double **, int n);
void str(double **, double **, double **, int n);
void prod(double **, double **, double **, int n);

double e_time(void);
double t1, t2, t3 = 0, elapsed;

int main(int argc, char *argv[])
{
    int i, j, n, size;
    double **A, **B, **C;

    if (argc <= 1) {
        fprintf(stderr, "error: argument which indicates the problem size is required.\n");
        exit(1);
    }
    n = atoi(argv[1]);
    size = n + 1;

    /* Get the memory space for arrays A, B, C. */
    A = (double **)malloc(size * sizeof(double *));
    B = (double **)malloc(size * sizeof(double *));
    C = (double **)malloc(size * sizeof(double *));
    for (i = 1; i <= size; i++) {
        *(A + i) = (double *)malloc(size * sizeof(double));
        *(B + i) = (double *)malloc(size * sizeof(double));
        *(C + i) = (double *)malloc(size * sizeof(double));
    }


    input(A, B, n);
    //reset vectors
    for(i = 1; i < size; i++) {
        for(j = 1; j < size; j++) {
            C[i][j] = 0.0;
        }
    }
#ifdef DEBUG
    printf("A = \n");
    output(A, n);
    printf("B = \n");
    output(B, n);
#endif

    t1 = e_time();
    str(A, B, C, n);
    prod(A, B, C, n);
    t2 = e_time();
    //elapsed = (t2 -t1) + t3;
    elapsed  = t2 -t1;
    printf("t1:%f, t2:%f, elapsed:%f\n", t1, t2, elapsed);
    printf("Performance : %f MFLOPS\n", (2.0*size*size*size)/elapsed/1000000.0);

#ifdef DEBUG
    printf("C = \n");
    output(C, n);
#endif

    for (i = 1; i <= size; i++) {
        free(*(A+i));
        free(*(B+i));
        free(*(C+i));
    }
    free(A);
    free(B);
    free(C);

    return 0;
}

void input(double **A, double **B, int n)
{
    int i, j;

    for (i = 1; i <= n; i++) {
        for (j = 1; j <= n; j++) {
            A[i][j] = (double)rand()/RAND_MAX;
            B[i][j] = (double)rand()/RAND_MAX;
        }
    }
}

void output(double **O, int n)
{
    int i, j;

    for (i = 1; i <= n; i++) {
        for (j = 1; j <= n; j++) 
            printf("%3.6f ", O[i][j]);
        printf("\n");
    }
}

void str(double **A, double **B, double **C, int n)
{
    int i, j, mid, size;
    double **X, **Y, **M1, **M2, **M3, **M4, **M5, **M6, **M7;
    double t4, t5;

    if (n <= BASE) {
        t4 = e_time();
        prod(A, B, C, n);
        t5 = e_time();
        t3 += (t5 - t4);
        return;
    }
    mid = n / 2;

    /* Get the memory space for arrays X, Y, M1--M7 here. */
    size = n + 1;
    X  = (double **)malloc(size * sizeof(double *));
    Y  = (double **)malloc(size * sizeof(double *));
    M1 = (double **)malloc(size * sizeof(double *));
    M2 = (double **)malloc(size * sizeof(double *));
    M3 = (double **)malloc(size * sizeof(double *));
    M4 = (double **)malloc(size * sizeof(double *));
    M5 = (double **)malloc(size * sizeof(double *));
    M6 = (double **)malloc(size * sizeof(double *));
    M7 = (double **)malloc(size * sizeof(double *));
    for (i = 1; i <= mid; i++) {
        *(X+i)  = (double *)malloc(size*sizeof(double));
        *(Y+i)  = (double *)malloc(size*sizeof(double));
        *(M1+i) = (double *)malloc(size*sizeof(double));
        *(M2+i) = (double *)malloc(size*sizeof(double));
        *(M3+i) = (double *)malloc(size*sizeof(double));
        *(M4+i) = (double *)malloc(size*sizeof(double));
        *(M5+i) = (double *)malloc(size*sizeof(double));
        *(M6+i) = (double *)malloc(size*sizeof(double));
        *(M7+i) = (double *)malloc(size*sizeof(double));
    }

    t4 = e_time();
    /* Compute M1 */
    for (i = 1; i <= mid; i++)
        for (j = 1; j <= mid; j++) {
            X[i][j] = A[i][mid+j] - A[mid+i][mid+j];
            Y[i][j] = B[mid+i][j] + B[mid+i][mid+j];
        }
    str(X,Y,M1,mid);


    /* Compute M2--M7 */
    // M2
    for (i = 1; i <= mid; i++) {
        for (j = 1; j <= mid; j++) {
            X[i][j] = A[i][j] + A[mid+i][mid+j];
            Y[i][j] = B[i][j] + B[mid+i][mid+j];
        }
    }    
    str(X,Y,M2,mid);

    // M3
    for (i = 1; i <= mid; i++) {
        for (j = 1; j <= mid; j++) {
            X[i][j] = A[i][j] - A[mid+i][j];
            Y[i][j] = B[i][j] + B[i][mid+j];
        }
    }    
    str(X,Y,M3,mid);

    // M4
    for (i = 1; i <= mid; i++) {
        for (j = 1; j <= mid; j++) {
            X[i][j] = A[i][j] + A[i][mid+j];
            Y[i][j] = B[mid+i][mid+j];
        }
    }    
    str(X,Y,M4,mid);

    // M5
    for (i = 1; i <= mid; i++) {
        for (j = 1; j <= mid; j++) {
            X[i][j] = A[i][j];
            Y[i][j] = B[i][mid+j] - B[mid+i][mid+j];
        }
    }    
    str(X,Y,M5,mid);

    // M6
    for (i = 1; i <= mid; i++) {
        for (j = 1; j <= mid; j++) {
            X[i][j] = A[mid+i][mid+j];
            Y[i][j] = B[mid+i][j] - B[i][j];
        }
    }    
    str(X,Y,M6,mid);

    // M7
    for (i = 1; i <= mid; i++) {
        for (j = 1; j <= mid; j++) {
            X[i][j] = A[mid+i][j] + A[mid+i][mid+j];
            Y[i][j] = B[i][j];
        }
    }    
    str(X,Y,M7,mid);


    for (i = 1; i <= mid; i++) {
        for (j = 1; j <= mid; j++) {
            C[i][j] = M1[i][j]+M2[i][j]-M4[i][j]+M6[i][j]; /* C11 */
            C[i][mid+j] = M4[i][j]+M5[i][j]; /* C12 */
            C[mid+i][j] = M6[i][j]+M7[i][j]; /* C21 */
            C[mid+i][mid+j] = M2[i][j]-M3[i][j]+M5[i][j]-M7[i][j]; /* C22 */
        }
    }



    /* Free the memory space for X,Y,M1--M7 here */
    for (i = 1; i <= mid; i++) {
        free(*(X+i));   free(*(Y+i));  free(*(M1+i));
        free(*(M2+i));  free(*(M3+i)); free(*(M4+i));
        free(*(M5+i));  free(*(M6+i)); free(*(M7+i));
    } 
    t5 = e_time();
    t3 += t5 - t4;
    free(X);  free(Y);  free(M1);
    free(M2); free(M3); free(M4);
    free(M5); free(M6); free(M7);

    return;
}

/* The conventional O(n^3) matrix multiplication algorithm
   which is used to compute the product of two 2 by 2 matrices. */
void prod(double **A, double **B, double **C, int n)
{
    int i, j, k;

    for (i = 1; i <= n; i++) 
        for (j = 1; j <= n; j++) {
            C[i][j] = 0;
            for (k = 1; k <= n; k++) {
                C[i][j] = C[i][j] + A[i][k] * B[k][j]; 

            }
        }
}


    /*************************************************************/
    /* functions related to time measurement                     */
    /*************************************************************/
    double e_time(void) {

        static struct timeval now;
        static struct timezone tz;

        gettimeofday(&now, &tz);
        return (double)(now.tv_sec  + now.tv_usec/1000000.0);
    }
  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 3

checkベストアンサー

0

メモリの確保/解放を高速化する手法としては、最初にまとめて確保して領域を使い回す方法(メモリプール)がありますが、まじめに実装しようとすると手間がかかります。
実のところ、malloc/freeの中身はメモリプールの仕組みそのものなので、それなりに高速です。
ただ、汎用故のオーバーヘッドもあるので、とにかく高速に処理したいというのなら、自前で工夫するしかないですね。

ということで、ちょっとコードを追ってみたところ、どうやら最初に引数で渡したnの値で再起呼び出しの深さも確定するようです。これなら、あらかじめ再帰呼び出し分も含めて全てのメモリ領域を一括で確保し、領域を使い回すということができそうです。2次元配列をさらに配列で管理して、そのインデックスをスタックポインタのようにstr関数の入り口で上げて出口で下げるイメージです。

あと、コードに気になる点が一つ。
main関数内のメモリを確保/解放するところで
for (i = 1; i <= size; i++)
とやっていますが、正しくは
for (i = 1; i < size; i++)
ですね。

追記です。
2重ポインタから固定サイズの2次元配列で書き換え、実現することは可能でしょうか。 
サイズがコマンドライン引数で渡されているので、固定サイズというのは無理な気がします。上記説明は、str関数の中でmalloc/freeが繰り返し実行されないようにするための一つの方法です。

なんだか面白そうだったので、実際に試してしまいました。
main関数に入る前のところに以下のコードを追加。
struct Arrays
{
    double **X, **Y, **M1, **M2, **M3, **M4, **M5, **M6, **M7;
};

Arrays    *arrays;
int        arrays_size;
int        index = 0;

main関数のA,B,Cをmallocしているところの後ろあたりに以下のコードを追加。
    arrays_size = 20;
    arrays = (Arrays *)malloc(sizeof(Arrays) * arrays_size);
    for(int ai = 0; ai < arrays_size; ai++)
    {
        Arrays *ptr = &arrays[ai];
        ptr->X = (double **)malloc(size * sizeof(double *));
        ptr->Y = (double **)malloc(size * sizeof(double *));
        ptr->M1 = (double **)malloc(size * sizeof(double *));
        ptr->M2 = (double **)malloc(size * sizeof(double *));
        ptr->M3 = (double **)malloc(size * sizeof(double *));
        ptr->M4 = (double **)malloc(size * sizeof(double *));
        ptr->M5 = (double **)malloc(size * sizeof(double *));
        ptr->M6 = (double **)malloc(size * sizeof(double *));
        ptr->M7 = (double **)malloc(size * sizeof(double *));
        for(i = 1; i < size; i++)
        {
            ptr->X[i] = (double *)malloc(size*sizeof(double));
            ptr->Y[i] = (double *)malloc(size*sizeof(double));
            ptr->M1[i] = (double *)malloc(size*sizeof(double));
            ptr->M2[i] = (double *)malloc(size*sizeof(double));
            ptr->M3[i] = (double *)malloc(size*sizeof(double));
            ptr->M4[i] = (double *)malloc(size*sizeof(double));
            ptr->M5[i] = (double *)malloc(size*sizeof(double));
            ptr->M6[i] = (double *)malloc(size*sizeof(double));
            ptr->M7[i] = (double *)malloc(size*sizeof(double));
        }
    }
arrays_size = 20;は本来であればちゃんと計算すべきでしょうけど、手抜きしました。nが100万くらいまではいけると思います。
あと、main関数から出るところでfreeしてあげてください。

str関数のif (n <= BASE)の処理の後ろ(mid = n / 2;の後ろあたり)に以下のコードを追加。
    Arrays *ptr = &arrays[index];
    index++;        // これ重要
    X = ptr->X;
    Y = ptr->Y;
    M1 = ptr->M1;
    M2 = ptr->M2;
    M3 = ptr->M3;
    M4 = ptr->M4;
    M5 = ptr->M5;
    M6 = ptr->M6;
    M7 = ptr->M7;

str関数の最後
    index--;    // これ重要
    return;
}

str関数内のmallocとfreeの部分は削除してしまってください。

手元の環境でコマンドライン引数を500で試してみましたが、修正前は約2秒で、修正後は約0.5秒でした。効果はあるようです。


投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2015/11/11 14:08 編集

    保留

    キャンセル

  • 2015/11/11 14:47


    ”サイズがコマンドライン引数で渡されているので、固定サイズというのは無理な気がします。” - この時点で気づくべきですよね…
    また、構造体で纏めてそれ自体の配列で管理する等、プログラム初心者の自分にとって大変参考になります。
    最終的にこの方法で実装して結果を最速にすることが出来ましたので、ベストアンサーとさせていただきます。
    ご回答、追記ありがとうございました。

    キャンセル

0

こんにちは。

argv[1]が決まれば、str()の最後のパラメータnの最大値も決まると思います。
malloc/freeの回数を減らすだけであれば、その最大値分のメモリをmain()の頭で現状と同じ構造で獲得すればX[i][j]形式のままアクセスできますよ。
ちょっと安易ですが、X, Y, M1,...をグローバル変数とすれば改造量は最小限で済みます。

他に、あまりお薦めではないですが、str()の最後のパラメータのnの最大値を決めれば通常の2次元配列を使ってできると思います。
ただ、記述にちょっと自信がないのと、たぶんお望みの方法ではないと思うので詳細は省略します。


【追記】
見落としてました。str()を再帰呼び出ししてますね。私の回答ではうまくいきません。
2次元配列は動的獲得できなかったと思います。
一次元配列にしてmalloc/free回数を減らすくらいしか思いつきません。
お役にたてず申し訳ない。


【更に追記】
良いサイトがありました。
> C言語で2次元配列を動的に割り当てる4つの方法

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2015/11/11 14:04 編集

    ご掲載頂いたサイトの配列自体へのポインタでのやり方が参考になりました、最終的には固定サイズの配列ではなく、こちらのような方法で記述していきます。
    ご回答、追記の方ありがとうございした。

    キャンセル

0

malloc/freeがボトルネックになっているのであれば、アロケータを自前で実装しては如何でしょうか?シングルスレッドかつ確保/解放がstack的(sbrk的)用法で事足りるので、相応の処理性能向上が期待できると思います。サンプルコードを添付します。
 - 既存コードのmallocをmyMallocに、freeをmyFreeに置換
 - myFreeの呼び出しをmyMallocの逆順に変更
 - 初期化/終了処理としてinitializeMyHeapとfinalizeMyHeapの呼び出しを追加
 - initializeMyHeapのiPoolSizeは最大確保サイズ[B], iMaxBlockは最大確保ブロック数です.
 
typedef struct _MyHeap {
    void* pPool;
    int iPoolSize;
    int* pBlockSize;
    int iMaxBlock;
    void* pCurrent;
    int iCurrentBlock;
} MyHeap;

static MyHeap* pThis = NULL;

///////////////////////////////////////////////////////////////////////////////
void initializeMyHeap(int iPoolSize, int iMaxBlock) {
    assert(pThis == NULL);
    assert(iPoolSize > 0);
    assert(iMaxBlock > 0);

    pThis = malloc(sizeof(MyHeap));
    pThis->iPoolSize = iPoolSize;
    pThis->iMaxBlock = iMaxBlock;
    pThis->pPool = malloc(iPoolSize);
    pThis->pBlockSize = malloc(sizeof(int) * iMaxBlock);
    pThis->pCurrent = pThis->pPool;
    pThis->iCurrentBlock = 0;
}

///////////////////////////////////////////////////////////////////////////////
void finalizeMyHeap() {
    assert(pThis != NULL);

    free(pThis->pPool);
    free(pThis->pBlockSize);
    free(pThis);
}

///////////////////////////////////////////////////////////////////////////////
void* myMalloc(int size) {
    void* rv;

    assert(size > 0);
    assert((char*)pThis->pCurrent - (char*)pThis->pPool + size < pThis->iPoolSize);
    assert(pThis->iCurrentBlock < pThis->iMaxBlock);

    rv = pThis->pCurrent;
    (char*)pThis->pCurrent += size;
    pThis->pBlockSize[pThis->iCurrentBlock] = size;
    ++pThis->iCurrentBlock;

    return rv;
}

///////////////////////////////////////////////////////////////////////////////
void myFree(void* p) {
    assert(p != NULL);
    assert(pThis->iCurrentBlock > 0);
    assert((char*)pThis->pCurrent - pThis->pBlockSize[pThis->iCurrentBlock-1] == p);

    --pThis->iCurrentBlock;
    (char*)pThis->pCurrent -= pThis->pBlockSize[pThis->iCurrentBlock];
}

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2015/11/11 14:17

    メモリアロケータ自体の改善は全く考えになかったので参考になりました、
    mallocを使う別の方法と併用して改善していきます。
    ご回答、コード記述有り難うございます。

    キャンセル

同じタグがついた質問を見る

  • C

    4131questions

    C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。