C******************************************************************************
C
C               Example 15 : Omp_Jacobi.f
C
C   Objective : Jacobi method to solve Ax = b matrix system of 
C               linear equations. 
C               This example demonstrates the use of
C               PARALLEL section and DO directive
C               
C   Input     : Real Symmetric Positive definite Matrix Matrix_A and 
C               vector Input_B. Read files (mdatjac.inp) for Matrix A
C	             and (vdatjac.inp) for Vector b
C
C   Output    : The solution of  Ax=b if the method converges and 
C               the number of iterations if the method fails to converge
C  
C   Necessary Condition : Number of threads should be less than or 
C                         equal to number of processors of a SMP node.
C
C   Note:       Matrix_A should be strongly diagonally dominant in
C               order to insure convergence.
C
C******************************************************************************

      PROGRAM JacobiMethod

      include 'jacconst.h'

      real     Matrix_A(0:11, 0:11)
      real     x(0:11)
      real     Input_B(0:11)
      integer  n
      real     Tolerance
      integer  max_iter
      integer  converged
      integer  Jacobi

      n = SIZE  
      Tolerance = EPSILON  
      max_iter = MAX_ITERATIONS

      call Read_matrix(Matrix_A, n)
      call Read_vector(Input_B, n)

      converged = Jacobi(Matrix_A, x, Input_B, n, Tolerance, max_iter)

      if (converged .EQ. 1 ) then
          call Print_vector('The solution is', x, n)
      else
          print *, 'Failed to converge in ', max_iter,
     +                   ' iterations'
      endif
      end


C *******************************************************************  
C  Return 1 if iteration converged, 0 otherwise   
C  MATRIX_T is just a 2-dimensional array         
      integer function Jacobi(Matrix_A, x, Input_B, n, Tolerance, 
     +                        max_iter)
      real     Matrix_A(0:n-1, 0:n-1)          
      real     x(0:n-1)        
      real     Input_B(0:n-1)        
      integer  n          
      real     Tolerance        
      integer    max_iter   
      integer    i, j
      integer    Iteration
      real  retval
      real  x_old(0:11)

      real Distance 

C  Initialize x   

C$OMP PARALLEL

C$OMP DO
      do 100 i = 0 ,n-1
          x(i) = Input_B(i)
 100  continue
C$OMP END DO
C$OMP END PARALLEL


      Iteration = 0
      retval = 99 
      do while ((Iteration .LT. max_iter) .AND.
     +          (retval .GE. Tolerance))
          Iteration = Iteration + 1

C$OMP PARALLEL

C$OMP DO
          do 200 i = 0 ,n-1
              x_old(i) = x(i)
 200      continue
C$OMP END DO


C$OMP DO 
          do 300 i = 0, n-1
              x(i) = Input_B(i)
              do 400 j = 0 , i-1
                  x(i) = x(i) - Matrix_A(i, j)*x_old(j)
 400          continue
              do 500 j = i+1, n-1
                  x(i) = x(i) - Matrix_A(i, j)*x_old(j)
 500          continue
              x(i) = x(i)/Matrix_A(i,i)
 300      continue
C$OMP END DO 
          retval = Distance(x,x_old,n)

C$OMP END PARALLEL

      end do


      if (retval .LT. Tolerance) then
          Jacobi = 1
      else
          Jacobi = 0
      endif

      return
      end 


C *******************************************************************  
      real function Distance(x, y, n)
      real x(0:n-1)
      real y(0:n-1)
      integer n
      integer i
      real sum

      sum = 0.0

C$OMP PARALLEL DO REDUCTION(+:sum) 
      do 100 i = 0 , n - 1
          sum = sum + ((x(i) - y(i))*(x(i) - y(i)))
 100  continue
      Distance =  sqrt(sum)
      return
      end


C *******************************************************************  
      subroutine Read_matrix(Matrix_A, n)
      real    Matrix_A(0:n-1, 0:n-1)        
      integer n        
      integer i, j

C	 .......Read the Matrix From Input file ......
	   open(unit=12, file = './data/mdatjac.inp')
           read(12,*) NoofRows, NoofCols
           do i = 0,NoofRows-1
              read(12,*) (Matrix_A(i,j),j=0,NoofCols-1)
           enddo
	   close(12)
      call Print_matrix('Input Matrix is', Matrix_A, NoofRows, NoofCols)     


      return
      end


C *******************************************************************  
      subroutine Read_vector(Input_B, n)
      real  Input_B(0:n-1)      
      integer    n        
      integer i

C    ....... Read vector from input file .......
	   open(unit=13, file = './data/vdatjac.inp')
          read(13,*) NoofRows
          read(13,*) (Input_B(i), i=0,NoofRows-1)
	   close(13)
      call Print_vector('Input Vector is', Input_B,  NoofRows)

      return
      end   


C *******************************************************************  
      subroutine Print_matrix(title, Matrix_A, Nrows, Ncols)
      character *15    title
      real    Matrix_A(0:Nrows-1, 0:Ncols-1)        
      integer Nrows, Ncols        
      integer i, j

      print *, title
      do 100 i = 0  ,Nrows-1
          print *,( Matrix_A(i,j), j = 0, Ncols-1)
 100  continue
      return
      end

C***********************************************************************
      subroutine Print_vector(title, x, n)
      character *15  title   
      real  x(0:n-1)      
      integer    n        
      integer i

      print *, title
      print *, (x(i), i = 0, n-1)
      return
      end   


