! fix some PETSc versioning problems
#include "petscversion.h"

#ifdef _PETSc30ANDBEFORE
#  define _PETSc31P8ANDBEFORE
#endif

#ifdef _PETSc31P8ANDBEFORE
#  define _PETSc35BEFORE
#  define PetscBool PetscTruth
#  define MatCreateAIJ MatCreateMPIAIJ
#endif

#ifdef _PETSc35BEFORE
#  define PETSC_DEFAULT_REAL PETSC_DEFAULT_DOUBLE_PRECISION
#  define KSPSetOperators(ksp,Apetsc,Bpetsc,ierr) KSPSetOperators(ksp,Apetsc,Bpetsc,DIFFERENT_NONZERO_PATTERN,ierr)
#endif

#if PETSC_VERSION_GT(3,8,0)

module petsc_gtc
#include "petsc/finclude/petsc.h"
use precision, only: lk
use global_parameters, only: mype, gtcout, partd_comm
use petsc_array, only: newcomm, nproc_newcomm, myrank_newcomm
use petscvec
use petscmat
use petscpc
use petsc
implicit none
type petscData
  Vec ::  x,b,x_recv
  Mat :: Apetsc
  KSP ::  ksp
  PC  ::  pc
  IS  ::  is_from,is_to
  VecScatter ::    gather
end type petscData
! stored data for petsc poisson solver
type(petscData) :: petsc_ps, petsc_ps0
public petsc_init, petsc_final, &
       lapmat_pscreate_fperp, &
       lapmat_pssolver_fperp, petsc_ps, petsc_ps0
private

save

contains

!======================================================================
  
subroutine petsc_init
#include "petsc/finclude/petsc.h"
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscmat.h>
#include <petsc/finclude/petscksp.h>
use petscvec
use petscmat
use petscpc
use petsc
implicit none

  
    PetscErrorCode ierr
    integer ierror
  
  
    call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
    call MPI_Comm_dup(partd_comm,newcomm,ierror)
  
    call mpi_comm_size(newcomm,nproc_newcomm,ierror)
  
    call mpi_comm_rank(newcomm,myrank_newcomm,ierror)
  
    if(ierror /= 0) then
        write(*,*)mype,'main_petsc: newcomm error',ierror
        call MPI_ABORT(MPI_COMM_WORLD,1,ierror)
    endif
          
    return
end subroutine petsc_init
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
subroutine petsc_final(petsc_comm)

#include "petsc/finclude/petsc.h"
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscmat.h>
#include <petsc/finclude/petscksp.h>
  use petscvec
  use petscmat
  use petscpc
  use petsc
  implicit none
  type(petscData), intent(in) :: petsc_comm
  
    PetscErrorCode ierr
  
    call VecDestroy(petsc_comm%x,ierr)
    call VecDestroy(petsc_comm%b,ierr)
    call VecDestroy(petsc_comm%x_recv,ierr)
    call VecScatterDestroy(petsc_comm%gather,ierr)
    call ISDestroy(petsc_comm%is_from,ierr)
    call ISDestroy(petsc_comm%is_to,ierr)
    call MatDestroy(petsc_comm%Apetsc,ierr)
    call KSPDestroy(petsc_comm%ksp,ierr)
  
    call PetscFinalize(ierr)
  
    return
end subroutine petsc_final
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lapmat_pscreate_fperp(petsc_comm,mindexlap,matrixsize,lusera,luserp,lusers)

#include "petsc/finclude/petsc.h"
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscmat.h>
#include <petsc/finclude/petscksp.h>
    use global_parameters,only: istep
    use petscvec
    use petscmat
    use petscpc
    use petsc
    implicit none

    type(petscData), intent(in) :: petsc_comm
    integer,intent(in) :: mindexlap,matrixsize
    integer,dimension(:),intent(in)  :: luserp
    integer,dimension(-1:),intent(in)  :: lusers
    real(lk),dimension(:),intent(in) :: lusera

    double precision tol,value(mindexlap)
    PetscErrorCode ierr
    PetscInt n,i,j,ij,ilow,ihigh,nz
    PetscInt idx(1)
    PetscBool flg
    integer col(mindexlap)
  
    n=matrixsize
  
  ! PETSc 3.7 adds one additional argument at the beginning
    call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
    call VecCreate(newcomm,petsc_comm%x,ierr)
    call VecSetSizes(petsc_comm%x,PETSC_DECIDE,n,ierr)
    call VecSetFromOptions(petsc_comm%x,ierr)
    call VecDuplicate(petsc_comm%x,petsc_comm%b,ierr)
  
  !*
    call VecCreate(PETSC_COMM_SELF,petsc_comm%x_recv,ierr)
    call VecSetSizes(petsc_comm%x_recv,n,n,ierr)
    call VecSetType(petsc_comm%x_recv,VECSEQ,ierr)
    call VecSetFromOptions(petsc_comm%x_recv,ierr)
  
    call ISCreateStride(newcomm,n,0,1,petsc_comm%is_from,ierr)
    call ISCreateStride(PETSC_COMM_SELF,n,0,1,petsc_comm%is_to,ierr)
  
    call VecScatterCreate(petsc_comm%x,petsc_comm%is_from,petsc_comm%x_recv,&
                          petsc_comm%is_to,petsc_comm%gather,ierr)
  !*
  
  !     call MatCreate(newcomm,petsc_comm%Apetsc,ierr)
  !     CHKERRQ(ierr)
  !     call MatSetSizes(petsc_comm%Apetsc,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
  !     CHKERRQ(ierr)
  !     call MatSetType(petsc_comm%Apetsc,MATMPIAIJ,ierr)
  !     CHKERRQ(ierr)
  !*optimized matrix creation
    call MatCreateAIJ(newcomm,PETSC_DECIDE,PETSC_DECIDE,n,n, &
       mindexlap,PETSC_NULL_INTEGER,mindexlap,PETSC_NULL_INTEGER, &
       petsc_comm%Apetsc,ierr)
  
    call MatSetFromOptions(petsc_comm%Apetsc,ierr)
    call MatGetOwnershipRange(petsc_comm%Apetsc,ilow,ihigh,ierr)
  
    ij=lusers(ilow-1)
    do i=ilow,ihigh-1
        nz=lusers(i)-lusers(i-1)
        do j=1,nz
            ij=ij+1
            col(j) = luserp(ij)
            value(j) = lusera(ij)
        enddo
        idx(1)=i
        call MatSetValues(petsc_comm%Apetsc,1,idx,nz,col,value,INSERT_VALUES,ierr)
    enddo
  
    call MatAssemblyBegin(petsc_comm%Apetsc,MAT_FINAL_ASSEMBLY,ierr)
    call MatAssemblyEnd(petsc_comm%Apetsc,MAT_FINAL_ASSEMBLY,ierr)
  !     call MatTranspose(petsc_comm%Apetsc,Bpetsc1,ierr)
  !     call MatAXPY(Bpetsc1,1.d0,petsc_comm%Apetsc,DIFFERENT_NONZERO_PATTERN,ierr)
  !     call MatPartitioningSetAdjacency(mpart, Bpetsc1, ierr )
  !     call MatCopy(petsc_comm%Apetsc,Bpetsc1,DIFFERENT_NONZERO_PATTERN,ierr)
  
    call KSPCreate(newcomm,petsc_comm%ksp,ierr)
    call KSPSetOperators(petsc_comm%ksp,petsc_comm%Apetsc,petsc_comm%Apetsc,ierr)
  !     call KSPGetPC(petsc_comm%ksp,petsc_comm%pc,ierr)
  !     call PCSetType(petsc_comm%pc,PCJACOBI,ierr)
  ! 1.d-4 seems to be adequate
  !  tol=1.d-4
    tol=1.d-6
  !     tol=1.d-8
    call KSPSetTolerances(petsc_comm%ksp,tol,PETSC_DEFAULT_REAL, &
         PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
    call KSPSetFromOptions(petsc_comm%ksp,ierr)
  
    if(mype==0 .and. istep==1)then
        write(gtcout,*)'Finish create lap-petsc matrix'
        call flush(gtcout)
    endif
  
  
  
    return
end subroutine lapmat_pscreate_fperp
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lapmat_pssolver_fperp(petsc_comm, luserb,luserx)
#include "petsc/finclude/petsc.h"
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscmat.h>
#include <petsc/finclude/petscksp.h>
    use global_parameters,only:istep
    use petscvec
    use petscmat
    use petscpc
    use petsc
    implicit none

    type(petscData), intent(in) :: petsc_comm
    real(lk),dimension(0:),intent(in) :: luserb
    real(lk),dimension(0:),intent(inout) :: luserx

    PetscScalar temp
    PetscErrorCode ierr
    PetscInt i,ilow,ihigh,isize,ilo,ihi
  !  PetscInt isi
  !  PetscOffset idx
  !  PetscScalar vec_x(1)
    PetscScalar, pointer :: vec_x(:)
  
    call VecGetLocalSize(petsc_comm%b,isize,ierr)
    call VecGetOwnershipRange(petsc_comm%b,ilow,ihigh,ierr)
    do i=ilow,ihigh-1
        temp=luserb(i)
        call VecSetValue(petsc_comm%b,i,temp,INSERT_VALUES,ierr)
    enddo
    call VecAssemblyBegin(petsc_comm%b,ierr)
    call VecAssemblyEnd(petsc_comm%b,ierr)
  
  !!!!!!!!!!!!!!!!!!!diag petsc  !!!!!!!
  !!!!     if(istep==1 .and. mype==0)then 
  !!!   write(*,*)'checking vector b:'
  !   call VecView(petsc_comm%b,PETSC_VIEWER_STDOUT_SELF,ierr)
  !   CHKERRQ(ierr)
  !!!     endif
  !!!     if(istep==1 .and. mype==0)then
  !!!   write(*,*)'checking matrix a:'
  !!!   call MatView(petsc_comm%Apetsc,PETSC_VIEWER_STDOUT_SELF,ierr)
  !!!   CHKERRQ(ierr)
  !!     endif
  
    call KSPSetInitialGuessNonzero(petsc_comm%ksp,Petsc_True,ierr)
    call KSPSolve(petsc_comm%ksp,petsc_comm%b,petsc_comm%x,ierr)
  
  !* 
    call VecScatterBegin(petsc_comm%gather,petsc_comm%x,petsc_comm%x_recv,INSERT_VALUES, &
         SCATTER_FORWARD,ierr)
    call VecScatterEnd(petsc_comm%gather,petsc_comm%x,petsc_comm%x_recv,INSERT_VALUES, &
         SCATTER_FORWARD,ierr)
  !*
    call VecGetOwnershipRange(petsc_comm%x_recv,ilo,ihi,ierr)
  !  call VecGetLocalSize(petsc_comm%x_recv,isi,ierr)
  !  call VecGetArray(petsc_comm%x_recv,vec_x,idx,ierr)
    call VecGetArrayF90(petsc_comm%x_recv,vec_x,ierr)
    do i=ilo,ihi-1
        luserx(i)=vec_x(i-ilo+1)
    enddo
  !  call VecRestoreArray(petsc_comm%x,vec_x,idx,ierr)
    call VecRestoreArrayF90(petsc_comm%x_recv,vec_x,ierr)
  
    return
end subroutine lapmat_pssolver_fperp
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module petsc_gtc

! for petsc v3.7 and before
#else
  module petsc_gtc
    use precision, only: lk
    use global_parameters, only: mype, gtcout, partd_comm
    use petsc_array, only: newcomm, nproc_newcomm, myrank_newcomm
    implicit none
    
    public petsc_init, petsc_final, &
           lapmat_pscreate_fperp, &
           lapmat_pssolver_fperp, &
           lapmat_pscreate_fperp0, &
           lapmat_pssolver_fperp0
    private
    
    save
    
    contains
    
!======================================================================
  
subroutine petsc_init
  
#include "petsc/finclude/petsc.h90"
  
    PetscErrorCode ierr
    integer ierror
  
  
    call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
    call MPI_Comm_dup(partd_comm,newcomm,ierror)
  
    call mpi_comm_size(newcomm,nproc_newcomm,ierror)
  
    call mpi_comm_rank(newcomm,myrank_newcomm,ierror)
  
    if(ierror /= 0) then
        write(*,*)mype,'main_petsc: newcomm error',ierror
        call MPI_ABORT(MPI_COMM_WORLD,1,ierror)
    endif
          
    return
end subroutine petsc_init
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
subroutine petsc_final
  
#include "petsc/finclude/petsc.h90"
  
    Vec   x5,b5,x_recv5,x0,b0,x_recv0
    Mat   Apetsc5,Apetsc0
    KSP   ksp5,ksp0
    PC    pc5,pc0
    IS    is_from5,is_to5,is_from0,is_to0
    VecScatter     gather5,gather0
    common/petsc_common0/Apetsc0,x0,b0,ksp0,pc0,x_recv0,is_from0,is_to0,gather0 
    common/petsc_common5/Apetsc5,x5,b5,ksp5,pc5,x_recv5,is_from5,is_to5,gather5
    PetscErrorCode ierr
  
    call VecDestroy(x5,ierr)
    call VecDestroy(b5,ierr)
    call VecDestroy(x_recv5,ierr)
    call VecScatterDestroy(gather5,ierr)
    call ISDestroy(is_from5,ierr)
    call ISDestroy(is_to5,ierr)
    call MatDestroy(Apetsc5,ierr)
    call KSPDestroy(ksp5,ierr)

    call VecDestroy(x0,ierr)
    call VecDestroy(b0,ierr)
    call VecDestroy(x_recv0,ierr)
    call VecScatterDestroy(gather0,ierr)
    call ISDestroy(is_from0,ierr)
    call ISDestroy(is_to0,ierr)
    call MatDestroy(Apetsc0,ierr)
    call KSPDestroy(ksp0,ierr)

  
    call PetscFinalize(ierr)
  
    return
end subroutine petsc_final
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lapmat_pscreate_fperp(mindexlap,matrixsize,lusera,luserp,lusers)
    use global_parameters,only: istep
  
#include "petsc/finclude/petsc.h90"
  
    Vec   x,b,x_recv
    Mat   Apetsc
    KSP   ksp
    PC    pc
    IS    is_from,is_to
    VecScatter     gather
  
    common/petsc_common5/Apetsc,x,b,ksp,pc,x_recv,is_from,is_to,gather
  
    integer,intent(in) :: mindexlap,matrixsize
    integer,dimension(:),intent(in)  :: luserp
    integer,dimension(-1:),intent(in)  :: lusers
    real(lk),dimension(:),intent(in) :: lusera

    double precision tol,value(mindexlap)
    PetscErrorCode ierr
    PetscInt n,i,j,ij,ilow,ihigh,nz
    PetscInt idx(1)
    PetscBool flg
    integer col(mindexlap)
  
    n=matrixsize
  
  ! PETSc 3.7 adds one additional argument at the beginning
    call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
    call VecCreate(newcomm,x,ierr)
    call VecSetSizes(x,PETSC_DECIDE,n,ierr)
    call VecSetFromOptions(x,ierr)
    call VecDuplicate(x,b,ierr)
  
  !*
    call VecCreate(PETSC_COMM_SELF,x_recv,ierr)
    call VecSetSizes(x_recv,n,n,ierr)
    call VecSetType(x_recv,VECSEQ,ierr)
    call VecSetFromOptions(x_recv,ierr)
  
    call ISCreateStride(newcomm,n,0,1,is_from,ierr)
    call ISCreateStride(PETSC_COMM_SELF,n,0,1,is_to,ierr)
  
    call VecScatterCreate(x,is_from,x_recv,is_to,gather,ierr)
  !*
  
  !     call MatCreate(newcomm,Apetsc,ierr)
  !     CHKERRQ(ierr)
  !     call MatSetSizes(Apetsc,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
  !     CHKERRQ(ierr)
  !     call MatSetType(Apetsc,MATMPIAIJ,ierr)
  !     CHKERRQ(ierr)
  !*optimized matrix creation
    call MatCreateAIJ(newcomm,PETSC_DECIDE,PETSC_DECIDE,n,n, &
       mindexlap,PETSC_NULL_INTEGER,mindexlap,PETSC_NULL_INTEGER, &
       Apetsc,ierr)
  
    call MatSetFromOptions(Apetsc,ierr)
    call MatGetOwnershipRange(Apetsc,ilow,ihigh,ierr)
  
    ij=lusers(ilow-1)
    do i=ilow,ihigh-1
        nz=lusers(i)-lusers(i-1)
        do j=1,nz
            ij=ij+1
            col(j) = luserp(ij)
            value(j) = lusera(ij)
        enddo
        idx(1)=i
        call MatSetValues(Apetsc,1,idx,nz,col,value,INSERT_VALUES,ierr)
    enddo
  
    call MatAssemblyBegin(Apetsc,MAT_FINAL_ASSEMBLY,ierr)
    call MatAssemblyEnd(Apetsc,MAT_FINAL_ASSEMBLY,ierr)
  !     call MatTranspose(Apetsc,Bpetsc1,ierr)
  !     call MatAXPY(Bpetsc1,1.d0,Apetsc,DIFFERENT_NONZERO_PATTERN,ierr)
  !     call MatPartitioningSetAdjacency(mpart, Bpetsc1, ierr )
  !     call MatCopy(Apetsc,Bpetsc1,DIFFERENT_NONZERO_PATTERN,ierr)
  
    call KSPCreate(newcomm,ksp,ierr)
    call KSPSetOperators(ksp,Apetsc,Apetsc,ierr)
  !     call KSPGetPC(ksp,pc,ierr)
  !     call PCSetType(pc,PCJACOBI,ierr)
  ! 1.d-4 seems to be adequate
  !  tol=1.d-4
    tol=1.d-6
  !     tol=1.d-8
    call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL, &
         PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
    call KSPSetFromOptions(ksp,ierr)
  
    if(mype==0 .and. istep==1)then
        write(gtcout,*)'Finish create lap-petsc matrix'
        call flush(gtcout)
    endif
  
  
  
    return
end subroutine lapmat_pscreate_fperp
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lapmat_pssolver_fperp(luserb,luserx)
    use global_parameters,only:istep
  
#include "petsc/finclude/petsc.h90"
  
    Vec   x,b,x_recv
    Mat   Apetsc
    KSP   ksp
    PC    pc
    IS    is_from,is_to
    VecScatter     gather
  
    common/petsc_common5/Apetsc,x,b,ksp,pc,x_recv,is_from,is_to,gather
    real(lk),dimension(0:),intent(in) :: luserb
    real(lk),dimension(0:),intent(inout) :: luserx
    integer ierror

    PetscScalar temp
    PetscErrorCode ierr
    PetscInt i,ilow,ihigh,isize,ilo,ihi
  !  PetscInt isi
  !  PetscOffset idx
  !  PetscScalar vec_x(1)
    PetscScalar, pointer :: vec_x(:)
  
    call VecGetLocalSize(b,isize,ierr)
    call VecGetOwnershipRange(b,ilow,ihigh,ierr)

    do i=ilow,ihigh-1
        temp=luserb(i)
        call VecSetValue(b,i,temp,INSERT_VALUES,ierr)
    enddo

    call VecAssemblyBegin(b,ierr)
    call VecAssemblyEnd(b,ierr)

  !!!!!!!!!!!!!!!!!!!diag petsc  !!!!!!!
  !!!!     if(istep==1 .and. mype==0)then 
  !!!   write(*,*)'checking vector b:'
  !   call VecView(b,PETSC_VIEWER_STDOUT_SELF,ierr)
  !   CHKERRQ(ierr)
  !!!     endif
  !!!     if(istep==1 .and. mype==0)then
  !!!   write(*,*)'checking matrix a:'
  !!!   call MatView(Apetsc,PETSC_VIEWER_STDOUT_SELF,ierr)
  !!!   CHKERRQ(ierr)
  !!     endif
  
    call KSPSetInitialGuessNonzero(ksp,Petsc_True,ierr)
    call KSPSolve(ksp,b,x,ierr)
  !* 
    call VecScatterBegin(gather,x,x_recv,INSERT_VALUES, &
         SCATTER_FORWARD,ierr)
    call VecScatterEnd(gather,x,x_recv,INSERT_VALUES, &
         SCATTER_FORWARD,ierr)
  !*
    call VecGetOwnershipRange(x_recv,ilo,ihi,ierr)
  !  call VecGetLocalSize(x_recv,isi,ierr)
  !  call VecGetArray(x_recv,vec_x,idx,ierr)
    call VecGetArrayF90(x_recv,vec_x,ierr)
    do i=ilo,ihi-1
        luserx(i)=vec_x(i-ilo+1)
    enddo

  !  call VecRestoreArray(x,vec_x,idx,ierr)
    call VecRestoreArrayF90(x_recv,vec_x,ierr)
  
    return
end subroutine lapmat_pssolver_fperp
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lapmat_pscreate_fperp0(mindexlap,matrixsize,lusera,luserp,lusers)
    use global_parameters,only: istep
  
#include "petsc/finclude/petsc.h90"
  
    Vec   x,b,x_recv
    Mat   Apetsc
    KSP   ksp
    PC    pc
    IS    is_from,is_to
    VecScatter     gather
  
    common/petsc_common0/Apetsc,x,b,ksp,pc,x_recv,is_from,is_to,gather
  
    integer,intent(in) :: mindexlap,matrixsize
    integer,dimension(:),intent(in)  :: luserp
    integer,dimension(-1:),intent(in)  :: lusers
    real(lk),dimension(:),intent(in) :: lusera

    double precision tol,value(mindexlap)
    PetscErrorCode ierr
    PetscInt n,i,j,ij,ilow,ihigh,nz
    PetscInt idx(1)
    PetscBool flg
    integer col(mindexlap)
  
    n=matrixsize
  
  ! PETSc 3.7 adds one additional argument at the beginning
    call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
    call VecCreate(newcomm,x,ierr)
    call VecSetSizes(x,PETSC_DECIDE,n,ierr)
    call VecSetFromOptions(x,ierr)
    call VecDuplicate(x,b,ierr)
  
  !*
    call VecCreate(PETSC_COMM_SELF,x_recv,ierr)
    call VecSetSizes(x_recv,n,n,ierr)
    call VecSetType(x_recv,VECSEQ,ierr)
    call VecSetFromOptions(x_recv,ierr)
  
    call ISCreateStride(newcomm,n,0,1,is_from,ierr)
    call ISCreateStride(PETSC_COMM_SELF,n,0,1,is_to,ierr)
  
    call VecScatterCreate(x,is_from,x_recv,is_to,gather,ierr)
  !*
  
  !     call MatCreate(newcomm,Apetsc,ierr)
  !     CHKERRQ(ierr)
  !     call MatSetSizes(Apetsc,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
  !     CHKERRQ(ierr)
  !     call MatSetType(Apetsc,MATMPIAIJ,ierr)
  !     CHKERRQ(ierr)
  !*optimized matrix creation
    call MatCreateAIJ(newcomm,PETSC_DECIDE,PETSC_DECIDE,n,n, &
       mindexlap,PETSC_NULL_INTEGER,mindexlap,PETSC_NULL_INTEGER, &
       Apetsc,ierr)
  
    call MatSetFromOptions(Apetsc,ierr)
    call MatGetOwnershipRange(Apetsc,ilow,ihigh,ierr)
  
    ij=lusers(ilow-1)
    do i=ilow,ihigh-1
        nz=lusers(i)-lusers(i-1)
        do j=1,nz
            ij=ij+1
            col(j) = luserp(ij)
            value(j) = lusera(ij)
        enddo
        idx(1)=i
        call MatSetValues(Apetsc,1,idx,nz,col,value,INSERT_VALUES,ierr)
    enddo
  
    call MatAssemblyBegin(Apetsc,MAT_FINAL_ASSEMBLY,ierr)
    call MatAssemblyEnd(Apetsc,MAT_FINAL_ASSEMBLY,ierr)
  !     call MatTranspose(Apetsc,Bpetsc1,ierr)
  !     call MatAXPY(Bpetsc1,1.d0,Apetsc,DIFFERENT_NONZERO_PATTERN,ierr)
  !     call MatPartitioningSetAdjacency(mpart, Bpetsc1, ierr )
  !     call MatCopy(Apetsc,Bpetsc1,DIFFERENT_NONZERO_PATTERN,ierr)
  
    call KSPCreate(newcomm,ksp,ierr)
    call KSPSetOperators(ksp,Apetsc,Apetsc,ierr)
  !     call KSPGetPC(ksp,pc,ierr)
  !     call PCSetType(pc,PCJACOBI,ierr)
  ! 1.d-4 seems to be adequate
  !  tol=1.d-4
    tol=1.d-6
  !     tol=1.d-8
    call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL, &
         PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
    call KSPSetFromOptions(ksp,ierr)
  
    if(mype==0 .and. istep==1)then
        write(gtcout,*)'Finish create lap-petsc matrix'
        call flush(gtcout)
    endif
  
  
  
    return
end subroutine lapmat_pscreate_fperp0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lapmat_pssolver_fperp0(luserb,luserx)
    use global_parameters,only:istep
  
#include "petsc/finclude/petsc.h90"
  
    Vec   x,b,x_recv
    Mat   Apetsc
    KSP   ksp
    PC    pc
    IS    is_from,is_to
    VecScatter     gather
  
    common/petsc_common0/Apetsc,x,b,ksp,pc,x_recv,is_from,is_to,gather
    real(lk),dimension(0:),intent(in) :: luserb
    real(lk),dimension(0:),intent(inout) :: luserx

    PetscScalar temp
    PetscErrorCode ierr
    PetscInt i,ilow,ihigh,isize,ilo,ihi
  !  PetscInt isi
  !  PetscOffset idx
  !  PetscScalar vec_x(1)
    PetscScalar, pointer :: vec_x(:)
  
    call VecGetLocalSize(b,isize,ierr)
    call VecGetOwnershipRange(b,ilow,ihigh,ierr)
    do i=ilow,ihigh-1
        temp=luserb(i)
        call VecSetValue(b,i,temp,INSERT_VALUES,ierr)
    enddo
    call VecAssemblyBegin(b,ierr)
    call VecAssemblyEnd(b,ierr)
  
  !!!!!!!!!!!!!!!!!!!diag petsc  !!!!!!!
  !!!!     if(istep==1 .and. mype==0)then 
  !!!   write(*,*)'checking vector b:'
  !   call VecView(b,PETSC_VIEWER_STDOUT_SELF,ierr)
  !   CHKERRQ(ierr)
  !!!     endif
  !!!     if(istep==1 .and. mype==0)then
  !!!   write(*,*)'checking matrix a:'
  !!!   call MatView(Apetsc,PETSC_VIEWER_STDOUT_SELF,ierr)
  !!!   CHKERRQ(ierr)
  !!     endif
  
    call KSPSetInitialGuessNonzero(ksp,Petsc_True,ierr)
    call KSPSolve(ksp,b,x,ierr)
  
  !* 
    call VecScatterBegin(gather,x,x_recv,INSERT_VALUES, &
         SCATTER_FORWARD,ierr)
    call VecScatterEnd(gather,x,x_recv,INSERT_VALUES, &
         SCATTER_FORWARD,ierr)
  !*
    call VecGetOwnershipRange(x_recv,ilo,ihi,ierr)
  !  call VecGetLocalSize(x_recv,isi,ierr)
  !  call VecGetArray(x_recv,vec_x,idx,ierr)
    call VecGetArrayF90(x_recv,vec_x,ierr)
    do i=ilo,ihi-1
        luserx(i)=vec_x(i-ilo+1)
    enddo
  !  call VecRestoreArray(x,vec_x,idx,ierr)
    call VecRestoreArrayF90(x_recv,vec_x,ierr)
  
    return
end subroutine lapmat_pssolver_fperp0
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

end module petsc_gtc
#endif

