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:       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_ARRAY,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_ARRAY,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*/