在fortran中填充(对称)距离矩阵的最有效的cpu方式

Question

我正在研究在有限元问题中实现周期性边界条件，我想将边界a上的节点与边界B上的节点配对，给定一个向量trans，该向量将边界a与边界B直线相连。边界a上的节点在列表g1中给出;在B中它们是g2。节点坐标在mesh%nodes(:,nodenum)中查找。

我决定通过为每个节点创建一个距离矩阵来做到这一点，我意识到这不是最有效的方法，说实话，我不希望通过优化这个算法来节省大量的时间。这个问题比较学术性。

我知道Fortran以列为主顺序存储，另一方面，数组将是对称的，当数组完成时，我想对它进行列切片以找到最近的节点。所以问题是如何填充它?

这是我天真的尝试。

subroutine autopair_nodes_in_groups(mesh, g1, g2, pairs, trans)
type(meshdata)                              :: mesh


integer(kind=sp)                            :: i,j
integer(kind=sp),dimension(:)               :: g1,g2
integer(kind=sp),dimension(:,:)             :: pairs

real(kind=dp)                               :: trans(3) !xyz translate

real(kind=dp)                               :: dist_mat(size(g1),size(g2)) 
real(kind=dp)                               :: p1(3), p2(3)

dist_mat = -1.0_wp

! make a distance matrix
do j=1,size(g2)
p2 = mesh%nodes(1:3,g2(j))-trans
do i=1,j
p1 = mesh%nodes(1:3,g1(i))
dist_mat(i,j) = norm2(p1-p2)              !equivalent to norm2(n1pos+trans-n2pos)
if (i.ne.j) dist_mat(j,i) = dist_mat(i,j) !fill symmetry
end do
end do

! Remainder of routine to find nearest nodes


end subroutine autopair_nodes_in_groups

据我所知，这个问题在内存访问方面是有效的，直到一个数组对称。

Answer 1

要进行快速的最近邻搜索，您应该实现一个搜索复杂度为O(log(N))的树结构，而不是查看所有点对点距离O(N^2)。

无论如何，关于对称矩阵处理，你将有:

! Storage size of a symmetric matrix
elemental integer function sym_size(N)
integer, intent(in) :: N
sym_size = (N*(N+1))/2
end function sym_size
! Compute pointer to an element in the symmetric matrix array
elemental integer function sym_ptr(N,i,j) 
integer, intent(in) :: N,i,j
integer :: irow,jcol
! Column-major storage
irow = merge(i,j,i>=j)
jcol = merge(j,i,i>=j)
! Skip previous columns
sym_ptr = N*(jcol-1)-((jcol-1)*(jcol-2))/2
! Locate in current column
sym_ptr = sym_ptr + (irow-jcol+1)
end function sym_ptr

然后做你的工作:

N = size(g2)
allocate(sym_dist_mat(sym_size(N)))
do j=1,size(g2)
p2 = mesh%nodes(1:3,g2(j))-trans
do i=j,size(g2)
p1 = mesh%nodes(1:3,g1(i))
sym_dist_mat(sym_ptr(N,i,j)) = norm2(p1-p2)  
end do
end do

然后minloc函数应该看起来像这样(未经测试):

! minloc-style
pure function symmetric_minloc(N,matrix,dim) result(loc_min)
integer, intent(in) :: N
real(8), intent(in) :: matrix(N*(N+1)/2)
integer, intent(in) :: dim
real(8) :: dim_min(N),min_column
integer :: loc_min(N)
select case (dim)
! Row/column does not matter: it's symmetric!
case (1,2)
dim_min = huge(dim_row)
loc_min = -1
ptr = 0
do j=1,N
! Diagonal element m(j,j)
ptr=ptr+1
min_column = matrix(ptr)
if (min_column<=dim_min(j)) then
loc_min(j) = j
dim_min(j) = min_column
end if
! Lower-diagonal part at column j
do i=j+1,N
ptr=ptr+1
! Applies to both this column,
if (matrix(ptr)<=dim_min(j)) then
loc_min(j) = i
dim_min(j) = matrix(ptr)
end if
! And the i-th column
if (matrix(ptr)<=dim_min(i)) then
loc_min(i) = j
dim_min(i) = matrix(ptr)
end if
end do
end do
case default
! Invalid dimension
loc_min = -1
end select
end function symmetric_minloc