OpenMP code not giving any result in Intel Fortran

Question

I wrote a serial code for the conjugate gradient method and tried to parallelize it with OpenMP ( my platform is Intel cluster )

When I use the serial code I am getting output as follows :-

Test data (matrix and right hand side) :
  5.00000  0.00000 -1.00000 -1.00000  1.00000
  0.00000  5.00000 -1.00000 -1.00000  1.00000
 -1.00000 -1.00000  5.00000 -1.00000  1.00000
 -1.00000 -1.00000 -1.00000  5.00000  1.00000
  Solution for the linear system:
  0.37500  0.37500  0.43750  0.43750

But when I am using Openmp for a particular do loop I am not getting the solution . The output is only this :

Test data (matrix and right hand side) :
  5.00000  0.00000 -1.00000 -1.00000  1.00000
  0.00000  5.00000 -1.00000 -1.00000  1.00000
 -1.00000 -1.00000  5.00000 -1.00000  1.00000
 -1.00000 -1.00000 -1.00000  5.00000  1.00000

There are no errors during compiling. Can anybody help ?

I have copy-pasted the OpenMP code.

  module cg
     contains
     subroutine cg_method (n,a,y,x)
     use omp_lib
     implicit none
     integer :: n, k , i , j
     integer,parameter :: sze = 24 , nt = 4
     double precision :: tol = 2.d-16
     double precision :: a(0:sze,0:sze),x(0:sze),y(0:sze)
     double precision :: d(0:n-1),g(0:n-1),auxD(0:n-1),alpha,beta
     double precision :: num,den,aux1,aux2,dist,xnorm

    !$call omp_set_num_threads(nt)
    ! start with x at origin

     do i = 0,n-1
        x(i) = 0.d0
     enddo

    ! initialize d and g
    ! d = -g = -(a*x-y) = y as x = 0

     do i = 0,n-1
       aux1 = y(i)
       d(i) = aux1
       g(i) = -aux1
     enddo

    ! perform at most n steps of CG algo

     do k = 0,n

    ! compute new alpha
    ! alpha = -(d(transp)*g)/(d(transp)*(a*d))

       num = 0.d0
       den = 0.d0
       do i = 0,n-1
          num = num + d(i)*g(i)
          aux1 = 0.d0
         do j = 0 , i-1
           aux1 = aux1 + a(j,i)*d(j)
         enddo
            do j = i,n-1
               aux1 = aux1 + a(i,j)*d(j)
            enddo
          auxD(i) = aux1
          den = den + d(i)*aux1
        enddo


     alpha = -num/den

    ! compute the norm of x and alpha*d and find a new x
    !x = x + alpha*d , then check if x is close in order to stop the process
    !before n complete steps

     xnorm = 0.d0
     dist = 0.d0

     do i = 0,n-1
        aux1 = x(i)
        xnorm = xnorm + aux1*aux1
        aux2 = alpha*d(i)
        dist = dist + aux2*aux2
        x(i) = aux1 + aux2
     enddo

    !compute new g : g + alpha*(a*d)

     do i = 0,n-1
        g(i) = g(i) + alpha*auxD(i)
     enddo

    ! compute new beta :
    ! beta = (g(transp)*(a*d))/(d(transp)*(a*d))

     num = 0.d0
     do i = 0,n-1
        num = num + g(i)*auxD(i)
     enddo

     beta = num/den

    !compute new d : d = -g + beta*d
    !$omp parallel default(none) shared(beta,d,g) private(i,n)
    !$omp do
     do i = 0,n-1
        d(i) = -g(i) + beta*d(i)
     enddo
    !$omp enddo
    !$omp end parallel
     enddo !k loop

     end subroutine cg_method 
     end module cg

     program test

     use cg
     use omp_lib
     integer,parameter :: sze = 24
     double precision :: a(0:sze,0:sze), y(0:sze),x(0:sze)
     integer :: n , i , j

     n = 4

     !define a matrix
     a(0,0)=5.d0;a(0,2)=-1.d0;a(0,3)=-1.d0
     a(1,1)=5.d0;a(1,2)=-1.d0;a(1,3)=-1.d0
     a(2,2)=5.d0;a(2,3)=-1.d0;a(3,3)=5.d0

    !define b vector
     y(0)=1.d0;y(1)=1.d0;y(2)=1.d0;y(3)=1.d0

     print*,'Test data (matrix and right hand side) :'
     do i = 0,n-1
       write(*,100) (a(j,i),j=0,i-1),(a(i,j),j=i,n-1),y(i)
     enddo

     call cg_method(n,a,y,x) !perform CG method
     print *,' Solution for the linear system:'
        write(*,100)  (x(i),i=0,n-1)
    100 format(10F9.5)
     end program test

OpenMP code not giving any result in Intel Fortran

Answers (1)

Related Questions