module petsc_stuff
use f90_kind

! Include and use statemnents for PETSC

#include "finclude/petscvecdef.h"
#include "finclude/petscmatdef.h"

use petscvec
use petscmat
implicit none

! Data type for a matrix alone

logical :: petsc_initialized = .false.     ! PetSC has been initialized

type petsc_matrix_type
  Mat :: A_PetSC                           ! The matrix
end type petsc_matrix_type

contains


subroutine init_petsc
use f90_kind
implicit none

PetscErrorCode :: ierr

call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

end subroutine init_petsc



subroutine stop_petsc
use f90_kind
implicit none

PetscErrorCode :: ierr

call PetscFinalize(ierr)

end subroutine stop_petsc



subroutine init_petsc_matrix(petsc_matrix,no_rows,no_cols,no_entries_row)
use f90_kind
implicit none

type(petsc_matrix_type), intent(inout) :: petsc_matrix
integer, intent(in) :: no_rows
integer, intent(in) :: no_cols
integer, intent(in) :: no_entries_row

integer, dimension(:), allocatable :: entries_per_rows
PetscErrorCode :: ierr

allocate(entries_per_rows(1:no_rows))
entries_per_rows = no_entries_row

call MatCreateSeqAIJ(PETSC_COMM_SELF,no_rows,no_cols, /*no_entries_row*/ 0,entries_per_rows,petsc_matrix,ierr)
call MatZeroEntries(petsc_matrix,ierr)

deallocate(entries_per_rows)

end subroutine init_petsc_matrix




subroutine finalize_petsc_matrix(petsc_matrix)
! Turns the matrix into its final form, i.e. no coefficients can be
! added after calling this subroutine.
use f90_kind
implicit none

type(petsc_matrix_type), intent(inout) :: petsc_matrix

PetscErrorCode :: ierr

call MatAssemblyBegin(petsc_matrix%A_PetSC,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyEnd(petsc_matrix%A_PetSC,MAT_FINAL_ASSEMBLY,ierr)

end subroutine finalize_petsc_matrix



subroutine get_info(petsc_matrix)
use f90_kind
implicit none

type(petsc_matrix_type), intent(inout) :: petsc_matrix
!integer :: m,n
PetscInt :: m,n

call MatGetSize(petsc_matrix%A_PetSC,m,n)
print *,'no of rows',m
print *,'no of cols',n

end subroutine get_info



subroutine addcoef_petsc_matrix(petsc_matrix,row,col,coef)
use f90_kind
implicit none

Mat :: petsc_matrix
integer, intent(in) :: row, col
real(dp), intent(in) :: coef
PetscInt :: m,n

PetscErrorCode :: ierr

call MatSetValue(petsc_matrix,row-1,col-1,coef,ADD_VALUES,ierr)

end subroutine addcoef_petsc_matrix



subroutine mat_prod(A,B,C)
! calc C = A * B
use f90_kind
implicit none

Mat :: A
Mat :: B
Mat :: C

call MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT_DOUBLE_PRECISION,C)

end subroutine mat_prod

end module petsc_stuff