Actual source code: ex44f.F90
1: program main ! Solves the linear system J x = f
2: #include finclude/petscdef.h
3: use petscksp; use petscda
4: Vec x,f; Mat J; DA da; KSP ksp; PetscErrorCode ierr
5: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
7: call DACreate1d(MPI_COMM_WORLD,DA_NONPERIODIC,8,1,1,PETSC_NULL_INTEGER,da,ierr)
8: call DACreateGlobalVector(da,x,ierr); call VecDuplicate(x,f,ierr)
9: call DAGetMatrix(da,MATAIJ,J,ierr)
11: call ComputeRHS(da,f,ierr)
12: call ComputeMatrix(da,J,ierr)
14: call KSPCreate(MPI_COMM_WORLD,ksp,ierr)
15: call KSPSetOperators(ksp,J,J,SAME_NONZERO_PATTERN,ierr)
16: call KSPSetFromOptions(ksp,ierr)
17: call KSPSolve(ksp,f,x,ierr)
19: call MatDestroy(J,ierr); call VecDestroy(x,ierr); call VecDestroy(f,ierr)
20: call KSPDestroy(ksp,ierr); call DADestroy(da,ierr)
21: call PetscFinalize(ierr)
22: end
23: subroutine ComputeRHS(da,x,ierr)
24: #include finclude/petscdef.h
25: use petscda
27: call DAGetInfo(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,ierr)
28: call DAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
29: hx = 1.d0/(mx-1)
30: call VecGetArrayF90(x,xx,ierr)
31: do i=xs,xs+xm-1
32: xx(i) = i*hx
33: enddo
34: call VecRestoreArrayF90(x,xx,ierr)
35: return
36: end
37: subroutine ComputeMatrix(da,J,ierr)
38: #include finclude/petscdef.h
39: use petscda
41: call DAGetInfo(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,ierr)
42: call DAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
43: hx = 1.d0/(mx-1)
44: do i=xs,xs+xm-1
45: if ((i .eq. 0) .or. (i .eq. mx-1)) then
46: call MatSetValue(J,i,i,1d0,INSERT_VALUES,ierr)
47: else
48: call MatSetValue(J,i,i-1,-hx,INSERT_VALUES,ierr)
49: call MatSetValue(J,i,i+1,-hx,INSERT_VALUES,ierr)
50: call MatSetValue(J,i,i,2*hx,INSERT_VALUES,ierr)
51: endif
52: enddo
53: call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr); call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr)
54: return
55: end