Actual source code: ex16f90.F90
1: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
3: ! Copyright (c) 2002-2013, 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: ! QEP 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 <finclude/slepcqepdef.h>
34: use slepcqep
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 <finclude/petsc.h>
42: !#include <finclude/slepc.h>
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: ! Declarations
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: !
48: ! Variables:
49: ! M,C,K problem matrices
50: ! solver quadratic eigenproblem solver context
52: #if defined(PETSC_USE_FORTRAN_DATATYPES)
53: type(Mat) M, C, K
54: type(QEP) solver
55: #else
56: Mat M, C, K
57: QEP solver
58: #endif
59: QEPType tname
60: PetscReal tol
61: PetscInt N, nx, ny, i, j, Istart, Iend, II
62: PetscInt nev, maxit
63: PetscMPIInt rank
64: PetscErrorCode ierr
65: PetscBool flg
66: PetscScalar one, mone, four
68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: ! Beginning of program
70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
73: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
74: nx = 10
75: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',nx,flg,ierr)
76: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-m',ny,flg,ierr)
77: if (.not. flg) then
78: ny = nx
79: endif
80: N = nx*ny
81: if (rank .eq. 0) then
82: write(*,100) N, nx, ny
83: endif
84: 100 format (/'Quadratic Eigenproblem, N=',I6,' (',I4,'x',I4,' grid)')
86: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: ! Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
88: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: ! ** K is the 2-D Laplacian
91: call MatCreate(PETSC_COMM_WORLD,K,ierr)
92: call MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr)
93: call MatSetFromOptions(K,ierr)
94: call MatSetUp(K,ierr)
95: call MatGetOwnershipRange(K,Istart,Iend,ierr)
96: mone = -1.0
97: four = 4.0
98: do II=Istart,Iend-1
99: i = II/nx
100: j = II-i*nx
101: if (i .gt. 0) then
102: call MatSetValue(K,II,II-nx,mone,INSERT_VALUES,ierr)
103: endif
104: if (i .lt. ny-1) then
105: call MatSetValue(K,II,II+nx,mone,INSERT_VALUES,ierr)
106: endif
107: if (j .gt. 0) then
108: call MatSetValue(K,II,II-1,mone,INSERT_VALUES,ierr)
109: endif
110: if (j .lt. nx-1) then
111: call MatSetValue(K,II,II+1,mone,INSERT_VALUES,ierr)
112: endif
113: call MatSetValue(K,II,II,four,INSERT_VALUES,ierr)
114: end do
115: call MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY,ierr)
116: call MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY,ierr)
118: ! ** C is the zero matrix
119: call MatCreate(PETSC_COMM_WORLD,C,ierr)
120: call MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr)
121: call MatSetFromOptions(C,ierr)
122: call MatSetUp(C,ierr)
123: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
124: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
126: ! ** M is the identity matrix
127: call MatCreate(PETSC_COMM_WORLD,M,ierr)
128: call MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr)
129: call MatSetFromOptions(M,ierr)
130: call MatSetUp(M,ierr)
131: call MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY,ierr)
132: call MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY,ierr)
133: one = 1.0
134: call MatShift(M,one,ierr)
136: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: ! Create the eigensolver and set various options
138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: ! ** Create eigensolver context
141: call QEPCreate(PETSC_COMM_WORLD,solver,ierr)
143: ! ** Set matrices and problem type
144: call QEPSetOperators(solver,M,C,K,ierr)
145: call QEPSetProblemType(solver,QEP_GENERAL,ierr)
147: ! ** Set solver parameters at runtime
148: call QEPSetFromOptions(solver,ierr)
150: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: ! Solve the eigensystem
152: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: call QEPSolve(solver,ierr)
156: ! ** Optional: Get some information from the solver and display it
157: call QEPGetType(solver,tname,ierr)
158: if (rank .eq. 0) then
159: write(*,120) tname
160: endif
161: 120 format (' Solution method: ',A)
162: call QEPGetDimensions(solver,nev,PETSC_NULL_INTEGER, &
163: & PETSC_NULL_INTEGER,ierr)
164: if (rank .eq. 0) then
165: write(*,130) nev
166: endif
167: 130 format (' Number of requested eigenvalues:',I4)
168: call QEPGetTolerances(solver,tol,maxit,ierr)
169: if (rank .eq. 0) then
170: write(*,140) tol, maxit
171: endif
172: 140 format (' Stopping condition: tol=',1P,E10.4,', maxit=',I4)
174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: ! Display solution and clean up
176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: call QEPPrintSolution(solver,PETSC_NULL_OBJECT,ierr)
179: call QEPDestroy(solver,ierr)
180: call MatDestroy(K,ierr)
181: call MatDestroy(C,ierr)
182: call MatDestroy(M,ierr)
183: call SlepcFinalize(ierr)
184: end