Actual source code: ex16f90.F90
slepc-3.11.2 2019-07-30
1: !
2: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
4: ! Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain
5: !
6: ! This file is part of SLEPc.
7: ! SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: !
10: ! Program usage: mpiexec -n <np> ./ex16f90 [-help] [-n <n>] [-m <m>] [SLEPc opts]
11: !
12: ! Description: Simple example that solves a quadratic eigensystem with the
13: ! PEP object. This is the Fortran90 equivalent to ex16.c
14: !
15: ! The command line options are:
16: ! -n <n>, where <n> = number of grid subdivisions in x dimension
17: ! -m <m>, where <m> = number of grid subdivisions in y dimension
18: !
19: ! ----------------------------------------------------------------------
20: !
21: program main
22: #include <slepc/finclude/slepcpep.h>
23: use slepcpep
24: implicit none
26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
27: ! Declarations
28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29: !
30: ! Variables:
31: ! M,C,K problem matrices
32: ! pep polynomial eigenproblem solver context
34: Mat M, C, K, A(3)
35: PEP pep
36: PEPType tname
37: PetscInt N, nx, ny, i, j, Istart, Iend, II
38: PetscInt nev, ithree
39: PetscMPIInt rank
40: PetscErrorCode ierr
41: PetscBool flg, terse
42: PetscScalar mone, two, four, val
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: ! Beginning of program
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
49: if (ierr .ne. 0) then
50: print*,'SlepcInitialize failed'
51: stop
52: endif
53: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr);CHKERRA(ierr)
54: nx = 10
55: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',nx,flg,ierr);CHKERRA(ierr)
56: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',ny,flg,ierr);CHKERRA(ierr)
57: if (.not. flg) then
58: ny = nx
59: endif
60: N = nx*ny
61: if (rank .eq. 0) then
62: write(*,100) N, nx, ny
63: endif
64: 100 format (/'Quadratic Eigenproblem, N=',I6,' (',I4,'x',I4,' grid)')
66: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: ! Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: ! ** K is the 2-D Laplacian
71: call MatCreate(PETSC_COMM_WORLD,K,ierr);CHKERRA(ierr)
72: call MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr);CHKERRA(ierr)
73: call MatSetFromOptions(K,ierr);CHKERRA(ierr)
74: call MatSetUp(K,ierr);CHKERRA(ierr)
75: call MatGetOwnershipRange(K,Istart,Iend,ierr);CHKERRA(ierr)
76: mone = -1.0
77: four = 4.0
78: do II=Istart,Iend-1
79: i = II/nx
80: j = II-i*nx
81: if (i .gt. 0) then
82: call MatSetValue(K,II,II-nx,mone,INSERT_VALUES,ierr);CHKERRA(ierr)
83: endif
84: if (i .lt. ny-1) then
85: call MatSetValue(K,II,II+nx,mone,INSERT_VALUES,ierr);CHKERRA(ierr)
86: endif
87: if (j .gt. 0) then
88: call MatSetValue(K,II,II-1,mone,INSERT_VALUES,ierr);CHKERRA(ierr)
89: endif
90: if (j .lt. nx-1) then
91: call MatSetValue(K,II,II+1,mone,INSERT_VALUES,ierr);CHKERRA(ierr)
92: endif
93: call MatSetValue(K,II,II,four,INSERT_VALUES,ierr);CHKERRA(ierr)
94: end do
95: call MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
96: call MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
98: ! ** C is the 1-D Laplacian on horizontal lines
99: call MatCreate(PETSC_COMM_WORLD,C,ierr);CHKERRA(ierr)
100: call MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr);CHKERRA(ierr)
101: call MatSetFromOptions(C,ierr);CHKERRA(ierr)
102: call MatSetUp(C,ierr);CHKERRA(ierr)
103: call MatGetOwnershipRange(C,Istart,Iend,ierr);CHKERRA(ierr)
104: two = 2.0
105: do II=Istart,Iend-1
106: i = II/nx
107: j = II-i*nx
108: if (j .gt. 0) then
109: call MatSetValue(C,II,II-1,mone,INSERT_VALUES,ierr);CHKERRA(ierr)
110: endif
111: if (j .lt. nx-1) then
112: call MatSetValue(C,II,II+1,mone,INSERT_VALUES,ierr);CHKERRA(ierr)
113: endif
114: call MatSetValue(C,II,II,two,INSERT_VALUES,ierr);CHKERRA(ierr)
115: end do
116: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
117: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
119: ! ** M is a diagonal matrix
120: call MatCreate(PETSC_COMM_WORLD,M,ierr);CHKERRA(ierr)
121: call MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr);CHKERRA(ierr)
122: call MatSetFromOptions(M,ierr);CHKERRA(ierr)
123: call MatSetUp(M,ierr);CHKERRA(ierr)
124: call MatGetOwnershipRange(M,Istart,Iend,ierr);CHKERRA(ierr)
125: do II=Istart,Iend-1
126: val = II+1
127: call MatSetValue(M,II,II,val,INSERT_VALUES,ierr);CHKERRA(ierr)
128: end do
129: call MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
130: call MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
132: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: ! Create the eigensolver and set various options
134: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: ! ** Create eigensolver context
137: call PEPCreate(PETSC_COMM_WORLD,pep,ierr);CHKERRA(ierr)
139: ! ** Set matrices and problem type
140: A(1) = K
141: A(2) = C
142: A(3) = M
143: ithree = 3
144: call PEPSetOperators(pep,ithree,A,ierr);CHKERRA(ierr)
145: call PEPSetProblemType(pep,PEP_GENERAL,ierr);CHKERRA(ierr)
147: ! ** Set solver parameters at runtime
148: call PEPSetFromOptions(pep,ierr);CHKERRA(ierr)
150: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: ! Solve the eigensystem
152: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: call PEPSolve(pep,ierr);CHKERRA(ierr)
156: ! ** Optional: Get some information from the solver and display it
157: call PEPGetType(pep,tname,ierr);CHKERRA(ierr)
158: if (rank .eq. 0) then
159: write(*,120) tname
160: endif
161: 120 format (' Solution method: ',A)
162: call PEPGetDimensions(pep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
163: if (rank .eq. 0) then
164: write(*,130) nev
165: endif
166: 130 format (' Number of requested eigenvalues:',I4)
168: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169: ! Display solution and clean up
170: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: ! ** show detailed info unless -terse option is given by user
173: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr);CHKERRA(ierr)
174: if (terse) then
175: call PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_NULL_VIEWER,ierr);CHKERRA(ierr)
176: else
177: call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr);CHKERRA(ierr)
178: call PEPReasonView(pep,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
179: call PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
180: call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
181: endif
182: call PEPDestroy(pep,ierr);CHKERRA(ierr)
183: call MatDestroy(K,ierr);CHKERRA(ierr)
184: call MatDestroy(C,ierr);CHKERRA(ierr)
185: call MatDestroy(M,ierr);CHKERRA(ierr)
186: call SlepcFinalize(ierr)
187: end
189: !/*TEST
190: !
191: ! test:
192: ! suffix: 1
193: ! args: -pep_nev 4 -pep_ncv 19 -terse
194: ! requires: !complex
195: !
196: !TEST*/