# Program aborts for no obvious reason with openMP

Hello, I am having a program using the following subroutine:

!=========================================================================================================

implicit none

! solve a system of equations [A]{x}={b}, using the conjugate gradient method

integer n, nrows ! n: number of values in the column representation of sparse matrix, nrows: number of rows in system of equations
real8 Kffv(n),bvec(nrows), xvec(nrows)
real
8 rvec(nrows)
real8 normv1,normv2, normv0

integer rowIndex(nrows+1)
integer Icol(n)
real
8 Kdiag(nrows)
integer Idiag(nrows)
real8 pvec1(nrows), vec1(nrows)
real
8 alpha1,beta1,denom1

integer i1,i2, iiter, ierr

ierr = 0
xvec = 0.d0
pvec1 = bvec
normv1 = 0.d0

do iiter=1,1000

if(iiter.eq.1) then
rvec = bvec
pvec1 = bvec

normv1 = dot_product(bvec,bvec)

normv2 = normv1
normv0 = normv1
else
beta1 = normv1/normv2
!write(101,) ‘step 1…’
!\$OMP PARALLEL DO
do i1 = 1,nrows
pvec1(i1) = rvec(i1) + beta1
pvec1(i1)
end do !i1

endif
! begin by calculating vec1, which is the product [A]{p}
vec1 = 0.d0
!write(101,*) ‘step 2…’
!\$OMP PARALLEL DO
do i1 = 1,nrows
do i2 = rowIndex(i1), rowIndex(i1+1) - rowIndex(1)
vec1(i1) = vec1(i1) + Kffv(i2)*pvec1(Icol(i2))
end do !i2
end do !i1

denom1 = dot_product(vec1,vec1)

alpha1 = normv1/denom1
!write(101,) ‘step 3…’
!\$OMP PARALLEL DO
do i1 = 1,nrows
xvec(i1) = xvec(i1) + alpha1
pvec1(i1)
rvec(i1) = rvec(i1) - alpha1*vec1(i1)
end do !i1

normv2 = normv1

normv1 = dot_product(rvec,rvec)

if(normv1.le.1.d-24normv0) then
exit
elseif(normv0.le.1.d-8.and.normv1.le.1.d-38) then
exit
endif
write(501,
) iiter,normv1
end do! iiter

if(iiter.ge.1000) ierr = 1

bvec = xvec ! finally, when we return from the subroutine, we replace the RHS vector bvec by the solution vector xvec

!write(501,) ‘iiter=’,iiter
!write(101,
) normv1,normv0
!stop
return