Line data Source code
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[] = "Eigenvalue problem associated with a Markov model of a random walk on a triangular grid. "
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 the initial vector.\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 : PetscErrorCode MatMarkovModel(PetscInt m,Mat A);
23 :
24 9 : int main(int argc,char **argv)
25 : {
26 9 : Vec v0; /* initial vector */
27 9 : Mat A; /* operator matrix */
28 9 : EPS eps; /* eigenproblem solver context */
29 9 : EPSType type;
30 9 : EPSStop stop;
31 9 : PetscReal thres;
32 9 : PetscInt N,m=15,nev;
33 9 : PetscMPIInt rank;
34 9 : PetscBool terse;
35 :
36 9 : PetscFunctionBeginUser;
37 9 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
38 :
39 9 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
40 9 : N = m*(m+1)/2;
41 9 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",N,m));
42 :
43 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44 : Compute the operator matrix that defines the eigensystem, Ax=kx
45 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
46 :
47 9 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
48 9 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
49 9 : PetscCall(MatSetFromOptions(A));
50 9 : PetscCall(MatMarkovModel(m,A));
51 :
52 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53 : Create the eigensolver and set various options
54 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
55 :
56 : /*
57 : Create eigensolver context
58 : */
59 9 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
60 :
61 : /*
62 : Set operators. In this case, it is a standard eigenvalue problem
63 : */
64 9 : PetscCall(EPSSetOperators(eps,A,NULL));
65 9 : PetscCall(EPSSetProblemType(eps,EPS_NHEP));
66 :
67 : /*
68 : Set solver parameters at runtime
69 : */
70 9 : PetscCall(EPSSetFromOptions(eps));
71 :
72 : /*
73 : Set the initial vector. This is optional, if not done the initial
74 : vector is set to random values
75 : */
76 9 : PetscCall(MatCreateVecs(A,&v0,NULL));
77 9 : PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
78 9 : if (!rank) {
79 5 : PetscCall(VecSetValue(v0,0,1.0,INSERT_VALUES));
80 5 : PetscCall(VecSetValue(v0,1,1.0,INSERT_VALUES));
81 5 : PetscCall(VecSetValue(v0,2,1.0,INSERT_VALUES));
82 : }
83 9 : PetscCall(VecAssemblyBegin(v0));
84 9 : PetscCall(VecAssemblyEnd(v0));
85 9 : PetscCall(EPSSetInitialSpace(eps,1,&v0));
86 :
87 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88 : Solve the eigensystem
89 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90 :
91 9 : PetscCall(EPSSolve(eps));
92 :
93 : /*
94 : Optional: Get some information from the solver and display it
95 : */
96 9 : PetscCall(EPSGetType(eps,&type));
97 9 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
98 9 : PetscCall(EPSGetStoppingTest(eps,&stop));
99 9 : if (stop!=EPS_STOP_THRESHOLD) {
100 8 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
101 8 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
102 : } else {
103 1 : PetscCall(EPSGetThreshold(eps,&thres,NULL));
104 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using threshold: %.4g\n",(double)thres));
105 : }
106 :
107 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108 : Display solution and clean up
109 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110 :
111 : /* show detailed info unless -terse option is given by user */
112 9 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
113 9 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
114 : else {
115 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
116 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
117 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
118 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
119 : }
120 9 : PetscCall(EPSDestroy(&eps));
121 9 : PetscCall(MatDestroy(&A));
122 9 : PetscCall(VecDestroy(&v0));
123 9 : PetscCall(SlepcFinalize());
124 : return 0;
125 : }
126 :
127 : /*
128 : Matrix generator for a Markov model of a random walk on a triangular grid.
129 :
130 : This subroutine generates a test matrix that models a random walk on a
131 : triangular grid. This test example was used by G. W. Stewart ["{SRRIT} - a
132 : FORTRAN subroutine to calculate the dominant invariant subspaces of a real
133 : matrix", Tech. report. TR-514, University of Maryland (1978).] and in a few
134 : papers on eigenvalue problems by Y. Saad [see e.g. LAA, vol. 34, pp. 269-295
135 : (1980) ]. These matrices provide reasonably easy test problems for eigenvalue
136 : algorithms. The transpose of the matrix is stochastic and so it is known
137 : that one is an exact eigenvalue. One seeks the eigenvector of the transpose
138 : associated with the eigenvalue unity. The problem is to calculate the steady
139 : state probability distribution of the system, which is the eigevector
140 : associated with the eigenvalue one and scaled in such a way that the sum all
141 : the components is equal to one.
142 :
143 : Note: the code will actually compute the transpose of the stochastic matrix
144 : that contains the transition probabilities.
145 : */
146 9 : PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
147 : {
148 9 : const PetscReal cst = 0.5/(PetscReal)(m-1);
149 9 : PetscReal pd,pu;
150 9 : PetscInt Istart,Iend,i,j,jmax,ix=0;
151 :
152 9 : PetscFunctionBeginUser;
153 9 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
154 144 : for (i=1;i<=m;i++) {
155 135 : jmax = m-i+1;
156 1215 : for (j=1;j<=jmax;j++) {
157 1080 : ix = ix + 1;
158 1080 : if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */
159 600 : if (j!=jmax) {
160 525 : pd = cst*(PetscReal)(i+j-1);
161 : /* north */
162 525 : if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES));
163 455 : else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES));
164 : /* east */
165 525 : if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES));
166 455 : else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES));
167 : }
168 : /* south */
169 600 : pu = 0.5 - cst*(PetscReal)(i+j-3);
170 600 : if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES));
171 : /* west */
172 1080 : if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES));
173 : }
174 : }
175 9 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
176 9 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
177 9 : PetscFunctionReturn(PETSC_SUCCESS);
178 : }
179 :
180 : /*TEST
181 :
182 : test:
183 : suffix: 1
184 : nsize: 2
185 : args: -eps_largest_real -eps_nev 4 -eps_two_sided {{0 1}} -eps_krylovschur_locking {{0 1}} -ds_parallel synchronized -terse
186 : filter: sed -e "s/90424/90423/" | sed -e "s/85715/85714/"
187 :
188 : test:
189 : suffix: 2
190 : args: -eps_threshold_relative .9 -eps_ncv 10 -terse
191 : filter: sed -e "s/-0.85714/0.85714/" | sed -e "s/90424/90423/" | sed -e "s/-1.00000, 1.00000/1.00000, -1.00000/" | sed -e "s/-0.97137, 0.97137/0.97137, -0.97137/" | sed -e "s/-0.90423, 0.90423/0.90423, -0.90423/"
192 :
193 : TEST*/
|
