Reference passing is changing the values of a matrix

Question

I'm trying to make a code in Fortran90 (compiling with ifort) in which I multiply two matrices. I'm writing code for this because one of the matrices is sparse, so you can do the multiplication without allocating memory for the whole matrix.

I have two subroutines. The first one, multiplies the sparse matrix (k in diagonal, and l at both sides of the diagonal) and a vector (b). The result is passed to the main function trough pointer r. I chose using subroutines and no functions because with this, I don't need to allocate memory again inside the subroutine.

    subroutine matSparXVec(k, l, b, r)

    implicit none
    real, intent(in)        ::k, l
    real, intent(in), dimension(:)  ::b
    real, intent(out), dimension(:) ::r
    integer             ::ierr, i

    r = (/&
        k*b(1)+l*b(2), &
        (l*b(i-1)+k*b(i)+l*b(i+1),i=2,size(b)-1),& 
        l*b(size(b)-1)+k*b(size(b))&
        /)

end subroutine matSparXVec

The second subroutine uses the first one to multiply a sparse matrix (k, l) with another matrix (B):

subroutine matSparXMat(k, l, B, R)
    implicit none
    real, intent(in)            ::k, l
    real, intent(in), dimension(:,:)    ::B
    real, intent(out), dimension(:,:)   ::R

    call  matsparXVec(k, l, B(1, :), R(1, :))
    R(1, :) = R(1, :) + l * B(1, :)

    do i = 2, (size(R)-2)
        call matsparXVec(k,l,B(i,:),R(i,:))
        R(i,:)=R(i-1,:)+R(i,:)
    enddo 

    call matsparXVec(k,l,B(size(B),:),R(size(R),:))
    R(size(R),:) = R(size(R),:) + l * B(size(B)-1,:)


end subroutine matSparXmat

Now the issue is that, in some place of the subroutine matSparXmat, I'm modifying the data pointed by B(:,:). In example, with the code:

    implicit none
    real, dimension(:,:), allocatable   ::BB, RR
    integer                 ::i, j, ierr, n, m
    real, parameter             ::k = 4.0, l = 1.0

    n=3 !Dimension vector
    m=3 !Dimension del segundo orden

    allocate(RR(n,m), stat=ierr)
    allocate(BB(n,m), stat=ierr)

    forall(i = 1:size(BB(:,1)), j=1:size(BB(1,:)))
        BB(i,j)=i+j
        RR(i,j) = 0 
    endforall

    do i=1,size(BB(:,1))
        print *, BB(:,i)
    enddo

    call matSparXMat(k, l, BB, RR)

    do i=1,size(BB(:,1))
        print *, BB(:,i)
    enddo 

I get the output:

 2.000000       3.000000       4.000000    
   3.000000       4.000000       5.000000    
   4.000000       5.000000       6.000000  

  4.6526092E+33   3.000000      1.9366391E+31
   3.000000       4.000000       5.000000    
   4.000000       5.000000       6.000000 

In which you can see that the values of BB has been modified.

Answer 1

Looks like it went wrong here:

do i = 2, (size(R)-2)                      !<---- here
   call matsparXVec(k,l,B(i,:),R(i,:))
   R(i,:)=R(i-1,:)+R(i,:)
enddo 

size(R)=9 but size(R(:,1))=3. However, using gfortran I do get the correct output with an error

*** glibc detected *** ./sparse_mults: free(): invalid next size (fast): 0x00000000014e4010 ***

But when I use ifort 2013, I get the value 1.93...E+31 that you have in the top right corner, but otherwise it is correct. Not sure what is going on, but if I come up with something I'll let you know.