4

このトピックに関連する同様のスレッドを検索しようとしましたが、誰もバンド システムを気にしていないようです (...)。

C コードから LAPACK/ScaLAPACK を使用してバンド行列を解くことに興味があります。まず、ScaLAPACK で何かを試みる前に、LAPACK で逐次的なソリューションを実現したいと考えています。

問題: 両方の言語の行優先/列優先の違いが、ソリューション プロセスに影響しているようです。私が解決しようとしているシステムは次のとおりです。

線形システム

次のコードは、そのマトリックスを、ここで指定された 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');

私が言ったように、私の解決策が得られるので、行列を翻訳する方法が間違っているのではないかと心配しています。

[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、因数分解が完了したことを意味しますがU(1,1) = 0、A = LU.

4

2 に答える 2

1

わかりました、「解決済み」と呼ぶためにこれに答えます。

これらのファイルは、機能するコードを表しています。誰かが同じ問題を抱えている場合に備えて、これを利用できるようにしています。 C の LAPACK を使用して、帯状の連立方程式を解きます。

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

そして、実装に関する追加の提案があれば、喜んで歓迎します!

私の次のステップは、ScaLAPACK で解くためにコア行列を配布することです。

ありがとう！

于 2013-03-21T17:02:39.593 に答える