GCC Code Coverage Report


Directory: ./
File: src/mfn/tutorials/ex23.c
Date: 2025-10-04 04:19:13
Exec Total Coverage
Lines: 70 70 100.0%
Functions: 2 2 100.0%
Branches: 199 302 65.9%

Line Branch Exec Source
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
11 static char help[] = "Computes exp(t*A)*v for a matrix associated with a Markov model.\n\n"
12 "The command line options are:\n"
13 " -t <t>, where <t> = time parameter (multiplies the matrix).\n"
14 " -m <m>, where <m> = number of grid subdivisions in each dimension.\n\n"
15 "To draw the solution run with -mfn_view_solution draw -draw_pause -1\n\n";
16
17 #include <slepcmfn.h>
18
19 /*
20 User-defined routines
21 */
22 PetscErrorCode MatMarkovModel(PetscInt m,Mat A);
23
24 10 int main(int argc,char **argv)
25 {
26 10 Mat A; /* problem matrix */
27 10 MFN mfn;
28 10 FN f;
29 10 PetscReal tol,norm;
30 10 PetscScalar t=2.0;
31 10 Vec v,y;
32 10 PetscInt N,m=15,ncv,maxit,its;
33 10 MFNConvergedReason reason;
34
35
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
10 PetscFunctionBeginUser;
36
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
37
38
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
39
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-t",&t,NULL));
40 10 N = m*(m+1)/2;
41
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov y=exp(t*A)*e_1, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",N,m));
42
43 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44 Compute the transition probability matrix, A
45 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
46
47
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
48
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
49
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatSetFromOptions(A));
50
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatMarkovModel(m,A));
51
52 /* set v = e_1 */
53
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatCreateVecs(A,NULL,&y));
54
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatCreateVecs(A,NULL,&v));
55
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(VecSetValue(v,0,1.0,INSERT_VALUES));
56
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(VecAssemblyBegin(v));
57
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(VecAssemblyEnd(v));
58
59 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60 Create the solver and set various options
61 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
62 /*
63 Create matrix function solver context
64 */
65
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
66
67 /*
68 Set operator matrix, the function to compute, and other options
69 */
70
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MFNSetOperator(mfn,A));
71
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MFNGetFN(mfn,&f));
72
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(FNSetType(f,FNEXP));
73
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(FNSetScale(f,t,1.0));
74
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MFNSetTolerances(mfn,1e-07,PETSC_CURRENT));
75
76 /*
77 Set solver parameters at runtime
78 */
79
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MFNSetFromOptions(mfn));
80
81 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82 Solve the problem, y=exp(t*A)*v
83 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84
85
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MFNSolve(mfn,v,y));
86
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MFNGetConvergedReason(mfn,&reason));
87
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
10 PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge");
88
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(VecNorm(y,NORM_2,&norm));
89
90 /*
91 Optional: Get some information from the solver and display it
92 */
93
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MFNGetIterationNumber(mfn,&its));
94
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",its));
95
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MFNGetDimensions(mfn,&ncv));
96
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Subspace dimension: %" PetscInt_FMT "\n",ncv));
97
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MFNGetTolerances(mfn,&tol,&maxit));
98
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit));
99
100 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101 Display solution and clean up
102 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computed vector at time t=%.4g has norm %g\n\n",(double)PetscRealPart(t),(double)norm));
104
105 /*
106 Free work space
107 */
108
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MFNDestroy(&mfn));
109
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatDestroy(&A));
110
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(VecDestroy(&v));
111
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(VecDestroy(&y));
112
3/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
10 PetscCall(SlepcFinalize());
113 return 0;
114 }
115
116 /*
117 Matrix generator for a Markov model of a random walk on a triangular grid.
118 See ex5.c for additional details.
119 */
120 10 PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
121 {
122 10 const PetscReal cst = 0.5/(PetscReal)(m-1);
123 10 PetscReal pd,pu;
124 10 PetscInt Istart,Iend,i,j,jmax,ix=0;
125
126
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
10 PetscFunctionBeginUser;
127
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
128
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
160 for (i=1;i<=m;i++) {
129 150 jmax = m-i+1;
130
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1350 for (j=1;j<=jmax;j++) {
131 1200 ix = ix + 1;
132
2/4
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
1200 if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */
133
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
1200 if (j!=jmax) {
134 1050 pd = cst*(PetscReal)(i+j-1);
135 /* north */
136
6/8
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
1050 if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES));
137
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
910 else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES));
138 /* east */
139
6/8
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
1050 if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES));
140
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
910 else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES));
141 }
142 /* south */
143 1200 pu = 0.5 - cst*(PetscReal)(i+j-3);
144
6/8
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
1200 if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES));
145 /* west */
146
6/8
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
1200 if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES));
147 }
148 }
149
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
150
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
10 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
151
5/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
2 PetscFunctionReturn(PETSC_SUCCESS);
152 }
153
154 /*TEST
155
156 test:
157 suffix: 1
158 args: -mfn_ncv 6
159
160 TEST*/
161