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: return
74: end
76: subroutine ComputeMatrix(da,J,ierr)
77: use petscdm
78: use petscmat
79: implicit none
80: Mat J
81: DM da
82: PetscErrorCode ierr
83: PetscInt xs,xm,i,mx
84: PetscScalar hx,one
86: one = 1.0
87: 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))
88: PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
89: hx = 1.0_PETSC_REAL_KIND/(mx-1)
90: do i=xs,xs+xm-1
91: if ((i .eq. 0) .or. (i .eq. mx-1)) then
92: PetscCall(MatSetValue(J,i,i,one,INSERT_VALUES,ierr))
93: else
94: PetscCall(MatSetValue(J,i,i-1,-hx,INSERT_VALUES,ierr))
95: PetscCall(MatSetValue(J,i,i+1,-hx,INSERT_VALUES,ierr))
96: PetscCall(MatSetValue(J,i,i,2*hx,INSERT_VALUES,ierr))
97: endif
98: enddo
99: PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr))
100: PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr))
101: return
102: end
104: !/*TEST
105: !
106: ! test:
107: ! args: -ksp_converged_reason
108: !
109: !TEST*/