Actual source code: ex16f90.F90
1: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
3: ! Copyright (c) 2002-2011, 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
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*K+k*X+M)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 MatGetOwnershipRange(K,Istart,Iend,ierr)
94: do II=Istart,Iend-1
95: i = II/nx
96: j = II-i*nx
97: if (i .gt. 0) then
98: call MatSetValue(K,II,II-nx,-1.D0,INSERT_VALUES,ierr)
99: endif
100: if (i .lt. ny-1) then
101: call MatSetValue(K,II,II+nx,-1.D0,INSERT_VALUES,ierr)
102: endif
103: if (j .gt. 0) then
104: call MatSetValue(K,II,II-1,-1.D0,INSERT_VALUES,ierr)
105: endif
106: if (j .lt. nx-1) then
107: call MatSetValue(K,II,II+1,-1.D0,INSERT_VALUES,ierr)
108: endif
109: call MatSetValue(K,II,II,4.D0,INSERT_VALUES,ierr)
110: end do
111: call MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY,ierr)
112: call MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY,ierr)
114: ! ** C is the zero matrix
115: call MatCreate(PETSC_COMM_WORLD,C,ierr)
116: call MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr)
117: call MatSetFromOptions(C,ierr)
118: call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
119: call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
121: ! ** M is the identity matrix
122: call MatCreate(PETSC_COMM_WORLD,M,ierr)
123: call MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr)
124: call MatSetFromOptions(M,ierr)
125: call MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY,ierr)
126: call MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY,ierr)
127: call MatShift(M,1.D0,ierr)
129: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: ! Create the eigensolver and set various options
131: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: ! ** Create eigensolver context
134: call QEPCreate(PETSC_COMM_WORLD,solver,ierr)
136: ! ** Set matrices and problem type
137: call QEPSetOperators(solver,M,C,K,ierr)
138: call QEPSetProblemType(solver,QEP_GENERAL,ierr)
140: ! ** Set solver parameters at runtime
141: call QEPSetFromOptions(solver,ierr)
143: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: ! Solve the eigensystem
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: call QEPSolve(solver,ierr)
148:
149: ! ** Optional: Get some information from the solver and display it
150: call QEPGetType(solver,tname,ierr)
151: if (rank .eq. 0) then
152: write(*,120) tname
153: endif
154: 120 format (' Solution method: ',A)
155: call QEPGetDimensions(solver,nev,PETSC_NULL_INTEGER, &
156: & PETSC_NULL_INTEGER,ierr)
157: if (rank .eq. 0) then
158: write(*,130) nev
159: endif
160: 130 format (' Number of requested eigenvalues:',I4)
161: call QEPGetTolerances(solver,tol,maxit,ierr)
162: if (rank .eq. 0) then
163: write(*,140) tol, maxit
164: endif
165: 140 format (' Stopping condition: tol=',1P,E10.4,', maxit=',I4)
167: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: ! Display solution and clean up
169: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: call QEPPrintSolution(solver,PETSC_NULL_OBJECT,ierr)
172: call QEPDestroy(solver,ierr)
173: call MatDestroy(K,ierr)
174: call MatDestroy(C,ierr)
175: call MatDestroy(M,ierr)
176: call SlepcFinalize(ierr)
177: end