2

我需要求解一个二维泊松方程，即 for AX=B 中的方程组，其中 A 是 n×n 矩阵，B 是 n×1 向量。作为 2D Poisson 问题的离散化矩阵，我知道只有 5 个对角线不会为空。Lapack 不提供解决此特定问题的函数，但它具有求解带状矩阵方程组的函数，即 DGBTRF（用于 LU 因式分解）和 DGBTRS。现在，5条对角线是：主对角线，主对角线上方和下方的第一条对角线以及主对角线上方和下方的两条对角线，m对角线wrt。在阅读了有关带存储的 lapack 文档后，我了解到我必须创建一个 (3*m+1)-by-n 矩阵来以带存储格式存储 A，我们将这个矩阵称为 AB。现在的问题：

1) dgbtrs 和 dgbtrs_ 有什么区别？英特尔 MKL 两者都提供，但我不明白为什么

2) dgbtrf 要求带存储矩阵是一个数组。我应该按行还是按列线性化 AB？

3）这是调用这两个函数的正确方法吗？

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 administrativa。只需调用dgbtrs（原因是在某些架构上，fortran 例程名称附加了下划线，而在其他情况下则没有，并且名称可能是大写或小写 - 英特尔 MKL#define是正确的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。lambda 是一个很小的数字，以避免 k=0 的分歧。5 对角方程是泊松方程的带限逼近，它通过有限差分逼近微分算子。

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

于 2013-08-28T19:56:21.200 回答