GCC Code Coverage Report

Directory: ./
File: src/eps/tutorials/ex18.c
Date: 2026-02-22 03:58:10
Exec Total Coverage
Lines: 76 82 92.7%
Functions: 3 3 100.0%
Branches: 175 296 59.1%
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[] = "Solves the same problem as in ex5, but with a user-defined sorting criterion.\n\n"
12 "It is a standard nonsymmetric eigenproblem with real eigenvalues and the rightmost eigenvalue is known to be 1.\n"
13 "This example illustrates how the user can set a custom spectrum selection.\n\n"
14 "The command line options are:\n"
15 " -m <m>, where <m> = number of grid subdivisions in each dimension.\n\n";
16
17 #include <slepceps.h>
18
19 /*
20 User-defined routines
21 */
22
23 PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx);
24 PetscErrorCode MatMarkovModel(PetscInt m,Mat A);
25
26 8 int main(int argc,char **argv)
27 {
28 8 Mat A; /* operator matrix */
29 8 EPS eps; /* eigenproblem solver context */
30 8 EPSType type;
31 8 PetscScalar target=0.5;
32 8 PetscInt N,m=15,nev;
33 8 PetscBool terse;
34 8 PetscViewer viewer;
35 8 char str[50];
36
37
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 8 PetscFunctionBeginUser;
38
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
39
40
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
41 8 N = m*(m+1)/2;
42
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n",N,m));
43
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-target",&target,NULL));
44
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(SlepcSNPrintfScalar(str,sizeof(str),target,PETSC_FALSE));
45
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Searching closest eigenvalues to the right of %s.\n\n",str));
46
47 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48 Compute the operator matrix that defines the eigensystem, Ax=kx
49 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50
51
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
52
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
53
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(MatSetFromOptions(A));
54
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(MatMarkovModel(m,A));
55
56 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57 Create the eigensolver and set various options
58 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59
60 /*
61 Create eigensolver context
62 */
63
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
64
65 /*
66 Set operators. In this case, it is a standard eigenvalue problem
67 */
68
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(EPSSetOperators(eps,A,NULL));
69
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(EPSSetProblemType(eps,EPS_NHEP));
70
71 /*
72 Set the custom comparing routine in order to obtain the eigenvalues
73 closest to the target on the right only
74 */
75
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(EPSSetEigenvalueComparison(eps,MyEigenSort,&target));
76
77 /*
78 Set solver parameters at runtime
79 */
80
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(EPSSetFromOptions(eps));
81
82 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83 Solve the eigensystem
84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85
86
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(EPSSolve(eps));
87
88 /*
89 Optional: Get some information from the solver and display it
90 */
91
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(EPSGetType(eps,&type));
92
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
93
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
94
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
95
96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97 Display solution and clean up
98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99
100 /* show detailed info unless -terse option is given by user */
101
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
102
5/8
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 8 if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
103 else {
104 ✗ PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
105 ✗ PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
106 ✗ PetscCall(EPSConvergedReasonView(eps,viewer));
107 ✗ PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,viewer));
108 ✗ PetscCall(PetscViewerPopFormat(viewer));
109 }
110
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(EPSDestroy(&eps));
111
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(MatDestroy(&A));
112
2/6
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 8 PetscCall(SlepcFinalize());
113 return 0;
114 }
115
116 /*
117 Matrix generator for a Markov model of a random walk on a triangular grid.
118
119 This subroutine generates a test matrix that models a random walk on a
120 triangular grid. This test example was used by G. W. Stewart ["{SRRIT} - a
121 FORTRAN subroutine to calculate the dominant invariant subspaces of a real
122 matrix", Tech. report. TR-514, University of Maryland (1978).] and in a few
123 papers on eigenvalue problems by Y. Saad [see e.g. LAA, vol. 34, pp. 269-295
124 (1980) ]. These matrices provide reasonably easy test problems for eigenvalue
125 algorithms. The transpose of the matrix is stochastic and so it is known
126 that one is an exact eigenvalue. One seeks the eigenvector of the transpose
127 associated with the eigenvalue unity. The problem is to calculate the steady
128 state probability distribution of the system, which is the eigevector
129 associated with the eigenvalue one and scaled in such a way that the sum all
130 the components is equal to one.
131
132 Note: the code will actually compute the transpose of the stochastic matrix
133 that contains the transition probabilities.
134 */
135 8 PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
136 {
137 8 const PetscReal cst = 0.5/(PetscReal)(m-1);
138 8 PetscReal pd,pu;
139 8 PetscInt Istart,Iend,i,j,jmax,ix=0;
140
141
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 8 PetscFunctionBeginUser;
142
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
143
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
 128 for (i=1;i<=m;i++) {
144 120 jmax = m-i+1;
145
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
 1080 for (j=1;j<=jmax;j++) {
146 960 ix = ix + 1;
147
2/4
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
 960 if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */
148
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
 960 if (j!=jmax) {
149 840 pd = cst*(PetscReal)(i+j-1);
150 /* north */
151
6/8
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 840 if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES));
152
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 728 else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES));
153 /* east */
154
6/8
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 840 if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES));
155
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 728 else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES));
156 }
157 /* south */
158 960 pu = 0.5 - cst*(PetscReal)(i+j-3);
159
6/8
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 960 if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES));
160 /* west */
161
6/8
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 960 if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES));
162 }
163 }
164
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
165
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 8 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
166
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);
167 }
168
169 /*
170 Function for user-defined eigenvalue ordering criterion.
171
172 Given two eigenvalues ar+i*ai and br+i*bi, the subroutine must choose
173 one of them as the preferred one according to the criterion.
174 In this example, the preferred value is the one closest to the target,
175 but on the right side.
176 */
177 54964 PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx)
178 {
179 54964 PetscScalar target = *(PetscScalar*)ctx;
180 54964 PetscReal da,db;
181 54964 PetscBool aisright,bisright;
182
183
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 54964 PetscFunctionBeginUser;
184
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
 54964 if (PetscRealPart(target) < PetscRealPart(ar)) aisright = PETSC_TRUE;
185 30448 else aisright = PETSC_FALSE;
186
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
 54964 if (PetscRealPart(target) < PetscRealPart(br)) bisright = PETSC_TRUE;
187 30904 else bisright = PETSC_FALSE;
188
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
 54964 if (aisright == bisright) {
189 /* both are on the same side of the target */
190 46932 da = SlepcAbsEigenvalue(ar-target,ai);
191 46932 db = SlepcAbsEigenvalue(br-target,bi);
192
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
 46932 if (da < db) *r = -1;
193
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
 40944 else if (da > db) *r = 1;
194 ✗ else *r = 0;
195
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
 8032 } else if (aisright && !bisright) *r = -1; /* 'a' is on the right */
196 3788 else *r = 1; /* 'b' is on the right */
197
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.
 54964 PetscFunctionReturn(PETSC_SUCCESS);
198 }
199
200 /*TEST
201
202 test:
203 suffix: 1
204 args: -eps_nev 4 -terse
205 requires: !single
206 filter: sed -e "s/[+-]0\.0*i//g"
207
208 TEST*/
209