Actual source code: ex23f90.F90

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: !
  2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
  4: !  Copyright (c) 2002-2021, 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> ./ex23f90 [-help] [-t <t>] [-m <m>] [SLEPc opts]
 11: !
 12: !  Description: Computes exp(t*A)*v for a matrix associated with a
 13: !  Markov model. This is the Fortran90 equivalent to ex23.c
 14: !
 15: !  The command line options are:
 16: !    -t <t>, where <t> = time parameter (multiplies the matrix)
 17: !    -m <m>, where <m> = number of grid subdivisions in each dimension
 18: !
 19: ! ----------------------------------------------------------------------
 20: !
 21:       program main
 22: #include <slepc/finclude/slepcmfn.h>
 23:       use slepcmfn
 24:       implicit none

 26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 27: !     Declarations
 28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 29: !
 30: !  Variables:
 31: !     A      problem matrix
 32: !     mfn    matrix function solver context

 34:       Mat            A
 35:       MFN            mfn
 36:       FN             f
 37:       PetscReal      tol, norm, cst, pd, pu
 38:       PetscScalar    t, z
 39:       Vec            v, y
 40:       PetscInt       N, m, ncv, maxit, its, ii, jj
 41:       PetscInt       i, j, jmax, ix, Istart, Iend
 42:       PetscMPIInt    rank
 43:       PetscErrorCode ierr
 44:       PetscBool      flg
 45:       MFNConvergedReason reason

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

 51:       call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
 52:       if (ierr .ne. 0) then
 53:         print*,'SlepcInitialize failed'
 54:         stop
 55:       endif
 56:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr);CHKERRA(ierr)
 57:       m = 15
 58:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr);CHKERRA(ierr)
 59:       t = 2.0
 60:       call PetscOptionsGetScalar(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-t',t,flg,ierr);CHKERRA(ierr)
 61:       N = m*(m+1)/2
 62:       if (rank .eq. 0) then
 63:         write(*,100) N, m
 64:       endif
 65:  100  format (/'Markov y=exp(t*A)*e_1, N=',I6,' (m=',I4,')')

 67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 68: !     Compute the transition probability matrix, A
 69: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 71:       call MatCreate(PETSC_COMM_WORLD,A,ierr);CHKERRA(ierr)
 72:       call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr);CHKERRA(ierr)
 73:       call MatSetFromOptions(A,ierr);CHKERRA(ierr)
 74:       call MatSetUp(A,ierr);CHKERRA(ierr)
 75:       call MatGetOwnershipRange(A,Istart,Iend,ierr);CHKERRA(ierr)
 76:       ix = 0
 77:       cst = 0.5/real(m-1)
 78:       do i=1,m
 79:         jmax = m-i+1
 80:         do j=1,jmax
 81:           ix = ix + 1
 82:           ii = ix - 1
 83:           if (ix-1.ge.Istart .and. ix.le.Iend) then
 84:             if (j.ne.jmax) then
 85:               pd = cst*(i+j-1)
 86:               !** north
 87:               if (i.eq.1) then
 88:                 z = 2.0*pd
 89:               else
 90:                 z = pd
 91:               end if
 92:               call MatSetValue(A,ii,ix,z,INSERT_VALUES,ierr);CHKERRA(ierr)
 93:               !** east
 94:               if (j.eq.1) then
 95:                 z = 2.0*pd
 96:               else
 97:                 z = pd
 98:               end if
 99:               jj = ix+jmax-1
100:               call MatSetValue(A,ii,jj,z,INSERT_VALUES,ierr);CHKERRA(ierr)
101:             end if

103:             !** south
104:             pu = 0.5 - cst*(i+j-3)
105:             z = pu
106:             if (j.gt.1) then
107:               jj = ix-2
108:               call MatSetValue(A,ii,jj,z,INSERT_VALUES,ierr);CHKERRA(ierr)
109:             end if
110:             !** west
111:             if (i.gt.1) then
112:               jj = ix-jmax-2
113:               call MatSetValue(A,ii,jj,z,INSERT_VALUES,ierr);CHKERRA(ierr)
114:             end if
115:           end if
116:         end do
117:       end do
118:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
119:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)

121: !     ** Set v = e_1
122:       call MatCreateVecs(A,y,v,ierr);CHKERRA(ierr)
123:       ii = 0
124:       z = 1.0
125:       call VecSetValue(v,ii,z,INSERT_VALUES,ierr);CHKERRA(ierr)
126:       call VecAssemblyBegin(v,ierr);CHKERRA(ierr)
127:       call VecAssemblyEnd(v,ierr);CHKERRA(ierr)

129: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: !     Create the solver and set various options
131: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

133: !     ** Create matrix function solver context
134:       call MFNCreate(PETSC_COMM_WORLD,mfn,ierr);CHKERRA(ierr)

136: !     ** Set operator matrix, the function to compute, and other options
137:       call MFNSetOperator(mfn,A,ierr);CHKERRA(ierr)
138:       call MFNGetFN(mfn,f,ierr);CHKERRA(ierr)
139:       call FNSetType(f,FNEXP,ierr);CHKERRA(ierr)
140:       z = 1.0
141:       call FNSetScale(f,t,z,ierr);CHKERRA(ierr)
142:       tol = 0.0000001
143:       call MFNSetTolerances(mfn,tol,PETSC_DEFAULT_INTEGER,ierr);CHKERRA(ierr)

145: !     ** Set solver parameters at runtime
146:       call MFNSetFromOptions(mfn,ierr);CHKERRA(ierr)

148: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: !     Solve the problem, y=exp(t*A)*v
150: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

152:       call MFNSolve(mfn,v,y,ierr);CHKERRA(ierr)
153:       call MFNGetConvergedReason(mfn,reason,ierr);CHKERRA(ierr)
154:       if (reason.lt.0) then; SETERRA(PETSC_COMM_WORLD,1,'Solver did not converge'); endif
155:       call VecNorm(y,NORM_2,norm,ierr);CHKERRA(ierr)

157: !     ** Optional: Get some information from the solver and display it
158:       call MFNGetIterationNumber(mfn,its,ierr);CHKERRA(ierr)
159:       if (rank .eq. 0) then
160:         write(*,120) its
161:       endif
162:  120  format (' Number of iterations of the method: ',I4)
163:       call MFNGetDimensions(mfn,ncv,ierr);CHKERRA(ierr)
164:       if (rank .eq. 0) then
165:         write(*,130) ncv
166:       endif
167:  130  format (' Subspace dimension:',I4)
168:       call MFNGetTolerances(mfn,tol,maxit,ierr);CHKERRA(ierr)
169:       if (rank .eq. 0) then
170:         write(*,140) tol,maxit
171:       endif
172:  140  format (' Stopping condition: tol=',f10.7,' maxit=',I4)

174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: !     Display solution and clean up
176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

178:       if (rank .eq. 0) then
179:         write(*,150) PetscRealPart(t),norm
180:       endif
181:  150  format (' Computed vector at time t=',f4.1,' has norm ',f8.5)

183:       call MFNDestroy(mfn,ierr);CHKERRA(ierr)
184:       call MatDestroy(A,ierr);CHKERRA(ierr)
185:       call VecDestroy(v,ierr);CHKERRA(ierr)
186:       call VecDestroy(y,ierr);CHKERRA(ierr)
187:       call SlepcFinalize(ierr)
188:       end

190: !/*TEST
191: !
192: !   test:
193: !      suffix: 1
194: !      args: -mfn_ncv 6
195: !
196: !TEST*/