Actual source code: test2f.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: !  Description: Simple example to test the NEP Fortran interface.
     11: !
     12: ! ----------------------------------------------------------------------
     13: !
     14:       program main
     15: #include <slepc/finclude/slepcnep.h>
     16:       use slepcnep
     17:       implicit none

     19: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     20: !     Declarations
     21: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     22:       Mat                A(3),B
     23:       FN                 f(3),g
     24:       NEP                nep
     25:       DS                 ds
     26:       RG                 rg
     27:       PetscReal          tol
     28:       PetscScalar        coeffs(2),tget,val
     29:       PetscInt           n,i,its,Istart,Iend
     30:       PetscInt           nev,ncv,mpd,nterm
     31:       PetscInt           nc,np
     32:       NEPWhich           which
     33:       NEPConvergedReason reason
     34:       NEPType            tname
     35:       NEPRefine          refine
     36:       NEPRefineScheme    rscheme
     37:       NEPConv            conv
     38:       NEPStop            stp
     39:       NEPProblemType     ptype
     40:       MatStructure       mstr
     41:       PetscMPIInt        rank
     42:       PetscErrorCode     ierr
     43:       PetscViewerAndFormat vf

     45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     46: !     Beginning of program
     47: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     49:       PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
     50:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
     51:       n = 20
     52:       if (rank .eq. 0) then
     53:         write(*,100) n
     54:       endif
     55:  100  format (/'Diagonal Nonlinear Eigenproblem, n =',I3,' (Fortran)')

     57: !     Matrices
     58:       PetscCallA(MatCreate(PETSC_COMM_WORLD,A(1),ierr))
     59:       PetscCallA(MatSetSizes(A(1),PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
     60:       PetscCallA(MatSetFromOptions(A(1),ierr))
     61:       PetscCallA(MatGetOwnershipRange(A(1),Istart,Iend,ierr))
     62:       do i=Istart,Iend-1
     63:         val = i+1
     64:         PetscCallA(MatSetValue(A(1),i,i,val,INSERT_VALUES,ierr))
     65:       enddo
     66:       PetscCallA(MatAssemblyBegin(A(1),MAT_FINAL_ASSEMBLY,ierr))
     67:       PetscCallA(MatAssemblyEnd(A(1),MAT_FINAL_ASSEMBLY,ierr))

     69:       PetscCallA(MatCreate(PETSC_COMM_WORLD,A(2),ierr))
     70:       PetscCallA(MatSetSizes(A(2),PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
     71:       PetscCallA(MatSetFromOptions(A(2),ierr))
     72:       PetscCallA(MatGetOwnershipRange(A(2),Istart,Iend,ierr))
     73:       do i=Istart,Iend-1
     74:         val = 1
     75:         PetscCallA(MatSetValue(A(2),i,i,val,INSERT_VALUES,ierr))
     76:       enddo
     77:       PetscCallA(MatAssemblyBegin(A(2),MAT_FINAL_ASSEMBLY,ierr))
     78:       PetscCallA(MatAssemblyEnd(A(2),MAT_FINAL_ASSEMBLY,ierr))

     80:       PetscCallA(MatCreate(PETSC_COMM_WORLD,A(3),ierr))
     81:       PetscCallA(MatSetSizes(A(3),PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
     82:       PetscCallA(MatSetFromOptions(A(3),ierr))
     83:       PetscCallA(MatGetOwnershipRange(A(3),Istart,Iend,ierr))
     84:       do i=Istart,Iend-1
     85:         val = real(n)/real(i+1)
     86:         PetscCallA(MatSetValue(A(3),i,i,val,INSERT_VALUES,ierr))
     87:       enddo
     88:       PetscCallA(MatAssemblyBegin(A(3),MAT_FINAL_ASSEMBLY,ierr))
     89:       PetscCallA(MatAssemblyEnd(A(3),MAT_FINAL_ASSEMBLY,ierr))

     91: !     Functions: f0=-lambda, f1=1.0, f2=sqrt(lambda)
     92:       PetscCallA(FNCreate(PETSC_COMM_WORLD,f(1),ierr))
     93:       PetscCallA(FNSetType(f(1),FNRATIONAL,ierr))
     94:       nc = 2
     95:       coeffs(1) = -1.0
     96:       coeffs(2) = 0.0
     97:       PetscCallA(FNRationalSetNumerator(f(1),nc,coeffs,ierr))

     99:       PetscCallA(FNCreate(PETSC_COMM_WORLD,f(2),ierr))
    100:       PetscCallA(FNSetType(f(2),FNRATIONAL,ierr))
    101:       nc = 1
    102:       coeffs(1) = 1.0
    103:       PetscCallA(FNRationalSetNumerator(f(2),nc,coeffs,ierr))

    105:       PetscCallA(FNCreate(PETSC_COMM_WORLD,f(3),ierr))
    106:       PetscCallA(FNSetType(f(3),FNSQRT,ierr))

    108: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    109: !     Create eigensolver and test interface functions
    110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    112:       PetscCallA(NEPCreate(PETSC_COMM_WORLD,nep,ierr))
    113:       nterm = 3
    114:       mstr = SAME_NONZERO_PATTERN
    115:       PetscCallA(NEPSetSplitOperator(nep,nterm,A,f,mstr,ierr))
    116:       PetscCallA(NEPGetSplitOperatorInfo(nep,nterm,mstr,ierr))
    117:       if (rank .eq. 0) then
    118:         write(*,110) nterm
    119:       endif
    120:  110  format (' Nonlinear function with ',I2,' terms')
    121:       i = 0
    122:       PetscCallA(NEPGetSplitOperatorTerm(nep,i,B,g,ierr))
    123:       PetscCallA(MatView(B,PETSC_NULL_VIEWER,ierr))
    124:       PetscCallA(FNView(g,PETSC_NULL_VIEWER,ierr))

    126:       PetscCallA(NEPSetType(nep,NEPRII,ierr))
    127:       PetscCallA(NEPGetType(nep,tname,ierr))
    128:       if (rank .eq. 0) then
    129:         write(*,120) tname
    130:       endif
    131:  120  format (' Type set to ',A)

    133:       PetscCallA(NEPGetProblemType(nep,ptype,ierr))
    134:       if (rank .eq. 0) then
    135:         write(*,130) ptype
    136:       endif
    137:  130  format (' Problem type before changing = ',I2)
    138:       PetscCallA(NEPSetProblemType(nep,NEP_RATIONAL,ierr))
    139:       PetscCallA(NEPGetProblemType(nep,ptype,ierr))
    140:       if (rank .eq. 0) then
    141:         write(*,140) ptype
    142:       endif
    143:  140  format (' ... changed to ',I2)

    145:       np = 1
    146:       tol = 1e-9
    147:       its = 2
    148:       refine = NEP_REFINE_SIMPLE
    149:       rscheme = NEP_REFINE_SCHEME_EXPLICIT
    150:       PetscCallA(NEPSetRefine(nep,refine,np,tol,its,rscheme,ierr))
    151:       PetscCallA(NEPGetRefine(nep,refine,np,tol,its,rscheme,ierr))
    152:       if (rank .eq. 0) then
    153:         write(*,190) refine,tol,its,rscheme
    154:       endif
    155:  190  format (' Refinement: ',I2,', tol=',F12.9,', its=',I2,', scheme=',I2)

    157:       tget = 1.1
    158:       PetscCallA(NEPSetTarget(nep,tget,ierr))
    159:       PetscCallA(NEPGetTarget(nep,tget,ierr))
    160:       PetscCallA(NEPSetWhichEigenpairs(nep,NEP_TARGET_MAGNITUDE,ierr))
    161:       PetscCallA(NEPGetWhichEigenpairs(nep,which,ierr))
    162:       if (rank .eq. 0) then
    163:         write(*,200) which,PetscRealPart(tget)
    164:       endif
    165:  200  format (' Which = ',I2,', target = ',F4.1)

    167:       nev = 1
    168:       ncv = 12
    169:       PetscCallA(NEPSetDimensions(nep,nev,ncv,PETSC_DETERMINE_INTEGER,ierr))
    170:       PetscCallA(NEPGetDimensions(nep,nev,ncv,mpd,ierr))
    171:       if (rank .eq. 0) then
    172:         write(*,210) nev,ncv,mpd
    173:       endif
    174:  210  format (' Dimensions: nev=',I2,', ncv=',I2,', mpd=',I2)

    176:       tol = 1.0e-6
    177:       its = 200
    178:       PetscCallA(NEPSetTolerances(nep,tol,its,ierr))
    179:       PetscCallA(NEPGetTolerances(nep,tol,its,ierr))
    180:       if (rank .eq. 0) then
    181:         write(*,220) tol,its
    182:       endif
    183:  220  format (' Tolerance =',F9.6,', max_its =',I4)

    185:       PetscCallA(NEPSetConvergenceTest(nep,NEP_CONV_ABS,ierr))
    186:       PetscCallA(NEPGetConvergenceTest(nep,conv,ierr))
    187:       PetscCallA(NEPSetStoppingTest(nep,NEP_STOP_BASIC,ierr))
    188:       PetscCallA(NEPGetStoppingTest(nep,stp,ierr))
    189:       if (rank .eq. 0) then
    190:         write(*,230) conv,stp
    191:       endif
    192:  230  format (' Convergence test =',I2,', stopping test =',I2)

    194:       PetscCallA(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,vf,ierr))
    195:       PetscCallA(NEPMonitorSet(nep,NEPMONITORFIRST,vf,PetscViewerAndFormatDestroy,ierr))
    196:       PetscCallA(NEPMonitorConvergedCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,PETSC_NULL_VEC,vf,ierr))
    197:       PetscCallA(NEPMonitorSet(nep,NEPMONITORCONVERGED,vf,NEPMonitorConvergedDestroy,ierr))
    198:       PetscCallA(NEPMonitorCancel(nep,ierr))

    200:       PetscCallA(NEPGetDS(nep,ds,ierr))
    201:       PetscCallA(DSView(ds,PETSC_NULL_VIEWER,ierr))
    202:       PetscCallA(NEPSetFromOptions(nep,ierr))

    204:       PetscCallA(NEPGetRG(nep,rg,ierr))
    205:       PetscCallA(RGView(rg,PETSC_NULL_VIEWER,ierr))

    207:       PetscCallA(NEPSolve(nep,ierr))
    208:       PetscCallA(NEPGetConvergedReason(nep,reason,ierr))
    209:       if (rank .eq. 0) then
    210:         write(*,240) reason
    211:       endif
    212:  240  format (' Finished - converged reason =',I2)

    214: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    215: !     Display solution and clean up
    216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    217:       PetscCallA(NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
    218:       PetscCallA(NEPDestroy(nep,ierr))
    219:       PetscCallA(MatDestroy(A(1),ierr))
    220:       PetscCallA(MatDestroy(A(2),ierr))
    221:       PetscCallA(MatDestroy(A(3),ierr))
    222:       PetscCallA(FNDestroy(f(1),ierr))
    223:       PetscCallA(FNDestroy(f(2),ierr))
    224:       PetscCallA(FNDestroy(f(3),ierr))

    226:       PetscCallA(SlepcFinalize(ierr))
    227:       end

    229: !/*TEST
    230: !
    231: !   test:
    232: !      suffix: 1
    233: !      requires: !single
    234: !
    235: !TEST*/