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:       program main
 17: #include <slepc/finclude/slepcfn.h>
 18:       use slepcfn
 19:       implicit none

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

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

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

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

 40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 41: !     Polynomial p(x)
 42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43:       na = 5
 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,na,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:       na = 0
 61:       nb = 3
 62:       q(1) = -3.1
 63:       q(2) = 1.1
 64:       q(3) = 1.0
 65:       PetscCallA(FNSetType(fn,FNRATIONAL,ierr))
 66:       PetscCallA(FNRationalSetNumerator(fn,na,PETSC_NULL_SCALAR,ierr))
 67:       PetscCallA(FNRationalSetDenominator(fn,nb,q,ierr))
 68:       PetscCallA(FNView(fn,PETSC_NULL_VIEWER,ierr))
 69:       x = 2.2
 70:       PetscCallA(FNEvaluateFunction(fn,x,y,ierr))
 71:       PetscCallA(FNEvaluateDerivative(fn,x,yp,ierr))
 72:       call PrintInfo(x,y,yp)

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

 96:       PetscCallA(FNRationalGetNumerator(fn,n,pp,ierr))
 97:       if (rank .eq. 0) then
 98:         write(*,100) 'Numerator',(PetscRealPart(pp(i)),i=1,n)
 99:       end if
100:       PetscCallA(FNRationalGetDenominator(fn,n,qq,ierr))
101:       if (rank .eq. 0) then
102:         write(*,100) 'Denominator',(PetscRealPart(qq(i)),i=1,n)
103:       end if
104:  100  format (A15,10F6.1)

106: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: !     Constant
108: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109:       na = 1
110:       nb = 0
111:       five = 5.0
112:       PetscCallA(FNSetType(fn,FNRATIONAL,ierr))
113:       PetscCallA(FNRationalSetNumerator(fn,na,five,ierr))
114:       PetscCallA(FNRationalSetDenominator(fn,nb,PETSC_NULL_SCALAR,ierr))
115:       PetscCallA(FNView(fn,PETSC_NULL_VIEWER,ierr))
116:       x = 2.2
117:       PetscCallA(FNEvaluateFunction(fn,x,y,ierr))
118:       PetscCallA(FNEvaluateDerivative(fn,x,yp,ierr))
119:       call PrintInfo(x,y,yp)

121: !     *** Clean up
122:       PetscCallA(FNDestroy(fn,ierr))
123:       PetscCallA(SlepcFinalize(ierr))
124:       end

126: ! -----------------------------------------------------------------

128:       subroutine PrintInfo(x,y,yp)
129: #include <slepc/finclude/slepcfn.h>
130:       use slepcfn
131:       implicit none
132:       PetscScalar    x,y,yp
133:       PetscReal      re,im
134:       PetscMPIInt    rank
135:       PetscErrorCode ierr

137:       PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
138:       if (rank .eq. 0) then
139:         re = PetscRealPart(y)
140:         im = PetscImaginaryPart(y)
141:         if (abs(im).lt.1.d-10) then
142:           write(*,110) 'f', PetscRealPart(x), re
143:         else
144:           write(*,120) 'f', PetscRealPart(x), re, im
145:         endif
146:         re = PetscRealPart(yp)
147:         im = PetscImaginaryPart(yp)
148:         if (abs(im).lt.1.d-10) then
149:           write(*,110) 'f''', PetscRealPart(x), re
150:         else
151:           write(*,120) 'f''', PetscRealPart(x), re, im
152:         endif
153:       endif
154:  110  format (A2,'(',F4.1,') = ',F10.5)
155:  120  format (A2,'(',F4.1,') = ',F10.5,SP,F9.5,'i')

157:       end

159: !/*TEST
160: !
161: !   test:
162: !      suffix: 1
163: !      nsize: 1
164: !      requires: !single
165: !
166: !TEST*/