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*/