*     -*-Fortran-*- Set Emacs Fortran mode explicitly
***********************************************************************
*     
*     $Id: driver.cpp,v 1.7 1995/05/16 16:15:09 dweber Exp dweber $
*     /*
*     $Log: driver.cpp,v $
# Revision 1.7  1995/05/16  16:15:09  dweber
# Final revision for project completion.  Notes:
# 	#ifdef CMF is back to #ifdef CMFORTRAN
#      No equivalent to CMF_descriptors.  Am looking into pointers.
#      Method #6 (and 16) appear to produce incorrect results.
#
# Revision 1.6  1995/04/18  01:02:09  dweber
# Isolated CMF & HPF specific features with #ifdef statements.
# Tried to update all HPF directives to be HPF specific - not CMF holdovers.
# Deleted the CMF select(var) - case(val) statements, and am using
#         if-then-elseif statements.
# Created architecture types to match the architecture.  This will need
#         some work as HPF has no concept similar to the CMF_architecture
#         functions.
# Since HPF doesn't define a type like CMF_sizeof_descriptor, I had to
#         hardcode the array layouts for HPF.  This may change.
#
# Revision 1.5  1995/03/03  06:05:52  jtowns
# Run all 6 algorithms.
# Run the left preconditioned algorithms.
# Only run 3 times each in the interest of time.
#     
# Revision 1.4  1995/02/22  22:59:32  jtowns
# Forgot to remove and enddo...
#     
# Revision 1.3  1995/02/22  22:48:47  jtowns
# Let the include file set iprecon and remove if-block for calling cmstabdrv.
#     
# Revision 1.2  1995/02/22  22:20:07  jtowns
# Fix includes to be directives handled by new Makefile with
# preprocessing.
#     
# Revision 1.1  1995/02/11  20:23:15  jtowns
# Initial revision
#     
*     */
*     
***********************************************************************
*     
      subroutine driver (ac,an,as,ae,aw,at,ab,u,rhs,
     .     im,ilower,iupper,jm,jlower,jupper,km,klower,kupper)
*     
***********************************************************************
*     
      implicit none

      integer i

#include "driver.h"
#ifdef CMFORTRAN
#include "cm/CMF_defs.h"
#endif

      integer im,iupper,ilower,jm,jupper,jlower,km,klower,kupper
#ifdef CMFORTRAN
      real, array(im,jm,km) :: ac,an,as,ae,aw,at,ab,u,rhs
cmf$ layout ac(:news,:news,:news)
cmf$ layout an(:news,:news,:news),as(:news,:news,:news)
cmf$ layout aw(:news,:news,:news),ae(:news,:news,:news)
cmf$ layout at(:news,:news,:news),ab(:news,:news,:news)
cmf$ layout u(:news,:news,:news),rhs(:news,:news,:news)
#elif HPF
      real, DIMENSION(im,jm,km) :: ac,an,as,ae,aw,at,ab,u,rhs
CHPF$ DISTRIBUTE ac(BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE an(BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE as(BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE aw(BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ae(BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE at(BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ab(BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE u(BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE rhs(BLOCK, BLOCK, BLOCK)
#endif
      
#ifdef CMFORTRAN
      real, array(7,iupper-ilower+1,jupper-jlower+1,kupper-klower+1)::
     .     a,q,p
cmf$ layout a(:serial,:news,:news,:news)
cmf$ layout q(:serial,:news,:news,:news)
cmf$ layout p(:serial,:news,:news,:news)
#elif HPF
      real, DIMENSION(7,iupper-ilower+1,jupper-jlower+1,
     .                kupper-klower+1)::a,q,p
CHPF$ DISTRIBUTE a(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE q(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE p(*, BLOCK, BLOCK, BLOCK)
#endif

#ifdef CMFORTRAN
      integer, array(3,1,1,1) :: ia,iq,ip
cmf$ layout ia(:serial,:news,:news,:news)
cmf$ layout iq(:serial,:news,:news,:news)
cmf$ layout ip(:serial,:news,:news,:news)
#elif HPF
      integer, DIMENSION(3,1,1,1) :: ia,iq,ip
CHPF$ DISTRIBUTE ia(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE iq(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE ip(*, BLOCK, BLOCK, BLOCK)
#endif

#ifdef CMFORTRAN
      real, array(1,iupper-ilower+1,jupper-jlower+1,kupper-klower+1)::
     .     x,b
cmf$ layout x(:serial,:news,:news,:news)
cmf$ layout b(:serial,:news,:news,:news)
#elif HPF
      real, DIMENSION(1,iupper-ilower+1,jupper-jlower+1,
     .     kupper-klower+1)::x,b
CHPF$ DISTRIBUTE x(*, BLOCK, BLOCK, BLOCK)
CHPF$ DISTRIBUTE b(*, BLOCK, BLOCK, BLOCK)
#endif

      integer iparm(20), ier
      real rparm(20)
#ifdef CMFORTRAN
cmf$ layout iparm(:serial)
cmf$ layout rparm(:serial)
#elif HPF
CHPF$ DISTRIBUTE iparm(*)
CHPF$ DISTRIBUTE rparm(*)
#endif

      external usrmv,neumann_left,preconr,split_left,split_right

      integer iscale,order,m,mm
      integer arch
      real flops,mflops,pe_peak,peak,efficiency,etime

      real matvec_ops,preconl_ops,preconr_ops

      character out*15
*      character*10 out(6)
*      data out(1) /"out.32"/,
*     .     out(2) /"out.64"/,
*     .     out(3) /"out.128"/,
*     .     out(4) /"out.256"/,
*     .     out(5) /"out.512"/

*     
***********************************************************************
*     
      
*     Find out what type of CM we are on
#ifdef CMFORTRAN
      cm5vu = CMF_cm5_vu
      cm5sparc = CMF_cm5_sparc
      cm2slice = CMF_cm2_slicewise
      cm2paris = CMF_cm2_paris
      cm200sl = CMF_cm200_slicewise
      cm200ps = CMF_cm200_paris
      arch = CMF_architecture()
      nodes = CMF_number_of_processors()
#elif HPF
      cm5vu = 0
      cm5sparc = 0
      cm2slice = 0
      cm2paris = 0
      cm200sl = 0
      cm200ps = 0
      arch = hpfarch
      nodes = NUMBER_OF_PROCESSORS()
#endif


*     Open the output file, name based on conditions
#ifdef CMFORTRAN
      call outname(out, method, iprecon, nodes/4)
#elif HPF
      call outname(out, method, iprecon, nodes)
#endif
      open(unit=10, file=out, status="UNKNOWN", position="APPEND")

      write(10,*)
      write(10,*)'>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
      write(10,*)'idim = ',im,'     jdim = ',jm,'     kdim = ',km
      write(10,*)'ilower = ',ilower,'     iupper = ',iupper
      write(10,*)'jlower = ',jlower,'     jupper = ',jupper
      write(10,*)'klower = ',klower,'     kupper = ',kupper
      call flush(10)

*     Set up the coefficient array a
      a(1,:,:,:)=ac(ilower:iupper,jlower:jupper,klower:kupper)
      a(2,:,:,:)=an(ilower:iupper,jlower:jupper,klower:kupper)
      a(3,:,:,:)=as(ilower:iupper,jlower:jupper,klower:kupper)
      a(4,:,:,:)=ae(ilower:iupper,jlower:jupper,klower:kupper)
      a(5,:,:,:)=aw(ilower:iupper,jlower:jupper,klower:kupper)
      a(6,:,:,:)=at(ilower:iupper,jlower:jupper,klower:kupper)
      a(7,:,:,:)=ab(ilower:iupper,jlower:jupper,klower:kupper)

      ia(1,1,1,1)=iupper-ilower+1
      ia(2,1,1,1)=jupper-jlower+1
      ia(3,1,1,1)=kupper-klower+1

*     Set up initial guess x
      x(1,:,:,:)=u(ilower:iupper,jlower:jupper,klower:kupper)

*     Set up right hand side b
      b(1,:,:,:)=rhs(ilower:iupper,jlower:jupper,klower:kupper)

*     Find out what type of CM we are on
      write(10,*)

      if (arch .eq. cm5vu) then
         write(10,*)'Architecture          : CM-5 with Vector Units'
         write(10,*)
         write(10,*)'Number of nodes       : ',nodes/4
         write(10,*)'Number of Vector Units: ',nodes
         pe_peak = 32.
         peak= pe_peak*nodes

      elseif (arch .eq. cm5sparc) then
         write(10,*)'Architecture    : CM-5 with SPARC-only'
         write(10,*)
         write(10,*)'Number of nodes : ',nodes
         pe_peak = 5.
         peak= pe_peak*nodes

      elseif (arch .eq. cm2slice) then
         write(10,*)
     .        'Architecture                  : CM-2 Slicewise Model'
         write(10,*)
         write(10,*)'Number of Floating Point Units: ',nodes
         pe_peak = 7.
         peak= pe_peak*nodes

      elseif (arch .eq. cm2paris) then
         write(10,*)'Architecture                   : CM-2 Paris Model'
         write(10,*)
         write(10,*)'Number of bit-serial processors: ',nodes
         pe_peak = 7.
         peak= pe_peak*nodes

      elseif (arch .eq. cm200sl) then
         write(10,*)'Architecture    : CM-200 Slicewise Model'
         write(10,*)
         write(10,*)'Number of Floating Point Units: ',nodes
         pe_peak = 7.
         peak= pe_peak*nodes

      elseif (arch .eq. cm200ps) then
         write(10,*)'Architecture    : CM-200 Paris Model'
         write(10,*)
         write(10,*)'Number of bit-serial processors: ',nodes
         pe_peak = 7.
         peak= pe_peak*nodes

      elseif (arch .eq. hpfarch) then
         write(10,*)'Architecture    : HPF Processors'
         write(10,*)
         write(10,*)'Number of processors: ', nodes
         pe_peak = 300.
         peak = pe_peak*nodes

      else
         write(10,*)'Using unsupported machine type'

      endif

      write(10,*)
      write(10,*)

      write(10,1235)'method','iters','Elapsed','Busy',
     .     'FLOPS','MFLOPS','Eff %'
 1234 format(x,i6,i7,f10.2,f9.2,f8.1,f9.1,f8.1)
 1235 format(x,a6,a7,a10  ,a9  ,a8  ,a9  ,a8)

      call flush(10)

      write(10,*)

*     Apply diagonal scaling preconditioner if we can
      iscale=1
      if (any(a(1,:,:,:).eq.0.)) iscale=0
      if (iprecon.eq.3) iscale=0
      if (iscale.eq.1) then
         b(1,:,:,:)=b(1,:,:,:)/a(1,:,:,:)
         a(2,:,:,:)=a(2,:,:,:)/a(1,:,:,:)
         a(3,:,:,:)=a(3,:,:,:)/a(1,:,:,:)
         a(4,:,:,:)=a(4,:,:,:)/a(1,:,:,:)
         a(5,:,:,:)=a(5,:,:,:)/a(1,:,:,:)
         a(6,:,:,:)=a(6,:,:,:)/a(1,:,:,:)
         a(7,:,:,:)=a(7,:,:,:)/a(1,:,:,:)
         a(1,:,:,:)=1.
      endif

      q= a
      iq= ia
      p= a
      ip= ia

      do i = 1, tstitr

         x=0.

*     Get default rparm and iparm values, then set what we want
         call dfault (rparm,iparm)

         iparm(1) = ntest
         iparm(2) = (im-2)*(jm-2)*(km-2)/2
         iparm(2) = 1500
         iparm(4) = level
         iparm(5) = iounit
         iparm(6) = iprecon
         iparm(7) = method
         iparm(8) = irstrt
         iparm(9) = maxrstrt
         iparm(11) = timer

         rparm(1) = zeta

*     Open the output file
         open (unit=iounit,file='cmstab.out',status='UNKNOWN')

*     Call the CMStab solver
         if (iprecon.ne.3) then
            call cmstabdrv (usrmv,neumann_left,preconr,a,ia,x,b,
     .           q,iq,p,ip,rparm,iparm,ier)
*     else
*     call cmstabdrv (usrmv,split_left,split_right,a,ia,x,b,
*     .                q,iq,p,ip,rparm,iparm,ier)
         endif

         close (unit=iounit)

*     Calculate the number of operations and other neat stuff
*     This is correct for the provided routines using any of the 
*     methods with either no preconditioning or left preconditioning
         etime= rparm(4)
         iters= iparm(3)
         matvec_ops= 13
         if (iprecon .eq. 0) then
            preconl_ops= 0
            preconr_ops= 0
         elseif (iprecon .eq. 1) then
            preconl_ops= 12.*3.
            preconr_ops= 0
         elseif (iprecon .eq. 2) then
            preconl_ops= 0
            preconr_ops= 0
         endif
      
         order = (iupper-ilower+1)*(jupper-jlower+1)*
     .        (kupper-klower+1)

         call flop_count(mflops,flops,peak,efficiency,
     .        order,method,ntest,iprecon,iters,etime,
     .        matvec_ops,preconl_ops,preconr_ops)
*     
************************************************************************

         if (ier.ne.0) then
            write(10,1234)method+10*iprecon,iparm(3),
     .           rparm(3),rparm(4),
     .           flops*1.e-9,mflops,efficiency
            write(10,*) ' Encountered Error:',ier
            call flush(10)
            exit
         endif

 1111    write(10,1234)method+10*iprecon,iparm(3),rparm(3),
     .        rparm(4),
     .        flops*1.e-9,mflops,efficiency

         call flush(10)
      enddo

      u(ilower:iupper,jlower:jupper,klower:kupper)=x(1,:,:,:)

      return
      end
*     
***********************************************************************
*     



*     
***********************************************************************
*     
*     outname
*     generates an output file name.  The old method named the file
* by the number of nodes used.  Here the name is based on the nodes
* used and the method chosen.
*
      subroutine outname(fname, method, iprecon, nodes)
      implicit none
      character fname*(*)
      integer method, iprecon, nodes
      integer pcmeth
      character nds*3, meth*2

      pcmeth = 10 * iprecon + method
      if (pcmeth .lt. 10) then
         write(meth, 'A1,I1') "0", pcmeth
      else
         write(meth, 'I2') pcmeth
      endif

      if (nodes .lt. 10) then
         write(nds, 'A2,I1') "00", nodes
      elseif (nodes .lt. 100) then
         write(nds, 'A1,I2') "0", nodes
      else
         write(nds, 'I3') nodes
      endif

      write(fname, 'A4,A3,A1,A2') "out.", nds, ".", meth

      return
      end

