*******************************************************************************
* Sample driver for AZTEC package.  The software is tested by setting up a    *
* system of equations and a right hand side and then solving the system of    *
* equations using AZTECs iterative solvers.                                   *
*                                                                             *
* This program create an MSR matrix corresponding to a 7pt discrete           * 
* approximation to the 3D Poisson operator on an n x n x n square.            *
* n must be <= 50.                                                            *
*******************************************************************************
       program main
       include 'mpif.h'
       include 'aztecfn.h'
       integer   n
       common    /global/ n
C 
       double precision b(0:125000),x(0:125000),tmp(0:125000)
C
       integer    i
       external function add_row_7pt
C
       integer proc_config(0:AZ_PROC_SIZE), options(0:AZ_OPTIONS_SIZE)
       double precision params(0:AZ_PARAMS_SIZE)
       integer data_org(0:125000)
       double precision status(0:AZ_STATUS_SIZE)
       integer update(0:125000), external(0:125000)
       integer update_index(0:125000), extern_index(0:125000)
       integer bindx(0:875000)
       double  precision val(0:875000),s
       integer N_update,ierror,recb(0:63),disp(0:63)
       CHARACTER*8 met(0:4)
C 
       CALL MPI_Init(ierror)      
C
       met(0) = 'cg'
       met(1) = 'gmres'
       met(2) = 'cgs'
       met(3) = 'tfqmr'
       met(4) = 'bicgstab'
C
C      get number of processors and the name of this processor
C
       call AZ_processor_info(proc_config)
       IAM = proc_config(AZ_node)
       NPROCS = proc_config (AZ_N_procs)
C
       if (IAM.eq.0) then
      write(6,1)
    1 format(' '//
     $       '  Victor N. Datsyuk'/
     $       '  High Performance Computers Department'/
     $       '  Rostov State University'/
     $       '  Rostov-on-Don, phone 28-54-66'//
     $       '  email: root@rsuss1.rnd.runnet.ru'/
     $       '  Internet: http://rsuss1.rnd.runnet.ru'/)
      write(6,2)
    2 format(' '/
     $ '  Aztec example driver'//
     $ '  This program create an MSR matrix corresponding to'/
     $ '  a 7pt discrete approximation to the 3D Poisson operator'/ 
     $ '  on an n x n x n square and solve it various methods.'/    
     $ '  n must be <= 50.'//)                                                    

       write (6,*) 'input number grid point on each direction n = '
       read(*,*) n
       endif
c       call AZ_broadcast(n,4,proc_config,AZ_PACK)
c       call AZ_broadcast(NULL,0,proc_config,AZ_SEND)
c
       call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
       nz = n*n*n 
        DO 777 K = 0,4
C
C      Define paritioning:matrix rows in ascending order assigned
C      to this node
C
       call AZ_read_update(N_update,update,proc_config,nz,1,0)
C
C      create the matrix: each processor creates only rows
C      appearing in update[] (using global col. numbers).
C
       bindx(0) = N_update+1
       do 250 i = 0, N_update-1
       call add_row_7pt(update(i),i,val,bindx)
250    continue
c
c      create the right hand side of the system  under rule:
c      b(i) = sum(A(i,j)*j), asuming that x(i) = i ;
c       
       do 341 i = 0, N_update-1
          in = bindx(i)
          ik = bindx(i+1) - 1
          s = val(i)*(update(i)+1)
c          s = val(i)
       do 342 j = in,ik
          s = s + val(j)*(bindx(j)+1)
c          s = s + val(j)
 342   continue
          tmp(i) = s      
 341   continue
c  
C      convert matrix to a local distributed matrix */
C
       call AZ_transform(proc_config,external,bindx,val,update,
     $                   update_index,extern_index,data_org,
     $                   N_update,NULL,NULL,NULL,NULL,AZ_MSR_MATRIX)
C

C      initialize AZTEC options
C
       call AZ_defaults(options, params)
       options(AZ_solver) = K
       options(4) = 1
       options(AZ_output) = AZ_last
       options(AZ_max_iter) = 2000
       params(AZ_tol) = 1.d-12

       if (IAM.eq.0) then
       write(6,780)
  780  format (16x,11(2H**))
       write (6,*) '               Method solve is =',met(k)
       write(6,780)
       write (6,*) '               Dimension matrices =',nz,
     +'   NPROCS = ', NPROCS
       endif       
C
C      Set rhs (delta function at grid center) and initialize guess
C
       do 350 i = 0, N_update-1
          x(update_index(i)) = 0.0
          b(update_index(i)) = tmp(i)    
 350    continue    
C
C      solve the system of equations using b  as the right hand side
C
       call AZ_solve(x,b, options, params, NULL,bindx,NULL,NULL,
     $               NULL,val, data_org, status, proc_config)
c
      DO 415 I = 0,N_update - 1
 415  tmp(i) = x(update_index(i))
c
c     gather all parts of solve on all processors
c
      CALL MPI_ALLgather(N_update,1,MPI_INTEGER,recb,1,
     *MPI_INTEGER,MPI_COMM_WORLD,ierror)
      CALL MPI_BARRIER(MPI_COMM_WORLD,ierror)
      DISP(0) = 0
      DO 404 I = 1,NPROCS-1
 404  DISP(I) = DISP(I-1) + RECB(I-1)
      CALL MPI_ALLgatherv(tmp,N_update,MPI_DOUBLE_PRECISION,x,recb,
     *disp,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,ierror)
      CALL MPI_BARRIER(MPI_COMM_WORLD,ierror)
       nz1 = nz-1
       if (IAM.eq.0) then
       do 348 i = 0,nz1,nz1
         s = x(i) 
       write (6,956) i,s           
 348    continue
  956  format (2x,'i=',I12,2x,'x(i)=',F18.8)
       endif 
 777   continue
       CALL MPI_FINALIZE(IERROR)
       stop
       end

C*********************************************************************
C*********************************************************************
C
       subroutine add_row_7pt(row,location,val,bindx)
C
       integer row,location,bindx(0:*)
       double precision val(0:*)
       integer   n
       common    /global/ n
C
C
C  Parameters:
C     row          == global row number of the new row to be added.
C     location     == local row where diagonal of the new row will be stored.
C     val,bindx    == (see user's guide). On output, val[] and bindx[]
C                     are appended such that the new row has been added.
C
       integer k
C
C      check neighbors in each direction and add nonzero if neighbor exits
C
       n3 = n*n*n - 1
       k = bindx(location)
       bindx(k)  = row + 1
       if (bindx(k).le.n3) then
          val(k) = -1.
          k = k + 1
       endif
       bindx(k)  = row - 1
       if (bindx(k) .ge. 0) then
          val(k) = -1.
           k = k + 1
       endif
       bindx(k)  = row + n
       if (bindx(k).le.n3) then
          val(k) = -1.
          k = k + 1
       endif
       bindx(k)  = row - n
       if (bindx(k).ge. 0) then
          val(k) = -1.
          k = k + 1
       endif
       bindx(k)  = row + n*n
       if (bindx(k) .le. n3) then
          val(k) = -1.
          k = k + 1
       endif
       bindx(k)  = row - n*n
       if (bindx(k) .ge. 0) then
          val(k) = -1.
          k = k + 1
       endif
       bindx(location+1) = k
       val(location)   = 6.0
       return
       end