2

2D ポアソン方程式、つまり AX=B の方程式系を解く必要があります。ここで、A は n 行 n 列の行列で、B は n 行 1 列のベクトルです。A は 2D ポアソン問題の離散化行列であるため、5 つの対角線のみが null にならないことがわかっています。Lapack には、この特定の問題を解決する関数はありませんが、帯行列連立方程式を解く関数、すなわち DGBTRF (LU 因数分解用) と DGBTRS があります。ここで、5 つの対角線は、主対角線、主対角線の上下の最初の対角線、および主対角線に対して m 対角線分だけ上下にある 2 つの対角線です。バンド ストレージに関する lapack ドキュメントを読んだ後、A をバンド ストレージ形式で格納するには (3*m+1) 行 n 列の行列を作成する必要があることを知りました。この行列を AB と呼びましょう。今質問：

1) dgbtrs と dgbtrs_ の違いは何ですか? インテル® MKL は両方を提供しますが、その理由がわかりません

2) dgbtrf では、バンド ストレージ マトリックスが配列である必要があります。行または列で AB を線形化する必要がありますか?

3) これは 2 つの関数を呼び出す正しい方法ですか?

int n, m;
double *AB;
/*... fill n, m, AB, with appropriate numbers */

int *pivots;
int nrows = 3 * m + 1, info, rhs = 1;
dgbtrf_(&n, &n, &m, &m, AB, &nrows, pivots, &info);
char trans = 'N';
dgbtrs_(&trans, &n, &m, &m, &rhs, AB, &nrows, pivots, B, &n, &info);
4

2 に答える 2

0

  1. また、 と も提供DGBTRSしDGBTRS_ます。これらは、気にする必要のない fortran 管理者です。呼び出すだけdgbtrsです (一部のアーキテクチャでは、fortran ルーチン名にアンダースコアが追加されている場合とそうでない場合があり、名前が大文字または小文字#defineのいずれかである場合があります。インテル® MKL はdgbtrs.

  2. LAPACK ルーチンは、列の主要な行列 (つまり、Fortran スタイル) を想定しています。列を次々に格納します。使用しなければならないバンド ストレージは難しくありません: http://www.netlib.org/lapack/lug/node124.html。

  3. 私には良いように思えますが、事前に小さな問題で試してみてください (ちなみに、常に良い考えです)。また、必ずゼロ以外を処理してくださいinfo(これがエラーの報告方法です)。

より良いスタイルはMKL_INT、plain の代わりに使用することですint。これは、正しい型への typedef です (一部のアーキテクチャでは異なる場合があります)。

pivotsまた、を呼び出す前に、必ずメモリを割り当ててdgbtrfください。

于 2012-06-06T19:45:52.930 に答える
0

これは論外かもしれません。しかし、ポアソン方程式の場合は、FFT ベースのソリューションの方がはるかに高速です。潜在的なフィールドの 2D FFT を実行し、-(k^2+lambda^2) で割ってから IFFT を実行します。ラムダは、k=0 の発散を避けるための小さな数です。5 対角方程式は、微分演算子を有限差分で近似する、ポアソン方程式の帯域制限近似です。

http://en.wikipedia.org/wiki/Screened_Poisson_equation

于 2013-08-28T19:56:21.200 に答える