PETSc

subroutine sparse_2_petsc( m_in , m_out )
#include <petsc/finclude/petscksp.h>
 
    use system_env
    use petscksp
    use cylindricalRZ, only: mgridcore,MaxVertexValency,&
                             grid_num_tri
    use global_parameters, only: partd_comm
 
    !-----------------------------
 
    implicit none
 
    real, intent(in) :: m_in( 9*grid_num_tri , 3 )
    Mat, intent(out) :: m_out
 
    PetscErrorCode  ierr
    PetscInt    ::  m,i1,i2,i,j,k,rstart,rend,M1,N1
    PetscScalar ::  ss
    PetscInt,pointer,dimension(:) ::  dnnz,onnz,idx
    Integer     :: nn
 
    !------------------------------------
 
    m = mgridcore
 
    allocate( dnnz(0:m-1) )
    allocate( onnz(0:m-1) )
    allocate( idx(1) )
 
    !------------------------------------
 
    do i=0,m-1
     dnnz(i) = 2
     onnz(i) = MaxVertexValency
    end do
 
    !------------------------------------
 
    call MatCreateAIJ( partd_comm , PETSC_DETERMINE , &
            PETSC_DETERMINE , m , m , PETSC_DECIDE , &
            dnnz , PETSC_DECIDE , onnz , m_out , ierr )
 
    call MatSetFromOptions( m_out , ierr )
    call MatSetUp( m_out , ierr )
 
    !-------------------------------------
 
    deallocate(dnnz)
    deallocate(onnz)
 
    !----------------------------------------------------------------------
 
    call MatSetOption( m_out , MAT_NEW_NONZERO_ALLOCATION_ERR , &
           PETSC_FALSE , ierr )
 
    call MatGetOwnershipRange( m_out , rstart , rend , ierr )
    call MatGetSize( m_out , M1 , N1 , ierr )
 
    !---------------------------------------------------------------------- 
 
                                                                 !
    nn = int( m_in(1,1) )                                        !
                                                                 !
    do i = 2 , (nn+1)                                            !
                                                                 !
       i1 = int( m_in(i,1) )                                     !
       i2 = int( m_in(i,2) )                                     !
       ss = m_in(i,3)                                            !
                                                                 !
       if( ((i1-1)>=rstart) .AND. ((i1-1)<rend) )  then          !
          call MatSetValue(m_out,i1-1,i2-1,ss,ADD_VALUES,ierr)   !
       endif                                                     !
                                                                 !
    enddo                                                        !
                                                                 !
    !-----------------------------------------------------------------------
 
    call MatAssemblyBegin(m_out,MAT_FINAL_ASSEMBLY,ierr)
    call MatAssemblyEnd(m_out,MAT_FINAL_ASSEMBLY,ierr)
 
    !--------------------------------
 
end subroutine sparse_2_petsc

!===========================================================================
 
subroutine petsc_2_solver( m_in , solver )
#include <petsc/finclude/petscksp.h>
 
   use system_env
   use petscksp
   use global_parameters, only: partd_comm
 
   !--------------------------------
 
   implicit none
 
   Mat, intent(in) :: m_in
   KSP, intent(out) :: solver
 
   PetscErrorCode  ierr
   double precision :: tol
 
   !---------------------------------
 
   call KSPCreate( partd_comm , solver , ierr )
   call KSPSetOperators( solver , m_in , m_in , ierr )
 
   tol = 0.0000001
 
   call KSPSetTolerances( solver , tol , PETSC_DEFAULT_REAL , &
       PETSC_DEFAULT_REAL , PETSC_DEFAULT_INTEGER , ierr )
 
   !----------------------------------
 
end subroutine petsc_2_solver
 
!===============================================================