Search code examples
loopsfortransubtraction

Subtract a vector from a matrix (row-wise)


Suppose I have a matrix A in fortran which is (n,m), and a vector B which is (1,m). I want to subtract the vector B from all rows of A without using a loop.

As of now I have only been able to do it with a loop:

PROGRAM Subtract
IMPLICIT NONE

REAL, DIMENSION(250,5) :: A
INTEGER, DIMENSION(1,5) :: B
INTEGER :: i

B(1,1) = 1
B(1,2) = 2
B(1,3) = 3
B(1,4) = 4
B(1,5) = 5

CALL RANDOM_NUMBER(A)

do i=1,250
    A(i,:) = A(i,:) - B(1,:)
end do


end program

But this is very inefficient. E.g. in matlab one can do it in one line using the function reptmat. Any suggestion on better way to do this?


Solution

  • You can use spread for this:

    A = A - spread( B(1,:), 1, 250 )
    

    Please note that Fortran is column-major, so B(1,:) is in general not contiguous in memory, and a temporary array is created. [ It is in your case since you have only one column - but still worth mentioning. ]

    In the same way, looping over the first index of A is inefficient. It would probably speed things up a lot if you transposed your matrices. Then, a loop solution might even be faster than the one using spread. (But this depends on the compiler. )