*============================================================================*
*
*  Program Application illustrates the use of libpjac calls.
*
*  Note that in this example, half of the PEs in the current partition
*  are selected to participate in the solution phase.
*
*  All processors in the partition must call the RSPJAC routine.
*
*============================================================================*

	program Application

	implicit none

*----------------------------------------------------------------------------*
*  Allocate variables.
*----------------------------------------------------------------------------*

	include 'common.inc'

	logical*8  vop, pop

	character*1  op1, op2

	integer*8  c, n, n_c, my_pid, case, MAX, IMAX

	real*8	ZRO

	intrinsic  MY_PE

	my_pid = MY_PE()

*----------------------------------------------------------------------------*
*  Print problem sizing data.
*----------------------------------------------------------------------------*

	c = IMAX(1,N$PES/2)

	if (my_pid .eq. mgr) then

	  print *
	  print *, 'psize: ', N$PES, ' PE(s)'
	  print *, 'using: ', c,     ' PE(s)'
	  print *
	  print *, 'input: n, case, ZRO, MAX, vop, pop'
	  read  *,  n, case, ZRO, MAX, op1, op2
	  print *

	  vop = .FALSE.
	  pop = .FALSE.

	  if (op1 .eq. 'y' .or. op1 .eq. 'Y') vop = .TRUE.
	  if (op2 .eq. 'y' .or. op2 .eq. 'Y') pop = .TRUE.

	  if (case .eq. 0) case = 3

	  n_c = n/c

	endif

	call BARRIER()

	if (my_pid .eq. mgr) then

	  print *, '----------------------------------'
	  print *, 'n       = ', n
	  print *, 'case    = ', case
	  print *, 'c       = ', c
	  print *, 'n/c     = ', n_c
	  print *, 'ZRO     = ', ZRO
	  print *, 'MAX     = ', MAX
	  print *, '----------------------------------'

	else

	  call SHMEM_GET (n,    n,    1, mgr)
	  call SHMEM_GET (case, case, 1, mgr)
	  call SHMEM_GET (n_c,  n_c,  1, mgr)
	  call SHMEM_GET (ZRO,  ZRO,  1, mgr)
	  call SHMEM_GET (MAX,  MAX,  1, mgr)

	  call SHMEM_GET (vop,  vop,  1, mgr)
	  call SHMEM_GET (pop,  pop,  1, mgr)

	endif

	call Application_Slave (n, case, n_c, c,
     ^                          ZRO, MAX, vop, pop)
	end

*============================================================================*
*  Subroutine Application_Slave allocates the local arrays on all PEs,
*  then opens the global data file and reads the proper matrix section
*  into the local array.
*============================================================================*

	subroutine Application_Slave (n, case, n_c, c,
     ^                                ZRO, MAX, vop, pop)
	implicit none

*----------------------------------------------------------------------------*
*  Allocate local variables.
*----------------------------------------------------------------------------*

	include 'common.inc'

	logical*8  vop, pop

	integer*8  my_pid, n, case, n_c, c, MAX,
     ^		   i, ier, err, abort

	real*8  SECONDR, timer(0:31), ZRO, A, d

	pointer (A_ptr,A)
	pointer (d_ptr,d)

	intrinsic  MY_PE

	my_pid = MY_PE()

	abort = 0

	pSync = -1

*----------------------------------------------------------------------------*
*  Generate a dataset using all PEs in the partition.
*----------------------------------------------------------------------------*

	timer(0) = SECONDR()

c	if (my_pid .eq. mgr)
c     ^	  call MAKMAT (n, n, case, 'pjac.dat', 9, ier)
c
c	call BARRIER()

*----------------------------------------------------------------------------*
*  Allocate local arrays using a subset of the partition equal to N$PES/2,
*  then get local data from global data file.
*----------------------------------------------------------------------------*

	timer(1) = SECONDR()

	call SHPALLOC (A_ptr, n*n/c, err, abort)
	if (err .ne. 0) print *, 'SHPALLOC error ', err

	call SHPALLOC (d_ptr, n/c,   err, abort)
	if (err .ne. 0) print *, 'SHPALLOC error ', err

	if (my_pid .lt. c)
     ^	  call FGET (A, n*n/c, 'pjac.dat', 10, my_pid+1, ier)

	call BARRIER()

*----------------------------------------------------------------------------*
*  Call real, symmetric parallel Jacobian solver routine.  Note that all
*  PEs in the partition must call RSPJAC even though all PEs might not
*  participate in the solution phase.
*----------------------------------------------------------------------------*

	timer(2) = SECONDR()

	call RSPJAC (A, n, n, c, d, ZRO, MAX, vop, ier)

	if (ier .ne. 0) then
	  if (my_pid .eq. mgr) print *, 'ier=', ier
	endif

*----------------------------------------------------------------------------*
*  Output data then reset PE count.
*----------------------------------------------------------------------------*

	timer(3) = SECONDR()

	if (my_pid .lt. c) then

	  call FPUT (d, n/c, 'evals.dat', 11, my_pid+1, ier)
	  if (ier .ne. 0) print *, 'FPUT error ', ier

	  if (vop) then
	    call FPUT (A, n*n/c, 'evecs.dat', 12, my_pid+1, ier)
	    if (ier .ne. 0) print *, 'FPUT error ', ier
	  endif

	endif

	call BARRIER()

*----------------------------------------------------------------------------*
*  Show solution data.
*----------------------------------------------------------------------------*

	timer(4) = SECONDR()

	if (my_pid .eq. mgr .and. pop) then

	  print *, 'eigenvalues:'
	  call SHOMAT (n/N$PES, n/N$PES, N$PES, 'evals.dat', 9, ier)

	  if (vop) then
	    print *, 'eigenvectors:'
	    call SHOMAT (n, n, n, 'evecs.dat', 9, ier)
	  endif

	endif

	timer(5) = SECONDR()

*----------------------------------------------------------------------------*
*  Print timings, if indicated, and exit.
*----------------------------------------------------------------------------*

	if (my_pid .eq. mgr) call Print_Timing_Data (timer)

	end


*============================================================================*
*  Subroutine Print_Timing_Data writes a report to STDOUT.
*============================================================================*

	subroutine Print_Timing_Data (timer)

	implicit none

*----------------------------------------------------------------------------*
*  Declare variables.
*----------------------------------------------------------------------------*

	real*8  timer(0:31)
	
*----------------------------------------------------------------------------*
*  Print timing data.
*----------------------------------------------------------------------------*

	print *
	write (*,100) ' makmat        ', timer(1)-timer(0)
	write (*,100) ' fget          ', timer(2)-timer(1)
	write (*,100) ' rspjac        ', timer(3)-timer(2)
	write (*,100) ' fput          ', timer(4)-timer(3)
	write (*,100) ' shomat        ', timer(5)-timer(4)
	write (*,100) ' TOTAL         ', timer(5)-timer(0)
	print *

100	format (A16, F12.4, 's')

*----------------------------------------------------------------------------*
*  Return to caller.
*----------------------------------------------------------------------------*

	return
	end