如何从矩阵中删除几列答案

【问题标题】：How to remove several columns from a matrix如何从矩阵中删除几列
【发布时间】：2018-01-21 04:06:54
【问题描述】：

Integer :: NBE,ierr,RN,i,j
Real(kind=8), allocatable :: AA1(:,:),AA2(:,:)

NBE=40 
RN=3*NBE-2
Allocate(AA1(3*NBE,3*NBE),AA2(3*NBE,RN),stat=ierr)

If (ierr .ne. 0) Then
    print *, 'allocate steps failed 1'
    pause
End If

Do i=1,3*NBE
  Do j=1,3*NBE
     AA1(i,j)=1
  End Do
End Do

我想从矩阵AA1 中删除第 97 列和第 113 列，然后这个矩阵变为AA2。我就是想知道有没有Fortran的命令可以实现这个操作？

【问题讨论】：

标签： matrix fortran

【解决方案1】：

Alexander Vogt's answer 给出了使用向量下标 来选择要包含的数组元素的概念。该答案使用

构造向量下标数组

[(i,i=1,96),(i=98,112),(i=114,3*NBE)]

有些人可能会考虑

AA2 = AA1(:,[(i,i=1,96),(i=98,112),(i=114,3*NBE)])

阅读不清晰。可以使用“临时”索引向量

integer selected_columns(RN)
selected_columns = [(i,i=1,96),(i=98,112),(i=114,3*NBE)]
AA2 = AA1(:,selected_columns)

但这并不能解决数组构造函数不好的问题，尤其是在更复杂的情况下。相反，我们可以创建一个掩码并使用我们常用的技术：

logical column_wanted(3*NBE)
integer, allocatable :: selected_columns(:)

! Create a mask of whether a column is wanted
column_wanted = .TRUE.
column_wanted([97,113]) = .FALSE.

! Create a list of indexes of wanted columns
selected_columns = PACK([(i,i=1,3*NBE)],column_wanted)
AA2 = AA1(:,selected_columns)

【讨论】：

好答案。但是，我在 PACK 命令中发现了一个错误。应该是 PACK([(i,i=1,3*NBE)],column_wanted)。
谢谢，我混淆了哪个数组的大小。在前面的答案中也有一个错误。现在更正（我希望）。

【解决方案2】：

这是一个简单的衬里：

AA2 = AA1(:,[(i,i=1,96),(i=98,112),(i=114,3*NBE)])

解释：

（内部）为索引[1,...,96,98,...,112,114,...,3*NBE]构造一个临时数组
（外部）复制矩阵，只考虑索引数组中的列

好的，我屈服于@IanBush...更简单的是做三个专门的任务：

AA2(:,1:96)   = AA1(:,1:96)
AA2(:,97:111) = AA1(:,98:112)
AA2(:,112:)   = AA1(:,114:)

【讨论】：

简单！？我肯定会把它写成 3 行，每个“块”一个
@IanBush Mhm，好吧......虽然我更喜欢让我的代码尽可能简洁。我更新了我的答案以反映您的建议。

【解决方案3】：

我家里没有fortran编译器，所以无法测试。但我会这样做：

i = 0
DO j = 1, 3*NBE
    IF (j == 97 .OR. j == 113) CYCLE
    i = i + 1
    AA2(:, i) = AA1(:, j)
END DO

CYCLE 命令意味着不应再执行循环的其余部分，而应开始下一次迭代。因此，i 不会增加，所以当 j=96 时，然后是 i=96，当 j=98 时，然后是 i=97，当 j=114 时，然后是 i=112。

多说几句：由于 Fortran 的内存布局，您希望以最快的速度循环第一个索引，依此类推。因此，如果您将代码更改为：

Do j=1,3*NBE    ! Outer loop over second index
  Do i=1,3*NBE  ! Inner loop over first index
    AA1(i,j)=1
  End Do
End Do

（当然，只需AA1(:,:) = 1 或AA1 = 1 就可以更轻松地完成如此简单的初始化。

【讨论】：