! This program is a F90/HPF code that performs the Jacobi iteration
! to solve Laplace's equation.  The square (64X64) grid
! is masked so that the boundary condition of Phi = 1.0 remains
! constant while the interior points are allowed to vary.
! The WHERE construct and the CSHIFT intrinsic handle the 
! duties of determining which points should be updated and
! the manner of updating them, in this case the Jacobi iteration.
! To see if we have the "right" answer, we test the final field
! to see if it satisfies the original equation, i.e.  if
! (delta Phi in x direction) + (delta Phi in y direction) = 0.0.
! This is done using CSHIFTS as well.  "test" is the largest
! value of the sum, and it should be small if the "right"
! answers are being output.


		Program Laplace
		Implicit none
    	Real Phiold(64,64),Phinew(64,64)  
	    Logical Hide(64,64)
    	Real Tol,average,test
    	Integer i,j
!HPF$ PROCESSORS PROC(2,2)
!HPF$ ALIGN with Phiold :: Phinew
!HPF$ DISTRIBUTE Phiold (BLOCK,BLOCK) onto PROC

	    Hide = .false.
	    Hide(2:63,2:63) = .true.
    	Phiold = 1.0

	    Where(Hide .eq. .true.)
	     	Phiold=0.0
    	End Where
	
	    Phinew=Phiold
		tol=1.0

	   DO WHILE(Tol .gt. 1.0e-5)
	    	Where(Hide .eq. .true.) 
          	  Phinew= 0.25*(CSHIFT(Phiold,dim=1,shift=1) +
     $			CSHIFT(Phiold,dim=1,shift=-1)  + 
     $ 				 CSHIFT(Phiold,dim=2,shift=1) + 
     $        CSHIFT(Phiold,dim=2,shift=-1))
		  END WHERE
		Tol=Maxval(Abs(Phinew-Phiold),mask=Hide)
		Phiold=Phinew
	  END DO
	    average=Sum(Phinew)/64.0/64.0
     	test=Maxval(Abs(CSHIFT(Phinew,dim=2,shift=1)+
     $  CSHIFT(Phinew,dim=1,shift=1) - 
     $    2.0*Phinew))
	    Write(6,*) "If test is small, we are right."
		Write(6,*) "Here is test:",test

	END