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 | 50 | int main(int argc,char **argv) | |
27 | { | ||
28 | 50 | Vec v0; /* initial vector */ | |
29 | 50 | Mat A; /* operator matrix */ | |
30 | 50 | EPS eps; /* eigenproblem solver context */ | |
31 | 50 | ST st; /* spectral transformation associated */ | |
32 | 50 | PetscReal tol=PETSC_SMALL; | |
33 | 50 | PetscScalar target=0.5; | |
34 | 50 | PetscInt N,m=15,nev; | |
35 | 50 | char str[50]; | |
36 | |||
37 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
50 | 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.
|
50 | 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.
|
50 | PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); |
41 | 50 | 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.
|
50 | 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.
|
50 | 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.
|
50 | 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.
|
50 | 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.
|
50 | 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.
|
50 | 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.
|
50 | 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.
|
50 | 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.
|
50 | 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.
|
50 | 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.
|
50 | PetscCall(EPSSetProblemType(eps,EPS_NHEP)); |
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.
|
50 | PetscCall(EPSSetTolerances(eps,tol,PETSC_CURRENT)); |
71 | |||
72 | /* | ||
73 | Set the custom comparing routine in order to obtain the eigenvalues | ||
74 | closest to the target on the right only | ||
75 | */ | ||
76 |
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.
|
50 | PetscCall(EPSSetEigenvalueComparison(eps,MyEigenSort,&target)); |
77 | |||
78 | /* | ||
79 | Set solver parameters at runtime | ||
80 | */ | ||
81 |
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.
|
50 | PetscCall(EPSSetFromOptions(eps)); |
82 | |||
83 | /* | ||
84 | Set the preconditioner based on A - target * I | ||
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.
|
50 | PetscCall(EPSGetST(eps,&st)); |
87 |
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.
|
50 | PetscCall(STSetShift(st,target)); |
88 | |||
89 | /* | ||
90 | Set the initial vector. This is optional, if not done the initial | ||
91 | vector is set to random values | ||
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.
|
50 | PetscCall(MatCreateVecs(A,&v0,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.
|
50 | PetscCall(VecSet(v0,1.0)); |
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.
|
50 | PetscCall(EPSSetInitialSpace(eps,1,&v0)); |
96 | |||
97 | /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | ||
98 | Solve the eigensystem | ||
99 | - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ | ||
100 | |||
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.
|
50 | PetscCall(EPSSolve(eps)); |
102 |
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.
|
50 | PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL)); |
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.
|
50 | PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev)); |
104 | |||
105 | /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | ||
106 | Display solution and clean up | ||
107 | - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ | ||
108 | |||
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.
|
50 | PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL)); |
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.
|
50 | 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.
|
50 | PetscCall(MatDestroy(&A)); |
112 |
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.
|
50 | PetscCall(VecDestroy(&v0)); |
113 |
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.
|
50 | PetscCall(SlepcFinalize()); |
114 | return 0; | ||
115 | } | ||
116 | |||
117 | 50 | PetscErrorCode MatMarkovModel(PetscInt m,Mat A) | |
118 | { | ||
119 | 50 | const PetscReal cst = 0.5/(PetscReal)(m-1); | |
120 | 50 | PetscReal pd,pu; | |
121 | 50 | PetscInt Istart,Iend,i,j,jmax,ix=0; | |
122 | |||
123 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
50 | PetscFunctionBeginUser; |
124 |
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.
|
50 | PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); |
125 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
800 | for (i=1;i<=m;i++) { |
126 | 750 | jmax = m-i+1; | |
127 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
6750 | for (j=1;j<=jmax;j++) { |
128 | 6000 | ix = ix + 1; | |
129 |
2/4✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
|
6000 | if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */ |
130 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
6000 | if (j!=jmax) { |
131 | 5250 | pd = cst*(PetscReal)(i+j-1); | |
132 | /* north */ | ||
133 |
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.
|
5250 | if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES)); |
134 |
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.
|
4550 | else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES)); |
135 | /* east */ | ||
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.
|
5250 | if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,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.
|
4550 | else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES)); |
138 | } | ||
139 | /* south */ | ||
140 | 6000 | pu = 0.5 - cst*(PetscReal)(i+j-3); | |
141 |
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.
|
6000 | if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES)); |
142 | /* west */ | ||
143 |
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.
|
6000 | if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES)); |
144 | } | ||
145 | } | ||
146 |
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.
|
50 | PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); |
147 |
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.
|
50 | PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); |
148 |
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.
|
10 | PetscFunctionReturn(PETSC_SUCCESS); |
149 | } | ||
150 | |||
151 | /* | ||
152 | Function for user-defined eigenvalue ordering criterion. | ||
153 | |||
154 | Given two eigenvalues ar+i*ai and br+i*bi, the subroutine must choose | ||
155 | one of them as the preferred one according to the criterion. | ||
156 | In this example, the preferred value is the one closest to the target, | ||
157 | but on the right side. | ||
158 | */ | ||
159 | 204175 | PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx) | |
160 | { | ||
161 | 204175 | PetscScalar target = *(PetscScalar*)ctx; | |
162 | 204175 | PetscReal da,db; | |
163 | 204175 | PetscBool aisright,bisright; | |
164 | |||
165 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
204175 | PetscFunctionBeginUser; |
166 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
204175 | if (PetscRealPart(target) < PetscRealPart(ar)) aisright = PETSC_TRUE; |
167 | 121736 | else aisright = PETSC_FALSE; | |
168 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
204175 | if (PetscRealPart(target) < PetscRealPart(br)) bisright = PETSC_TRUE; |
169 | 176033 | else bisright = PETSC_FALSE; | |
170 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
204175 | if (aisright == bisright) { |
171 | /* both are on the same side of the target */ | ||
172 | 140998 | da = SlepcAbsEigenvalue(ar-target,ai); | |
173 | 140998 | db = SlepcAbsEigenvalue(br-target,bi); | |
174 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
140998 | if (da < db) *r = -1; |
175 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 5 times.
|
38629 | else if (da > db) *r = 1; |
176 | 23 | else *r = 0; | |
177 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
63177 | } else if (aisright && !bisright) *r = -1; /* 'a' is on the right */ |
178 | 4440 | else *r = 1; /* 'b' is on the right */ | |
179 |
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.
|
204175 | PetscFunctionReturn(PETSC_SUCCESS); |
180 | } | ||
181 | |||
182 | /*TEST | ||
183 | |||
184 | testset: | ||
185 | args: -eps_nev 4 | ||
186 | requires: !single | ||
187 | output_file: output/test11_1.out | ||
188 | test: | ||
189 | suffix: 1 | ||
190 | args: -eps_type {{krylovschur arnoldi lapack}} -st_type sinvert | ||
191 | test: | ||
192 | suffix: 1_ks_cayley | ||
193 | args: -st_type cayley -st_cayley_antishift 1 | ||
194 | |||
195 | test: | ||
196 | suffix: 2 | ||
197 | args: -target 0.77 -eps_type gd -eps_nev 4 -eps_tol 1e-7 -eps_gd_krylov_start -eps_gd_blocksize 3 | ||
198 | requires: double | ||
199 | filter: sed -e "s/[+-]0\.0*i//g" | ||
200 | |||
201 | TEST*/ | ||
202 |