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: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
16: if (ierr .ne. 0) then
17: print*,'Unable to initialize PETSc'
18: stop
19: endif
20: call DMDACreate1d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,eight,one,one,PETSC_NULL_INTEGER,da,ierr);CHKERRA(ierr)
21: call DMSetFromOptions(da,ierr)
22: call DMSetUp(da,ierr)
23: call DMCreateGlobalVector(da,x,ierr);CHKERRA(ierr)
24: call VecDuplicate(x,f,ierr);CHKERRA(ierr)
25: call DMSetMatType(da,MATAIJ,ierr);CHKERRA(ierr)
26: call DMCreateMatrix(da,J,ierr);CHKERRA(ierr)
28: call ComputeRHS(da,f,ierr);CHKERRA(ierr)
29: call ComputeMatrix(da,J,ierr);CHKERRA(ierr)
31: call KSPCreate(MPI_COMM_WORLD,ksp,ierr);CHKERRA(ierr)
32: call KSPSetOperators(ksp,J,J,ierr);CHKERRA(ierr)
33: call KSPSetFromOptions(ksp,ierr);CHKERRA(ierr)
34: call KSPSolve(ksp,f,x,ierr);CHKERRA(ierr)
36: call MatDestroy(J,ierr);CHKERRA(ierr)
37: call VecDestroy(x,ierr);CHKERRA(ierr)
38: call VecDestroy(f,ierr);CHKERRA(ierr)
39: call KSPDestroy(ksp,ierr);CHKERRA(ierr)
40: call DMDestroy(da,ierr);CHKERRA(ierr)
41: call PetscFinalize(ierr)
42: end
44: ! AVX512 crashes without this..
45: block data init
46: implicit none
47: PetscScalar sd
48: common /cb/ sd
49: data sd /0/
50: end
51: subroutine knl_workaround(xx)
52: implicit none
53: PetscScalar xx
54: PetscScalar sd
55: common /cb/ sd
56: sd = sd+xx
57: end
59: subroutine ComputeRHS(da,x,ierr)
60: use petscdmda
61: 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: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
69: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
70: & PETSC_NULL_INTEGER,ierr);PetscCall(ierr)
71: call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);PetscCall(ierr)
72: hx = 1.0_PETSC_REAL_KIND/(mx-1)
73: call DMDAVecGetArrayF90(da,x,xx,ierr);PetscCall(ierr)
74: do i=xs,xs+xm-1
75: call knl_workaround(xx(i))
76: xx(i) = i*hx
77: enddo
78: call DMDAVecRestoreArrayF90(da,x,xx,ierr);PetscCall(ierr)
79: return
80: end
82: subroutine ComputeMatrix(da,J,ierr)
83: use petscdm
84: use petscmat
85: implicit none
86: Mat J
87: DM da
88: PetscErrorCode ierr
89: PetscInt xs,xm,i,mx
90: PetscScalar hx,one
92: one = 1.0
93: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
94: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
95: & PETSC_NULL_INTEGER,ierr);PetscCall(ierr)
96: call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);PetscCall(ierr)
97: hx = 1.0_PETSC_REAL_KIND/(mx-1)
98: do i=xs,xs+xm-1
99: if ((i .eq. 0) .or. (i .eq. mx-1)) then
100: call MatSetValue(J,i,i,one,INSERT_VALUES,ierr);PetscCall(ierr)
101: else
102: call MatSetValue(J,i,i-1,-hx,INSERT_VALUES,ierr);PetscCall(ierr)
103: call MatSetValue(J,i,i+1,-hx,INSERT_VALUES,ierr);PetscCall(ierr)
104: call MatSetValue(J,i,i,2*hx,INSERT_VALUES,ierr);PetscCall(ierr)
105: endif
106: enddo
107: call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr);PetscCall(ierr)
108: call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr);PetscCall(ierr)
109: return
110: end
112: !/*TEST
113: !
114: ! test:
115: ! args: -ksp_converged_reason
116: !
117: !TEST*/