Actual source code: ex16f90.F90
slepc-3.20.1 2023-11-27
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> ./ex16f90 [-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: 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: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
49: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
50: nx = 10
51: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',nx,flg,ierr))
52: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',ny,flg,ierr))
53: if (.not. flg) then
54: ny = nx
55: endif
56: N = nx*ny
57: if (rank .eq. 0) then
58: write(*,100) N, nx, ny
59: endif
60: 100 format (/'Quadratic Eigenproblem, N=',I6,' (',I4,'x',I4,' grid)')
62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: ! Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
64: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: ! ** K is the 2-D Laplacian
67: PetscCallA(MatCreate(PETSC_COMM_WORLD,K,ierr))
68: PetscCallA(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr))
69: PetscCallA(MatSetFromOptions(K,ierr))
70: PetscCallA(MatSetUp(K,ierr))
71: PetscCallA(MatGetOwnershipRange(K,Istart,Iend,ierr))
72: mone = -1.0
73: four = 4.0
74: do II=Istart,Iend-1
75: i = II/nx
76: j = II-i*nx
77: if (i .gt. 0) then
78: PetscCallA(MatSetValue(K,II,II-nx,mone,INSERT_VALUES,ierr))
79: endif
80: if (i .lt. ny-1) then
81: PetscCallA(MatSetValue(K,II,II+nx,mone,INSERT_VALUES,ierr))
82: endif
83: if (j .gt. 0) then
84: PetscCallA(MatSetValue(K,II,II-1,mone,INSERT_VALUES,ierr))
85: endif
86: if (j .lt. nx-1) then
87: PetscCallA(MatSetValue(K,II,II+1,mone,INSERT_VALUES,ierr))
88: endif
89: PetscCallA(MatSetValue(K,II,II,four,INSERT_VALUES,ierr))
90: end do
91: PetscCallA(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY,ierr))
92: PetscCallA(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY,ierr))
94: ! ** C is the 1-D Laplacian on horizontal lines
95: PetscCallA(MatCreate(PETSC_COMM_WORLD,C,ierr))
96: PetscCallA(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr))
97: PetscCallA(MatSetFromOptions(C,ierr))
98: PetscCallA(MatSetUp(C,ierr))
99: PetscCallA(MatGetOwnershipRange(C,Istart,Iend,ierr))
100: two = 2.0
101: do II=Istart,Iend-1
102: i = II/nx
103: j = II-i*nx
104: if (j .gt. 0) then
105: PetscCallA(MatSetValue(C,II,II-1,mone,INSERT_VALUES,ierr))
106: endif
107: if (j .lt. nx-1) then
108: PetscCallA(MatSetValue(C,II,II+1,mone,INSERT_VALUES,ierr))
109: endif
110: PetscCallA(MatSetValue(C,II,II,two,INSERT_VALUES,ierr))
111: end do
112: PetscCallA(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr))
113: PetscCallA(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr))
115: ! ** M is a diagonal matrix
116: PetscCallA(MatCreate(PETSC_COMM_WORLD,M,ierr))
117: PetscCallA(MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr))
118: PetscCallA(MatSetFromOptions(M,ierr))
119: PetscCallA(MatSetUp(M,ierr))
120: PetscCallA(MatGetOwnershipRange(M,Istart,Iend,ierr))
121: do II=Istart,Iend-1
122: val = II+1
123: PetscCallA(MatSetValue(M,II,II,val,INSERT_VALUES,ierr))
124: end do
125: PetscCallA(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY,ierr))
126: PetscCallA(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY,ierr))
128: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: ! Create the eigensolver and set various options
130: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: ! ** Create eigensolver context
133: PetscCallA(PEPCreate(PETSC_COMM_WORLD,pep,ierr))
135: ! ** Set matrices and problem type
136: A(1) = K
137: A(2) = C
138: A(3) = M
139: ithree = 3
140: PetscCallA(PEPSetOperators(pep,ithree,A,ierr))
141: PetscCallA(PEPSetProblemType(pep,PEP_GENERAL,ierr))
143: ! ** Set solver parameters at runtime
144: PetscCallA(PEPSetFromOptions(pep,ierr))
146: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: ! Solve the eigensystem
148: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: PetscCallA(PEPSolve(pep,ierr))
152: ! ** Optional: Get some information from the solver and display it
153: PetscCallA(PEPGetType(pep,tname,ierr))
154: if (rank .eq. 0) then
155: write(*,120) tname
156: endif
157: 120 format (' Solution method: ',A)
158: PetscCallA(PEPGetDimensions(pep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
159: if (rank .eq. 0) then
160: write(*,130) nev
161: endif
162: 130 format (' Number of requested eigenvalues:',I4)
164: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: ! Display solution and clean up
166: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: ! ** show detailed info unless -terse option is given by user
169: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr))
170: if (terse) then
171: PetscCallA(PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_NULL_VIEWER,ierr))
172: else
173: PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr))
174: PetscCallA(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD,ierr))
175: PetscCallA(PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD,ierr))
176: PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr))
177: endif
178: PetscCallA(PEPDestroy(pep,ierr))
179: PetscCallA(MatDestroy(K,ierr))
180: PetscCallA(MatDestroy(C,ierr))
181: PetscCallA(MatDestroy(M,ierr))
182: PetscCallA(SlepcFinalize(ierr))
183: end
185: !/*TEST
186: !
187: ! test:
188: ! suffix: 1
189: ! args: -pep_nev 4 -pep_ncv 19 -terse
190: ! requires: !complex
191: !
192: !TEST*/