ニュートン法による陰関数描画の並列計算

###前提・実現したいこと
陰関数描画をニュートン法を用いて行い、それを並列計算させるという課題に取り組んでいます。実際にプログラムを書いて実行してみたところ、描画をするための"ans_0"～"ans_3"までは作られて図形をプロットすることはできました。しかし、実行時間が表示されず実行後に次のようなメッセージが出てしまいます。
一体このメッセージは何を表していて、どう改善すべきなのでしょうか。
よろしくお願いいたします。
なお、この問題で用いている陰関数はハートを書くもので
(x^2+y^2-1)^3-(x^2)*(y^2)=0
です。

###発生している問題・エラーメッセージ

*** glibc detected *** ./pfigure: double free or corruption (out): 0x00007f7a3c0008c0 ***
中止

###該当のソースコード

C言語
1#include<stdio.h>
2#include<math.h>
3#include<omp.h>
4#define DIV 100
5#define MAXITER 10
6#define ZERO 0.0
7#define ONE 1.0
8#define TWO 2.0
9#define THREE 3.0
10#define EPS 1.0E-15
11double fxy(double x,double y){
12  return (x*x+y*y-ONE)*(x*x+y*y-ONE)*(x*x+y*y-ONE)-x*x*y*y*y;
13}
14double dfx(double x,double y){
15  return THREE*(x*x+y*y-ONE)*(x*x+y*y-ONE)*(TWO*x)-TWO*y*y*y*x;
16}
17double dfy(double x,double y){
18  return THREE*(x*x+y*y-ONE)*(x*x+y*y-ONE)*(TWO*y)-THREE*x*x*y*y;
19}
20double newtonx(double x,double y){
21  int i;
22  double x2;
23  for(i=0;i<MAXITER;i++){
24    x2=x-fxy(x,y)/dfx(x,y);
25    if(fabs(x2-x)<=EPS){
26      return x2;
27    }
28    if(fabs(x2)>TWO){
29      return x2;
30    }
31    x=x2;
32  }
33  return -THREE;
34}
35double newtony(double x,double y){
36  int i;
37  double y2;
38  for(i=0;i<MAXITER;i++){
39    y2=y-fxy(x,y)/dfy(x,y);
40    if(fabs(y2-y)<=EPS){
41      return y2;
42    }
43    if(fabs(y2)>TWO){
44      return y2;
45    }
46    y=y2;
47  }
48  return -THREE;
49}
50int main(){
51  int i,j;
52  double unit,x,y;
53  FILE *fp;
54  char filename[10];
55  unit=TWO/DIV;
56  double t1,t2;
57
58  t1=omp_get_wtime();
59
60#pragma omp parallel private(x,y,i,j,filename)
61  {
62    sprintf(filename,"ans_%d",omp_get_thread_num());
63    fp=fopen(filename,"w");
64
65#pragma omp for
66    for(i=-DIV;i<=DIV;i++){
67      for(j=-DIV;j<=DIV;j++){
68	x=newtonx(unit*i,unit*j);
69	if(fabs(x)<=TWO)
70	  fprintf(fp,"%f %f\n",x,unit*j);
71	y=newtony(unit*i,unit*j);
72	if(fabs(y)<=TWO)
73	  fprintf(fp,"%f %f\n",unit*i,y);
74      }
75    }
76
77    fclose(fp);
78  }
79  t2=omp_get_wtime();
80  printf("TIME=%10.5g\n",t2-t1);
81
82  return 0;
83}
84

LouiS0616

2017/06/28 07:26 編集

とりあえず、#define ZERO 0.0 みたいのはやめた方がいいです。ひどいコードの一例として非常に有名です。

行動規範の内容に同意します

回答1件

ベストアンサー

すみません、私はOpenMPは使ったことがないのですが...
誰も反応しないようですので、推測でものを言わせていただきます。

エラーメッセージを見るに、メモリの二重解放が発生しているように思えます。
コード内で明示的にいじっているのはファイルポインタしかないので、試しに
fclose(fp)の一行をコメントアウトしてみてください。

それで動きますようなら、私の推測がおそらく正しいことになります。
個人的には、次の部分があやしいような気がします。

C
1sprintf(filename,"ans_%d",omp_get_thread_num());
2fp=fopen(filename,"w");

スレッドのID毎にfilenameを変更しているのに、同じファイルポインタに入れてよいのでしょうか？

こうは書けませんか？

C
1#pragma omp parallel private(x,y,i,j,filename)
2{
3  FILE *fp;
4  sprintf(filename,"ans_%d",omp_get_thread_num());
5  fp=fopen(filename,"w");
6  // 中略
7  fclose(fp);
8}