Actual source code: ex16f.F90
1: !
2: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
4: ! Copyright (c) 2002-, 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> ./ex16f [-help] [-n <n>] [-m <m>] [SLEPc opts]
11: !
12: ! Description: Simple example that solves a quadratic eigensystem with PEP.
13: ! 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: #include <slepc/finclude/slepcpep.h>
22: program ex16f
23: use slepcpep
24: implicit none
26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
27: ! Declarations
28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30: Mat :: M, C, K, A(3) ! problem matrices
31: PEP :: pep ! polynomial eigenproblem solver context
32: PEPType :: tname
33: PetscInt :: N, nx, ny, i, j, Istart, Iend, II
34: PetscInt :: nev
35: PetscMPIInt :: rank
36: PetscErrorCode :: ierr
37: PetscBool :: flg, terse
38: PetscScalar :: val
39: PetscScalar, parameter :: two = 2.0, four = 4.0, mone = -1.0
41: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: ! Beginning of program
43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
46: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
47: nx = 10
48: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', nx, flg, ierr))
49: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', ny, flg, ierr))
50: if (.not. flg) then
51: ny = nx
52: end if
53: N = nx*ny
54: if (rank == 0) then
55: write (*, '(/a,i6,a,i4,a,i4,a)') 'Quadratic Eigenproblem, N=', N, ' (', nx, 'x', ny, ' grid)'
56: end if
58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: ! Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: ! ** K is the 2-D Laplacian
63: PetscCallA(MatCreate(PETSC_COMM_WORLD, K, ierr))
64: PetscCallA(MatSetSizes(K, PETSC_DECIDE, PETSC_DECIDE, N, N, ierr))
65: PetscCallA(MatSetFromOptions(K, ierr))
66: PetscCallA(MatGetOwnershipRange(K, Istart, Iend, ierr))
67: do II = Istart, Iend - 1
68: i = II/nx
69: j = II - i*nx
70: if (i > 0) then
71: PetscCallA(MatSetValue(K, II, II - nx, mone, INSERT_VALUES, ierr))
72: end if
73: if (i < ny - 1) then
74: PetscCallA(MatSetValue(K, II, II + nx, mone, INSERT_VALUES, ierr))
75: end if
76: if (j > 0) then
77: PetscCallA(MatSetValue(K, II, II - 1, mone, INSERT_VALUES, ierr))
78: end if
79: if (j < nx - 1) then
80: PetscCallA(MatSetValue(K, II, II + 1, mone, INSERT_VALUES, ierr))
81: end if
82: PetscCallA(MatSetValue(K, II, II, four, INSERT_VALUES, ierr))
83: end do
84: PetscCallA(MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY, ierr))
85: PetscCallA(MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY, ierr))
87: ! ** C is the 1-D Laplacian on horizontal lines
88: PetscCallA(MatCreate(PETSC_COMM_WORLD, C, ierr))
89: PetscCallA(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N, ierr))
90: PetscCallA(MatSetFromOptions(C, ierr))
91: PetscCallA(MatGetOwnershipRange(C, Istart, Iend, ierr))
92: do II = Istart, Iend - 1
93: i = II/nx
94: j = II - i*nx
95: if (j > 0) then
96: PetscCallA(MatSetValue(C, II, II - 1, mone, INSERT_VALUES, ierr))
97: end if
98: if (j < nx - 1) then
99: PetscCallA(MatSetValue(C, II, II + 1, mone, INSERT_VALUES, ierr))
100: end if
101: PetscCallA(MatSetValue(C, II, II, two, INSERT_VALUES, ierr))
102: end do
103: PetscCallA(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY, ierr))
104: PetscCallA(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY, ierr))
106: ! ** M is a diagonal matrix
107: PetscCallA(MatCreate(PETSC_COMM_WORLD, M, ierr))
108: PetscCallA(MatSetSizes(M, PETSC_DECIDE, PETSC_DECIDE, N, N, ierr))
109: PetscCallA(MatSetFromOptions(M, ierr))
110: PetscCallA(MatGetOwnershipRange(M, Istart, Iend, ierr))
111: do II = Istart, Iend - 1
112: val = II + 1
113: PetscCallA(MatSetValue(M, II, II, val, INSERT_VALUES, ierr))
114: end do
115: PetscCallA(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY, ierr))
116: PetscCallA(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY, ierr))
118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: ! Create the eigensolver and set various options
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: ! ** Create eigensolver context
123: PetscCallA(PEPCreate(PETSC_COMM_WORLD, pep, ierr))
125: ! ** Set matrices and problem type
126: A(1) = K
127: A(2) = C
128: A(3) = M
129: PetscCallA(PEPSetOperators(pep, 3_PETSC_INT_KIND, A, ierr))
130: PetscCallA(PEPSetProblemType(pep, PEP_GENERAL, ierr))
132: ! ** Set solver parameters at runtime
133: PetscCallA(PEPSetFromOptions(pep, ierr))
135: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: ! Solve the eigensystem
137: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: PetscCallA(PEPSolve(pep, ierr))
141: ! ** Optional: Get some information from the solver and display it
142: PetscCallA(PEPGetType(pep, tname, ierr))
143: if (rank == 0) then
144: write (*, '(a,a)') ' Solution method: ', tname
145: end if
146: PetscCallA(PEPGetDimensions(pep, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
147: if (rank == 0) then
148: write (*, '(a,i4)') ' Number of requested eigenvalues:', nev
149: end if
151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: ! Display solution and clean up
153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: ! ** show detailed info unless -terse option is given by user
156: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
157: if (terse) then
158: PetscCallA(PEPErrorView(pep, PEP_ERROR_BACKWARD, PETSC_NULL_VIEWER, ierr))
159: else
160: PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
161: PetscCallA(PEPConvergedReasonView(pep, PETSC_VIEWER_STDOUT_WORLD, ierr))
162: PetscCallA(PEPErrorView(pep, PEP_ERROR_BACKWARD, PETSC_VIEWER_STDOUT_WORLD, ierr))
163: PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
164: end if
165: PetscCallA(PEPDestroy(pep, ierr))
166: PetscCallA(MatDestroy(K, ierr))
167: PetscCallA(MatDestroy(C, ierr))
168: PetscCallA(MatDestroy(M, ierr))
169: PetscCallA(SlepcFinalize(ierr))
170: end program ex16f
172: !/*TEST
173: !
174: ! test:
175: ! suffix: 1
176: ! args: -pep_nev 4 -pep_ncv 19 -terse
177: ! requires: !complex
178: ! test:
179: ! suffix: 2
180: ! args: -pep_nev 4 -pep_smallest_real -terse
181: ! output_file: output/ex16f_1.out
182: ! requires: !complex
183: !
184: !TEST*/