Actual source code: test1f.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> ./test1f [-help]
     11: !
     12: !  Description: Test rational function in Fortran.
     13: !
     14: ! ----------------------------------------------------------------------
     15: !
     16: #include <slepc/finclude/slepcfn.h>
     17: program test1f
     18:   use slepcfn
     19:   implicit none

     21: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     22: ! Declarations
     23: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     25:   FN             :: fn
     26:   PetscInt       :: i, n
     27:   PetscMPIInt    :: rank
     28:   PetscErrorCode :: ierr
     29:   PetscScalar    :: x, y, yp, p(10), q(10)
     30:   PetscScalar    :: pp(10), qq(10), tau, eta
     31:   PetscScalar, parameter :: five = 5.0

     33: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     34: ! Beginning of program
     35: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     37:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
     38:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
     39:   PetscCallA(FNCreate(PETSC_COMM_WORLD, fn, ierr))

     41: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     42: ! Polynomial p(x)
     43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     44:   p(1) = -3.1
     45:   p(2) = 1.1
     46:   p(3) = 1.0
     47:   p(4) = -2.0
     48:   p(5) = 3.5
     49:   PetscCallA(FNSetType(fn, FNRATIONAL, ierr))
     50:   PetscCallA(FNRationalSetNumerator(fn, 5_PETSC_INT_KIND, p, ierr))
     51:   PetscCallA(FNView(fn, PETSC_NULL_VIEWER, ierr))
     52:   x = 2.2
     53:   PetscCallA(FNEvaluateFunction(fn, x, y, ierr))
     54:   PetscCallA(FNEvaluateDerivative(fn, x, yp, ierr))
     55:   call PrintInfo(x, y, yp)

     57: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     58: ! Inverse of polynomial 1/q(x)
     59: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     60:   q(1) = -3.1
     61:   q(2) = 1.1
     62:   q(3) = 1.0
     63:   PetscCallA(FNSetType(fn, FNRATIONAL, ierr))
     64:   PetscCallA(FNRationalSetNumerator(fn, 0_PETSC_INT_KIND, PETSC_NULL_SCALAR_ARRAY, ierr))
     65:   PetscCallA(FNRationalSetDenominator(fn, 3_PETSC_INT_KIND, q, ierr))
     66:   PetscCallA(FNView(fn, PETSC_NULL_VIEWER, ierr))
     67:   x = 2.2
     68:   PetscCallA(FNEvaluateFunction(fn, x, y, ierr))
     69:   PetscCallA(FNEvaluateDerivative(fn, x, yp, ierr))
     70:   call PrintInfo(x, y, yp)

     72: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     73: ! Rational p(x)/q(x)
     74: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     75:   p(1) = 1.1
     76:   p(2) = 1.1
     77:   q(1) = 1.0
     78:   q(2) = -2.0
     79:   q(3) = 3.5
     80:   PetscCallA(FNSetType(fn, FNRATIONAL, ierr))
     81:   PetscCallA(FNRationalSetNumerator(fn, 2_PETSC_INT_KIND, p, ierr))
     82:   PetscCallA(FNRationalSetDenominator(fn, 3_PETSC_INT_KIND, q, ierr))
     83:   tau = 1.2
     84:   eta = 0.5
     85:   PetscCallA(FNSetScale(fn, tau, eta, ierr))
     86:   PetscCallA(FNView(fn, PETSC_NULL_VIEWER, ierr))
     87:   x = 2.2
     88:   PetscCallA(FNEvaluateFunction(fn, x, y, ierr))
     89:   PetscCallA(FNEvaluateDerivative(fn, x, yp, ierr))
     90:   call PrintInfo(x, y, yp)

     92:   PetscCallA(FNRationalGetNumerator(fn, n, pp, ierr))
     93:   if (rank == 0) then
     94:     write (*, '(a15,10f6.1)') 'Numerator', (PetscRealPart(pp(i)), i=1, n)
     95:   end if
     96:   PetscCallA(FNRationalGetDenominator(fn, n, qq, ierr))
     97:   if (rank == 0) then
     98:     write (*, '(a15,10f6.1)') 'Denominator', (PetscRealPart(qq(i)), i=1, n)
     99:   end if

    101: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    102: ! Constant
    103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    104:   PetscCallA(FNSetType(fn, FNRATIONAL, ierr))
    105:   PetscCallA(FNRationalSetNumerator(fn, 1_PETSC_INT_KIND, [five], ierr))
    106:   PetscCallA(FNRationalSetDenominator(fn, 0_PETSC_INT_KIND, PETSC_NULL_SCALAR_ARRAY, ierr))
    107:   PetscCallA(FNView(fn, PETSC_NULL_VIEWER, ierr))
    108:   x = 2.2
    109:   PetscCallA(FNEvaluateFunction(fn, x, y, ierr))
    110:   PetscCallA(FNEvaluateDerivative(fn, x, yp, ierr))
    111:   call PrintInfo(x, y, yp)

    113: ! *** Clean up
    114:   PetscCallA(FNDestroy(fn, ierr))
    115:   PetscCallA(SlepcFinalize(ierr))

    117: contains

    119:   subroutine PrintInfo(x, y, yp)
    120:     use slepcfn
    121:     implicit none
    122:     PetscScalar    :: x, y, yp
    123:     PetscReal      :: re, im
    124:     PetscMPIInt    :: rank
    125:     PetscErrorCode :: ierr

    127:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
    128:     if (rank == 0) then
    129:       re = PetscRealPart(y)
    130:       im = PetscImaginaryPart(y)
    131:       if (abs(im) < 1.d-10) then
    132:         write (*, '(a3,f3.1,a,f10.5)') 'f(', PetscRealPart(x), ') = ', re
    133:       else
    134:         write (*, '(a3,f3.1,a,f10.5,sp,f9.5,a)') 'f(', PetscRealPart(x), ') = ', re, im, 'i'
    135:       end if
    136:       re = PetscRealPart(yp)
    137:       im = PetscImaginaryPart(yp)
    138:       if (abs(im) < 1.d-10) then
    139:         write (*, '(a3,f3.1,a,f10.5)') 'f''(', PetscRealPart(x), ') = ', re
    140:       else
    141:         write (*, '(a3,f3.1,a,f10.5,sp,f9.5,a)') 'f''(', PetscRealPart(x), ') = ', re, im, 'i'
    142:       end if
    143:     end if

    145:   end subroutine

    147: end program test1f

    149: !/*TEST
    150: !
    151: !   test:
    152: !      suffix: 1
    153: !      nsize: 1
    154: !      requires: !single
    155: !
    156: !TEST*/