Actual source code: ex44f.F90

petsc-main 2021-04-20
Report Typos and Errors
  1:       program main              !   Solves the linear system  J x = f
  2: #include <petsc/finclude/petsc.h>
  3:       use petscmpi  ! or mpi or mpi_f08
  4:       use petscksp
  5:       implicit none
  6:       Vec x,f
  7:       Mat J
  8:       DM da
  9:       KSP ksp
 10:       PetscErrorCode ierr
 11:       PetscInt eight,one

 13:       eight = 8
 14:       one = 1
 15:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 16:       if (ierr .ne. 0) then
 17:         print*,'Unable to initialize PETSc'
 18:         stop
 19:       endif
 20:       call DMDACreate1d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,eight,one,one,PETSC_NULL_INTEGER,da,ierr);CHKERRA(ierr)
 21:       call DMSetFromOptions(da,ierr)
 22:       call DMSetUp(da,ierr)
 23:       call DMCreateGlobalVector(da,x,ierr);CHKERRA(ierr)
 24:       call VecDuplicate(x,f,ierr);CHKERRA(ierr)
 25:       call DMSetMatType(da,MATAIJ,ierr);CHKERRA(ierr)
 26:       call DMCreateMatrix(da,J,ierr);CHKERRA(ierr)

 28:       call ComputeRHS(da,f,ierr);CHKERRA(ierr)
 29:       call ComputeMatrix(da,J,ierr);CHKERRA(ierr)

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

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

 44: ! AVX512 crashes without this..
 45:       block data init
 46:       implicit none
 47:       PetscScalar sd
 48:       common /cb/ sd
 49:       data sd /0/
 50:       end
 51:       subroutine knl_workarround(xx)
 52:       implicit none
 53:       PetscScalar xx
 54:       PetscScalar sd
 55:       common /cb/ sd
 56:       sd = sd+xx
 57:       end

 59:       subroutine  ComputeRHS(da,x,ierr)
 60:       use petscdmda
 61:       implicit none
 62:       DM da
 63:       Vec x
 64:       PetscErrorCode ierr
 65:       PetscInt xs,xm,i,mx
 66:       PetscScalar hx
 67:       PetscScalar, pointer :: xx(:)
 68:       call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
 69:      &                 PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
 70:      &                 PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
 71:       call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
 72:       hx     = 1.0_PETSC_REAL_KIND/(mx-1)
 73:       call DMDAVecGetArrayF90(da,x,xx,ierr);CHKERRQ(ierr)
 74:       do i=xs,xs+xm-1
 75:         call knl_workarround(xx(i))
 76:         xx(i) = i*hx
 77:       enddo
 78:       call DMDAVecRestoreArrayF90(da,x,xx,ierr);CHKERRQ(ierr)
 79:       return
 80:       end

 82:       subroutine ComputeMatrix(da,J,ierr)
 83:       use petscdm
 84:       use petscmat
 85:       implicit none
 86:       Mat J
 87:       DM da
 88:       PetscErrorCode ierr
 89:       PetscInt xs,xm,i,mx
 90:       PetscScalar hx,one

 92:       one = 1.0
 93:       call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
 94:      &                 PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,   &
 95:      &                 PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
 96:       call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
 97:       hx     = 1.0_PETSC_REAL_KIND/(mx-1)
 98:       do i=xs,xs+xm-1
 99:         if ((i .eq. 0) .or. (i .eq. mx-1)) then
100:           call MatSetValue(J,i,i,one,INSERT_VALUES,ierr);CHKERRQ(ierr)
101:         else
102:           call MatSetValue(J,i,i-1,-hx,INSERT_VALUES,ierr);CHKERRQ(ierr)
103:           call MatSetValue(J,i,i+1,-hx,INSERT_VALUES,ierr);CHKERRQ(ierr)
104:           call MatSetValue(J,i,i,2*hx,INSERT_VALUES,ierr);CHKERRQ(ierr)
105:         endif
106:       enddo
107:       call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
108:       call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
109:       return
110:       end

112: !/*TEST
113: !
114: !   test:
115: !      args: -ksp_converged_reason
116: !
117: !TEST*/