Actual source code: ex23f.F90
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: #include <slepc/finclude/slepcmfn.h>
22: program ex23f
23: use slepcmfn
24: implicit none
26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
27: ! Declarations
28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30: Mat :: A ! problem matrix
31: MFN :: mfn ! matrix function solver context
32: FN :: f
33: PetscReal :: tol, norm, cst, pd, pu
34: PetscScalar :: t, z
35: Vec :: v, y
36: PetscInt :: N, m, ncv, maxit, its, ii, jj
37: PetscInt :: i, j, jmax, ix, Istart, Iend
38: PetscMPIInt :: rank
39: PetscErrorCode :: ierr
40: PetscBool :: flg
41: MFNConvergedReason :: reason
42: PetscScalar, parameter :: one = 1.0
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: ! Beginning of program
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
49: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
50: m = 15
51: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
52: t = 2.0
53: PetscCallA(PetscOptionsGetScalar(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-t', t, flg, ierr))
54: N = m*(m + 1)/2
55: if (rank == 0) then
56: write (*, '(/a,i6,a,i4,a)') 'Markov y=exp(t*A)*e_1, N=', N, ' (m=', m, ')'
57: end if
59: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: ! Compute the transition probability matrix, A
61: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
64: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N, ierr))
65: PetscCallA(MatSetFromOptions(A, ierr))
66: PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
67: ix = 0
68: cst = 0.5/real(m - 1, PETSC_REAL_KIND)
69: do i = 1, m
70: jmax = m - i + 1
71: do j = 1, jmax
72: ix = ix + 1
73: ii = ix - 1
74: if (ix - 1 >= Istart .and. ix <= Iend) then
75: if (j /= jmax) then
76: pd = cst*(i + j - 1)
77: !** north
78: if (i == 1) then
79: z = 2.0*pd
80: else
81: z = pd
82: end if
83: PetscCallA(MatSetValue(A, ii, ix, z, INSERT_VALUES, ierr))
84: !** east
85: if (j == 1) then
86: z = 2.0*pd
87: else
88: z = pd
89: end if
90: jj = ix + jmax - 1
91: PetscCallA(MatSetValue(A, ii, jj, z, INSERT_VALUES, ierr))
92: end if
94: !** south
95: pu = 0.5 - cst*(i + j - 3)
96: z = pu
97: if (j > 1) then
98: jj = ix - 2
99: PetscCallA(MatSetValue(A, ii, jj, z, INSERT_VALUES, ierr))
100: end if
101: !** west
102: if (i > 1) then
103: jj = ix - jmax - 2
104: PetscCallA(MatSetValue(A, ii, jj, z, INSERT_VALUES, ierr))
105: end if
106: end if
107: end do
108: end do
109: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
110: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
112: ! ** Set v = e_1
113: PetscCallA(MatCreateVecs(A, y, v, ierr))
114: ii = 0
115: z = 1.0
116: PetscCallA(VecSetValue(v, ii, z, INSERT_VALUES, ierr))
117: PetscCallA(VecAssemblyBegin(v, ierr))
118: PetscCallA(VecAssemblyEnd(v, ierr))
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: ! Create the solver and set various options
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: ! ** Create matrix function solver context
125: PetscCallA(MFNCreate(PETSC_COMM_WORLD, mfn, ierr))
127: ! ** Set operator matrix, the function to compute, and other options
128: PetscCallA(MFNSetOperator(mfn, A, ierr))
129: PetscCallA(MFNGetFN(mfn, f, ierr))
130: PetscCallA(FNSetType(f, FNEXP, ierr))
131: PetscCallA(FNSetScale(f, t, one, ierr))
132: tol = 0.0000001
133: PetscCallA(MFNSetTolerances(mfn, tol, PETSC_CURRENT_INTEGER, ierr))
135: ! ** Set solver parameters at runtime
136: PetscCallA(MFNSetFromOptions(mfn, ierr))
138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: ! Solve the problem, y=exp(t*A)*v
140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: PetscCallA(MFNSolve(mfn, v, y, ierr))
143: PetscCallA(MFNGetConvergedReason(mfn, reason, ierr))
144: PetscCheckA(reason%v >= 0, PETSC_COMM_WORLD, PETSC_ERR_NOT_CONVERGED, 'Solver did not converge')
145: PetscCallA(VecNorm(y, NORM_2, norm, ierr))
147: ! ** Optional: Get some information from the solver and display it
148: PetscCallA(MFNGetIterationNumber(mfn, its, ierr))
149: if (rank == 0) then
150: write (*, '(a,i4)') ' Number of iterations of the method: ', its
151: end if
152: PetscCallA(MFNGetDimensions(mfn, ncv, ierr))
153: if (rank == 0) then
154: write (*, '(a,i4)') ' Subspace dimension:', ncv
155: end if
156: PetscCallA(MFNGetTolerances(mfn, tol, maxit, ierr))
157: if (rank == 0) then
158: write (*, '(a,f10.7,a,i4)') ' Stopping condition: tol=', tol, ' maxit=', maxit
159: end if
161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: ! Display solution and clean up
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: if (rank == 0) then
166: write (*, '(a,f4.1,a,f8.5)') ' Computed vector at time t=', PetscRealPart(t), ' has norm ', norm
167: end if
169: PetscCallA(MFNDestroy(mfn, ierr))
170: PetscCallA(MatDestroy(A, ierr))
171: PetscCallA(VecDestroy(v, ierr))
172: PetscCallA(VecDestroy(y, ierr))
173: PetscCallA(SlepcFinalize(ierr))
174: end program ex23f
176: !/*TEST
177: !
178: ! test:
179: ! suffix: 1
180: ! args: -mfn_ncv 6
181: !
182: !TEST*/