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