Actual source code: ex23f90.F90
slepc-3.18.3 2023-03-24
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> ./ex23f90 [-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(MatSetUp(A,ierr))
71: PetscCallA(MatGetOwnershipRange(A,Istart,Iend,ierr))
72: ix = 0
73: cst = 0.5/real(m-1)
74: do i=1,m
75: jmax = m-i+1
76: do j=1,jmax
77: ix = ix + 1
78: ii = ix - 1
79: if (ix-1.ge.Istart .and. ix.le.Iend) then
80: if (j.ne.jmax) then
81: pd = cst*(i+j-1)
82: !** north
83: if (i.eq.1) then
84: z = 2.0*pd
85: else
86: z = pd
87: end if
88: PetscCallA(MatSetValue(A,ii,ix,z,INSERT_VALUES,ierr))
89: !** east
90: if (j.eq.1) then
91: z = 2.0*pd
92: else
93: z = pd
94: end if
95: jj = ix+jmax-1
96: PetscCallA(MatSetValue(A,ii,jj,z,INSERT_VALUES,ierr))
97: end if
99: !** south
100: pu = 0.5 - cst*(i+j-3)
101: z = pu
102: if (j.gt.1) then
103: jj = ix-2
104: PetscCallA(MatSetValue(A,ii,jj,z,INSERT_VALUES,ierr))
105: end if
106: !** west
107: if (i.gt.1) then
108: jj = ix-jmax-2
109: PetscCallA(MatSetValue(A,ii,jj,z,INSERT_VALUES,ierr))
110: end if
111: end if
112: end do
113: end do
114: PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
115: PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))
117: ! ** Set v = e_1
118: PetscCallA(MatCreateVecs(A,y,v,ierr))
119: ii = 0
120: z = 1.0
121: PetscCallA(VecSetValue(v,ii,z,INSERT_VALUES,ierr))
122: PetscCallA(VecAssemblyBegin(v,ierr))
123: PetscCallA(VecAssemblyEnd(v,ierr))
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: ! Create the solver and set various options
127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: ! ** Create matrix function solver context
130: PetscCallA(MFNCreate(PETSC_COMM_WORLD,mfn,ierr))
132: ! ** Set operator matrix, the function to compute, and other options
133: PetscCallA(MFNSetOperator(mfn,A,ierr))
134: PetscCallA(MFNGetFN(mfn,f,ierr))
135: PetscCallA(FNSetType(f,FNEXP,ierr))
136: z = 1.0
137: PetscCallA(FNSetScale(f,t,z,ierr))
138: tol = 0.0000001
139: PetscCallA(MFNSetTolerances(mfn,tol,PETSC_DEFAULT_INTEGER,ierr))
141: ! ** Set solver parameters at runtime
142: PetscCallA(MFNSetFromOptions(mfn,ierr))
144: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: ! Solve the problem, y=exp(t*A)*v
146: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: PetscCallA(MFNSolve(mfn,v,y,ierr))
149: PetscCallA(MFNGetConvergedReason(mfn,reason,ierr))
150: if (reason.lt.0) then; SETERRA(PETSC_COMM_WORLD,1,'Solver did not converge'); endif
151: PetscCallA(VecNorm(y,NORM_2,norm,ierr))
153: ! ** Optional: Get some information from the solver and display it
154: PetscCallA(MFNGetIterationNumber(mfn,its,ierr))
155: if (rank .eq. 0) then
156: write(*,120) its
157: endif
158: 120 format (' Number of iterations of the method: ',I4)
159: PetscCallA(MFNGetDimensions(mfn,ncv,ierr))
160: if (rank .eq. 0) then
161: write(*,130) ncv
162: endif
163: 130 format (' Subspace dimension:',I4)
164: PetscCallA(MFNGetTolerances(mfn,tol,maxit,ierr))
165: if (rank .eq. 0) then
166: write(*,140) tol,maxit
167: endif
168: 140 format (' Stopping condition: tol=',f10.7,' maxit=',I4)
170: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: ! Display solution and clean up
172: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: if (rank .eq. 0) then
175: write(*,150) PetscRealPart(t),norm
176: endif
177: 150 format (' Computed vector at time t=',f4.1,' has norm ',f8.5)
179: PetscCallA(MFNDestroy(mfn,ierr))
180: PetscCallA(MatDestroy(A,ierr))
181: PetscCallA(VecDestroy(v,ierr))
182: PetscCallA(VecDestroy(y,ierr))
183: PetscCallA(SlepcFinalize(ierr))
184: end
186: !/*TEST
187: !
188: ! test:
189: ! suffix: 1
190: ! args: -mfn_ncv 6
191: !
192: !TEST*/