*******************************************************************************
* 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 'az_aztecf.h'
       
       integer   n 
       common    /global/ n
C 
C
       integer i,j,IAM,NPROCS,ierr,nz,K,in,ik,nz1,nzz,no
       integer proc_config(0:AZ_PROC_SIZE), options(0:AZ_OPTIONS_SIZE)
       double precision params(0:AZ_PARAMS_SIZE)
       double precision status(0:AZ_STATUS_SIZE)
       double precision s,time,total
       integer N_update,ierror,recb(0:63),disp(0:63)
       CHARACTER*8 met(0:4)

       integer, dimension(:), allocatable:: data_org, update
       integer, dimension(:), allocatable:: external, update_index
       integer, dimension(:), allocatable:: extern_index, bindx       
       real(KIND=8), dimension(:), allocatable:: b,x,tmp,val

 
       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.'//)                                                    
    
       open(10, FILE='size.dat', FORM='FORMATTED')
       read(10,*) 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
       nzz = nz/NPROCS + 1


       allocate( data_org(0:nzz) )
       allocate( update(0:nzz) ) 
       allocate( external(0:nzz) )
       allocate( update_index(0:nzz) )
       allocate( extern_index(0:nzz) )
       allocate( bindx(0:nzz*7) )
       allocate( b(0:nzz) )
       allocate( x(0:nz) )
       allocate( tmp(0:nzz) )
       allocate( val(0:nzz*7) )

       total = 0.0 
        DO 777 K = 0,4
        if(k.eq.1) goto 777
       time = MPI_Wtime()

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
c       options(AZ_precond) = AZ_dom_decomp
c       options(AZ_subdomain_solve) = AZ_ilut
       options(4) = 1
       options(AZ_output) = AZ_none
       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 
       time = MPI_Wtime() - time
       total = total + time
       if(iam.eq.0) write(*,435) time
 435   FORMAT(2x,'TIME CALCULATION (SEC.) = ',F12.4)
       CALL FLUSH(6)
 777   continue
       if(iam.eq.0) write(*,436) total
 436   FORMAT(/,2x,' TOTAL TIME CALCULATION (SEC.) = ',F12.4)       
       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,n3
       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