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

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

新規登録して質問してみよう
ただいま回答率
85.50%
C

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

Q&A

解決済

3回答

5055閲覧

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

Takabatake_Y

総合スコア13

C

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

0グッド

1クリップ

投稿2015/11/10 09:27

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

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

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

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

###ソースコード

C

1#include <stdio.h> 2#include <math.h> 3#include <stdlib.h> 4#include <sys/time.h> 5#include <sys/resource.h> 6 7#define BASE 2 8 9void input(double **, double **, int n); 10void output(double **, int n); 11void str(double **, double **, double **, int n); 12void prod(double **, double **, double **, int n); 13 14double e_time(void); 15double t1, t2, t3 = 0, elapsed; 16 17int main(int argc, char *argv[]) 18{ 19 int i, j, n, size; 20 double **A, **B, **C; 21 22 if (argc <= 1) { 23 fprintf(stderr, "error: argument which indicates the problem size is required.\n"); 24 exit(1); 25 } 26 n = atoi(argv[1]); 27 size = n + 1; 28 29 /* Get the memory space for arrays A, B, C. */ 30 A = (double **)malloc(size * sizeof(double *)); 31 B = (double **)malloc(size * sizeof(double *)); 32 C = (double **)malloc(size * sizeof(double *)); 33 for (i = 1; i <= size; i++) { 34 *(A + i) = (double *)malloc(size * sizeof(double)); 35 *(B + i) = (double *)malloc(size * sizeof(double)); 36 *(C + i) = (double *)malloc(size * sizeof(double)); 37 } 38 39 40 input(A, B, n); 41 //reset vectors 42 for(i = 1; i < size; i++) { 43 for(j = 1; j < size; j++) { 44 C[i][j] = 0.0; 45 } 46 } 47#ifdef DEBUG 48 printf("A = \n"); 49 output(A, n); 50 printf("B = \n"); 51 output(B, n); 52#endif 53 54 t1 = e_time(); 55 str(A, B, C, n); 56 prod(A, B, C, n); 57 t2 = e_time(); 58 //elapsed = (t2 -t1) + t3; 59 elapsed = t2 -t1; 60 printf("t1:%f, t2:%f, elapsed:%f\n", t1, t2, elapsed); 61 printf("Performance : %f MFLOPS\n", (2.0*size*size*size)/elapsed/1000000.0); 62 63#ifdef DEBUG 64 printf("C = \n"); 65 output(C, n); 66#endif 67 68 for (i = 1; i <= size; i++) { 69 free(*(A+i)); 70 free(*(B+i)); 71 free(*(C+i)); 72 } 73 free(A); 74 free(B); 75 free(C); 76 77 return 0; 78} 79 80void input(double **A, double **B, int n) 81{ 82 int i, j; 83 84 for (i = 1; i <= n; i++) { 85 for (j = 1; j <= n; j++) { 86 A[i][j] = (double)rand()/RAND_MAX; 87 B[i][j] = (double)rand()/RAND_MAX; 88 } 89 } 90} 91 92void output(double **O, int n) 93{ 94 int i, j; 95 96 for (i = 1; i <= n; i++) { 97 for (j = 1; j <= n; j++) 98 printf("%3.6f ", O[i][j]); 99 printf("\n"); 100 } 101} 102 103void str(double **A, double **B, double **C, int n) 104{ 105 int i, j, mid, size; 106 double **X, **Y, **M1, **M2, **M3, **M4, **M5, **M6, **M7; 107 double t4, t5; 108 109 if (n <= BASE) { 110 t4 = e_time(); 111 prod(A, B, C, n); 112 t5 = e_time(); 113 t3 += (t5 - t4); 114 return; 115 } 116 mid = n / 2; 117 118 /* Get the memory space for arrays X, Y, M1--M7 here. */ 119 size = n + 1; 120 X = (double **)malloc(size * sizeof(double *)); 121 Y = (double **)malloc(size * sizeof(double *)); 122 M1 = (double **)malloc(size * sizeof(double *)); 123 M2 = (double **)malloc(size * sizeof(double *)); 124 M3 = (double **)malloc(size * sizeof(double *)); 125 M4 = (double **)malloc(size * sizeof(double *)); 126 M5 = (double **)malloc(size * sizeof(double *)); 127 M6 = (double **)malloc(size * sizeof(double *)); 128 M7 = (double **)malloc(size * sizeof(double *)); 129 for (i = 1; i <= mid; i++) { 130 *(X+i) = (double *)malloc(size*sizeof(double)); 131 *(Y+i) = (double *)malloc(size*sizeof(double)); 132 *(M1+i) = (double *)malloc(size*sizeof(double)); 133 *(M2+i) = (double *)malloc(size*sizeof(double)); 134 *(M3+i) = (double *)malloc(size*sizeof(double)); 135 *(M4+i) = (double *)malloc(size*sizeof(double)); 136 *(M5+i) = (double *)malloc(size*sizeof(double)); 137 *(M6+i) = (double *)malloc(size*sizeof(double)); 138 *(M7+i) = (double *)malloc(size*sizeof(double)); 139 } 140 141 t4 = e_time(); 142 /* Compute M1 */ 143 for (i = 1; i <= mid; i++) 144 for (j = 1; j <= mid; j++) { 145 X[i][j] = A[i][mid+j] - A[mid+i][mid+j]; 146 Y[i][j] = B[mid+i][j] + B[mid+i][mid+j]; 147 } 148 str(X,Y,M1,mid); 149 150 151 /* Compute M2--M7 */ 152 // M2 153 for (i = 1; i <= mid; i++) { 154 for (j = 1; j <= mid; j++) { 155 X[i][j] = A[i][j] + A[mid+i][mid+j]; 156 Y[i][j] = B[i][j] + B[mid+i][mid+j]; 157 } 158 } 159 str(X,Y,M2,mid); 160 161 // M3 162 for (i = 1; i <= mid; i++) { 163 for (j = 1; j <= mid; j++) { 164 X[i][j] = A[i][j] - A[mid+i][j]; 165 Y[i][j] = B[i][j] + B[i][mid+j]; 166 } 167 } 168 str(X,Y,M3,mid); 169 170 // M4 171 for (i = 1; i <= mid; i++) { 172 for (j = 1; j <= mid; j++) { 173 X[i][j] = A[i][j] + A[i][mid+j]; 174 Y[i][j] = B[mid+i][mid+j]; 175 } 176 } 177 str(X,Y,M4,mid); 178 179 // M5 180 for (i = 1; i <= mid; i++) { 181 for (j = 1; j <= mid; j++) { 182 X[i][j] = A[i][j]; 183 Y[i][j] = B[i][mid+j] - B[mid+i][mid+j]; 184 } 185 } 186 str(X,Y,M5,mid); 187 188 // M6 189 for (i = 1; i <= mid; i++) { 190 for (j = 1; j <= mid; j++) { 191 X[i][j] = A[mid+i][mid+j]; 192 Y[i][j] = B[mid+i][j] - B[i][j]; 193 } 194 } 195 str(X,Y,M6,mid); 196 197 // M7 198 for (i = 1; i <= mid; i++) { 199 for (j = 1; j <= mid; j++) { 200 X[i][j] = A[mid+i][j] + A[mid+i][mid+j]; 201 Y[i][j] = B[i][j]; 202 } 203 } 204 str(X,Y,M7,mid); 205 206 207 for (i = 1; i <= mid; i++) { 208 for (j = 1; j <= mid; j++) { 209 C[i][j] = M1[i][j]+M2[i][j]-M4[i][j]+M6[i][j]; /* C11 */ 210 C[i][mid+j] = M4[i][j]+M5[i][j]; /* C12 */ 211 C[mid+i][j] = M6[i][j]+M7[i][j]; /* C21 */ 212 C[mid+i][mid+j] = M2[i][j]-M3[i][j]+M5[i][j]-M7[i][j]; /* C22 */ 213 } 214 } 215 216 217 218 /* Free the memory space for X,Y,M1--M7 here */ 219 for (i = 1; i <= mid; i++) { 220 free(*(X+i)); free(*(Y+i)); free(*(M1+i)); 221 free(*(M2+i)); free(*(M3+i)); free(*(M4+i)); 222 free(*(M5+i)); free(*(M6+i)); free(*(M7+i)); 223 } 224 t5 = e_time(); 225 t3 += t5 - t4; 226 free(X); free(Y); free(M1); 227 free(M2); free(M3); free(M4); 228 free(M5); free(M6); free(M7); 229 230 return; 231} 232 233/* The conventional O(n^3) matrix multiplication algorithm 234 which is used to compute the product of two 2 by 2 matrices. */ 235void prod(double **A, double **B, double **C, int n) 236{ 237 int i, j, k; 238 239 for (i = 1; i <= n; i++) 240 for (j = 1; j <= n; j++) { 241 C[i][j] = 0; 242 for (k = 1; k <= n; k++) { 243 C[i][j] = C[i][j] + A[i][k] * B[k][j]; 244 245 } 246 } 247} 248 249 250 /*************************************************************/ 251 /* functions related to time measurement */ 252 /*************************************************************/ 253 double e_time(void) { 254 255 static struct timeval now; 256 static struct timezone tz; 257 258 gettimeofday(&now, &tz); 259 return (double)(now.tv_sec + now.tv_usec/1000000.0); 260 }

気になる質問をクリップする

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

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

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答3

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関数に入る前のところに以下のコードを追加。

C

1struct Arrays 2{ 3 double **X, **Y, **M1, **M2, **M3, **M4, **M5, **M6, **M7; 4}; 5 6Arrays *arrays; 7int arrays_size; 8int index = 0;

main関数のA,B,Cをmallocしているところの後ろあたりに以下のコードを追加。

C

1 arrays_size = 20; 2 arrays = (Arrays *)malloc(sizeof(Arrays) * arrays_size); 3 for(int ai = 0; ai < arrays_size; ai++) 4 { 5 Arrays *ptr = &arrays[ai]; 6 ptr->X = (double **)malloc(size * sizeof(double *)); 7 ptr->Y = (double **)malloc(size * sizeof(double *)); 8 ptr->M1 = (double **)malloc(size * sizeof(double *)); 9 ptr->M2 = (double **)malloc(size * sizeof(double *)); 10 ptr->M3 = (double **)malloc(size * sizeof(double *)); 11 ptr->M4 = (double **)malloc(size * sizeof(double *)); 12 ptr->M5 = (double **)malloc(size * sizeof(double *)); 13 ptr->M6 = (double **)malloc(size * sizeof(double *)); 14 ptr->M7 = (double **)malloc(size * sizeof(double *)); 15 for(i = 1; i < size; i++) 16 { 17 ptr->X[i] = (double *)malloc(size*sizeof(double)); 18 ptr->Y[i] = (double *)malloc(size*sizeof(double)); 19 ptr->M1[i] = (double *)malloc(size*sizeof(double)); 20 ptr->M2[i] = (double *)malloc(size*sizeof(double)); 21 ptr->M3[i] = (double *)malloc(size*sizeof(double)); 22 ptr->M4[i] = (double *)malloc(size*sizeof(double)); 23 ptr->M5[i] = (double *)malloc(size*sizeof(double)); 24 ptr->M6[i] = (double *)malloc(size*sizeof(double)); 25 ptr->M7[i] = (double *)malloc(size*sizeof(double)); 26 } 27 }

arrays_size = 20;は本来であればちゃんと計算すべきでしょうけど、手抜きしました。nが100万くらいまではいけると思います。
あと、main関数から出るところでfreeしてあげてください。

str関数のif (n <= BASE)の処理の後ろ(mid = n / 2;の後ろあたり)に以下のコードを追加。

C

1 Arrays *ptr = &arrays[index]; 2 index++; // これ重要 3 X = ptr->X; 4 Y = ptr->Y; 5 M1 = ptr->M1; 6 M2 = ptr->M2; 7 M3 = ptr->M3; 8 M4 = ptr->M4; 9 M5 = ptr->M5; 10 M6 = ptr->M6; 11 M7 = ptr->M7;

str関数の最後

C

1 index--; // これ重要 2 return; 3}

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

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

投稿2015/11/11 03:20

編集2015/11/11 05:06
catsforepaw

総合スコア5938

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

Takabatake_Y

2015/11/11 05:47

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

0

malloc/freeがボトルネックになっているのであれば、アロケータを自前で実装しては如何でしょうか?シングルスレッドかつ確保/解放がstack的(sbrk的)用法で事足りるので、相応の処理性能向上が期待できると思います。サンプルコードを添付します。

  • 既存コードのmallocをmyMallocに、freeをmyFreeに置換
  • myFreeの呼び出しをmyMallocの逆順に変更
  • 初期化/終了処理としてinitializeMyHeapとfinalizeMyHeapの呼び出しを追加
  • initializeMyHeapのiPoolSizeは最大確保サイズ[B], iMaxBlockは最大確保ブロック数です.

lang

1typedef struct _MyHeap { 2 void* pPool; 3 int iPoolSize; 4 int* pBlockSize; 5 int iMaxBlock; 6 void* pCurrent; 7 int iCurrentBlock; 8} MyHeap; 9 10static MyHeap* pThis = NULL; 11 12/////////////////////////////////////////////////////////////////////////////// 13void initializeMyHeap(int iPoolSize, int iMaxBlock) { 14 assert(pThis == NULL); 15 assert(iPoolSize > 0); 16 assert(iMaxBlock > 0); 17 18 pThis = malloc(sizeof(MyHeap)); 19 pThis->iPoolSize = iPoolSize; 20 pThis->iMaxBlock = iMaxBlock; 21 pThis->pPool = malloc(iPoolSize); 22 pThis->pBlockSize = malloc(sizeof(int) * iMaxBlock); 23 pThis->pCurrent = pThis->pPool; 24 pThis->iCurrentBlock = 0; 25} 26 27/////////////////////////////////////////////////////////////////////////////// 28void finalizeMyHeap() { 29 assert(pThis != NULL); 30 31 free(pThis->pPool); 32 free(pThis->pBlockSize); 33 free(pThis); 34} 35 36/////////////////////////////////////////////////////////////////////////////// 37void* myMalloc(int size) { 38 void* rv; 39 40 assert(size > 0); 41 assert((char*)pThis->pCurrent - (char*)pThis->pPool + size < pThis->iPoolSize); 42 assert(pThis->iCurrentBlock < pThis->iMaxBlock); 43 44 rv = pThis->pCurrent; 45 (char*)pThis->pCurrent += size; 46 pThis->pBlockSize[pThis->iCurrentBlock] = size; 47 ++pThis->iCurrentBlock; 48 49 return rv; 50} 51 52/////////////////////////////////////////////////////////////////////////////// 53void myFree(void* p) { 54 assert(p != NULL); 55 assert(pThis->iCurrentBlock > 0); 56 assert((char*)pThis->pCurrent - pThis->pBlockSize[pThis->iCurrentBlock-1] == p); 57 58 --pThis->iCurrentBlock; 59 (char*)pThis->pCurrent -= pThis->pBlockSize[pThis->iCurrentBlock]; 60} 61

投稿2015/11/11 01:32

編集2015/11/11 01:43
退会済みユーザー

退会済みユーザー

総合スコア0

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

Takabatake_Y

2015/11/11 05:17

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

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/10 09:58

編集2015/11/10 13:49
Chironian

総合スコア23272

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

Takabatake_Y

2015/11/11 05:20 編集

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問