Fortran - 逻辑索引

Question

假设我有一个矩阵A 是(m x n) 和一个向量B 这是(m x 1)。这个向量B 是一个由 0 和 1 组成的向量。

同样让标量 s 成为 B 中元素的总和。

我想创建一个矩阵C，它是s x n，对应于等于1的B的行，以及一个向量D，它是s x 1，这些元素在A中的位置.

Take as an example: 

      A = [1, 2, 3; 
           4, 5, 6; 
           7, 8, 9; 
           10, 11, 12; 
           13, 14, 15 ]

        B = [0; 1; 0; 1; 1]

    Then:
C = [ 4, 5, 6; 
     10, 11, 12; 
     13, 14, 15 ]
and
D = [2; 
     4
     5]

到目前为止，我的 fortran 代码如下所示：

PROGRAM test1
IMPLICIT NONE

 REAL, DIMENSION(5,3) :: A
INTEGER, DIMENSION(5,1) :: B = 0
INTEGER :: i, j, k

k = 1

!Create A matrix
do i=1,5
    do j=1,3
        A(i,j) = k
        k = k+1
    end do
end do

!Create B matrix
B(2,1) = 1
B(4,1) = 1
B(5,1) = 1


end program

在 matlab 中，我可以通过以下方式简单地创建 C：C = A(logical(B),:), and D 以类似的方式。

如何在 fortran 中做到这一点，避免循环？我查看了where 和forall 语句，但它们并不是我想要的。

Answer 1

正如@francescalus 所建议的，PACK 内在函数是 Matlab 逻辑切片的 Fortran 等价物，但只是部分。 Matlab 逻辑索引有两种类型：

• 整个数组逻辑索引：掩码必须与数组具有相同的形状，并且返回值的秩为 1（Matlab 中的向量）。之所以如此，是因为真值的位置是任意的，因此您不能保证基本上在矩阵中戳洞的结果是矩形的。 这就是 PACK 在 Fortran 中所做的事情。

  c = a(a < 1); % Matlab: a(m,n) -> c(s)
  C = PACK(A, A < 1) ! Fortran: A(m,n) -> C(s)
  

• 一维逻辑索引：掩码必须是一维的，用于在数组的一维内进行选择。其他维度可以选择整体或带有自己的索引。 这就是你想要的。

  b = a(:,1) > 0; % a(m,n) gives logical b(m)
  c = a(b,:); % a(m,n) selected with b(m) -> c(s,n)
  

但是，PACK 不直接承认这种用法。 @high-performance-mark 的解决方案向您展示了如何执行此壮举：SPREAD 基本上是 repmat，因此他的解决方案在 Matlab 中将如下所示：

b = a(:,1) > 0; % a(m,n) gives logical b(m)
bMat = repmat(b, 1, size(a,2)); % n copies of b(m) -> bMat(m,n)
c = reshape(a(bMat), [sum(b), size(a,2)]); % a(m,n) -> c(s,n)

因为a(bMat) 不是一个矩阵，而是一个大小为 s*n 的向量，因为使用了全矩阵选择形式，所以需要进行最终整形。另一个答案的 Fortran 代码正是这样做的：

c = RESHAPE( &
        PACK(a, & ! The matrix we are selecting from
            SPREAD(b==1, ! The == 1 is because you were using an INTEGER array, not LOGICAL
                   dim=2, ncopies=SIZE(a,2)) & ! This is the repmat-like code
        ), & ! The result of this PACK is cVec(s*n)
        [COUNT(b==1),SIZE(a,2)]) ! Reshape into (s,n)

分配给 D 的代码非常相似，但我们没有使用 A，而是从动态生成的数组中选择包含从 1 到 A 的第一维长度（即与 B) 的大小相同。这次不需要对结果进行任何 SPREAD 或 RESHAPE 掩码，因为源数组是一维的。

顺便说一句，Fortran 确实直接支持整数向量索引，因此您可以使用代码首先生成 D 向量（带有 B 的真实元素的索引），然后 然后 做：

C = A(D,:) ! Works the same in Fortran!

为自己节省传播和重塑。

Answer 2

对于需要基于某个向量 u 修剪矩阵 K 的问题（例如，在 FEA 中反转刚度矩阵时），您不能只使用 Ktrim = PACK(K,u/=0)，因为 K 是二维的。我发现这样的东西对于删除与向量为零的矩阵相对应的元素很有用：

n_trim = COUNT(u/=0)
NonZero(1:n_trim) = PACK([1:n], u/=0)
Ktrim = K(NonZero(1:n_trim),NonZero(1:n_trim))

其中 U 是一个向量，其零值应从矩阵 K 中移除，并将结果存储到 Ktrim。 K 的大小为 n*n，U 的大小为 n*1，Ktrim 的大小为 n_trim*n_trim。

用一些“反转”函数invKtrim(1:n_trim,1:n_trim) = invert(Ktrim)反转后，您可以像这样“取消修剪”：

invK = 0
invK(NonZero,NonZero) = invKtrim(1:n_trim,1:n_trim)