MODULE poisson
  IMPLICIT NONE
  INTEGER, PARAMETER:: p = SELECTED_REAL_KIND(12,50)
CONTAINS
SUBROUTINE init(a, b, f, n)                                                     
  IMPLICIT NONE
  INTEGER :: n, i, j
  REAL(KIND=p), DIMENSION(0:n+1,0:n+1), INTENT(OUT) :: &
      a, b, f                   

  !HPF$ DISTRIBUTE a(CYCLIC,*)
  !HPF$ ALIGN b(:,:) WITH a(:,:)
  !HPF$ ALIGN f(:,:) WITH a(:,:)

  a = 0.0_p
  b = 0.0_p
  f = 0.0_p
  !HPF$ INDEPENDENT, NEW(i)
  DO j = 1, n
    !HPF$ INDEPENDENT
    DO i = 1, n
      f(j,i) = EXP(-(REAL(i)/(n+1)-0.5_p)**2 - &
          (REAL(j)/(n+1)-0.5_p)**2)
    END DO
  END DO
END SUBROUTINE init
SUBROUTINE sweep(a, f, b, n)                                                    
  IMPLICIT NONE
  INTEGER :: n, i, j
  REAL(KIND=p), DIMENSION(0:n+1,0:n+1), INTENT(IN) :: a, f                       
  REAL(KIND=p), DIMENSION(0:n+1,0:n+1), INTENT(OUT) :: b                         
  REAL(KIND=p) :: h                                                              

  !HPF$ INHERIT a
  !HPF$ INHERIT b
  !HPF$ INHERIT f

  h = 1.0/REAL(n+1)
  FORALL (i = 1:n, j = 1:n)                                                     
    b(i,j) = 0.25_p*(a(i-1,j) + a(i,j+1) + a(i,j-1) + &
        a(i+1,j)) - h*h*f(i,j)            
  END FORALL
END SUBROUTINE sweep
FUNCTION diff(a, b, n) RESULT(d)                                                
  IMPLICIT NONE
  INTEGER :: n, i, j
  REAL(KIND=p), DIMENSION(0:n+1,0:n+1), INTENT(IN) :: a, b                       
  REAL(KIND=p), DIMENSION(0:n+1,0:n+1) :: c                                      
  REAL(KIND=p) :: d                                                              

  !HPF$ DISTRIBUTE a *(BLOCK,BLOCK)
  !HPF$ DISTRIBUTE b *(BLOCK,BLOCK)
  !HPF$ ALIGN c(:,:) WITH a(:,:)

  c = 0.0_p
  FORALL (i = 1:n, j = 1:n)                                                     
    c(i,j) = (a(i,j)-b(i,j))**2                                                 
  END FORALL
  d = SUM(c)                                                                    
END FUNCTION diff
END MODULE POISSON

PROGRAM poistesti
  USE poisson, ONLY : p, init, sweep, diff
  IMPLICIT NONE
  INTEGER, PARAMETER :: n = 160
  REAL(KIND=p) :: diff, dnorm = 1.0, tol = 1.0e-6
  REAL(KIND=p), DIMENSION(0:n+1,0:n+1) :: a, b, f
  INTEGER :: i = 0, maxiter = 500

  !HPF$ DISTRIBUTE a(BLOCK,BLOCK)
  !HPF$ ALIGN b(:,:) WITH a(:,:)
  !HPF$ ALIGN f(:,:) WITH a(:,:)

  CALL init(a,b,f,n)

  DO WHILE (dnorm > tol .AND. i <= maxiter)                                     
    CALL sweep(a,f,b,n)
    CALL sweep(b,f,a,n)
    dnorm = diff(a,b,n)
    IF (MOD(i,50) == 0) WRITE (*,*) 'Erotus: ', dnorm                           
    i = i + 1                                                                   
  END DO
END PROGRAM poistesti