Actual source code: ex45f.F
1: program main
2: implicit none
4: #include <finclude/petscsys.h>
5: #include <finclude/petscvec.h>
6: #include <finclude/petscmat.h>
7: #include <finclude/petscpc.h>
8: #include <finclude/petscksp.h>
9: #include <finclude/petscdm.h>
10: #include <finclude/petscdmda.h>
12: PetscInt is,js,iw,jw
13: PetscErrorCode ierr
14: KSP ksp
15: DM dm
16: external ComputeRHS,ComputeMatrix,ComputeInitialGuess
18: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
19: call KSPCreate(MPI_COMM_WORLD,ksp,ierr)
20: call DMDACreate2D(MPI_COMM_WORLD, DMDA_BOUNDARY_NONE, &
21: & DMDA_BOUNDARY_NONE, DMDA_STENCIL_STAR, 3,3, PETSC_DECIDE, &
22: & PETSC_DECIDE,1,1, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, &
23: & dm, ierr)
24: call DMSetInitialGuess(dm,ComputeInitialGuess,ierr)
25: call DMSetFunction(dm,ComputeRHS,ierr)
26: call DMSetJacobian(dm,ComputeMatrix,ierr)
27: call KSPSetDM(ksp,dm,ierr)
28: call DMDAGetCorners(dm,is,js,PETSC_NULL_INTEGER,iw,jw, &
29: & PETSC_NULL_INTEGER,ierr)
30: call KSPSetFromOptions(ksp,ierr)
31: call KSPSetUp(ksp,ierr)
32: call KSPSolve(ksp,PETSC_NULL_OBJECT,PETSC_NULL_OBJECT,ierr)
33: call KSPDestroy(ksp,ierr)
34: call DMDestroy(dm,ierr)
35: call PetscFinalize(ierr)
36: return
37: end
40: subroutine ComputeInitialGuess(dm,b,ierr)
41: implicit none
42: PetscErrorCode ierr
43: DM dm
44: Vec b
45: PetscScalar h
46: h=0.0
47: call VecSet(b,h,ierr)
48: end subroutine
50: subroutine ComputeRHS(dm,x,b,ierr)
51: implicit none
52: PetscErrorCode ierr
53: DM dm
54: Vec x,b
55: PetscScalar h
56: h=1.0
57: call VecSet(b,h,ierr)
58: end subroutine
60: subroutine ComputeMatrix(dm,x,A,B,str,ierr)
61: implicit none
62: #include <finclude/petscvec.h>
63: #include <finclude/petscmat.h>
64: PetscErrorCode ierr
65: DM dm
66: Vec x
67: Mat A,B
68: MatStructure str
69: PetscInt rstart,rend,i
70: PetscScalar h
71: h = 1.0
72: call MatGetOwnershipRange(B,rstart,rend,ierr)
73: do 10, i=rstart,rend-1
74: call MatSetValues(B,1,i,1,i,h,INSERT_VALUES,ierr)
75: 10 continue
76: call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)
77: call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)
78: if ( A .ne. B) then
79: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
80: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
81: endif
82: end subroutine