###前提・実現したいこと
陰関数描画をニュートン法を用いて行い、それを並列計算させるという課題に取り組んでいます。実際にプログラムを書いて実行してみたところ、描画をするための"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
とりあえず、#define ZERO 0.0 みたいのはやめた方がいいです。ひどいコードの一例として非常に有名です。
回答1件
あなたの回答
tips
プレビュー