Actual source code: test17f.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: Test Fortran interface of spectrum-slicing Krylov-Schur.
     11: !
     12: ! ----------------------------------------------------------------------
     13: !
     14:       program main
     15: #include <slepc/finclude/slepceps.h>
     16:       use slepceps
     17:       implicit none

     19: #define MAXSUB 16
     20: #define MAXSHI 16

     22: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     23: !     Declarations
     24: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     25:       Mat            A,B,As,Bs,Au
     26:       EPS            eps
     27:       ST             st
     28:       KSP            ksp
     29:       PC             pc
     30:       Vec            v
     31:       PetscScalar    value
     32:       PetscInt       n,m,i,j,k,Istart,Iend
     33:       PetscInt       nev,ncv,mpd,nval
     34:       PetscInt       row,col,nloc,nlocs,mlocs
     35:       PetscInt       II,npart,inertias(MAXSHI)
     36:       PetscBool      flg,lock
     37:       PetscMPIInt    nprc,rank
     38:       PetscReal      int0,int1,keep,subint(MAXSUB)
     39:       PetscReal      shifts(MAXSHI)
     40:       PetscScalar    eval,one,mone,zero
     41:       PetscErrorCode ierr
     42:       MPI_Comm       comm

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

     48:       PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
     49:       if (ierr .ne. 0) then
     50:         print*,'SlepcInitialize failed'
     51:         stop
     52:       endif
     53:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,nprc,ierr))
     54:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
     55:       n = 35
     56:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))
     57:       m = n*n
     58:       if (rank .eq. 0) then
     59:         write(*,100) n
     60:       endif
     61:  100  format (/'Spectrum-slicing test, n =',I3,' (Fortran)'/)

     63:       PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
     64:       PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr))
     65:       PetscCallA(MatSetFromOptions(A,ierr))
     66:       PetscCallA(MatCreate(PETSC_COMM_WORLD,B,ierr))
     67:       PetscCallA(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr))
     68:       PetscCallA(MatSetFromOptions(B,ierr))
     69:       PetscCallA(MatGetOwnershipRange(A,Istart,Iend,ierr))
     70:       do II=Istart,Iend-1
     71:         i = II/n
     72:         j = II-i*n
     73:         value = -1.0
     74:         row = II
     75:         if (i>0) then
     76:           col = II-n
     77:           PetscCallA(MatSetValue(A,row,col,value,INSERT_VALUES,ierr))
     78:         endif
     79:         if (i<n-1) then
     80:           col = II+n
     81:           PetscCallA(MatSetValue(A,row,col,value,INSERT_VALUES,ierr))
     82:         endif
     83:         if (j>0) then
     84:           col = II-1
     85:           PetscCallA(MatSetValue(A,row,col,value,INSERT_VALUES,ierr))
     86:         endif
     87:         if (j<n-1) then
     88:           col = II+1
     89:           PetscCallA(MatSetValue(A,row,col,value,INSERT_VALUES,ierr))
     90:         endif
     91:         col = II
     92:         value = 4.0
     93:         PetscCallA(MatSetValue(A,row,col,value,INSERT_VALUES,ierr))
     94:         value = 2.0
     95:         PetscCallA(MatSetValue(B,row,col,value,INSERT_VALUES,ierr))
     96:       enddo
     97:       if (Istart .eq. 0) then
     98:         row = 0
     99:         col = 0
    100:         value = 6.0
    101:         PetscCallA(MatSetValue(B,row,col,value,INSERT_VALUES,ierr))
    102:         row = 0
    103:         col = 1
    104:         value = -1.0
    105:         PetscCallA(MatSetValue(B,row,col,value,INSERT_VALUES,ierr))
    106:         row = 1
    107:         col = 0
    108:         value = -1.0
    109:         PetscCallA(MatSetValue(B,row,col,value,INSERT_VALUES,ierr))
    110:         row = 1
    111:         col = 1
    112:         value = 1.0
    113:         PetscCallA(MatSetValue(B,row,col,value,INSERT_VALUES,ierr))
    114:       endif
    115:       PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
    116:       PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))
    117:       PetscCallA(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
    118:       PetscCallA(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))

    120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    121: !     Create eigensolver and set various options
    122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    124:       PetscCallA(EPSCreate(PETSC_COMM_WORLD,eps,ierr))
    125:       PetscCallA(EPSSetOperators(eps,A,B,ierr))
    126:       PetscCallA(EPSSetProblemType(eps,EPS_GHEP,ierr))
    127:       PetscCallA(EPSSetType(eps,EPSKRYLOVSCHUR,ierr))

    129: !     Set interval and other settings for spectrum slicing

    131:       PetscCallA(EPSSetWhichEigenpairs(eps,EPS_ALL,ierr))
    132:       int0 = 1.1
    133:       int1 = 1.3
    134:       PetscCallA(EPSSetInterval(eps,int0,int1,ierr))
    135:       PetscCallA(EPSGetST(eps,st,ierr))
    136:       PetscCallA(STSetType(st,STSINVERT,ierr))
    137:       if (nprc>0) then
    138:         npart = nprc
    139:         PetscCallA(EPSKrylovSchurSetPartitions(eps,npart,ierr))
    140:       endif
    141:       PetscCallA(EPSKrylovSchurGetKSP(eps,ksp,ierr))
    142:       PetscCallA(KSPGetPC(ksp,pc,ierr))
    143:       PetscCallA(KSPSetType(ksp,KSPPREONLY,ierr))
    144:       PetscCallA(PCSetType(pc,PCCHOLESKY,ierr))

    146: !     Test interface functions of Krylov-Schur solver

    148:       PetscCallA(EPSKrylovSchurGetRestart(eps,keep,ierr))
    149:       if (rank .eq. 0) then
    150:         write(*,110) keep
    151:       endif
    152:  110  format (' Restart parameter before changing = ',f7.4)
    153:       keep = 0.4
    154:       PetscCallA(EPSKrylovSchurSetRestart(eps,keep,ierr))
    155:       PetscCallA(EPSKrylovSchurGetRestart(eps,keep,ierr))
    156:       if (rank .eq. 0) then
    157:         write(*,120) keep
    158:       endif
    159:  120  format (' ... changed to ',f7.4)

    161:       PetscCallA(EPSKrylovSchurGetLocking(eps,lock,ierr))
    162:       if (rank .eq. 0) then
    163:         write(*,130) lock
    164:       endif
    165:  130  format (' Locking flag before changing = ',L4)
    166:       PetscCallA(EPSKrylovSchurSetLocking(eps,PETSC_FALSE,ierr))
    167:       PetscCallA(EPSKrylovSchurGetLocking(eps,lock,ierr))
    168:       if (rank .eq. 0) then
    169:         write(*,140) lock
    170:       endif
    171:  140  format (' ... changed to ',L4)

    173:       PetscCallA(EPSKrylovSchurGetDimensions(eps,nev,ncv,mpd,ierr))
    174:       if (rank .eq. 0) then
    175:         write(*,150) nev,ncv,mpd
    176:       endif
    177:  150  format (' Sub-solve dimensions before changing: nev=',I2,', ncv=',I2,', mpd=',I2)
    178:       nev = 30
    179:       ncv = 60
    180:       mpd = 60
    181:       PetscCallA(EPSKrylovSchurSetDimensions(eps,nev,ncv,mpd,ierr))
    182:       PetscCallA(EPSKrylovSchurGetDimensions(eps,nev,ncv,mpd,ierr))
    183:       if (rank .eq. 0) then
    184:         write(*,160) nev,ncv,mpd
    185:       endif
    186:  160  format (' ... changed to: nev=',I2,', ncv=',I2,', mpd=',I2)

    188:       if (nprc>0) then
    189:         PetscCallA(EPSKrylovSchurGetPartitions(eps,npart,ierr))
    190:         if (rank .eq. 0) then
    191:           write(*,170) npart
    192:         endif
    193:  170    format (' Using ',I2,' partitions')
    194:         if (npart>MAXSUB) then; SETERRA(PETSC_COMM_SELF,1,'Too many subintervals'); endif

    196:         subint(1) = int0
    197:         subint(npart+1) = int1
    198:         do i=2,npart
    199:           subint(i) = int0+(i-1)*(int1-int0)/npart
    200:         enddo
    201:         PetscCallA(EPSKrylovSchurSetSubintervals(eps,subint,ierr))
    202:         PetscCallA(EPSKrylovSchurGetSubintervals(eps,subint,ierr))
    203:         if (rank .eq. 0) then
    204:           write(*,*) 'Using sub-interval separations ='
    205:           do i=2,npart
    206:             write(*,180) subint(i)
    207:           enddo
    208:         endif
    209:  180    format (f7.4)
    210:       endif

    212:       PetscCallA(EPSSetFromOptions(eps,ierr))

    214: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    215: !     Compute all eigenvalues in interval and display info
    216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    218:       PetscCallA(EPSSetUp(eps,ierr))
    219:       PetscCallA(EPSKrylovSchurGetInertias(eps,k,PETSC_NULL_REAL_ARRAY,PETSC_NULL_INTEGER_ARRAY,ierr))
    220:       if (k>MAXSHI) then; SETERRA(PETSC_COMM_SELF,1,'Too many shifts'); endif
    221:       PetscCallA(EPSKrylovSchurGetInertias(eps,k,shifts,inertias,ierr))
    222:       if (rank .eq. 0) then
    223:         write(*,*) 'Inertias after EPSSetUp:'
    224:         do i=1,k
    225:           write(*,185) shifts(i),inertias(i)
    226:         enddo
    227:       endif
    228:  185  format (' .. ',f4.1,' (',I3,')')

    230:       PetscCallA(EPSSolve(eps,ierr))
    231:       PetscCallA(EPSGetDimensions(eps,nev,ncv,mpd,ierr))
    232:       PetscCallA(EPSGetInterval(eps,int0,int1,ierr))
    233:       if (rank .eq. 0) then
    234:         write(*,190) nev,int0,int1
    235:       endif
    236:  190  format (' Found ',I2,' eigenvalues in interval [',f7.4,',',f7.4,']')

    238:       if (nprc>0) then
    239:         PetscCallA(EPSKrylovSchurGetSubcommInfo(eps,k,nval,v,ierr))
    240:         if (rank .eq. 0) then
    241:           write(*,200) rank,k,nval
    242:           do i=0,nval-1
    243:             PetscCallA(EPSKrylovSchurGetSubcommPairs(eps,i,eval,v,ierr))
    244:             write(*,210) PetscRealPart(eval)
    245:           enddo
    246:         endif
    247:  200    format (' Process ',I2,' has worked in sub-interval ',I2,', containing ',I2,' eigenvalues')
    248:  210    format (f7.4)
    249:         PetscCallA(VecDestroy(v,ierr))

    251:         PetscCallA(EPSKrylovSchurGetSubcommMats(eps,As,Bs,ierr))
    252:         PetscCallA(MatGetLocalSize(A,nloc,PETSC_NULL_INTEGER,ierr))
    253:         PetscCallA(MatGetLocalSize(As,nlocs,mlocs,ierr))
    254:         if (rank .eq. 0) then
    255:           write(*,220) rank,nloc,nlocs
    256:         endif
    257:  220    format (' Process ',I2,' owns ',I5,', rows of the global',' matrices, and ',I5,' rows in the subcommunicator')

    259: !       modify A on subcommunicators
    260:         PetscCallA(PetscObjectGetComm(As,comm,ierr))
    261:         PetscCallA(MatCreate(comm,Au,ierr))
    262:         PetscCallA(MatSetSizes(Au,nlocs,mlocs,m,m,ierr))
    263:         PetscCallA(MatSetFromOptions(Au,ierr))
    264:         PetscCallA(MatGetOwnershipRange(Au,Istart,Iend,ierr))
    265:         do II=Istart,Iend-1
    266:           value = 0.5
    267:           PetscCallA(MatSetValue(Au,II,II,value,INSERT_VALUES,ierr))
    268:         end do
    269:         PetscCallA(MatAssemblyBegin(Au,MAT_FINAL_ASSEMBLY,ierr))
    270:         PetscCallA(MatAssemblyEnd(Au,MAT_FINAL_ASSEMBLY,ierr))
    271:         one = 1.0
    272:         mone = -1.0
    273:         zero = 0.0
    274:         PetscCallA(EPSKrylovSchurUpdateSubcommMats(eps,one,mone,Au,zero,zero,PETSC_NULL_MAT,DIFFERENT_NONZERO_PATTERN,PETSC_TRUE,ierr))
    275:         PetscCallA(MatDestroy(Au,ierr))
    276:       endif

    278:       PetscCallA(EPSDestroy(eps,ierr))
    279:       PetscCallA(MatDestroy(A,ierr))
    280:       PetscCallA(MatDestroy(B,ierr))

    282:       PetscCallA(SlepcFinalize(ierr))
    283:       end

    285: !/*TEST
    286: !
    287: !   test:
    288: !      suffix: 1
    289: !      nsize: 2
    290: !
    291: !TEST*/