2010-11-18 12 views
1

Mathematicaから非常に複雑な行列(〜1.3MBのプレーンテキスト)をFortranプログラムで使用するために出力しようとしています。私がこれを行うと(Spliceを介して)、変数に数値が与えられると、結果の行列は〜2%だけオフになります。これは、正確にゼロである固有値が必要であり、固有ベクトルの構成が正確である必要があるため、問題である。Mathematica CForm/FortranFormのバグ?

精度、正しい変数、適切な対角化コードなどの通常のデューデリジェンスをすべて実行しました。そのような大きな行列やMathematicaがFortranFormの出力を乱してしまったことが起こりました。

私はMathematicaに行列のCFormを与え、それを試しました。 FortranFormマトリックスと同じ(マシンの精度内)でした。

誰かがこの種の問題に遭遇しましたか?何が原因かもしれないか考えていますか?これを理解するためにMathematicaでフォーマットされたFortranコードの25000行を調べる必要があります。

EDIT:問題のマトリックスは複雑で、大きくはありません。それはわずか6x6ですが、各要素は三角関数、対数、さまざまなルーツとパワーを含めて、代数的に非常に乱雑です。

本発明のマトリックスの(1,1)元素のPlaintext,C codeおよびFortran codeである。正当なパラメータ値は、0 <ラムダ、カッパ、Y *** < 1; 100と1000

+0

小さな行列でエラーを再現できますか? –

+0

エクスポートする象徴的な行列であることを正しく理解していますか?そして、Fortran/Cで評価し、数値結果をMathematicaにインポートすると、Mathematicaの記号式を評価するのと比べて、2%のエントリーが削除されますか? – Janus

+0

@HighPerformanceMark:いいえ、はるかに単純な行列を使ってテストしました。それらはすべて機械精度で再現されています。 – Timo

答えて

1

最初の要素の式を見ると、これは有限精度浮動小数点演算のためにほぼ確実に問題になります。私にはこの問題に取り組む2つの方法があるようです。 1つ目は、Cコードの精度を上げることです。多分、マクロプリプロセッサを使って(あるいは単に見つけて置き換えて)、Cコード内のすべての宣言をdoubleからlong doubleに変更することができます。使用するコンパイラに応じて、64ビットから80または128になります(長いdoubleを標準doubleとして扱うMicrosoft's compilerを使用しない限り)。あなたが使用しているCコンパイラを投稿すると、そのルートに行くためにどのような変更が必要かを調べるのに役立ちます。最後にMathematicaやdoubleだけを使っている別のCプログラムに渡す必要がある場合は、doubleにキャストすることができます。

また、Cコンパイラは異なる浮動小数点モードを持つことができます。モードが正確か厳格か高速かを確認する価値があります。高速または高速計算の場合は、丸めの問題に寄与している可能性があります。

もう1つの方法は、浮動小数点精度を最適化するためにMathematicaを使用して各要素の条件を並べ替えることです。一般的に(a * b)* c =(c * b)* a浮動小数点演算の土地でこれは当てはまりません。私が聞いたことが最善のヒューリスティックは、次のとおりです。

  1. 追加:あなたは維持するように、値が互いに非常に接近している差し引いた値を回避または遅延:最大に
  2. 減算を丸め誤差を最小限にするために、最小から用語を追加します可能な限り多くの有効桁数
  3. 桁数/除算:非常に小さい桁数で掛けたり、有効桁数の理由で非常に大きな数値で除算したりしないでください。

実際には、これはこれらの式を因数分解することがあなたの最善の策であり、特に問題のある用語です。 a1〜= a2の場合、(a1-a2)(b-d)式は(a1-a2)* b-(a1-a2)* dよりも精度/有効桁数が大きくなります。

+0

しかし、痛みでしたが、表現をリファクタリングするのに役立ちました。 – Timo

1

それは以下のような簡単なものになるだろうとの間に他のすべて:かなり異なる出力で例えばMathematicaは数値評価について非常に賢いので、これら二つの式

N[(10^18+100)-10^18] 
N[10^18+100]-N[10^18] 

resul

Out[3]= 100. 
Out[4]= 128. 

CまたはFortranではベビーシッターを受けず、2行目と同じように評価されます。
Mathematicaで1つの行列エントリを見て、C/Fortranと同じように評価して、それがあなたの見た目と一致するかどうかを調べることはできますか?

EDIT:
中間の精度に問題のための簡単なテストが

In[10]:= Compile[{a,b},a+100-b][10^18,10^18] 
Out[10]= 128. 

はあなたの行列のコンパイル済みのバージョンがと一致するかどうかを確認することができるように、すべての精密検査を低下するCompileを使用しているようですC/Fortran?

+0

私たちの行列の要素はせいぜいO(10^6)であり、個々の変数はO(1)-O(1000)なので、単精度で十分ですし、確かに倍精度でもFortranで使用します。また、Mathematicaの精度を人為的に4桁に減らし、依然として元の結果の0.01%以内に収まっています。 – Timo

+0

さて、N [expr、4]で10^-4以下の誤差が出るのはボックスが言うことです:正しく理解すれば、N [、k]は結果をk桁の精度に正確に保障します。私は確かに、より多くの作業精度でもN [、4]より精度が低くなる状況を考えることができます。 – Janus

+0

@Timo:Mathematicaを愚かなものにするためにあらゆる種類の巧妙なコードを試した後、私はコンパイルが中間精度が問題かどうかをテストする簡単な方法であることに気付きました。上記の編集を参照してください。 – Janus