稀疏 BLAS 矩阵行向量积覆盖循环索引答案

【问题标题】：Sparse BLAS matrix-row vector product overwriting loop index稀疏 BLAS 矩阵行向量积覆盖循环索引
【发布时间】：2014-07-18 16:40:09
【问题描述】：

我正在使用从http://math.nist.gov/~KRemington/fspblas/ 下载的 NIST Sparse BLAS v0.5 矩阵-矩阵乘法例程，一次将矩阵乘以列向量。在特定点调用例程后（根据矩阵而变化），它包含的 do 循环索引将被覆盖。我认为这可能与我将其作为上三角矩阵报告给例程有关，但不明白为什么它应该对单行向量产生影响。我已经使用此例程乘以整个（相同）矩阵而没有任何问题。重现问题的伪代码（使其结束第四个 [外部] 循环）：

program row_vec_mult_test   
integer lwork,ind1,ind2,ind(4),ind_int,i,j,row_ind,val_ind
integer descra(9),H_col(:),H_pntrb(1),H_pntre(1),H_col_tmp(4)
double precision in(10),out(10),work(5),H_val(:),
 1 out_tmp(1),H(10,10),H_val_tmp(4)

allocatable H_col,H_val

descra(1) = 1
descra(2) = 2
descra(3) = 0
descra(4) = 1
descra(5) = 0
lwork = 5

do i = 1,10
 in(i) = i
end do

H = 0.0d0
H(1,6) = 3
H(2,8) = 6
H(4,5) = 2
H(7,9) = 1
H(8,9) = 7
H(9,10) = 5

do i = 1,10
 do j = 1,10
  if (i==j) then
   H(i,j) = i
  else
   H(j,i) = H(i,j)
  end if
 end do
end do

do ind1 = 1,10
 row_ind = 1
 val_ind = 1
 H_pntrb(1) = 1
 H_col_tmp = 0
 H_val_tmp = 0.0d0
 do ind2 = 1,10
  if(H(ind1,ind2) /= 0.0d0) then
   H_val_tmp(val_ind) = H(ind1,ind2)
   H_col_tmp(val_ind) = ind2
   val_ind = val_ind + 1
   row_ind = row_ind + 1
   H_pntre(1) = row_ind
  else
  end if
 end do
 where(H_val_tmp /= 0)
  ind = 1
 elsewhere
  ind = 0
 end where
 ind_int = sum(ind)
 if (allocated(H_val)) deallocate(H_val)
 allocate(H_val(ind_int))
 H_val = H_val_tmp(:ind_int)
 if (allocated(H_col)) deallocate(H_col)
 allocate(H_col(ind_int))
 H_col = H_col_tmp(:ind_int)
 !print *,
 !print *, ind1
 !print *, int(H_val)
 !print *, H_col
 !print *, H_pntrb
 !print *, H_pntre
 call dcsrmm(0,1,1,10,1.0d0,descra,H_val,H_col,H_pntrb,
 1  H_pntre,in,10,0.0d0,out_tmp,1,work,lwork)
   out(ind1) = out_tmp(1)
  end do

  do i = 1,10
   print *, out(i)
  end doroutine to

end

然后输出：

At line 81 of file csr_row_mult_test.f
Fortran runtime error: Array reference out of bounds for array 'out', upper bound of dimension 1 exceeded (1074790400 > 10)

Backtrace for this error:
+ in the main program
+ /lib64/libc.so.6(__libc_start_main+0xfd) [0x3c5cc1ed1d]

在 RHEL6 上编译，gfortran 4.4.7，-fimplicit-none -Wall -Wtabs -fbounds-check -fbacktrace -O5

附：虽然与我之前的问题相似，但它们实际上是不同的程序，并且该问题的解决方案都没有解决这个问题。

【问题讨论】：

您从哪里获得 NIST 稀疏 BLAS，什么版本？如果您自己构建它，您还应该使用这些边界检查标志来编译库。
将信息添加到原始帖子中。但是，在检查版本时环顾四周，我发现了一个更新的版本，我将尝试使用。

标签： matrix fortran sparse-matrix blas

【解决方案1】：

问题是我对对称矩形矩阵的规范。我没有尝试将其更改为指定非对称矩阵，尽管这可能有效。在我的问题中引用的 Sparse BLAS 实现中，有一个向量-向量乘法例程 ddoti，因此修改矩阵的行以表示单个向量并有效地使用它可以消除该问题。

顺便说一句，在该库的后续实现中包含更详细的错误报告，这就是我确定错误来源的方式。

【讨论】：