Actual source code: ex23f90.F90
slepc-3.11.2 2019-07-30
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> ./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*/