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_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*/