Actual source code: test3f.F

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: !
  2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
  4: !  Copyright (c) 2002-2019, 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> ./test3f [-help] [-n <n>] [all SLEPc options]
 11: !
 12: !  Description: square root of the 2-D Laplacian.
 13: !
 14: !  The command line options are:
 15: !    -n <n>, where <n> = matrix rows and columns
 16: !
 17: ! ----------------------------------------------------------------------
 18: !
 19:       program main
 20: #include <slepc/finclude/slepcmfn.h>
 21:       use slepcmfn
 22:       implicit none

 24: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 25: !     Declarations
 26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 27: !
 28:       Mat                A, B
 29:       MFN                mfn
 30:       FN                 f
 31:       MFNConvergedReason reason;
 32:       Vec                v, y
 33:       PetscInt           Nt, n, i, j, II
 34:       PetscInt           Istart, maxit, ncv
 35:       PetscInt           col, its, Iend
 36:       PetscScalar        val
 37:       PetscReal          tol, norm
 38:       PetscMPIInt        rank
 39:       PetscErrorCode     ierr
 40:       PetscBool          flg

 42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43: !     Beginning of program
 44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 46:       call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
 47:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 48:       n = 4
 49:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,  &
 50:      &                        '-n',n,flg,ierr)
 51:       Nt = n*n

 53:       if (rank .eq. 0) then
 54:         write(*,100) n
 55:       endif
 56:  100  format (/'nSquare root of Laplacian, n=',I3,' (Fortran)')

 58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 59: !     Compute the discrete 2-D Laplacian
 60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 62:       call MatCreate(PETSC_COMM_WORLD,A,ierr)
 63:       call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,Nt,Nt,ierr)
 64:       call MatSetFromOptions(A,ierr)
 65:       call MatSetUp(A,ierr)

 67:       call MatGetOwnershipRange(A,Istart,Iend,ierr)
 68:       do II=Istart,Iend-1
 69:         i = II/n
 70:         j = II-i*n
 71:         val = -1.0
 72:         if (i .gt. 0) then
 73:           col = II-n
 74:           call MatSetValue(A,II,col,val,INSERT_VALUES,ierr)
 75:         end if
 76:         if (i .lt. n-1) then
 77:           col = II+n
 78:           call MatSetValue(A,II,col,val,INSERT_VALUES,ierr)
 79:         end if
 80:         if (j .gt. 0) then
 81:           col = II-1
 82:           call MatSetValue(A,II,col,val,INSERT_VALUES,ierr)
 83:         end if
 84:         if (j .lt. n-1) then
 85:           col = II+1
 86:           call MatSetValue(A,II,col,val,INSERT_VALUES,ierr)
 87:         end if
 88:         val = 4.0
 89:         call MatSetValue(A,II,II,val,INSERT_VALUES,ierr)
 90:       enddo

 92:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
 93:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

 95:       call MatCreateVecs(A,PETSC_NULL_VEC,v,ierr)
 96:       i = 0
 97:       val = 1.0
 98:       call VecSetValue(v,i,val,INSERT_VALUES,ierr)
 99:       call VecAssemblyBegin(v,ierr)
100:       call VecAssemblyEnd(v,ierr)
101:       call VecDuplicate(v,y,ierr)

103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: !     Compute singular values
105: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

107:       call MFNCreate(PETSC_COMM_WORLD,mfn,ierr)
108:       call MFNSetOperator(mfn,A,ierr)
109:       call MFNGetFN(mfn,f,ierr)
110:       call FNSetType(f,FNSQRT,ierr)

112: !     ** test some interface functions
113:       call MFNGetOperator(mfn,B,ierr)
114:       call MatView(B,PETSC_VIEWER_STDOUT_WORLD,ierr)
115:       call MFNSetOptionsPrefix(mfn,'myprefix_',ierr)
116:       tol = 1e-4
117:       maxit = 500
118:       call MFNSetTolerances(mfn,tol,maxit,ierr)
119:       ncv = 6
120:       call MFNSetDimensions(mfn,ncv,ierr)
121:       call MFNSetErrorIfNotConverged(mfn,PETSC_TRUE,ierr)
122:       call MFNSetFromOptions(mfn,ierr)

124: !     ** query properties and print them
125:       call MFNGetTolerances(mfn,tol,maxit,ierr)
126:       if (rank .eq. 0) then
127:         write(*,110) tol,maxit
128:       endif
129:  110  format (/' Tolerance: ',F7.4,', maxit: ',I4)
130:       call MFNGetDimensions(mfn,ncv,ierr)
131:       if (rank .eq. 0) then
132:         write(*,120) ncv
133:       endif
134:  120  format (' Subspace dimension: ',I3)
135:       call MFNGetErrorIfNotConverged(mfn,flg,ierr)
136:       if (rank .eq. 0 .and. flg) then
137:         write(*,*) 'Erroring out if convergence fails'
138:       endif

140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: !     Call the solver
142: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

144:       call MFNSolve(mfn,v,y,ierr)
145:       call MFNGetConvergedReason(mfn,reason,ierr)
146:       if (rank .eq. 0) then
147:         write(*,130) reason
148:       endif
149:  130  format (' Converged reason:',I2)
150:       call MFNGetIterationNumber(mfn,its,ierr)
151: !     if (rank .eq. 0) then
152: !       write(*,140) its
153: !     endif
154: !140  format (' Number of iterations of the method:',I4)

156:       call VecNorm(y,NORM_2,norm,ierr)
157:       if (rank .eq. 0) then
158:         write(*,150) norm
159:       endif
160:  150  format (' sqrt(A)*v has norm ',F7.4)

162:       call MFNDestroy(mfn,ierr)
163:       call MatDestroy(A,ierr)
164:       call VecDestroy(v,ierr)
165:       call VecDestroy(y,ierr)

167:       call SlepcFinalize(ierr)
168:       end

170: !/*TEST
171: !
172: !   test:
173: !      suffix: 1
174: !      args: -info_exclude mfn -log_exclude mfn
175: !
176: !TEST*/