CODEWITHSUNDEEP
matrixfortran
Jeremy_Tamu
Jeremy_Tamu

Reputation: 755

How to remove several columns from a matrix

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

I want to remove columns 97 and 113 from the matrix AA1 and then this matrix becomes AA2. I just want to know if any command of Fortran can realize this operation?

Upvotes: 3

Views: 990

Answers (3)

francescalus
francescalus

Reputation: 32396

Alexander Vogt's answer gives the concept of using a vector subscript to select the elements of the array for inclusion. That answer constructs the vector subscript array using

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

Some may consider

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

to be less than clear in reading. One could use a "temporary" index vector

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

but that doesn't address the array constructor being not nice, especially in more complicated cases. Instead, we can create a mask and use our common techniques:

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)

Upvotes: 7

Alexander Vogt
Alexander Vogt

Reputation: 18108

Here's a simple one liner: 

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

Explanation:

  1. (Inner part) Construct a temporary array for the indices [1,...,96,98,...,112,114,...,3*NBE]

  2. (Outer part) Copy the matrix and only consider the columns in the index array

OK, I yield to @IanBush... Even simpler would be to do three dedicated assignments: 

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

Upvotes: 5

chw21
chw21

Reputation: 8140

I don't have a fortran compiler here at home, so I cant test it. But I'd do something line this:

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

The CYCLE command means that the rest of the loop should not be executed any more and the next iteration should start. Thereby, i will not get incremented, so when j=96, then i=96, when j=98 then i=97, and when j=114 then i=112.

A few more words: Due to Fortran's memory layout, you want to cycle over the first index the fastest, and so forth. So your code would run faster if you changed it to:

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

(Of course, such an easy initialisation can be done even easier with just AA1(:,:) = 1 of just AA1 = 1.

Upvotes: 1

Related Questions