Actual source code: ex44f.F90
petsc-3.14.6 2021-03-30
1: program main ! Solves the linear system J x = f
2: #include <petsc/finclude/petsc.h>
3: use petscksp
4: implicit none
5: Vec x,f
6: Mat J
7: DM da
8: KSP ksp
9: PetscErrorCode ierr
10: PetscInt eight,one
12: eight = 8
13: one = 1
14: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
15: if (ierr .ne. 0) then
16: print*,'Unable to initialize PETSc'
17: stop
18: endif
19: call DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,eight,one,one,PETSC_NULL_INTEGER,da,ierr);CHKERRA(ierr)
20: call DMSetFromOptions(da,ierr)
21: call DMSetUp(da,ierr)
22: call DMCreateGlobalVector(da,x,ierr);CHKERRA(ierr)
23: call VecDuplicate(x,f,ierr);CHKERRA(ierr)
24: call DMSetMatType(da,MATAIJ,ierr);CHKERRA(ierr)
25: call DMCreateMatrix(da,J,ierr);CHKERRA(ierr)
27: call ComputeRHS(da,f,ierr);CHKERRA(ierr)
28: call ComputeMatrix(da,J,ierr);CHKERRA(ierr)
30: call KSPCreate(MPI_COMM_WORLD,ksp,ierr);CHKERRA(ierr)
31: call KSPSetOperators(ksp,J,J,ierr);CHKERRA(ierr)
32: call KSPSetFromOptions(ksp,ierr);CHKERRA(ierr)
33: call KSPSolve(ksp,f,x,ierr);CHKERRA(ierr)
35: call MatDestroy(J,ierr);CHKERRA(ierr)
36: call VecDestroy(x,ierr);CHKERRA(ierr)
37: call VecDestroy(f,ierr);CHKERRA(ierr)
38: call KSPDestroy(ksp,ierr);CHKERRA(ierr)
39: call DMDestroy(da,ierr);CHKERRA(ierr)
40: call 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_workarround(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
61: DM da
62: Vec x
63: PetscErrorCode ierr
64: PetscInt xs,xm,i,mx
65: PetscScalar hx
66: PetscScalar, pointer :: xx(:)
67: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
68: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
69: & PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
70: call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
71: hx = 1.0_PETSC_REAL_KIND/(mx-1)
72: call DMDAVecGetArrayF90(da,x,xx,ierr);CHKERRQ(ierr)
73: do i=xs,xs+xm-1
74: call knl_workarround(xx(i))
75: xx(i) = i*hx
76: enddo
77: call DMDAVecRestoreArrayF90(da,x,xx,ierr);CHKERRQ(ierr)
78: return
79: end
81: subroutine ComputeMatrix(da,J,ierr)
82: use petscdm
83: use petscmat
84: implicit none
85: Mat J
86: DM da
87: PetscErrorCode ierr
88: PetscInt xs,xm,i,mx
89: PetscScalar hx,one
91: one = 1.0
92: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
93: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
94: & PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
95: call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
96: hx = 1.0_PETSC_REAL_KIND/(mx-1)
97: do i=xs,xs+xm-1
98: if ((i .eq. 0) .or. (i .eq. mx-1)) then
99: call MatSetValue(J,i,i,one,INSERT_VALUES,ierr);CHKERRQ(ierr)
100: else
101: call MatSetValue(J,i,i-1,-hx,INSERT_VALUES,ierr);CHKERRQ(ierr)
102: call MatSetValue(J,i,i+1,-hx,INSERT_VALUES,ierr);CHKERRQ(ierr)
103: call MatSetValue(J,i,i,2*hx,INSERT_VALUES,ierr);CHKERRQ(ierr)
104: endif
105: enddo
106: call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
107: call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
108: return
109: end
111: !/*TEST
112: !
113: ! test:
114: ! args: -ksp_converged_reason
115: !
116: !TEST*/