Actual source code: ex44f.F90

      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_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,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_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,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*/