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

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

subroutine conjGrad1(Kffv,bvec,n,nrows,rowIndex,Icol,Idiag,ierr)

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

real*8 Kffv(n),bvec(nrows), xvec(nrows)
real*8 rvec(nrows)

real

*8 normv1,normv2, normv0*

integer rowIndex(nrows+1)

integer Icol(n)

real8 Kdiag(nrows)

integer rowIndex(nrows+1)

integer Icol(n)

real

integer Idiag(nrows)

real

*8 pvec1(nrows), vec1(nrows)*

real8 alpha1,beta1,denom1

real

integer i1,i2, iiter, ierr

CALL OMP_SET_NUM_THREADS(2)

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-24*normv0) 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

end subroutine conjGrad1]

[/code]

==================================================

When I try to run my program, it aborts for no obvious reason, in one of the parallel regions. When I rerun the same exact program (without any change in the code) with 1 thread, it runs without any problem. What could I be doing wrong with my openMP directives?