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[] = "Solves the same problem as in ex5, with a user-defined stopping test."
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 stopping test function.\n\n"
14 : "The command line options are:\n"
15 : " -m <m>, where <m> = number of grid subdivisions in each dimension.\n"
16 : " -seconds <s>, where <s> = maximum time in seconds allowed for computation.\n\n";
17 :
18 : #include <slepceps.h>
19 : #include <petsctime.h>
20 :
21 : /*
22 : User-defined routines
23 : */
24 :
25 : PetscErrorCode MyStoppingTest(EPS,PetscInt,PetscInt,PetscInt,PetscInt,EPSConvergedReason*,void*);
26 : PetscErrorCode MatMarkovModel(PetscInt,Mat);
27 :
28 1 : int main(int argc,char **argv)
29 : {
30 1 : Mat A; /* operator matrix */
31 1 : EPS eps; /* eigenproblem solver context */
32 1 : PetscReal seconds=2.5; /* maximum time allowed for computation */
33 1 : PetscLogDouble deadline; /* time to abort computation */
34 1 : PetscInt N,m=15,nconv;
35 1 : PetscBool terse;
36 1 : PetscViewer viewer;
37 1 : EPSConvergedReason reason;
38 :
39 1 : PetscFunctionBeginUser;
40 1 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
41 :
42 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
43 1 : N = m*(m+1)/2;
44 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n",N,m));
45 1 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-seconds",&seconds,NULL));
46 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Maximum time for computation is set to %g seconds.\n\n",(double)seconds));
47 1 : deadline = seconds;
48 1 : PetscCall(PetscTimeAdd(&deadline));
49 :
50 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51 : Compute the operator matrix that defines the eigensystem, Ax=kx
52 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53 :
54 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
55 1 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
56 1 : PetscCall(MatSetFromOptions(A));
57 1 : PetscCall(MatMarkovModel(m,A));
58 :
59 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60 : Create the eigensolver and set various options
61 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
62 :
63 1 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
64 1 : PetscCall(EPSSetOperators(eps,A,NULL));
65 1 : PetscCall(EPSSetProblemType(eps,EPS_NHEP));
66 1 : PetscCall(EPSSetStoppingTestFunction(eps,MyStoppingTest,&deadline,NULL));
67 1 : PetscCall(EPSSetFromOptions(eps));
68 :
69 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70 : Solve the eigensystem
71 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72 :
73 1 : PetscCall(EPSSolve(eps));
74 :
75 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76 : Display solution and clean up
77 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78 :
79 : /* show detailed info unless -terse option is given by user */
80 1 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
81 1 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
82 : else {
83 1 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
84 1 : PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
85 1 : PetscCall(EPSGetConvergedReason(eps,&reason));
86 1 : if (reason!=EPS_CONVERGED_USER) {
87 0 : PetscCall(EPSConvergedReasonView(eps,viewer));
88 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,viewer));
89 : } else {
90 1 : PetscCall(EPSGetConverged(eps,&nconv));
91 1 : PetscCall(PetscViewerASCIIPrintf(viewer,"Eigensolve finished with %" PetscInt_FMT " converged eigenpairs; reason=%s\n",nconv,EPSConvergedReasons[reason]));
92 : }
93 1 : PetscCall(PetscViewerPopFormat(viewer));
94 : }
95 1 : PetscCall(EPSDestroy(&eps));
96 1 : PetscCall(MatDestroy(&A));
97 1 : PetscCall(SlepcFinalize());
98 : return 0;
99 : }
100 :
101 : /*
102 : Matrix generator for a Markov model of a random walk on a triangular grid.
103 :
104 : This subroutine generates a test matrix that models a random walk on a
105 : triangular grid. This test example was used by G. W. Stewart ["{SRRIT} - a
106 : FORTRAN subroutine to calculate the dominant invariant subspaces of a real
107 : matrix", Tech. report. TR-514, University of Maryland (1978).] and in a few
108 : papers on eigenvalue problems by Y. Saad [see e.g. LAA, vol. 34, pp. 269-295
109 : (1980) ]. These matrices provide reasonably easy test problems for eigenvalue
110 : algorithms. The transpose of the matrix is stochastic and so it is known
111 : that one is an exact eigenvalue. One seeks the eigenvector of the transpose
112 : associated with the eigenvalue unity. The problem is to calculate the steady
113 : state probability distribution of the system, which is the eigevector
114 : associated with the eigenvalue one and scaled in such a way that the sum all
115 : the components is equal to one.
116 :
117 : Note: the code will actually compute the transpose of the stochastic matrix
118 : that contains the transition probabilities.
119 : */
120 1 : PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
121 : {
122 1 : const PetscReal cst = 0.5/(PetscReal)(m-1);
123 1 : PetscReal pd,pu;
124 1 : PetscInt Istart,Iend,i,j,jmax,ix=0;
125 :
126 1 : PetscFunctionBeginUser;
127 1 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
128 351 : for (i=1;i<=m;i++) {
129 350 : jmax = m-i+1;
130 61775 : for (j=1;j<=jmax;j++) {
131 61425 : ix = ix + 1;
132 61425 : if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */
133 61425 : if (j!=jmax) {
134 61075 : pd = cst*(PetscReal)(i+j-1);
135 : /* north */
136 61075 : if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES));
137 60726 : else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES));
138 : /* east */
139 61075 : if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES));
140 60726 : else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES));
141 : }
142 : /* south */
143 61425 : pu = 0.5 - cst*(PetscReal)(i+j-3);
144 61425 : if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES));
145 : /* west */
146 61425 : if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES));
147 : }
148 : }
149 1 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
150 1 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
151 1 : PetscFunctionReturn(PETSC_SUCCESS);
152 : }
153 :
154 : /*
155 : Function for user-defined stopping test.
156 :
157 : Checks that the computing time has not exceeded the deadline.
158 : */
159 24 : PetscErrorCode MyStoppingTest(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
160 : {
161 24 : PetscLogDouble now,deadline = *(PetscLogDouble*)ctx;
162 :
163 24 : PetscFunctionBeginUser;
164 : /* check if usual termination conditions are met */
165 24 : PetscCall(EPSStoppingBasic(eps,its,max_it,nconv,nev,reason,NULL));
166 24 : if (*reason==EPS_CONVERGED_ITERATING) {
167 : /* check if deadline has expired */
168 24 : PetscCall(PetscTime(&now));
169 24 : if (now>deadline) *reason = EPS_CONVERGED_USER;
170 : }
171 24 : PetscFunctionReturn(PETSC_SUCCESS);
172 : }
173 :
174 : /*TEST
175 :
176 : test:
177 : suffix: 1
178 : args: -m 350 -seconds 0.6
179 : requires: !single
180 :
181 : TEST*/
|
