program main
#include <petsc/finclude/petscksp.h>
      use petscksp
      implicit none

      Vec     x,u,b
      Mat     A
      KSP    ksp
      PC     pc
      PetscInt i,j,II,JJ,m,n
      PetscInt Istart,Iend
      PetscInt one
      PetscErrorCode ierr
      PetscScalar  v

      call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
      if (ierr .ne. 0) then
        print*,'Unable to initialize PETSc'
        stop
      endif
      m      = 3
      n      = 3
      one    = 1

      call MatCreate(PETSC_COMM_WORLD,A,ierr)
      call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
      call MatSetFromOptions(A,ierr)
      call MatSetUp(A,ierr)

      call MatGetOwnershipRange(A,Istart,Iend,ierr)

      do 10, II=Istart,Iend-1
        v = -1.0
        i = II/n
        j = II - i*n
        if (i.gt.0) then
          JJ = II - n
          call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr)
        endif
        if (i.lt.m-1) then
          JJ = II + n
          call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr)
        endif
        if (j.gt.0) then
          JJ = II - 1
          call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr)
        endif
        if (j.lt.n-1) then
          JJ = II + 1
          call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr)
        endif
        v = 4.0
        call  MatSetValues(A,one,II,one,II,v,ADD_VALUES,ierr)
 10   continue

      call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
      call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

      call VecCreate(PETSC_COMM_WORLD,u,ierr)
      call VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
      call VecSetFromOptions(u,ierr)
      call VecDuplicate(u,b,ierr)
      call VecDuplicate(b,x,ierr)

      call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)

      call KSPSetOperators(ksp,A,A,ierr)

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !call KSPGetPC(ksp, pc)                                     ! This one does warn
      call KSPGetPC(ksp, pc, ierr)
      call PCSetType(pc, PCLU, ierr)
      !call PCSetType(pc, PCLU)                                   ! So does this (and it has a more-similar custom interface)
      call PCFactorSetMatSolverType(pc, "umfpack", ierr)
      !call PCFactorSetMatOrderingType(pc, MATORDERINGEXTERNAL)   ! no warning, segfault!
      call PCFactorSetMatOrderingType(pc, MATORDERINGEXTERNAL, ierr)
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      call KSPSetFromOptions(ksp,ierr)

      call KSPSolve(ksp,b,x,ierr)

      call VecDestroy(u,ierr)
      call VecDestroy(x,ierr)
      call VecDestroy(b,ierr)
      call MatDestroy(A,ierr)
      call KSPDestroy(ksp,ierr)

      call PetscFinalize(ierr)
      end