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

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

ただいまの
回答率

89.99%

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

解決済

回答 1

投稿

  • 評価
  • クリップ 0
  • VIEW 939

jackie687456

score 14

前提・実現したいこと

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

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

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

該当のソースコード

#include<stdio.h>
#include<math.h>
#include<omp.h>
#define DIV 100
#define MAXITER 10
#define ZERO 0.0
#define ONE 1.0
#define TWO 2.0
#define THREE 3.0
#define EPS 1.0E-15
double fxy(double x,double y){
  return (x*x+y*y-ONE)*(x*x+y*y-ONE)*(x*x+y*y-ONE)-x*x*y*y*y;
}
double dfx(double x,double y){
  return THREE*(x*x+y*y-ONE)*(x*x+y*y-ONE)*(TWO*x)-TWO*y*y*y*x;
}
double dfy(double x,double y){
  return THREE*(x*x+y*y-ONE)*(x*x+y*y-ONE)*(TWO*y)-THREE*x*x*y*y;
}
double newtonx(double x,double y){
  int i;
  double x2;
  for(i=0;i<MAXITER;i++){
    x2=x-fxy(x,y)/dfx(x,y);
    if(fabs(x2-x)<=EPS){
      return x2;
    }
    if(fabs(x2)>TWO){
      return x2;
    }
    x=x2;
  }
  return -THREE;
}
double newtony(double x,double y){
  int i;
  double y2;
  for(i=0;i<MAXITER;i++){
    y2=y-fxy(x,y)/dfy(x,y);
    if(fabs(y2-y)<=EPS){
      return y2;
    }
    if(fabs(y2)>TWO){
      return y2;
    }
    y=y2;
  }
  return -THREE;
}
int main(){
  int i,j;
  double unit,x,y;
  FILE *fp;
  char filename[10];
  unit=TWO/DIV;
  double t1,t2;

  t1=omp_get_wtime();

#pragma omp parallel private(x,y,i,j,filename)
  {
    sprintf(filename,"ans_%d",omp_get_thread_num());
    fp=fopen(filename,"w");

#pragma omp for
    for(i=-DIV;i<=DIV;i++){
      for(j=-DIV;j<=DIV;j++){
    x=newtonx(unit*i,unit*j);
    if(fabs(x)<=TWO)
      fprintf(fp,"%f %f\n",x,unit*j);
    y=newtony(unit*i,unit*j);
    if(fabs(y)<=TWO)
      fprintf(fp,"%f %f\n",unit*i,y);
      }
    }

    fclose(fp);
  }
  t2=omp_get_wtime();
  printf("TIME=%10.5g\n",t2-t1);

  return 0;
}
  • 気になる質問をクリップする

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

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

    クリップを取り消します

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

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

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

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

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

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

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

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

    質問の評価を下げる

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

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

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

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

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

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

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

    詳細な説明はこちら

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

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

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

質問への追記・修正、ベストアンサー選択の依頼

  • LouiS0616

    2017/06/28 15:42 編集

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

    キャンセル

回答 1

checkベストアンサー

+1

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

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

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

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

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

こうは書けませんか?

#pragma omp parallel private(x,y,i,j,filename)
{
  FILE *fp;
  sprintf(filename,"ans_%d",omp_get_thread_num());
  fp=fopen(filename,"w");
  // 中略
  fclose(fp);
}

最初に言ったように、OpenMPはいじったことがありませんので、
免罪符とは言いませんが、もしとんちんかんなことを言っていても大目に見ていただければと思います...

投稿

編集

  • 回答の評価を上げる

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

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

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

  • 回答の評価を下げる

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

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

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

  • 2017/06/30 23:10

    返信遅くなってすみません。わざわざご回答ありがとうございます。もう一度考え直してみます。もしこのように直してみては?という意見があれば確信がなくてもかまわないのでご助言いただけると嬉しいです。

    キャンセル

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

  • ただいまの回答率 89.99%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

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