Actual source code: ex16f.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> ./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: #include <slepc/finclude/slepcpep.h>
     22: program ex16f
     23:   use slepcpep
     24:   implicit none

     26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     27: ! Declarations
     28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     30:   Mat            :: M, C, K, A(3)   ! problem matrices
     31:   PEP            :: pep             ! polynomial eigenproblem solver context
     32:   PEPType        :: tname
     33:   PetscInt       :: N, nx, ny, i, j, Istart, Iend, II
     34:   PetscInt       :: nev
     35:   PetscMPIInt    :: rank
     36:   PetscErrorCode :: ierr
     37:   PetscBool      :: flg, terse
     38:   PetscScalar    :: val
     39:   PetscScalar, parameter :: two = 2.0, four = 4.0, mone = -1.0

     41: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     42: ! Beginning of program
     43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     45:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
     46:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
     47:   nx = 10
     48:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', nx, flg, ierr))
     49:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', ny, flg, ierr))
     50:   if (.not. flg) then
     51:     ny = nx
     52:   end if
     53:   N = nx*ny
     54:   if (rank == 0) then
     55:     write (*, '(/a,i6,a,i4,a,i4,a)') 'Quadratic Eigenproblem, N=', N, ' (', nx, 'x', ny, ' grid)'
     56:   end if

     58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     59: ! Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
     60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     62: ! ** K is the 2-D Laplacian
     63:   PetscCallA(MatCreate(PETSC_COMM_WORLD, K, ierr))
     64:   PetscCallA(MatSetSizes(K, PETSC_DECIDE, PETSC_DECIDE, N, N, ierr))
     65:   PetscCallA(MatSetFromOptions(K, ierr))
     66:   PetscCallA(MatGetOwnershipRange(K, Istart, Iend, ierr))
     67:   do II = Istart, Iend - 1
     68:     i = II/nx
     69:     j = II - i*nx
     70:     if (i > 0) then
     71:       PetscCallA(MatSetValue(K, II, II - nx, mone, INSERT_VALUES, ierr))
     72:     end if
     73:     if (i < ny - 1) then
     74:       PetscCallA(MatSetValue(K, II, II + nx, mone, INSERT_VALUES, ierr))
     75:     end if
     76:     if (j > 0) then
     77:       PetscCallA(MatSetValue(K, II, II - 1, mone, INSERT_VALUES, ierr))
     78:     end if
     79:     if (j < nx - 1) then
     80:       PetscCallA(MatSetValue(K, II, II + 1, mone, INSERT_VALUES, ierr))
     81:     end if
     82:     PetscCallA(MatSetValue(K, II, II, four, INSERT_VALUES, ierr))
     83:   end do
     84:   PetscCallA(MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY, ierr))
     85:   PetscCallA(MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY, ierr))

     87: ! ** C is the 1-D Laplacian on horizontal lines
     88:   PetscCallA(MatCreate(PETSC_COMM_WORLD, C, ierr))
     89:   PetscCallA(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N, ierr))
     90:   PetscCallA(MatSetFromOptions(C, ierr))
     91:   PetscCallA(MatGetOwnershipRange(C, Istart, Iend, ierr))
     92:   do II = Istart, Iend - 1
     93:     i = II/nx
     94:     j = II - i*nx
     95:     if (j > 0) then
     96:       PetscCallA(MatSetValue(C, II, II - 1, mone, INSERT_VALUES, ierr))
     97:     end if
     98:     if (j < nx - 1) then
     99:       PetscCallA(MatSetValue(C, II, II + 1, mone, INSERT_VALUES, ierr))
    100:     end if
    101:     PetscCallA(MatSetValue(C, II, II, two, INSERT_VALUES, ierr))
    102:   end do
    103:   PetscCallA(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY, ierr))
    104:   PetscCallA(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY, ierr))

    106: ! ** M is a diagonal matrix
    107:   PetscCallA(MatCreate(PETSC_COMM_WORLD, M, ierr))
    108:   PetscCallA(MatSetSizes(M, PETSC_DECIDE, PETSC_DECIDE, N, N, ierr))
    109:   PetscCallA(MatSetFromOptions(M, ierr))
    110:   PetscCallA(MatGetOwnershipRange(M, Istart, Iend, ierr))
    111:   do II = Istart, Iend - 1
    112:     val = II + 1
    113:     PetscCallA(MatSetValue(M, II, II, val, INSERT_VALUES, ierr))
    114:   end do
    115:   PetscCallA(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY, ierr))
    116:   PetscCallA(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY, ierr))

    118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    119: ! Create the eigensolver and set various options
    120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    122: ! ** Create eigensolver context
    123:   PetscCallA(PEPCreate(PETSC_COMM_WORLD, pep, ierr))

    125: ! ** Set matrices and problem type
    126:   A(1) = K
    127:   A(2) = C
    128:   A(3) = M
    129:   PetscCallA(PEPSetOperators(pep, 3_PETSC_INT_KIND, A, ierr))
    130:   PetscCallA(PEPSetProblemType(pep, PEP_GENERAL, ierr))

    132: ! ** Set solver parameters at runtime
    133:   PetscCallA(PEPSetFromOptions(pep, ierr))

    135: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    136: ! Solve the eigensystem
    137: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    139:   PetscCallA(PEPSolve(pep, ierr))

    141: ! ** Optional: Get some information from the solver and display it
    142:   PetscCallA(PEPGetType(pep, tname, ierr))
    143:   if (rank == 0) then
    144:     write (*, '(a,a)') ' Solution method: ', tname
    145:   end if
    146:   PetscCallA(PEPGetDimensions(pep, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
    147:   if (rank == 0) then
    148:     write (*, '(a,i4)') ' Number of requested eigenvalues:', nev
    149:   end if

    151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    152: ! Display solution and clean up
    153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    155: ! ** show detailed info unless -terse option is given by user
    156:   PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
    157:   if (terse) then
    158:     PetscCallA(PEPErrorView(pep, PEP_ERROR_BACKWARD, PETSC_NULL_VIEWER, ierr))
    159:   else
    160:     PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
    161:     PetscCallA(PEPConvergedReasonView(pep, PETSC_VIEWER_STDOUT_WORLD, ierr))
    162:     PetscCallA(PEPErrorView(pep, PEP_ERROR_BACKWARD, PETSC_VIEWER_STDOUT_WORLD, ierr))
    163:     PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
    164:   end if
    165:   PetscCallA(PEPDestroy(pep, ierr))
    166:   PetscCallA(MatDestroy(K, ierr))
    167:   PetscCallA(MatDestroy(C, ierr))
    168:   PetscCallA(MatDestroy(M, ierr))
    169:   PetscCallA(SlepcFinalize(ierr))
    170: end program ex16f

    172: !/*TEST
    173: !
    174: !   test:
    175: !      suffix: 1
    176: !      args: -pep_nev 4 -pep_ncv 19 -terse
    177: !      requires: !complex
    178: !   test:
    179: !      suffix: 2
    180: !      args: -pep_nev 4 -pep_smallest_real -terse
    181: !      output_file: output/ex16f_1.out
    182: !      requires: !complex
    183: !
    184: !TEST*/