Actual source code: test3f.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 PEP Fortran interface.
     11: !
     12: ! ----------------------------------------------------------------------
     13: !
     14:       program main
     15: #include <slepc/finclude/slepcpep.h>
     16:       use slepcpep
     17:       implicit none

     19: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     20: !     Declarations
     21: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     22:       Mat                A(3),B
     23:       PEP                pep
     24:       ST                 st
     25:       KSP                ksp
     26:       DS                 ds
     27:       PetscReal          tol,tolabs,alpha,lambda
     28:       PetscScalar        tget,val
     29:       PetscInt           n,i,its,Istart,Iend
     30:       PetscInt           nev,ncv,mpd,nmat,np
     31:       PEPWhich           which
     32:       PEPConvergedReason reason
     33:       PEPType            tname
     34:       PEPExtract         extr
     35:       PEPBasis           basis
     36:       PEPScale           scal
     37:       PEPRefine          refine
     38:       PEPRefineScheme    rscheme
     39:       PEPConv            conv
     40:       PEPStop            stp
     41:       PEPProblemType     ptype
     42:       PetscMPIInt        rank
     43:       PetscErrorCode     ierr
     44:       PetscViewerAndFormat vf

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

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

     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: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     92: !     Create eigensolver and test interface functions
     93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     95:       PetscCallA(PEPCreate(PETSC_COMM_WORLD,pep,ierr))
     96:       nmat = 3
     97:       PetscCallA(PEPSetOperators(pep,nmat,A,ierr))
     98:       PetscCallA(PEPGetNumMatrices(pep,nmat,ierr))
     99:       if (rank .eq. 0) then
    100:         write(*,110) nmat-1
    101:       endif
    102:  110  format (' Polynomial of degree ',I2)
    103:       i = 0
    104:       PetscCallA(PEPGetOperators(pep,i,B,ierr))
    105:       PetscCallA(MatView(B,PETSC_NULL_VIEWER,ierr))

    107:       PetscCallA(PEPSetType(pep,PEPTOAR,ierr))
    108:       PetscCallA(PEPGetType(pep,tname,ierr))
    109:       if (rank .eq. 0) then
    110:         write(*,120) tname
    111:       endif
    112:  120  format (' Type set to ',A)

    114:       PetscCallA(PEPGetProblemType(pep,ptype,ierr))
    115:       if (rank .eq. 0) then
    116:         write(*,130) ptype
    117:       endif
    118:  130  format (' Problem type before changing = ',I2)
    119:       PetscCallA(PEPSetProblemType(pep,PEP_HERMITIAN,ierr))
    120:       PetscCallA(PEPGetProblemType(pep,ptype,ierr))
    121:       if (rank .eq. 0) then
    122:         write(*,140) ptype
    123:       endif
    124:  140  format (' ... changed to ',I2)

    126:       PetscCallA(PEPGetExtract(pep,extr,ierr))
    127:       if (rank .eq. 0) then
    128:         write(*,150) extr
    129:       endif
    130:  150  format (' Extraction before changing = ',I2)
    131:       PetscCallA(PEPSetExtract(pep,PEP_EXTRACT_STRUCTURED,ierr))
    132:       PetscCallA(PEPGetExtract(pep,extr,ierr))
    133:       if (rank .eq. 0) then
    134:         write(*,160) extr
    135:       endif
    136:  160  format (' ... changed to ',I2)

    138:       alpha = .1
    139:       its = 5
    140:       lambda = 1.
    141:       scal = PEP_SCALE_SCALAR
    142:       PetscCallA(PEPSetScale(pep,scal,alpha,PETSC_NULL_VEC,PETSC_NULL_VEC,its,lambda,ierr))
    143:       PetscCallA(PEPGetScale(pep,scal,alpha,PETSC_NULL_VEC,PETSC_NULL_VEC,its,lambda,ierr))
    144:       if (rank .eq. 0) then
    145:         write(*,170) scal,alpha,its
    146:       endif
    147:  170  format (' Scaling: ',I2,', alpha=',F7.4,', its=',I2)

    149:       basis = PEP_BASIS_CHEBYSHEV1
    150:       PetscCallA(PEPSetBasis(pep,basis,ierr))
    151:       PetscCallA(PEPGetBasis(pep,basis,ierr))
    152:       if (rank .eq. 0) then
    153:         write(*,180) basis
    154:       endif
    155:  180  format (' Polynomial basis: ',I2)

    157:       np = 1
    158:       tol = 1e-9
    159:       its = 2
    160:       refine = PEP_REFINE_SIMPLE
    161:       rscheme = PEP_REFINE_SCHEME_SCHUR
    162:       PetscCallA(PEPSetRefine(pep,refine,np,tol,its,rscheme,ierr))
    163:       PetscCallA(PEPGetRefine(pep,refine,np,tol,its,rscheme,ierr))
    164:       if (rank .eq. 0) then
    165:         write(*,190) refine,tol,its,rscheme
    166:       endif
    167:  190  format (' Refinement: ',I2,', tol=',F7.4,', its=',I2,', schem=',I2)

    169:       tget = 4.8
    170:       PetscCallA(PEPSetTarget(pep,tget,ierr))
    171:       PetscCallA(PEPGetTarget(pep,tget,ierr))
    172:       PetscCallA(PEPSetWhichEigenpairs(pep,PEP_TARGET_MAGNITUDE,ierr))
    173:       PetscCallA(PEPGetWhichEigenpairs(pep,which,ierr))
    174:       if (rank .eq. 0) then
    175:         write(*,200) which,PetscRealPart(tget)
    176:       endif
    177:  200  format (' Which = ',I2,', target = ',F4.1)

    179:       nev = 4
    180:       PetscCallA(PEPSetDimensions(pep,nev,PETSC_DETERMINE_INTEGER,PETSC_DETERMINE_INTEGER,ierr))
    181:       PetscCallA(PEPGetDimensions(pep,nev,ncv,mpd,ierr))
    182:       if (rank .eq. 0) then
    183:         write(*,210) nev,ncv,mpd
    184:       endif
    185:  210  format (' Dimensions: nev=',I2,', ncv=',I2,', mpd=',I2)

    187:       tol = 2.2e-4
    188:       its = 200
    189:       PetscCallA(PEPSetTolerances(pep,tol,its,ierr))
    190:       PetscCallA(PEPGetTolerances(pep,tol,its,ierr))
    191:       if (rank .eq. 0) then
    192:         write(*,220) tol,its
    193:       endif
    194:  220  format (' Tolerance =',F8.5,', max_its =',I4)

    196:       PetscCallA(PEPSetConvergenceTest(pep,PEP_CONV_ABS,ierr))
    197:       PetscCallA(PEPGetConvergenceTest(pep,conv,ierr))
    198:       PetscCallA(PEPSetStoppingTest(pep,PEP_STOP_BASIC,ierr))
    199:       PetscCallA(PEPGetStoppingTest(pep,stp,ierr))
    200:       if (rank .eq. 0) then
    201:         write(*,230) conv,stp
    202:       endif
    203:  230  format (' Convergence test =',I2,', stopping test =',I2)

    205:       PetscCallA(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,vf,ierr))
    206:       PetscCallA(PEPMonitorSet(pep,PEPMONITORFIRST,vf,PetscViewerAndFormatDestroy,ierr))
    207:       PetscCallA(PEPMonitorConvergedCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,PETSC_NULL_VEC,vf,ierr))
    208:       PetscCallA(PEPMonitorSet(pep,PEPMONITORCONVERGED,vf,PEPMonitorConvergedDestroy,ierr))
    209:       PetscCallA(PEPMonitorCancel(pep,ierr))

    211:       PetscCallA(PEPGetST(pep,st,ierr))
    212:       PetscCallA(STGetKSP(st,ksp,ierr))
    213:       tol = 1.e-8
    214:       tolabs = 1.e-35
    215:       PetscCallA(KSPSetTolerances(ksp,tol,tolabs,PETSC_CURRENT_REAL,PETSC_CURRENT_INTEGER,ierr))
    216:       PetscCallA(STView(st,PETSC_NULL_VIEWER,ierr))
    217:       PetscCallA(PEPGetDS(pep,ds,ierr))
    218:       PetscCallA(DSView(ds,PETSC_NULL_VIEWER,ierr))

    220:       PetscCallA(PEPSetFromOptions(pep,ierr))
    221:       PetscCallA(PEPSolve(pep,ierr))
    222:       PetscCallA(PEPGetConvergedReason(pep,reason,ierr))
    223:       if (rank .eq. 0) then
    224:         write(*,240) reason
    225:       endif
    226:  240  format (' Finished - converged reason =',I2)

    228: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    229: !     Display solution and clean up
    230: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    231:       PetscCallA(PEPErrorView(pep,PEP_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
    232:       PetscCallA(PEPDestroy(pep,ierr))
    233:       PetscCallA(MatDestroy(A(1),ierr))
    234:       PetscCallA(MatDestroy(A(2),ierr))
    235:       PetscCallA(MatDestroy(A(3),ierr))

    237:       PetscCallA(SlepcFinalize(ierr))
    238:       end

    240: !/*TEST
    241: !
    242: !   test:
    243: !      suffix: 1
    244: !      args: -pep_tol 1e-6 -pep_ncv 22
    245: !      filter: sed -e "s/[+-]0\.0*i//g"
    246: !
    247: !TEST*/