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