Actual source code: test3f.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: !  Description: Simple example to test the PEP Fortran interface.
     11: !
     12: ! ----------------------------------------------------------------------
     13: !
     14: #include <slepc/finclude/slepcpep.h>
     15: program test3f
     16:   use slepcpep
     17:   implicit none

     19: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     20: ! Declarations
     21: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     23:   Mat                  :: A(3), B
     24:   PEP                  :: pep
     25:   ST                   :: st
     26:   KSP                  :: ksp
     27:   DS                   :: ds
     28:   PetscReal            :: tol, tolabs, alpha, lambda
     29:   PetscScalar          :: tget, val
     30:   PetscInt             :: n, i, its, Istart, Iend, nev, ncv, mpd, nmat, np
     31:   PEPWhich             :: which
     32:   PEPConvergedReason   :: reason
     33:   PEPType              :: tname
     34:   PEPExtract           :: extr
     35:   PEPBasis             :: basis
     36:   PEPScale             :: scal
     37:   PEPRefine            :: refine
     38:   PEPRefineScheme      :: rscheme
     39:   PEPConv              :: conv
     40:   PEPStop              :: stp
     41:   PEPProblemType       :: ptype
     42:   PetscMPIInt          :: rank
     43:   PetscErrorCode       :: ierr
     44:   PetscViewerAndFormat :: vf

     46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     47: ! Beginning of program
     48: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     50:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
     51:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
     52:   n = 20
     53:   if (rank == 0) then
     54:     write (*, '(/a,i3,a)') 'Diagonal Quadratic Eigenproblem, n =', n, ' (Fortran)'
     55:   end if

     57:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A(1), ierr))
     58:   PetscCallA(MatSetSizes(A(1), PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
     59:   PetscCallA(MatSetFromOptions(A(1), ierr))
     60:   PetscCallA(MatGetOwnershipRange(A(1), Istart, Iend, ierr))
     61:   do i = Istart, Iend - 1
     62:     val = i + 1
     63:     PetscCallA(MatSetValue(A(1), i, i, val, INSERT_VALUES, ierr))
     64:   end do
     65:   PetscCallA(MatAssemblyBegin(A(1), MAT_FINAL_ASSEMBLY, ierr))
     66:   PetscCallA(MatAssemblyEnd(A(1), MAT_FINAL_ASSEMBLY, ierr))

     68:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A(2), ierr))
     69:   PetscCallA(MatSetSizes(A(2), PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
     70:   PetscCallA(MatSetFromOptions(A(2), ierr))
     71:   PetscCallA(MatGetOwnershipRange(A(2), Istart, Iend, ierr))
     72:   do i = Istart, Iend - 1
     73:     val = 1
     74:     PetscCallA(MatSetValue(A(2), i, i, val, INSERT_VALUES, ierr))
     75:   end do
     76:   PetscCallA(MatAssemblyBegin(A(2), MAT_FINAL_ASSEMBLY, ierr))
     77:   PetscCallA(MatAssemblyEnd(A(2), MAT_FINAL_ASSEMBLY, ierr))

     79:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A(3), ierr))
     80:   PetscCallA(MatSetSizes(A(3), PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
     81:   PetscCallA(MatSetFromOptions(A(3), ierr))
     82:   PetscCallA(MatGetOwnershipRange(A(3), Istart, Iend, ierr))
     83:   do i = Istart, Iend - 1
     84:     val = real(n, PETSC_REAL_KIND)/real(i + 1, PETSC_REAL_KIND)
     85:     PetscCallA(MatSetValue(A(3), i, i, val, INSERT_VALUES, ierr))
     86:   end do
     87:   PetscCallA(MatAssemblyBegin(A(3), MAT_FINAL_ASSEMBLY, ierr))
     88:   PetscCallA(MatAssemblyEnd(A(3), MAT_FINAL_ASSEMBLY, ierr))

     90: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     91: ! Create eigensolver and test interface functions
     92: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     94:   PetscCallA(PEPCreate(PETSC_COMM_WORLD, pep, ierr))
     95:   nmat = 3
     96:   PetscCallA(PEPSetOperators(pep, nmat, A, ierr))
     97:   PetscCallA(PEPGetNumMatrices(pep, nmat, ierr))
     98:   if (rank == 0) then
     99:     write (*, '(a,i2)') ' Polynomial of degree ', nmat - 1
    100:   end if
    101:   i = 0
    102:   PetscCallA(PEPGetOperators(pep, i, B, ierr))
    103:   PetscCallA(MatView(B, PETSC_NULL_VIEWER, ierr))

    105:   PetscCallA(PEPSetType(pep, PEPTOAR, ierr))
    106:   PetscCallA(PEPGetType(pep, tname, ierr))
    107:   if (rank == 0) then
    108:     write (*, '(a,a)') ' Type set to ', tname
    109:   end if

    111:   PetscCallA(PEPGetProblemType(pep, ptype, ierr))
    112:   if (rank == 0) then
    113:     write (*, '(a,i2)') ' Problem type before changing = ', ptype
    114:   end if
    115:   PetscCallA(PEPSetProblemType(pep, PEP_HERMITIAN, ierr))
    116:   PetscCallA(PEPGetProblemType(pep, ptype, ierr))
    117:   if (rank == 0) then
    118:     write (*, '(a,i2)') ' ... changed to ', ptype
    119:   end if

    121:   PetscCallA(PEPGetExtract(pep, extr, ierr))
    122:   if (rank == 0) then
    123:     write (*, '(a,i2)') ' Extraction before changing = ', extr
    124:   end if
    125:   PetscCallA(PEPSetExtract(pep, PEP_EXTRACT_STRUCTURED, ierr))
    126:   PetscCallA(PEPGetExtract(pep, extr, ierr))
    127:   if (rank == 0) then
    128:     write (*, '(a,i2)') ' ... changed to ', extr
    129:   end if

    131:   alpha = .1
    132:   its = 5
    133:   lambda = 1.
    134:   scal = PEP_SCALE_SCALAR
    135:   PetscCallA(PEPSetScale(pep, scal, alpha, PETSC_NULL_VEC, PETSC_NULL_VEC, its, lambda, ierr))
    136:   PetscCallA(PEPGetScale(pep, scal, alpha, PETSC_NULL_VEC, PETSC_NULL_VEC, its, lambda, ierr))
    137:   if (rank == 0) then
    138:     write (*, '(a,i2,a,f7.4,a,i2)') ' Scaling: ', scal, ', alpha=', alpha, ', its=', its
    139:   end if

    141:   basis = PEP_BASIS_CHEBYSHEV1
    142:   PetscCallA(PEPSetBasis(pep, basis, ierr))
    143:   PetscCallA(PEPGetBasis(pep, basis, ierr))
    144:   if (rank == 0) then
    145:     write (*, '(a,i2)') ' Polynomial basis: ', basis
    146:   end if

    148:   np = 1
    149:   tol = 1e-9
    150:   its = 2
    151:   refine = PEP_REFINE_SIMPLE
    152:   rscheme = PEP_REFINE_SCHEME_SCHUR
    153:   PetscCallA(PEPSetRefine(pep, refine, np, tol, its, rscheme, ierr))
    154:   PetscCallA(PEPGetRefine(pep, refine, np, tol, its, rscheme, ierr))
    155:   if (rank == 0) then
    156:     write (*, '(a,i2,a,f7.4,a,i2,a,i2)') ' Refinement: ', refine, ', tol=', tol, ', its=', its, ', schem=', rscheme
    157:   end if

    159:   tget = 4.8
    160:   PetscCallA(PEPSetTarget(pep, tget, ierr))
    161:   PetscCallA(PEPGetTarget(pep, tget, ierr))
    162:   PetscCallA(PEPSetWhichEigenpairs(pep, PEP_TARGET_MAGNITUDE, ierr))
    163:   PetscCallA(PEPGetWhichEigenpairs(pep, which, ierr))
    164:   if (rank == 0) then
    165:     write (*, '(a,i2,a,f4.1)') ' Which = ', which, ', target = ', PetscRealPart(tget)
    166:   end if

    168:   nev = 4
    169:   PetscCallA(PEPSetDimensions(pep, nev, PETSC_DETERMINE_INTEGER, PETSC_DETERMINE_INTEGER, ierr))
    170:   PetscCallA(PEPGetDimensions(pep, nev, ncv, mpd, ierr))
    171:   if (rank == 0) then
    172:     write (*, '(a,i2,a,i2,a,i2)') ' Dimensions: nev=', nev, ', ncv=', ncv, ', mpd=', mpd
    173:   end if

    175:   tol = 2.2e-4
    176:   its = 200
    177:   PetscCallA(PEPSetTolerances(pep, tol, its, ierr))
    178:   PetscCallA(PEPGetTolerances(pep, tol, its, ierr))
    179:   if (rank == 0) then
    180:     write (*, '(a,f8.5,a,i4)') ' Tolerance =', tol, ', max_its =', its
    181:   end if

    183:   PetscCallA(PEPSetConvergenceTest(pep, PEP_CONV_ABS, ierr))
    184:   PetscCallA(PEPGetConvergenceTest(pep, conv, ierr))
    185:   PetscCallA(PEPSetStoppingTest(pep, PEP_STOP_BASIC, ierr))
    186:   PetscCallA(PEPGetStoppingTest(pep, stp, ierr))
    187:   if (rank == 0) then
    188:     write (*, '(a,i2,a,i2)') ' Convergence test =', conv, ', stopping test =', stp
    189:   end if

    191:   PetscCallA(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, vf, ierr))
    192:   PetscCallA(PEPMonitorSet(pep, PEPMONITORFIRST, vf, PetscViewerAndFormatDestroy, ierr))
    193:   PetscCallA(PEPMonitorConvergedCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, PETSC_NULL, vf, ierr))
    194:   PetscCallA(PEPMonitorSet(pep, PEPMONITORCONVERGED, vf, PetscViewerAndFormatDestroy, ierr))
    195:   PetscCallA(PEPMonitorCancel(pep, ierr))

    197:   PetscCallA(PEPGetST(pep, st, ierr))
    198:   PetscCallA(STGetKSP(st, ksp, ierr))
    199:   tol = 1.e-8
    200:   tolabs = 1.e-35
    201:   PetscCallA(KSPSetTolerances(ksp, tol, tolabs, PETSC_CURRENT_REAL, PETSC_CURRENT_INTEGER, ierr))
    202:   PetscCallA(STView(st, PETSC_NULL_VIEWER, ierr))
    203:   PetscCallA(PEPGetDS(pep, ds, ierr))
    204:   PetscCallA(DSView(ds, PETSC_NULL_VIEWER, ierr))

    206:   PetscCallA(PEPSetFromOptions(pep, ierr))
    207:   PetscCallA(PEPSolve(pep, ierr))
    208:   PetscCallA(PEPGetConvergedReason(pep, reason, ierr))
    209:   if (rank == 0) then
    210:     write (*, '(a,i2)') ' Finished - converged reason =', reason
    211:   end if

    213: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    214: ! Display solution and clean up
    215: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    216:   PetscCallA(PEPErrorView(pep, PEP_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
    217:   PetscCallA(PEPDestroy(pep, ierr))
    218:   PetscCallA(MatDestroy(A(1), ierr))
    219:   PetscCallA(MatDestroy(A(2), ierr))
    220:   PetscCallA(MatDestroy(A(3), ierr))

    222:   PetscCallA(SlepcFinalize(ierr))
    223: end program test3f

    225: !/*TEST
    226: !
    227: !   test:
    228: !      suffix: 1
    229: !      args: -pep_tol 1e-5 -pep_ncv 22
    230: !      filter: sed -e "s/[+-]0\.0*i//g"
    231: !
    232: !TEST*/