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 | 10 | int main(int argc,char **argv) | |
27 | { | ||
28 | 10 | Mat A; /* operator matrix */ | |
29 | 10 | EPS eps; /* eigenproblem solver context */ | |
30 | 10 | EPSType type; | |
31 | 10 | PetscScalar target=0.5; | |
32 | 10 | PetscInt N,m=15,nev; | |
33 | 10 | PetscBool terse; | |
34 | 10 | PetscViewer viewer; | |
35 | 10 | char str[50]; | |
36 | |||
37 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
10 | PetscFunctionBeginUser; |
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(SlepcInitialize(&argc,&argv,NULL,help)); |
39 | |||
40 |
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)); |
41 | 10 | N = m*(m+1)/2; | |
42 |
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 Model, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n",N,m)); |
43 |
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,"-target",&target,NULL)); |
44 |
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(SlepcSNPrintfScalar(str,sizeof(str),target,PETSC_FALSE)); |
45 |
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,"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 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)); |
52 |
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)); |
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(MatSetFromOptions(A)); |
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(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 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | 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 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(EPSSetOperators(eps,A,NULL)); |
69 |
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(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 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | 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 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(EPSSetFromOptions(eps)); |
81 | |||
82 | /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | ||
83 | Solve the eigensystem | ||
84 | - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ | ||
85 | |||
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(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 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(EPSGetType(eps,&type)); |
92 |
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," Solution method: %s\n\n",type)); |
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(EPSGetDimensions(eps,&nev,NULL,NULL)); |
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 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 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse)); |
102 |
5/8✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ 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.
|
10 | 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 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(EPSDestroy(&eps)); |
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(MatDestroy(&A)); |
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 | |||
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 | 10 | PetscErrorCode MatMarkovModel(PetscInt m,Mat A) | |
136 | { | ||
137 | 10 | const PetscReal cst = 0.5/(PetscReal)(m-1); | |
138 | 10 | PetscReal pd,pu; | |
139 | 10 | PetscInt Istart,Iend,i,j,jmax,ix=0; | |
140 | |||
141 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
10 | PetscFunctionBeginUser; |
142 |
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)); |
143 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
160 | for (i=1;i<=m;i++) { |
144 | 150 | jmax = m-i+1; | |
145 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
1350 | for (j=1;j<=jmax;j++) { |
146 | 1200 | ix = ix + 1; | |
147 |
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 */ |
148 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
1200 | if (j!=jmax) { |
149 | 1050 | pd = cst*(PetscReal)(i+j-1); | |
150 | /* north */ | ||
151 |
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)); |
152 |
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)); |
153 | /* east */ | ||
154 |
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)); |
155 |
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)); |
156 | } | ||
157 | /* south */ | ||
158 | 1200 | pu = 0.5 - cst*(PetscReal)(i+j-3); | |
159 |
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)); |
160 | /* west */ | ||
161 |
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)); |
162 | } | ||
163 | } | ||
164 |
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)); |
165 |
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)); |
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 | 68705 | PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx) | |
178 | { | ||
179 | 68705 | PetscScalar target = *(PetscScalar*)ctx; | |
180 | 68705 | PetscReal da,db; | |
181 | 68705 | PetscBool aisright,bisright; | |
182 | |||
183 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
68705 | PetscFunctionBeginUser; |
184 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
68705 | if (PetscRealPart(target) < PetscRealPart(ar)) aisright = PETSC_TRUE; |
185 | 38066 | else aisright = PETSC_FALSE; | |
186 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
68705 | if (PetscRealPart(target) < PetscRealPart(br)) bisright = PETSC_TRUE; |
187 | 38630 | else bisright = PETSC_FALSE; | |
188 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
68705 | if (aisright == bisright) { |
189 | /* both are on the same side of the target */ | ||
190 | 58671 | da = SlepcAbsEigenvalue(ar-target,ai); | |
191 | 58671 | db = SlepcAbsEigenvalue(br-target,bi); | |
192 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
58671 | if (da < db) *r = -1; |
193 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
51186 | else if (da > db) *r = 1; |
194 | ✗ | else *r = 0; | |
195 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
10034 | } else if (aisright && !bisright) *r = -1; /* 'a' is on the right */ |
196 | 4735 | 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.
|
68705 | 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 |