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[] = "Illustrates region filtering. "
12 : "Based on ex5.\n"
13 : "The command line options are:\n"
14 : " -m <m>, where <m> = number of grid subdivisions in each dimension.\n\n";
15 :
16 : #include <slepceps.h>
17 :
18 : /*
19 : User-defined routines
20 : */
21 : PetscErrorCode MatMarkovModel(PetscInt m,Mat A);
22 :
23 1 : int main(int argc,char **argv)
24 : {
25 1 : Mat A;
26 1 : EPS eps;
27 1 : ST st;
28 1 : RG rg;
29 1 : PetscReal radius,tol=PETSC_SMALL;
30 1 : PetscScalar target=0.5,kr,ki;
31 1 : PetscInt N,m=15,nev,i,nconv;
32 1 : PetscBool checkfile;
33 1 : char filename[PETSC_MAX_PATH_LEN];
34 1 : PetscViewer viewer;
35 1 : char str[50];
36 :
37 1 : PetscFunctionBeginUser;
38 1 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
39 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
40 1 : N = m*(m+1)/2;
41 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n",N,m));
42 1 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-target",&target,NULL));
43 1 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),target,PETSC_FALSE));
44 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Searching closest eigenvalues to the right of %s.\n\n",str));
45 :
46 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47 : Compute the operator matrix that defines the eigensystem, Ax=kx
48 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
49 :
50 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
51 1 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
52 1 : PetscCall(MatSetFromOptions(A));
53 1 : PetscCall(MatMarkovModel(m,A));
54 :
55 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56 : Create a standalone spectral transformation
57 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58 :
59 1 : PetscCall(STCreate(PETSC_COMM_WORLD,&st));
60 1 : PetscCall(STSetType(st,STSINVERT));
61 1 : PetscCall(STSetShift(st,target));
62 :
63 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64 : Create a region for filtering
65 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66 :
67 1 : PetscCall(RGCreate(PETSC_COMM_WORLD,&rg));
68 1 : PetscCall(RGSetType(rg,RGELLIPSE));
69 1 : radius = (1.0-PetscRealPart(target))/2.0;
70 1 : PetscCall(RGEllipseSetParameters(rg,target+radius,radius,0.1));
71 :
72 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73 : Create the eigensolver and set various options
74 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75 :
76 1 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
77 1 : PetscCall(EPSSetST(eps,st));
78 1 : PetscCall(EPSSetRG(eps,rg));
79 1 : PetscCall(EPSSetOperators(eps,A,NULL));
80 1 : PetscCall(EPSSetProblemType(eps,EPS_NHEP));
81 1 : PetscCall(EPSSetTolerances(eps,tol,PETSC_DEFAULT));
82 1 : PetscCall(EPSSetTarget(eps,target));
83 1 : PetscCall(EPSSetFromOptions(eps));
84 :
85 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86 : Solve the eigensystem and display solution
87 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88 :
89 1 : PetscCall(EPSSolve(eps));
90 1 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
91 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
92 1 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
93 :
94 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95 : Check file containing the eigenvalues
96 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97 1 : PetscCall(PetscOptionsGetString(NULL,NULL,"-checkfile",filename,sizeof(filename),&checkfile));
98 1 : if (checkfile) {
99 : #if defined(PETSC_HAVE_COMPLEX)
100 1 : PetscComplex *eigs,eval;
101 :
102 1 : PetscCall(EPSGetConverged(eps,&nconv));
103 1 : PetscCall(PetscMalloc1(nconv,&eigs));
104 1 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
105 1 : PetscCall(PetscViewerBinaryRead(viewer,eigs,nconv,NULL,PETSC_COMPLEX));
106 1 : PetscCall(PetscViewerDestroy(&viewer));
107 5 : for (i=0;i<nconv;i++) {
108 4 : PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,NULL,NULL));
109 : #if defined(PETSC_USE_COMPLEX)
110 4 : eval = kr;
111 : #else
112 : eval = PetscCMPLX(kr,ki);
113 : #endif
114 4 : PetscCheck(eval==eigs[i],PETSC_COMM_WORLD,PETSC_ERR_FILE_UNEXPECTED,"Eigenvalues in the file do not match");
115 : }
116 1 : PetscCall(PetscFree(eigs));
117 : #else
118 : SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"The -checkfile option requires C99 complex numbers");
119 : #endif
120 : }
121 :
122 1 : PetscCall(EPSDestroy(&eps));
123 1 : PetscCall(STDestroy(&st));
124 1 : PetscCall(RGDestroy(&rg));
125 1 : PetscCall(MatDestroy(&A));
126 1 : PetscCall(SlepcFinalize());
127 : return 0;
128 : }
129 :
130 1 : PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
131 : {
132 1 : const PetscReal cst = 0.5/(PetscReal)(m-1);
133 1 : PetscReal pd,pu;
134 1 : PetscInt Istart,Iend,i,j,jmax,ix=0;
135 :
136 1 : PetscFunctionBeginUser;
137 1 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
138 16 : for (i=1;i<=m;i++) {
139 15 : jmax = m-i+1;
140 135 : for (j=1;j<=jmax;j++) {
141 120 : ix = ix + 1;
142 120 : if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */
143 120 : if (j!=jmax) {
144 105 : pd = cst*(PetscReal)(i+j-1);
145 : /* north */
146 105 : if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES));
147 91 : else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES));
148 : /* east */
149 105 : if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES));
150 91 : else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES));
151 : }
152 : /* south */
153 120 : pu = 0.5 - cst*(PetscReal)(i+j-3);
154 120 : if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES));
155 : /* west */
156 120 : if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES));
157 : }
158 : }
159 1 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
160 1 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
161 1 : PetscFunctionReturn(PETSC_SUCCESS);
162 : }
163 :
164 : /*TEST
165 :
166 : test:
167 : suffix: 1
168 : args: -eps_nev 4 -eps_ncv 20 -eps_view_values binary:myvalues.bin -checkfile myvalues.bin
169 : output_file: output/test11_1.out
170 : requires: !single c99_complex
171 :
172 : TEST*/
|
