0

I am using the NIST Sparse BLAS v0.5 matrix-matrix multiplication routine, downloaded from http://math.nist.gov/~KRemington/fspblas/, to multiply a matrix one row at a time by a column vector. After calling the routine at a particular point (changes depending on the matrix), the do loop index it is contained within is being overwritten. I'm thinking it may have something to do with me reporting it as an upper triangular matrix to the routine, but don't see why it should make a difference for a single row vector. I've used this routine to multiply by the whole (same) matrix without any problem. Pseudocode reproducing the problem (makes it to end of fourth [outside] loop):

program row_vec_mult_test   
integer lwork,ind1,ind2,ind(4),ind_int,i,j,row_ind,val_ind
integer descra(9),H_col(:),H_pntrb(1),H_pntre(1),H_col_tmp(4)
double precision in(10),out(10),work(5),H_val(:),
 1 out_tmp(1),H(10,10),H_val_tmp(4)

allocatable H_col,H_val

descra(1) = 1
descra(2) = 2
descra(3) = 0
descra(4) = 1
descra(5) = 0
lwork = 5

do i = 1,10
 in(i) = i
end do

H = 0.0d0
H(1,6) = 3
H(2,8) = 6
H(4,5) = 2
H(7,9) = 1
H(8,9) = 7
H(9,10) = 5

do i = 1,10
 do j = 1,10
  if (i==j) then
   H(i,j) = i
  else
   H(j,i) = H(i,j)
  end if
 end do
end do

do ind1 = 1,10
 row_ind = 1
 val_ind = 1
 H_pntrb(1) = 1
 H_col_tmp = 0
 H_val_tmp = 0.0d0
 do ind2 = 1,10
  if(H(ind1,ind2) /= 0.0d0) then
   H_val_tmp(val_ind) = H(ind1,ind2)
   H_col_tmp(val_ind) = ind2
   val_ind = val_ind + 1
   row_ind = row_ind + 1
   H_pntre(1) = row_ind
  else
  end if
 end do
 where(H_val_tmp /= 0)
  ind = 1
 elsewhere
  ind = 0
 end where
 ind_int = sum(ind)
 if (allocated(H_val)) deallocate(H_val)
 allocate(H_val(ind_int))
 H_val = H_val_tmp(:ind_int)
 if (allocated(H_col)) deallocate(H_col)
 allocate(H_col(ind_int))
 H_col = H_col_tmp(:ind_int)
 !print *,
 !print *, ind1
 !print *, int(H_val)
 !print *, H_col
 !print *, H_pntrb
 !print *, H_pntre
 call dcsrmm(0,1,1,10,1.0d0,descra,H_val,H_col,H_pntrb,
 1  H_pntre,in,10,0.0d0,out_tmp,1,work,lwork)
   out(ind1) = out_tmp(1)
  end do

  do i = 1,10
   print *, out(i)
  end doroutine to

end

And output:

At line 81 of file csr_row_mult_test.f
Fortran runtime error: Array reference out of bounds for array 'out', upper bound of dimension 1 exceeded (1074790400 > 10)

Backtrace for this error:
+ in the main program
+ /lib64/libc.so.6(__libc_start_main+0xfd) [0x3c5cc1ed1d]

Compiled on RHEL6, gfortran 4.4.7, with -fimplicit-none -Wall -Wtabs -fbounds-check -fbacktrace -O5

P.S. Although similar to my previous question, they are in fact different programs, and none of the solutions to that problem solved this one.

alex_d
  • 111
  • 5
  • 1
    Where did you get NIST sparse BLAS, what version? If you built it yourself, you should also use those bounds checking flags for compiling the library. – steabert Jul 20 '14 at 07:16
  • Added the information to the original post. However, looking around a bit while checking the version I found a more recent release that I am going to attempt to use. – alex_d Jul 21 '14 at 13:50

1 Answers1

0

The problem was my specification of a symmetric rectangular matrix. I didn't try to change it to specify a non-symmetric matrix, although that may have worked. In the Sparse BLAS implementation referenced in my question there's a vector-vector multiplication routine ddoti, so modifying the row of the matrix to represent a single vector and using it effectively eliminated that problem.

As an aside, in the later implementations of the library include more verbose error reporting, which is how I determined the source of the error.

alex_d
  • 111
  • 5