Actual source code: ex23f.F90
slepc-3.21.1 2024-04-26
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: ! Program usage: mpiexec -n <np> ./ex23f [-help] [-t <t>] [-m <m>] [SLEPc opts]
11: !
12: ! Description: Computes exp(t*A)*v for a matrix from a Markov model.
13: ! 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: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
52: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
53: m = 15
54: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr))
55: t = 2.0
56: PetscCallA(PetscOptionsGetScalar(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-t',t,flg,ierr))
57: N = m*(m+1)/2
58: if (rank .eq. 0) then
59: write(*,100) N, m
60: endif
61: 100 format (/'Markov y=exp(t*A)*e_1, N=',I6,' (m=',I4,')')
63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: ! Compute the transition probability matrix, A
65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
68: PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr))
69: PetscCallA(MatSetFromOptions(A,ierr))
70: PetscCallA(MatGetOwnershipRange(A,Istart,Iend,ierr))
71: ix = 0
72: cst = 0.5/real(m-1)
73: do i=1,m
74: jmax = m-i+1
75: do j=1,jmax
76: ix = ix + 1
77: ii = ix - 1
78: if (ix-1.ge.Istart .and. ix.le.Iend) then
79: if (j.ne.jmax) then
80: pd = cst*(i+j-1)
81: !** north
82: if (i.eq.1) then
83: z = 2.0*pd
84: else
85: z = pd
86: end if
87: PetscCallA(MatSetValue(A,ii,ix,z,INSERT_VALUES,ierr))
88: !** east
89: if (j.eq.1) then
90: z = 2.0*pd
91: else
92: z = pd
93: end if
94: jj = ix+jmax-1
95: PetscCallA(MatSetValue(A,ii,jj,z,INSERT_VALUES,ierr))
96: end if
98: !** south
99: pu = 0.5 - cst*(i+j-3)
100: z = pu
101: if (j.gt.1) then
102: jj = ix-2
103: PetscCallA(MatSetValue(A,ii,jj,z,INSERT_VALUES,ierr))
104: end if
105: !** west
106: if (i.gt.1) then
107: jj = ix-jmax-2
108: PetscCallA(MatSetValue(A,ii,jj,z,INSERT_VALUES,ierr))
109: end if
110: end if
111: end do
112: end do
113: PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
114: PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))
116: ! ** Set v = e_1
117: PetscCallA(MatCreateVecs(A,y,v,ierr))
118: ii = 0
119: z = 1.0
120: PetscCallA(VecSetValue(v,ii,z,INSERT_VALUES,ierr))
121: PetscCallA(VecAssemblyBegin(v,ierr))
122: PetscCallA(VecAssemblyEnd(v,ierr))
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: ! Create the solver and set various options
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: ! ** Create matrix function solver context
129: PetscCallA(MFNCreate(PETSC_COMM_WORLD,mfn,ierr))
131: ! ** Set operator matrix, the function to compute, and other options
132: PetscCallA(MFNSetOperator(mfn,A,ierr))
133: PetscCallA(MFNGetFN(mfn,f,ierr))
134: PetscCallA(FNSetType(f,FNEXP,ierr))
135: z = 1.0
136: PetscCallA(FNSetScale(f,t,z,ierr))
137: tol = 0.0000001
138: PetscCallA(MFNSetTolerances(mfn,tol,PETSC_DEFAULT_INTEGER,ierr))
140: ! ** Set solver parameters at runtime
141: PetscCallA(MFNSetFromOptions(mfn,ierr))
143: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: ! Solve the problem, y=exp(t*A)*v
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: PetscCallA(MFNSolve(mfn,v,y,ierr))
148: PetscCallA(MFNGetConvergedReason(mfn,reason,ierr))
149: if (reason.lt.0) then; SETERRA(PETSC_COMM_WORLD,1,'Solver did not converge'); endif
150: PetscCallA(VecNorm(y,NORM_2,norm,ierr))
152: ! ** Optional: Get some information from the solver and display it
153: PetscCallA(MFNGetIterationNumber(mfn,its,ierr))
154: if (rank .eq. 0) then
155: write(*,120) its
156: endif
157: 120 format (' Number of iterations of the method: ',I4)
158: PetscCallA(MFNGetDimensions(mfn,ncv,ierr))
159: if (rank .eq. 0) then
160: write(*,130) ncv
161: endif
162: 130 format (' Subspace dimension:',I4)
163: PetscCallA(MFNGetTolerances(mfn,tol,maxit,ierr))
164: if (rank .eq. 0) then
165: write(*,140) tol,maxit
166: endif
167: 140 format (' Stopping condition: tol=',f10.7,' maxit=',I4)
169: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170: ! Display solution and clean up
171: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: if (rank .eq. 0) then
174: write(*,150) PetscRealPart(t),norm
175: endif
176: 150 format (' Computed vector at time t=',f4.1,' has norm ',f8.5)
178: PetscCallA(MFNDestroy(mfn,ierr))
179: PetscCallA(MatDestroy(A,ierr))
180: PetscCallA(VecDestroy(v,ierr))
181: PetscCallA(VecDestroy(y,ierr))
182: PetscCallA(SlepcFinalize(ierr))
183: end
185: !/*TEST
186: !
187: ! test:
188: ! suffix: 1
189: ! args: -mfn_ncv 6
190: !
191: !TEST*/