1:       program main              !   Solves the linear system  J x = f
  2: #include <petsc/finclude/petscksp.h>
  3: #include <petsc/finclude/petscdmda.h>
  4:       use petscmpi  ! or mpi or mpi_f08
  5:       use petscksp
  6:       use petscdmda
  7:       implicit none

  9:       Vec x,f
 10:       Mat J
 11:       DM da
 12:       KSP ksp
 13:       PetscErrorCode ierr
 14:       PetscInt eight,one

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

 27:       PetscCallA(ComputeRHS(da,f,ierr))
 28:       PetscCallA(ComputeMatrix(da,J,ierr))

 30:       PetscCallA(KSPCreate(MPI_COMM_WORLD,ksp,ierr))
 31:       PetscCallA(KSPSetOperators(ksp,J,J,ierr))
 32:       PetscCallA(KSPSetFromOptions(ksp,ierr))
 33:       PetscCallA(KSPSolve(ksp,f,x,ierr))

 35:       PetscCallA(MatDestroy(J,ierr))
 36:       PetscCallA(VecDestroy(x,ierr))
 37:       PetscCallA(VecDestroy(f,ierr))
 38:       PetscCallA(KSPDestroy(ksp,ierr))
 39:       PetscCallA(DMDestroy(da,ierr))
 40:       PetscCallA(PetscFinalize(ierr))
 41:       end

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

 58:       subroutine  ComputeRHS(da,x,ierr)
 59:       use petscdmda
 60:       implicit none

 62:       DM da
 63:       Vec x
 64:       PetscErrorCode ierr
 65:       PetscInt xs,xm,i,mx
 66:       PetscScalar hx
 67:       PetscScalar, pointer :: xx(:)
 68:       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_DMBOUNDARYTYPE,PETSC_NULL_DMBOUNDARYTYPE,PETSC_NULL_DMBOUNDARYTYPE,PETSC_NULL_DMDASTENCILTYPE,ierr))
 69:       PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
 70:       hx     = 1.0_PETSC_REAL_KIND/(mx-1)
 71:       PetscCall(DMDAVecGetArray(da,x,xx,ierr))
 72:       do i=xs,xs+xm-1
 73:         call knl_workaround(xx(i))
 74:         xx(i) = i*hx
 75:       enddo
 76:       PetscCall(DMDAVecRestoreArray(da,x,xx,ierr))
 77:       end

 79:       subroutine ComputeMatrix(da,J,ierr)
 80:       use petscdmda
 81:       use petscmat
 82:       implicit none

 84:       Mat J
 85:       DM da
 86:       PetscErrorCode ierr
 87:       PetscInt xs,xm,i,mx
 88:       PetscScalar hx,one

 90:       one = 1.0
 91:       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_DMBOUNDARYTYPE,PETSC_NULL_DMBOUNDARYTYPE,PETSC_NULL_DMBOUNDARYTYPE,PETSC_NULL_DMDASTENCILTYPE,ierr))
 92:       PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
 93:       hx     = 1.0_PETSC_REAL_KIND/(mx-1)
 94:       do i=xs,xs+xm-1
 95:         if ((i .eq. 0) .or. (i .eq. mx-1)) then
 96:           PetscCall(MatSetValue(J,i,i,one,INSERT_VALUES,ierr))
 97:         else
 98:           PetscCall(MatSetValue(J,i,i-1,-hx,INSERT_VALUES,ierr))
 99:           PetscCall(MatSetValue(J,i,i+1,-hx,INSERT_VALUES,ierr))
100:           PetscCall(MatSetValue(J,i,i,2*hx,INSERT_VALUES,ierr))
101:         endif
102:       enddo
103:       PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr))
104:       PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr))
105:       end

107: !/*TEST
108: !
109: !   test:
110: !      args: -ksp_converged_reason
111: !
112: !TEST*/