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 11 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(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,eight,one,one,PETSC_NULL_INTEGER,da,ierr);CHKERRA(ierr) 20 call DMSetFromOptions(da,ierr) 21 call DMSetUp(da,ierr) 22 call DMCreateGlobalVector(da,x,ierr);CHKERRA(ierr) 23 call VecDuplicate(x,f,ierr);CHKERRA(ierr) 24 call DMSetMatType(da,MATAIJ,ierr);CHKERRA(ierr) 25 call DMCreateMatrix(da,J,ierr);CHKERRA(ierr) 26 27 call ComputeRHS(da,f,ierr);CHKERRA(ierr) 28 call ComputeMatrix(da,J,ierr);CHKERRA(ierr) 29 30 call KSPCreate(MPI_COMM_WORLD,ksp,ierr);CHKERRA(ierr) 31 call KSPSetOperators(ksp,J,J,ierr);CHKERRA(ierr) 32 call KSPSetFromOptions(ksp,ierr);CHKERRA(ierr) 33 call KSPSolve(ksp,f,x,ierr);CHKERRA(ierr) 34 35 call MatDestroy(J,ierr);CHKERRA(ierr) 36 call VecDestroy(x,ierr);CHKERRA(ierr) 37 call VecDestroy(f,ierr);CHKERRA(ierr) 38 call KSPDestroy(ksp,ierr);CHKERRA(ierr) 39 call DMDestroy(da,ierr);CHKERRA(ierr) 40 call PetscFinalize(ierr) 41 end 42 43! AVX512 crashes without this.. 44 block data init 45 implicit none 46 PetscScalar sd 47 common /cb/ sd 48 data sd /0/ 49 end 50 subroutine knl_workarround(xx) 51 implicit none 52 PetscScalar xx 53 PetscScalar sd 54 common /cb/ sd 55 sd = sd+xx 56 end 57 58 subroutine ComputeRHS(da,x,ierr) 59 use petscdmda 60 implicit none 61 DM da 62 Vec x 63 PetscErrorCode ierr 64 PetscInt xs,xm,i,mx 65 PetscScalar hx 66 PetscScalar, pointer :: xx(:) 67 call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & 68 & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & 69 & PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr) 70 call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr) 71 hx = 1.0_PETSC_REAL_KIND/(mx-1) 72 call DMDAVecGetArrayF90(da,x,xx,ierr);CHKERRQ(ierr) 73 do i=xs,xs+xm-1 74 call knl_workarround(xx(i)) 75 xx(i) = i*hx 76 enddo 77 call DMDAVecRestoreArrayF90(da,x,xx,ierr);CHKERRQ(ierr) 78 return 79 end 80 81 subroutine ComputeMatrix(da,J,ierr) 82 use petscdm 83 use petscmat 84 implicit none 85 Mat J 86 DM da 87 PetscErrorCode ierr 88 PetscInt xs,xm,i,mx 89 PetscScalar hx,one 90 91 one = 1.0 92 call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & 93 & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & 94 & PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr) 95 call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr) 96 hx = 1.0_PETSC_REAL_KIND/(mx-1) 97 do i=xs,xs+xm-1 98 if ((i .eq. 0) .or. (i .eq. mx-1)) then 99 call MatSetValue(J,i,i,one,INSERT_VALUES,ierr);CHKERRQ(ierr) 100 else 101 call MatSetValue(J,i,i-1,-hx,INSERT_VALUES,ierr);CHKERRQ(ierr) 102 call MatSetValue(J,i,i+1,-hx,INSERT_VALUES,ierr);CHKERRQ(ierr) 103 call MatSetValue(J,i,i,2*hx,INSERT_VALUES,ierr);CHKERRQ(ierr) 104 endif 105 enddo 106 call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr) 107 call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr) 108 return 109 end 110 111!/*TEST 112! 113! test: 114! args: -ksp_converged_reason 115! 116!TEST*/ 117