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