Actual source code: ex44f.F90

petsc-master 2017-04-26
Report Typos and Errors
  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(MPI_COMM_WORLD,DM_BOUNDARY_NONE,eight,one,one,PETSC_NULL_INTEGER,da,ierr);CHKERRQ(ierr)
 20:       call DMSetFromOptions(da,ierr)
 21:       call DMSetUp(da,ierr)
 22:       call DMCreateGlobalVector(da,x,ierr);CHKERRQ(ierr)
 23:       call VecDuplicate(x,f,ierr);CHKERRQ(ierr)
 24:       call DMSetMatType(da,MATAIJ,ierr);CHKERRQ(ierr)
 25:       call DMCreateMatrix(da,J,ierr);CHKERRQ(ierr)

 27:       call ComputeRHS(da,f,ierr);CHKERRQ(ierr)
 28:       call ComputeMatrix(da,J,ierr);CHKERRQ(ierr)

 30:       call KSPCreate(MPI_COMM_WORLD,ksp,ierr);CHKERRQ(ierr)
 31:       call KSPSetOperators(ksp,J,J,ierr);CHKERRQ(ierr)
 32:       call KSPSetFromOptions(ksp,ierr);CHKERRQ(ierr)
 33:       call KSPSolve(ksp,f,x,ierr);CHKERRQ(ierr)

 35:       call MatDestroy(J,ierr);CHKERRQ(ierr)
 36:       call VecDestroy(x,ierr);CHKERRQ(ierr)
 37:       call VecDestroy(f,ierr);CHKERRQ(ierr)
 38:       call KSPDestroy(ksp,ierr);CHKERRQ(ierr)
 39:       call DMDestroy(da,ierr);CHKERRQ(ierr)
 40:       call PetscFinalize(ierr);CHKERRQ(ierr)
 41:       end

 43:       subroutine  ComputeRHS(da,x,ierr)
 44:       use petscdmda
 45:       implicit none
 46:       DM da
 47:       Vec x
 48:       PetscErrorCode ierr
 49:       PetscInt xs,xm,i,mx
 50:       PetscScalar hx
 51:       PetscScalar, pointer :: xx(:)
 52:       call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
 53:      &     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
 54:      &     PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
 55:       call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
 56:       hx     = 1.0_PETSC_REAL_KIND/(mx-1)
 57:       call DMDAVecGetArrayF90(da,x,xx,ierr);CHKERRQ(ierr)
 58:       do i=xs,xs+xm-1
 59:        xx(i) = i*hx
 60:       enddo
 61:       call DMDAVecRestoreArrayF90(da,x,xx,ierr);CHKERRQ(ierr)
 62:       return
 63:       end
 64: 
 65:       subroutine ComputeMatrix(da,J,ierr)
 66:       use petscdm
 67:       use petscmat
 68:       implicit none
 69:       Mat J
 70:       DM da
 71:       PetscErrorCode ierr
 72:       PetscInt xs,xm,i,mx
 73:       PetscScalar hx,one

 75:       one = 1.0
 76:       call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
 77:      &  PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,   &
 78:      &  PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
 79:       call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
 80:       hx     = 1.0_PETSC_REAL_KIND/(mx-1)
 81:       do i=xs,xs+xm-1
 82:         if ((i .eq. 0) .or. (i .eq. mx-1)) then
 83:           call MatSetValue(J,i,i,one,INSERT_VALUES,ierr);CHKERRQ(ierr)
 84:         else
 85:           call MatSetValue(J,i,i-1,-hx,INSERT_VALUES,ierr);CHKERRQ(ierr)
 86:           call MatSetValue(J,i,i+1,-hx,INSERT_VALUES,ierr);CHKERRQ(ierr)
 87:           call MatSetValue(J,i,i,2*hx,INSERT_VALUES,ierr);CHKERRQ(ierr)
 88:         endif
 89:       enddo
 90:       call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
 91:       call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
 92:       return
 93:       end