Actual source code: ex16f90.F90

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: !
  2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
  4: !  Copyright (c) 2002-2021, 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 the
 13: !  PEP object. 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:       call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
 49:       if (ierr .ne. 0) then
 50:         print*,'SlepcInitialize failed'
 51:         stop
 52:       endif
 53:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr);CHKERRA(ierr)
 54:       nx = 10
 55:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',nx,flg,ierr);CHKERRA(ierr)
 56:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',ny,flg,ierr);CHKERRA(ierr)
 57:       if (.not. flg) then
 58:         ny = nx
 59:       endif
 60:       N = nx*ny
 61:       if (rank .eq. 0) then
 62:         write(*,100) N, nx, ny
 63:       endif
 64:  100  format (/'Quadratic Eigenproblem, N=',I6,' (',I4,'x',I4,' grid)')

 66: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 67: !     Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
 68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 70: !     ** K is the 2-D Laplacian
 71:       call MatCreate(PETSC_COMM_WORLD,K,ierr);CHKERRA(ierr)
 72:       call MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr);CHKERRA(ierr)
 73:       call MatSetFromOptions(K,ierr);CHKERRA(ierr)
 74:       call MatSetUp(K,ierr);CHKERRA(ierr)
 75:       call MatGetOwnershipRange(K,Istart,Iend,ierr);CHKERRA(ierr)
 76:       mone = -1.0
 77:       four = 4.0
 78:       do II=Istart,Iend-1
 79:         i = II/nx
 80:         j = II-i*nx
 81:         if (i .gt. 0) then
 82:           call MatSetValue(K,II,II-nx,mone,INSERT_VALUES,ierr);CHKERRA(ierr)
 83:         endif
 84:         if (i .lt. ny-1) then
 85:           call MatSetValue(K,II,II+nx,mone,INSERT_VALUES,ierr);CHKERRA(ierr)
 86:         endif
 87:         if (j .gt. 0) then
 88:           call MatSetValue(K,II,II-1,mone,INSERT_VALUES,ierr);CHKERRA(ierr)
 89:         endif
 90:         if (j .lt. nx-1) then
 91:           call MatSetValue(K,II,II+1,mone,INSERT_VALUES,ierr);CHKERRA(ierr)
 92:         endif
 93:         call MatSetValue(K,II,II,four,INSERT_VALUES,ierr);CHKERRA(ierr)
 94:       end do
 95:       call MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
 96:       call MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)

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

119: !     ** M is a diagonal matrix
120:       call MatCreate(PETSC_COMM_WORLD,M,ierr);CHKERRA(ierr)
121:       call MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr);CHKERRA(ierr)
122:       call MatSetFromOptions(M,ierr);CHKERRA(ierr)
123:       call MatSetUp(M,ierr);CHKERRA(ierr)
124:       call MatGetOwnershipRange(M,Istart,Iend,ierr);CHKERRA(ierr)
125:       do II=Istart,Iend-1
126:         val = II+1
127:         call MatSetValue(M,II,II,val,INSERT_VALUES,ierr);CHKERRA(ierr)
128:       end do
129:       call MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
130:       call MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)

132: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: !     Create the eigensolver and set various options
134: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

136: !     ** Create eigensolver context
137:       call PEPCreate(PETSC_COMM_WORLD,pep,ierr);CHKERRA(ierr)

139: !     ** Set matrices and problem type
140:       A(1) = K
141:       A(2) = C
142:       A(3) = M
143:       ithree = 3
144:       call PEPSetOperators(pep,ithree,A,ierr);CHKERRA(ierr)
145:       call PEPSetProblemType(pep,PEP_GENERAL,ierr);CHKERRA(ierr)

147: !     ** Set solver parameters at runtime
148:       call PEPSetFromOptions(pep,ierr);CHKERRA(ierr)

150: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: !     Solve the eigensystem
152: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

154:       call PEPSolve(pep,ierr);CHKERRA(ierr)

156: !     ** Optional: Get some information from the solver and display it
157:       call PEPGetType(pep,tname,ierr);CHKERRA(ierr)
158:       if (rank .eq. 0) then
159:         write(*,120) tname
160:       endif
161:  120  format (' Solution method: ',A)
162:       call PEPGetDimensions(pep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
163:       if (rank .eq. 0) then
164:         write(*,130) nev
165:       endif
166:  130  format (' Number of requested eigenvalues:',I4)

168: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169: !     Display solution and clean up
170: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

172: !     ** show detailed info unless -terse option is given by user
173:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr);CHKERRA(ierr)
174:       if (terse) then
175:         call PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_NULL_VIEWER,ierr);CHKERRA(ierr)
176:       else
177:         call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr);CHKERRA(ierr)
178:         call PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
179:         call PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
180:         call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
181:       endif
182:       call PEPDestroy(pep,ierr);CHKERRA(ierr)
183:       call MatDestroy(K,ierr);CHKERRA(ierr)
184:       call MatDestroy(C,ierr);CHKERRA(ierr)
185:       call MatDestroy(M,ierr);CHKERRA(ierr)
186:       call SlepcFinalize(ierr)
187:       end

189: !/*TEST
190: !
191: !   test:
192: !      suffix: 1
193: !      args: -pep_nev 4 -pep_ncv 19 -terse
194: !      requires: !complex
195: !
196: !TEST*/