1:       program main              !   Solves the linear system  J x = f
  2: #include <petsc/finclude/petsc.h>
  3:       use petscmpi  ! or mpi or mpi_f08
  4:       use petscksp
  5:       implicit none
  6:       Vec x,f
  7:       Mat J
  8:       DM da
  9:       KSP ksp
 10:       PetscErrorCode ierr
 11:       PetscInt eight,one

 13:       eight = 8
 14:       one = 1
 15:       PetscCallA(PetscInitialize(ierr))
 16:       PetscCallA(DMDACreate1d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,eight,one,one,PETSC_NULL_INTEGER,da,ierr))
 17:       PetscCallA(DMSetFromOptions(da,ierr))
 18:       PetscCallA(DMSetUp(da,ierr))
 19:       PetscCallA(DMCreateGlobalVector(da,x,ierr))
 20:       PetscCallA(VecDuplicate(x,f,ierr))
 21:       PetscCallA(DMSetMatType(da,MATAIJ,ierr))
 22:       PetscCallA(DMCreateMatrix(da,J,ierr))

 24:       PetscCallA(ComputeRHS(da,f,ierr))
 25:       PetscCallA(ComputeMatrix(da,J,ierr))

 27:       PetscCallA(KSPCreate(MPI_COMM_WORLD,ksp,ierr))
 28:       PetscCallA(KSPSetOperators(ksp,J,J,ierr))
 29:       PetscCallA(KSPSetFromOptions(ksp,ierr))
 30:       PetscCallA(KSPSolve(ksp,f,x,ierr))

 32:       PetscCallA(MatDestroy(J,ierr))
 33:       PetscCallA(VecDestroy(x,ierr))
 34:       PetscCallA(VecDestroy(f,ierr))
 35:       PetscCallA(KSPDestroy(ksp,ierr))
 36:       PetscCallA(DMDestroy(da,ierr))
 37:       PetscCallA(PetscFinalize(ierr))
 38:       end

 40: ! AVX512 crashes without this..
 41:       block data init
 42:       implicit none
 43:       PetscScalar sd
 44:       common /cb/ sd
 45:       data sd /0/
 46:       end
 47:       subroutine knl_workaround(xx)
 48:       implicit none
 49:       PetscScalar xx
 50:       PetscScalar sd
 51:       common /cb/ sd
 52:       sd = sd+xx
 53:       end

 55:       subroutine  ComputeRHS(da,x,ierr)
 56:       use petscdmda
 57:       implicit none
 58:       DM da
 59:       Vec x
 60:       PetscErrorCode ierr
 61:       PetscInt xs,xm,i,mx
 62:       PetscScalar hx
 63:       PetscScalar, pointer :: xx(:)
 64:       PetscCall(DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
 65:       PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
 66:       hx     = 1.0_PETSC_REAL_KIND/(mx-1)
 67:       PetscCall(DMDAVecGetArrayF90(da,x,xx,ierr))
 68:       do i=xs,xs+xm-1
 69:         call knl_workaround(xx(i))
 70:         xx(i) = i*hx
 71:       enddo
 72:       PetscCall(DMDAVecRestoreArrayF90(da,x,xx,ierr))
 73:       end

 75:       subroutine ComputeMatrix(da,J,ierr)
 76:       use petscdm
 77:       use petscmat
 78:       implicit none
 79:       Mat J
 80:       DM da
 81:       PetscErrorCode ierr
 82:       PetscInt xs,xm,i,mx
 83:       PetscScalar hx,one

 85:       one = 1.0
 86:       PetscCall(DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
 87:       PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
 88:       hx     = 1.0_PETSC_REAL_KIND/(mx-1)
 89:       do i=xs,xs+xm-1
 90:         if ((i .eq. 0) .or. (i .eq. mx-1)) then
 91:           PetscCall(MatSetValue(J,i,i,one,INSERT_VALUES,ierr))
 92:         else
 93:           PetscCall(MatSetValue(J,i,i-1,-hx,INSERT_VALUES,ierr))
 94:           PetscCall(MatSetValue(J,i,i+1,-hx,INSERT_VALUES,ierr))
 95:           PetscCall(MatSetValue(J,i,i,2*hx,INSERT_VALUES,ierr))
 96:         endif
 97:       enddo
 98:       PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr))
 99:       PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr))
100:       end

102: !/*TEST
103: !
104: !   test:
105: !      args: -ksp_converged_reason
106: !
107: !TEST*/