我尝试搜索有关此主题的类似帖子,但似乎没有人关心带状系统(...)。
我有兴趣使用LAPACK / ScaLAPACK从C代码解决带状矩阵。首先,我想使用LAPACK实现顺序解决方案,然后再尝试使用ScaLAPACK。
问题:两种语言之间的行主/列主差异似乎影响了我的解决过程。这是我打算解决的系统:
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
。