使用LAPACK的DGBSV从C中解决一个带状系统问题

Question

我尝试搜索有关此主题的类似帖子，但似乎没有人关心带状系统（...）。

我有兴趣使用LAPACK / ScaLAPACK从C代码解决带状矩阵。首先，我想使用LAPACK实现顺序解决方案，然后再尝试使用ScaLAPACK。

问题：两种语言之间的行主/列主差异似乎影响了我的解决过程。这是我打算解决的系统：

以下代码将矩阵转换为LAPACK的带状数据结构，该结构在此处中指定。

  int rr = 6;     // Rank.
  int kl = 2;     // Number of lower diagonals.
  int ku = 1;     // Number of upper diagonals.
  int nrhs = 1;   // Number of RHS.
  double vals[36] = {0.0,   0.0,   0.0,   0.0,   0.0,   0.0,  // Req. ex. space.
                     0.0,   0.0,   0.0,   0.0,   0.0,   0.0,  // Req. ex. space.
                   666.0,   0.0,   0.0,   0.0,   0.0,  22.5,  // First up diag.
                     1.0, -50.0, -50.0, -50.0, -50.0,  -2.6,  // Main diagonal.
                    27.5,  27.5,  27.5,  27.5,   4.0, 666.0,  // First low diag.
                     0.0,   0.0,   0.0,  -1.0, 666.0, 666.0}; // 2nd low diag.

  int lda = rr;   // Leading dimension of the matrix.
  int ipiv[6];    // Information on pivoting array.
  double rhs[] = {1.0, 1.0, 1.0, 1.0, 1.0, 0.0};  // RHS.
  int ldb = lda;  // Leading dimension of the RHS.
  int info = 0;   // Evaluation variable for solution process.
  int ii;         // Iterator.
  int jj;         // Iterator.

  dgbsv_(&rr, &kl, &ku, &nrhs, vals, &lda, ipiv, rhs, &ldb, &info);

  printf("info = %d\n", info);
  for (ii = 0; ii < ldb; ii++) {
    printf("%f\n", rhs[ii]);
  }
  putchar('\n');

如我所说，我担心我的矩阵翻译方式不正确，因为它是按列主序排列的，而且还涉及到Fortran的索引方式，因为我的解决方案产生了以下结果:

[ejspeiro@node01 lapack-ex02]$ make runs
`pwd`/blogs < blogs.in
info = 1
1.000000
1.000000
1.000000
1.000000
1.000000
0.000000

从Fortran的返回值info = 1可以看出，分解已完成，但在LU分解中A = LU中的U(1,1) = 0。

- Eduardo

---

有没有办法让我把这个标记为已解决？我实际上得到了一个可运行的代码！ :D - Eduardo

Answer 1

正如您所注意到的那样，您的矩阵是以行主格式输入的。将其写成列主格式后，我们所看到的行将对应于列：

  double vals[36] = {0.0,   0.0,   0.0,   1.0,  27.5,   0.0,
                     0.0,   0.0,   0.0, -50.0,  27.5,   0.0,
                     0.0,   0.0,  22.5, -50.0,  27.5,   0.0,
                     0.0,   0.0,  22.5, -50.0,  27.5,  -1.0,
                     0.0,   0.0,  22.5, -50.0,   4.0,  0.0,
                     0.0,   0.0,  22.5,  -2.6,   0.0,  0.0};

- bossylobster

Answer 2

好的，我会按顺序回答这个问题，以便将其称为“已解决”。

这些文件代表着功能代码。我将其提供出来，以防有人遇到相同的问题：使用C中的LAPACK解决带状方程组。

1. https://dl.dropbox.com/u/5432016/banded/cmtk.h
2. https://dl.dropbox.com/u/5432016/banded/blogs.c

如果有人对实现方面有额外的建议，我将非常欢迎！

我的下一步是分发核心矩阵，以便使用ScaLAPACK进行求解。

谢谢！

- Eduardo

---

链接已失效。请勿链接将会失效的外部资源。 - Vladimir F Героям слава