求解带状矩阵方程组答案

【问题标题】：Solve a banded matrix system of equations求解带状矩阵方程组
【发布时间】：2012-06-06 09:29:45
【问题描述】：

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

【问题讨论】：

标签： c system lapack equations poisson

【解决方案1】：

它还提供DGBTRS 和DGBTRS_。这些是您不应该关心的fortran administrativa。只需调用dgbtrs（原因是在某些体系结构上，fortran 例程名称附加了下划线，而在其他情况下则没有，名称可能是大写或小写——英特尔 MKL #defines 是 dgbtrs 的正确名称）。
LAPACK 例程需要列主矩阵（即 Fortran 样式）：一个接一个地存储列。您必须使用的带状存储并不难：http://www.netlib.org/lapack/lug/node124.html。
这对我来说似乎不错，但请事先在小问题上尝试一下（顺便说一句，这总是一个好主意）。还要确保你处理非零info（这是报告错误的方式）。

更好的风格是使用MKL_INT 而不是普通的int，这是正确类型的typedef（在某些架构上可能不同）。

还要确保在调用dgbtrf 之前为pivots 分配内存。

【讨论】：

关于#2，这不会弄乱LU分解吗？我的意思是，没有选项可以告诉 dgbtrf AB 是行专业的。另外，我应该如何设置这两种情况的主要尺寸？
是的，您将解决 A'x = B。让我更正我对此的回答。对于前导维度，它们通常等于行数，除非您正在做一些时髦的事情（例如将子矩阵作为参数传递）。同样，解决 LAPACK 问题的最佳方法是尝试可以手工解决的小问题。

【解决方案2】：

这可能是题外话。但是对于泊松方程，基于 FFT 的解决方案要快得多。只需对势场进行 2D FFT，除以 -(k^2+lambda^2) 然后进行 IFFT。 lambda 是一个很小的数字，以避免 k=0 的分歧。 5对角方程是泊松方程的带限逼近，通过有限差分逼近微分算子。

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

【讨论】：