1

我想问你一个我在使用 Fortran 管理 CSR/CSC* 矩阵时经常遇到的问题。假设我们有一个具有 N 个实数值的向量 V。该向量之前已经分配了一定的大小。现在我们必须在它中间的索引 P 处添加一个值。蛮力代码将是：

allocate(tempV(N))
tempV=V
deallocate(V)
allocate(V(N+1))
V=(/tempV(1:P-1), newValue, tempV(P:N)/)
deallocate(tempV)

显然，如果做一次没问题，但重复上千次就不会那么有效了。每个我想插入的值，内存都会填满和清空 4 次。

我想知道解决这个问题的更好程序。您可以提出普通的 Fortran（首选），也可以通过 MKL/Lapack/Blas 等库提出一些解决方案。

附录：我可以用 RESHAPE 来做吗？通过这个定义（与我的 Fortran 书籍定义相同），我可以做类似的事情

REAL, DIMENSION(1:1) :: newPad = (/ newValue /)
V=RESHAPE(V, (/ N+1 /), PAD=newPad)

现在值已添加到 V 的末尾，所以我用

V=(/V(1:P-1), V(N+1:N+1), V(P:N) /)

这样，可以避免显式创建临时向量并丢失分配。

由于 RESHAPE 可以在库中进行并行化，它是否高效且可扩展？

* PS：为了清楚起见，CSR = 压缩稀疏行格式，CSC = 压缩稀疏列格式，更多信息在这里：

MKL 定义： http: //software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-9FCEB1C4-670D-4738-81D2-F378013412B0.htm

4

2 回答 2

3

move_alloc为此目的引入了Fortran 2003 子例程。它将分配状态、数组边界、动态类型、类型参数和值从源移动到目标，而不实际复制数据。源变量被释放。

用一个简短的例子修改你的代码move_alloc只需要一个复制操作：

allocate(tempV(N+1)) 
tempV(:P-1) = V(:P-1)
tempV(P)    = newValue
tempV(P+1:) = V(P:) 
call move_alloc(tempV, V)

当然，一次为多个项目执行此操作会减少分配开销，但这对您来说可能是不可能的。

编辑
至于指针的建议，如果您不能使用任何 F2003 功能，您可能可以使用有点像这样的子例程：

pure subroutine insert(arr, val, pos)
  real, pointer, intent(inout) :: arr(:)
  real, intent(in)             :: val
  integer, intent(in)          :: pos

  real, pointer :: temp(:)

  if(associated(arr)) then
    allocate(temp(lbound(arr,1):ubound(arr,1) + 1))
    ! ...or perhaps check/do something if pos is not within these bounds

    temp(:pos-1) = arr(:pos-1)
    temp(pos)    = val 
    temp(pos+1:) = arr(pos:)

    deallocate(arr)
    arr => temp
  endif
end subroutine insert

当然，您可以轻松地将其调整为适合您的目的或使其更通用。您可以将其与分配的指针变量一起使用：

real, pointer :: V(:)
! :
allocate(V(10))
V = 1.
! :
call insert(V, 3.141592, 5)
于 2013-01-30T16:45:44.953 回答
1

问题是，您是否需要在插入后立即使用 CSR 格式，或者您可以先以更易于插入的格式进行更多插入，然后再转换为 CSR。

你可以看看sparskit和它的链表存储格式，或者你可以为新元素创建一些链表，一旦你完成插入，你可以建立一个新的 CSR 格式的矩阵，节省大量的数据重组.

于 2013-01-31T11:38:43.893 回答