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