Actual source code: ex16f90.F90
slepc-3.6.3 2016-03-29
1: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
3: ! Copyright (c) 2002-2015, Universitat Politecnica de Valencia, Spain
4: !
5: ! This file is part of SLEPc.
6: !
7: ! SLEPc is free software: you can redistribute it and/or modify it under the
8: ! terms of version 3 of the GNU Lesser General Public License as published by
9: ! the Free Software Foundation.
10: !
11: ! SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
12: ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13: ! FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
14: ! more details.
15: !
16: ! You should have received a copy of the GNU Lesser General Public License
17: ! along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
18: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
19: !
20: ! Program usage: mpirun -np n ex16f90 [-help] [-n <n>] [-m <m>] [SLEPc opts]
21: !
22: ! Description: Simple example that solves a quadratic eigensystem with the
23: ! PEP object. This is the Fortran90 equivalent to ex16.c
24: !
25: ! The command line options are:
26: ! -n <n>, where <n> = number of grid subdivisions in x dimension
27: ! -m <m>, where <m> = number of grid subdivisions in y dimension
28: !
29: ! ----------------------------------------------------------------------
30: !
31: program main
33: #include <slepc/finclude/slepcpepdef.h>
34: use slepcpep
36: implicit none
38: ! For usage without modules, uncomment the following lines and remove
39: ! the previous lines between 'program main' and 'implicit none'
40: !
41: !#include <petsc/finclude/petsc.h>
42: !#include <slepc/finclude/slepc.h>
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: ! Declarations
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: !
48: ! Variables:
49: ! M,C,K problem matrices
50: ! pep polynomial eigenproblem solver context
52: #if defined(PETSC_USE_FORTRAN_DATATYPES)
53: type(Mat) M, C, K, A(3)
54: type(PEP) pep
55: #else
56: Mat M, C, K, A(3)
57: PEP pep
58: #endif
59: PEPType tname
60: PetscInt N, nx, ny, i, j, Istart, Iend, II
61: PetscInt nev, ithree
62: PetscMPIInt rank
63: PetscErrorCode ierr
64: PetscBool flg, terse
65: PetscScalar one, mone, four
67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: ! Beginning of program
69: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
72: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
73: nx = 10
74: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',nx,flg,ierr)
75: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-m',ny,flg,ierr)
76: if (.not. flg) then
77: ny = nx
78: endif
79: N = nx*ny
80: if (rank .eq. 0) then
81: write(*,100) N, nx, ny
82: endif
83: 100 format (/'Quadratic Eigenproblem, N=',I6,' (',I4,'x',I4,' grid)')
85: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: ! Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
87: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: ! ** K is the 2-D Laplacian
90: call MatCreate(PETSC_COMM_WORLD,K,ierr)
91: call MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr)
92: call MatSetFromOptions(K,ierr)
93: call MatSetUp(K,ierr)
94: call MatGetOwnershipRange(K,Istart,Iend,ierr)
95: mone = -1.0
96: four = 4.0
97: do II=Istart,Iend-1
98: i = II/nx
99: j = II-i*nx
100: if (i .gt. 0) then
101: call MatSetValue(K,II,II-nx,mone,INSERT_VALUES,ierr)
102: endif
103: if (i .lt. ny-1) then
104: call MatSetValue(K,II,II+nx,mone,INSERT_VALUES,ierr)
105: endif
106: if (j .gt. 0) then
107: call MatSetValue(K,II,II-1,mone,INSERT_VALUES,ierr)
108: endif
109: if (j .lt. nx-1) then
110: call MatSetValue(K,II,II+1,mone,INSERT_VALUES,ierr)
111: endif
112: call MatSetValue(K,II,II,four,INSERT_VALUES,ierr)
113: end do
114: call MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY,ierr)
115: call MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY,ierr)
117: ! ** C is the zero matrix
118: call MatCreate(PETSC_COMM_WORLD,C,ierr)
119: call MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr)
120: call MatSetFromOptions(C,ierr)
121: call MatSetUp(C,ierr)
122: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
123: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
125: ! ** M is the identity matrix
126: call MatCreate(PETSC_COMM_WORLD,M,ierr)
127: call MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr)
128: call MatSetFromOptions(M,ierr)
129: call MatSetUp(M,ierr)
130: call MatGetOwnershipRange(M,Istart,Iend,ierr)
131: one = 1.0
132: do II=Istart,Iend-1
133: call MatSetValue(M,II,II,one,INSERT_VALUES,ierr)
134: end do
135: call MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY,ierr)
136: call MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY,ierr)
138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: ! Create the eigensolver and set various options
140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: ! ** Create eigensolver context
143: call PEPCreate(PETSC_COMM_WORLD,pep,ierr)
145: ! ** Set matrices and problem type
146: A(1) = K
147: A(2) = C
148: A(3) = M
149: ithree = 3
150: call PEPSetOperators(pep,ithree,A,ierr)
151: call PEPSetProblemType(pep,PEP_GENERAL,ierr)
153: ! ** Set solver parameters at runtime
154: call PEPSetFromOptions(pep,ierr)
156: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: ! Solve the eigensystem
158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: call PEPSolve(pep,ierr)
162: ! ** Optional: Get some information from the solver and display it
163: call PEPGetType(pep,tname,ierr)
164: if (rank .eq. 0) then
165: write(*,120) tname
166: endif
167: 120 format (' Solution method: ',A)
168: call PEPGetDimensions(pep,nev,PETSC_NULL_INTEGER, &
169: & PETSC_NULL_INTEGER,ierr)
170: if (rank .eq. 0) then
171: write(*,130) nev
172: endif
173: 130 format (' Number of requested eigenvalues:',I4)
175: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: ! Display solution and clean up
177: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: ! ** show detailed info unless -terse option is given by user
180: call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-terse',terse,ierr)
181: if (terse) then
182: call PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_NULL_OBJECT,ierr)
183: else
184: call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, &
185: & PETSC_VIEWER_ASCII_INFO_DETAIL,ierr)
186: call PEPReasonView(pep,PETSC_VIEWER_STDOUT_WORLD,ierr)
187: call PEPErrorView(pep,PEP_ERROR_BACKWARD, &
188: & PETSC_VIEWER_STDOUT_WORLD,ierr)
189: call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr)
190: endif
191: call PEPDestroy(pep,ierr)
192: call MatDestroy(K,ierr)
193: call MatDestroy(C,ierr)
194: call MatDestroy(M,ierr)
195: call SlepcFinalize(ierr)
196: end