Actual source code: ex16f.F90

slepc-3.21.1 2024-04-26
Report Typos and Errors
```  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.
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:       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(MatGetOwnershipRange(K,Istart,Iend,ierr))
71:       mone = -1.0
72:       four = 4.0
73:       do II=Istart,Iend-1
74:         i = II/nx
75:         j = II-i*nx
76:         if (i .gt. 0) then
77:           PetscCallA(MatSetValue(K,II,II-nx,mone,INSERT_VALUES,ierr))
78:         endif
79:         if (i .lt. ny-1) then
80:           PetscCallA(MatSetValue(K,II,II+nx,mone,INSERT_VALUES,ierr))
81:         endif
82:         if (j .gt. 0) then
83:           PetscCallA(MatSetValue(K,II,II-1,mone,INSERT_VALUES,ierr))
84:         endif
85:         if (j .lt. nx-1) then
86:           PetscCallA(MatSetValue(K,II,II+1,mone,INSERT_VALUES,ierr))
87:         endif
88:         PetscCallA(MatSetValue(K,II,II,four,INSERT_VALUES,ierr))
89:       end do
90:       PetscCallA(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY,ierr))
91:       PetscCallA(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY,ierr))

93: !     ** C is the 1-D Laplacian on horizontal lines
94:       PetscCallA(MatCreate(PETSC_COMM_WORLD,C,ierr))
95:       PetscCallA(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr))
96:       PetscCallA(MatSetFromOptions(C,ierr))
97:       PetscCallA(MatGetOwnershipRange(C,Istart,Iend,ierr))
98:       two = 2.0
99:       do II=Istart,Iend-1
100:         i = II/nx
101:         j = II-i*nx
102:         if (j .gt. 0) then
103:           PetscCallA(MatSetValue(C,II,II-1,mone,INSERT_VALUES,ierr))
104:         endif
105:         if (j .lt. nx-1) then
106:           PetscCallA(MatSetValue(C,II,II+1,mone,INSERT_VALUES,ierr))
107:         endif
108:         PetscCallA(MatSetValue(C,II,II,two,INSERT_VALUES,ierr))
109:       end do
110:       PetscCallA(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr))
111:       PetscCallA(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr))

113: !     ** M is a diagonal matrix
114:       PetscCallA(MatCreate(PETSC_COMM_WORLD,M,ierr))
115:       PetscCallA(MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr))
116:       PetscCallA(MatSetFromOptions(M,ierr))
117:       PetscCallA(MatGetOwnershipRange(M,Istart,Iend,ierr))
118:       do II=Istart,Iend-1
119:         val = II+1
120:         PetscCallA(MatSetValue(M,II,II,val,INSERT_VALUES,ierr))
121:       end do
122:       PetscCallA(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY,ierr))
123:       PetscCallA(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY,ierr))

125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: !     Create the eigensolver and set various options
127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

129: !     ** Create eigensolver context
130:       PetscCallA(PEPCreate(PETSC_COMM_WORLD,pep,ierr))

132: !     ** Set matrices and problem type
133:       A(1) = K
134:       A(2) = C
135:       A(3) = M
136:       ithree = 3
137:       PetscCallA(PEPSetOperators(pep,ithree,A,ierr))
138:       PetscCallA(PEPSetProblemType(pep,PEP_GENERAL,ierr))

140: !     ** Set solver parameters at runtime
141:       PetscCallA(PEPSetFromOptions(pep,ierr))

143: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: !     Solve the eigensystem
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

147:       PetscCallA(PEPSolve(pep,ierr))

149: !     ** Optional: Get some information from the solver and display it
150:       PetscCallA(PEPGetType(pep,tname,ierr))
151:       if (rank .eq. 0) then
152:         write(*,120) tname
153:       endif
154:  120  format (' Solution method: ',A)
155:       PetscCallA(PEPGetDimensions(pep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
156:       if (rank .eq. 0) then
157:         write(*,130) nev
158:       endif
159:  130  format (' Number of requested eigenvalues:',I4)

161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: !     Display solution and clean up
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

165: !     ** show detailed info unless -terse option is given by user
166:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr))
167:       if (terse) then
168:         PetscCallA(PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_NULL_VIEWER,ierr))
169:       else
170:         PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr))
171:         PetscCallA(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD,ierr))
172:         PetscCallA(PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD,ierr))
173:         PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr))
174:       endif
175:       PetscCallA(PEPDestroy(pep,ierr))
176:       PetscCallA(MatDestroy(K,ierr))
177:       PetscCallA(MatDestroy(C,ierr))
178:       PetscCallA(MatDestroy(M,ierr))
179:       PetscCallA(SlepcFinalize(ierr))
180:       end

182: !/*TEST
183: !
184: !   test:
185: !      suffix: 1
186: !      args: -pep_nev 4 -pep_ncv 19 -terse
187: !      requires: !complex
188: !
189: !TEST*/
```