Actual source code: test15f.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> ./test15f [-help] [-n <n>] [all SLEPc options]
     11: !
     12: !  Description: Tests custom monitors from Fortran.
     13: !
     14: !  The command line options are:
     15: !    -n <n>, where <n> = number of grid points = matrix size
     16: !    -my_eps_monitor, activates the custom monitor
     17: !
     18: ! ----------------------------------------------------------------------
     19: !
     20: #include <slepc/finclude/slepceps.h>

     22: module test15fmodule
     23:   use slepceps
     24:   implicit none

     26: contains
     27:   ! --------------------------------------------------------------
     28:   ! MyEPSMonitor - This is a user-defined routine for monitoring
     29:   ! the EPS iterative solvers.
     30:   !
     31:   ! Input Parameters:
     32:   !   eps   - eigensolver context
     33:   !   its   - iteration number
     34:   !   nconv - number of converged eigenpairs
     35:   !   eigr  - real part of the eigenvalues
     36:   !   eigi  - imaginary part of the eigenvalues
     37:   !   errest- relative error estimates for each eigenpair
     38:   !   nest  - number of error estimates
     39:   !   dummy - optional user-defined monitor context (unused here)
     40:   !
     41:   subroutine MyEPSMonitor(eps, its, nconv, eigr, eigi, errest, nest, dummy, ierr)
     42:     use slepceps
     43:     implicit none

     45:     EPS            :: eps
     46:     PetscInt       :: its, nconv, nest, dummy
     47:     PetscScalar    :: eigr(*), eigi(*)
     48:     PetscReal      :: re, errest(*)
     49:     PetscMPIInt    :: rank
     50:     PetscErrorCode, intent(out) :: ierr

     52:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
     53:     if (its > 0 .and. rank == 0) then
     54:       re = PetscRealPart(eigr(nconv + 1))
     55:       write (6, '(i3,a,i2,a,f7.4,a,g10.3,a)') its, ' EPS nconv=', nconv, ' first unconverged value (error) ', re, ' (', errest(nconv + 1), ')'
     56:     end if
     57:     ierr = 0
     58:   end subroutine

     60: end module test15fmodule

     62: program test15f
     63:   use slepceps
     64:   use test15fmodule
     65:   implicit none

     67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     68: ! Declarations
     69: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     70: !
     71:   Mat            :: A     ! operator matrix
     72:   EPS            :: eps   ! eigenproblem solver context
     73:   EPSType        :: tname
     74:   PetscInt       :: n, i, Istart, Iend, nev
     75:   PetscInt       :: col(3)
     76:   PetscInt       :: i1, i2, i3
     77:   PetscMPIInt    :: rank
     78:   PetscErrorCode :: ierr
     79:   PetscBool      :: flg
     80:   PetscScalar    :: val(3)

     82: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     83: ! Beginning of program
     84: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     86:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
     87:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
     88:   n = 30
     89:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))

     91:   if (rank == 0) then
     92:     write (*, '(/a,i3,a)') '1-D Laplacian Eigenproblem, n =', n, ' (Fortran)'
     93:   end if

     95: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     96: ! Compute the operator matrix that defines the eigensystem, Ax=kx
     97: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     99:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
    100:   PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
    101:   PetscCallA(MatSetFromOptions(A, ierr))

    103:   i1 = 1
    104:   i2 = 2
    105:   i3 = 3
    106:   PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
    107:   if (Istart == 0) then
    108:     i = 0
    109:     col(1) = 0
    110:     col(2) = 1
    111:     val(1) = 2.0
    112:     val(2) = -1.0
    113:     PetscCallA(MatSetValues(A, i1, [i], i2, col, val, INSERT_VALUES, ierr))
    114:     Istart = Istart + 1
    115:   end if
    116:   if (Iend == n) then
    117:     i = n - 1
    118:     col(1) = n - 2
    119:     col(2) = n - 1
    120:     val(1) = -1.0
    121:     val(2) = 2.0
    122:     PetscCallA(MatSetValues(A, i1, [i], i2, col, val, INSERT_VALUES, ierr))
    123:     Iend = Iend - 1
    124:   end if
    125:   val(1) = -1.0
    126:   val(2) = 2.0
    127:   val(3) = -1.0
    128:   do i = Istart, Iend - 1
    129:     col(1) = i - 1
    130:     col(2) = i
    131:     col(3) = i + 1
    132:     PetscCallA(MatSetValues(A, i1, [i], i3, col, val, INSERT_VALUES, ierr))
    133:   end do

    135:   PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
    136:   PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

    138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    139: ! Create the eigensolver and display info
    140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    142: ! ** Create eigensolver context
    143:   PetscCallA(EPSCreate(PETSC_COMM_WORLD, eps, ierr))

    145: ! ** Set operators. In this case, it is a standard eigenvalue problem
    146:   PetscCallA(EPSSetOperators(eps, A, PETSC_NULL_MAT, ierr))
    147:   PetscCallA(EPSSetProblemType(eps, EPS_HEP, ierr))

    149: ! ** Set user-defined monitor
    150:   PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my_eps_monitor', flg, ierr))
    151:   if (flg) then
    152:     PetscCallA(EPSMonitorSet(eps, MyEPSMonitor, 0, PETSC_NULL_FUNCTION, ierr))
    153:   end if

    155: ! ** Set solver parameters at runtime
    156:   PetscCallA(EPSSetFromOptions(eps, ierr))

    158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    159: ! Solve the eigensystem
    160: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    162:   PetscCallA(EPSSolve(eps, ierr))

    164: ! ** Optional: Get some information from the solver and display it
    165:   PetscCallA(EPSGetType(eps, tname, ierr))
    166:   if (rank == 0) then
    167:     write (*, '(a,a)') ' Solution method: ', tname
    168:   end if
    169:   PetscCallA(EPSGetDimensions(eps, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
    170:   if (rank == 0) then
    171:     write (*, '(a,i2)') ' Number of requested eigenvalues:', nev
    172:   end if

    174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    175: ! Display solution and clean up
    176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    178:   PetscCallA(EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
    179:   PetscCallA(EPSDestroy(eps, ierr))
    180:   PetscCallA(MatDestroy(A, ierr))

    182:   PetscCallA(SlepcFinalize(ierr))
    183: end program test15f

    185: !/*TEST
    186: !
    187: !   test:
    188: !      suffix: 1
    189: !      args: -my_eps_monitor
    190: !      requires: double
    191: !
    192: !TEST*/