PROGRAM petscitest
      
      IMPLICIT NONE      

!  The following include statements are required for KSP Fortran programs:
!     petscsys.h       - base PETSc routines
!     petscvec.h    - vectors
!     petscmat.h    - matrices
!     petscpc.h     - preconditioners
!     petscksp.h    - Krylov subspace methods

#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscmat.h>
#include <finclude/petscpc.h>
#include <finclude/petscksp.h>

      KSP            :: solver
      Mat            :: A
      Vec            :: x,y,b
      PetscErrorCode :: ierr
      PetscMPIInt    :: rank,size
      PetscInt       :: N2
      PetscScalar    :: value
      PetscInt       :: nr,np,nt,npts,num_steps,ndiags,ai,i,j,k,M,N
      integer, dimension(:), allocatable :: veci
      real*8,  dimension(:), allocatable :: xvec
      real*8,  dimension(:), allocatable :: CSR_A
      integer, dimension(:), allocatable :: CSR_AI
      integer, dimension(:), allocatable :: CSR_AJ
      real ttmp1,ttmp2,flops
      real elapsed(2)
      real etime
     
      call PetscInitialize(PETSC_NULL_CHARACTER,ierr)      
      call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
      call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
      
      npts       = 25
      nr         = npts-2 
      np         = npts+2
      nt         = npts
      N          = 3*nr*np*nt
      IF (rank .eq. 0) THEN
          print*,'N: ',N
          print*,'cores: ',size
      END IF
      num_steps  = 2000 !Number of smvps to do in loop   
            
      value      = 1.0      
      
      print*,'MPI TEST:  My rank is:',rank
      
      !Create x, y, and b vectors:
      call VecCreate(PETSC_COMM_WORLD,x,ierr)
      call VecSetSizes(x,PETSC_DECIDE,N,ierr)
      !call VecSetType(x, VECMPI)
      call VecSetFromOptions(x,ierr)      
      
      call VecCreate(PETSC_COMM_WORLD,y,ierr)
      call VecSetSizes(y,PETSC_DECIDE,N,ierr)
      !call VecSetType(y, VECMPI)
      call VecSetFromOptions(y,ierr)    
      
      call VecCreate(PETSC_COMM_WORLD,b,ierr)
      call VecSetSizes(b,PETSC_DECIDE,N,ierr)
      !call VecSetType(b, VECMPI)
      call VecSetFromOptions(b,ierr)  
            
      !Initialize y vector to 0:      
      call VecSet(y,0.0,ierr)   
      
      !Initialize b and x vectors to (different) random values:      
      IF (rank .eq. 0) THEN
         allocate(xvec(N))
         allocate(veci(N))             
         !Create list of indices (0-based for petsc)
         DO i=0,N-1
            veci(i+1) = i
         END DO         
         call RANDOM_NUMBER(xvec)
         call VecSetValues(b,N,veci,xvec,INSERT_VALUES,ierr)         
         call RANDOM_NUMBER(xvec)         
         call VecSetValues(x,N,veci,xvec,INSERT_VALUES,ierr)                  
         deallocate(xvec)
         deallocate(veci)
      END IF
      
      !Assemble the x and b vectors across all cores:
      call VecAssemblyBegin(x,ierr);
      call VecAssemblyEnd(x,ierr);
      
      call VecAssemblyBegin(y,ierr);
      call VecAssemblyEnd(y,ierr);   
      
      call VecAssemblyBegin(b,ierr);
      call VecAssemblyEnd(b,ierr);     
            
      !call VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr)
      !call VecView(y,PETSC_VIEWER_STDOUT_WORLD,ierr)      
         
      !Create matrix:
          
      call MatCreate(PETSC_COMM_WORLD,A,ierr)      
      call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr)
      call MatSetType(A,MATMPIAIJ,ierr)        
      call MatSetFromOptions(A,ierr)
      !print*,'3nrt: ',3*nr*nt      
      i = 16
      IF(size .eq. 1) THEN
          j = 0
      ELSE
          j = 8
      END IF
          
      call MatMPIAIJSetPreallocation(A,i,PETSC_NULL_INTEGER,
     &                               j,PETSC_NULL_INTEGER,ierr)
            
      !Do not call this if using preallocation!
      !call MatSetUp(A,ierr)  
            
      call MatGetOwnershipRange(A,i,j,ierr)
      print*,'Rank ',rank,' has range ',i,' and ',j
          
      !Get MAS matrix in CSR format (random numbers for now):    
      IF (rank .eq. 0) THEN  
         call GET_RAND_MAS_MATRIX(CSR_A,CSR_AI,CSR_AJ,nr,nt,np,M) 
         print*,'Number of non-zero entries in matrix:',M 
         
         !Store matrix values one-by-one (inefficient:  better way
             !   more complicated - implement later)
       
         DO i=1,N
           !print*,'numofnonzerosinrowi:',CSR_AJ(i+1)-CSR_AJ(i)+1
            DO j=CSR_AJ(i)+1,CSR_AJ(i+1)
               call MatSetValue(A,i-1,CSR_AI(j),CSR_A(j),
     &                               INSERT_VALUES,ierr)      
            
            END DO
         END DO            
         print*,'Done setting matrix values...'     
      END IF
      
      !Assemble matrix A across all cores:
      call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
      print*,'between assembly'
      call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
          
      !call MatView(A,PETSC_VIEWER_DRAW_WORLD,ierr)
           
      
      !call csr Ax for validation/timing
      
      
      
      
     
      IF (rank .eq. 0) THEN         
         !Clean up CSR matrices:
         deallocate(CSR_A)
         deallocate(CSR_AI)
         deallocate(CSR_AJ)               
      END IF
      
           
      
      
          
      !Call matrix-vector-product y=Ax num-steps times:
      IF (rank .eq. 0) THEN
          ttmp1 = etime(elapsed)
      END IF
      DO i=1,num_steps
         call MatMult(A,x,y,ierr)
      END DO
      IF (rank .eq. 0) THEN
      ttmp2 = etime(elapsed)
      flops = ((2.0*M*num_steps)/(ttmp2-ttmp1))/1024/1024/1024 
      
      ttmp1 = 3.e9*(ttmp2-ttmp1)/(1.0*N*num_steps)
      print*,'PETSc y=Ax time:   ',ttmp1,'nsec/mp.'         
      print*,'PETSc y=Ax flops:  ',flops,'GFLOPS.'
      END IF
      
      
      !call VecView(y,PETSC_VIEWER_STDOUT_WORLD,ierr) 
          
      !Call CG with no preconditioner:    
      !Call CG with Jacobi PC:
      !Call CG with block-Jacobi PC:  
      !Call CG with ILU PC:
      !Set-up Shell Matrix:
      !Call CG?
      
      !Clean up:
      call MatDestroy(A,ierr)
      call VecDestroy(x,ierr)   
      call VecDestroy(y,ierr)
      call VecDestroy(b,ierr)
      
      call PetscFinalize(ierr)
c-----------------------------------------------------------------    
c-----------------------------------------------------------------   
c-----------------------------------------------------------------   
      contains
      
      
             !Generate MAS sparsity pattern of random numbers in DIA format:
      SUBROUTINE GET_RAND_MAS_MATRIX(CSR_A,CSR_AI,CSR_AJ,nr,nt,np,M)
         implicit none         
         
         real*8,    dimension(:,:),allocatable :: DIA_A
         integer,   dimension(:),  allocatable :: DIA_AI
         real*8,    dimension(:),  allocatable :: CSR_A
         integer,   dimension(:),  allocatable :: CSR_AI
         integer,   dimension(:),  allocatable :: CSR_AJ
         integer :: N,M,ndiags,ai,i,j,k,istart,iend,nr,nt,np  
         
         N = 3*nr*np*nt
         ndiags = 29
         allocate(DIA_AI(ndiags)) 
         !Set positions of diags:
         DIA_AI(1)  = -3*nr*nt
         DIA_AI(2)  = -3*nr*nt + 1
         DIA_AI(3)  = -3*nr*nt + 2
         DIA_AI(4)  = -3*(nr*nt-1) + 2
         DIA_AI(5)  = -3*(nr*nt-nr) + 1
         DIA_AI(6)  = -3*nr - 1
         DIA_AI(7)  = -3*nr 
         DIA_AI(8)  = -3*nr + 1
         DIA_AI(9)  = -3*(nr-1) + 1
         DIA_AI(10) = -5
         DIA_AI(11) = -4
         DIA_AI(12) = -3
         DIA_AI(13) = -2
         DIA_AI(14) = -1 
         DIA_AI(15) = 0 
         !Assume symmetric diag-structure (not values):
         DO i=(ndiags-ndiags/2),ndiags
           DIA_AI(i) = -DIA_AI(ndiags+1-i)
         END DO
         allocate(DIA_A(N,ndiags))     
         DIA_A(:,:) = -9999
         CALL RANDOM_SEED()                
         !MAS sparse pattern
         DO k=1,ndiags
            IF(DIA_AI(k) < 0) THEN 
               istart = -DIA_AI(k)+1
               iend   = N
            ELSE
                 istart = 1
                 iend   = N-DIA_AI(k)
              END IF 
              !Seven full diagonals:
              IF(k .eq. 15) THEN
                  call RANDOM_NUMBER(DIA_A(istart:iend,k))
              END IF
              
              IF((k .eq. 12) .or. (k .eq. 18)) THEN
                  call RANDOM_NUMBER(DIA_A(istart:iend,k))
              END IF
              
              IF((k .eq. 7) .or. (k .eq. 23)) THEN
                  call RANDOM_NUMBER(DIA_A(istart:iend,k))
              END IF
              
              IF((k .eq. 1) .or. (k .eq. 29)) THEN
                  call RANDOM_NUMBER(DIA_A(istart:iend,k))
              END IF 
              !Two fill-2, skip 1
              IF((k .eq. 14) .or. (k .eq. 16)) THEN    
                 i = istart
                DO
                 CALL RANDOM_NUMBER(DIA_A(i,k))
                 i = i+1
                 IF(i .gt. iend) EXIT
                 CALL RANDOM_NUMBER(DIA_A(i,k))
                 i = i+1
                 IF(i .gt. iend) EXIT
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT                 
                END DO 
              END IF
              !Six fill1-skip2
              IF((k .eq. 13) .or. (k .eq. 17)) THEN    
                 i = istart
                DO
                 CALL RANDOM_NUMBER(DIA_A(i,k))
                 i = i+1
                 IF(i .gt. iend) EXIT
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT                 
                END DO 
              END IF
              IF((k .eq. 11) .or. (k .eq. 19)) THEN    
                 i = istart
                DO
                 CALL RANDOM_NUMBER(DIA_A(i,k))
                 i = i+1
                 IF(i .gt. iend) EXIT
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT                 
                END DO 
              END IF
              IF((k .eq. 10) .or. (k .eq. 20)) THEN    
                 i = istart
                DO
                 CALL RANDOM_NUMBER(DIA_A(i,k))
                 i = i+1
                 IF(i .gt. iend) EXIT
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT                 
                END DO 
              END IF
              !Six skipinit-fill1-skip2
              IF((k .eq. 9) .or. (k .eq. 21)) THEN    
                 i = istart
                 DIA_A(i,k) = 0
                 i = i+1                 
                 DO
                   CALL RANDOM_NUMBER(DIA_A(i,k))
                   i = i+1
                  IF(i .gt. iend) EXIT
                   DIA_A(i,k) = 0
                   i = i+1
                   IF(i .gt. iend) EXIT
                   DIA_A(i,k) = 0
                   i = i+1
                   IF(i .gt. iend) EXIT                 
                END DO 
               END IF    
              IF((k .eq. 8) .or. (k .eq. 22)) THEN    
                 i = istart
                 DIA_A(i,k) = 0
                 i = i+1                 
                 DO
                   CALL RANDOM_NUMBER(DIA_A(i,k))
                   i = i+1
                   IF(i .gt. iend) EXIT
                   DIA_A(i,k) = 0
                   i = i+1
                   IF(i .gt. iend) EXIT
                   DIA_A(i,k) = 0
                   i = i+1
                   IF(i .gt. iend) EXIT                 
                 END DO 
              END IF    
              IF((k .eq. 6) .or. (k .eq. 24)) THEN    
                 i = istart
                 DIA_A(i,k) = 0
                 i = i+1                 
                 DO
                   CALL RANDOM_NUMBER(DIA_A(i,k))
                   i = i+1
                  IF(i .gt. iend) EXIT
                   DIA_A(i,k) = 0
                   i = i+1
                   IF(i .gt. iend) EXIT
                   DIA_A(i,k) = 0
                   i = i+1
                   IF(i .gt. iend) EXIT                 
                 END DO 
              END IF 
              !Eight skip2-fill1
              IF((k .eq. 5) .or. (k .eq. 25)) THEN    
                 i = istart
                DO
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT
                 CALL RANDOM_NUMBER(DIA_A(i,k))
                 i = i+1
                 IF(i .gt. iend) EXIT                 
                END DO 
              END IF    
              IF((k .eq. 4) .or. (k .eq. 26)) THEN    
                 i = istart
                DO
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT
                 CALL RANDOM_NUMBER(DIA_A(i,k))
                 i = i+1
                 IF(i .gt. iend) EXIT                 
                END DO 
              END IF    
              IF((k .eq. 3) .or. (k .eq. 27)) THEN    
                 i = istart
                DO
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT
                 CALL RANDOM_NUMBER(DIA_A(i,k))
                 i = i+1
                 IF(i .gt. iend) EXIT                 
                END DO 
              END IF    
              IF((k .eq. 2) .or. (k .eq. 28)) THEN    
                 i = istart
                DO
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT
                 DIA_A(i,k) = 0
                 i = i+1
                 IF(i .gt. iend) EXIT
                 CALL RANDOM_NUMBER(DIA_A(i,k))
                 i = i+1
                 IF(i .gt. iend) EXIT                 
                END DO 
              END IF               
          END DO    
      M = 0
      DO j=1,ndiags
         DO i=1,N
             IF((DIA_A(i,j) .ne. -9999) .and. (DIA_A(i,j) .ne. 0)) THEN
                 M = M+1
             END IF
         END DO
      END DO     
      allocate(CSR_A(M))
      allocate(CSR_AI(M))
      allocate(CSR_AJ(N+1))
      
      !Use zero-indexing for matrix indices for PETsc
      ai=1           
      CSR_AJ(1) = 0       
      DO i=1,N
         DO k=1,ndiags
            IF((DIA_A(i,k) .ne. 0) .and. (DIA_A(i,k) .ne. -9999)) THEN
               CSR_AI(ai) = i + DIA_AI(k) - 1
               CSR_A(ai)  = DIA_A(i,k)            
               ai = ai+1               
            END IF
         END DO
         CSR_AJ(i+1) = ai-1          
      END DO   
      deallocate(DIA_A,DIA_AI)  
           
      END SUBROUTINE
      END PROGRAM