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[] = "Computes the smallest nonzero eigenvalue of the Laplacian of a graph.\n\n"
12 : "This example illustrates EPSSetDeflationSpace(). The example graph corresponds to a "
13 : "2-D regular mesh. The command line options are:\n"
14 : " -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
15 : " -m <m>, where <m> = number of grid subdivisions in y dimension.\n\n";
16 :
17 : #include <slepceps.h>
18 :
19 2 : int main (int argc,char **argv)
20 : {
21 2 : EPS eps; /* eigenproblem solver context */
22 2 : Mat A; /* operator matrix */
23 2 : Vec x;
24 2 : EPSType type;
25 2 : PetscInt N,n=10,m,i,j,II,Istart,Iend,nev;
26 2 : PetscScalar w;
27 2 : PetscBool flag,terse;
28 :
29 2 : PetscFunctionBeginUser;
30 2 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
31 :
32 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
33 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
34 2 : if (!flag) m=n;
35 2 : N = n*m;
36 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nFiedler vector of a 2-D regular mesh, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
37 :
38 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39 : Compute the operator matrix that defines the eigensystem, Ax=kx
40 : In this example, A = L(G), where L is the Laplacian of graph G, i.e.
41 : Lii = degree of node i, Lij = -1 if edge (i,j) exists in G
42 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
43 :
44 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
45 2 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
46 2 : PetscCall(MatSetFromOptions(A));
47 :
48 2 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
49 202 : for (II=Istart;II<Iend;II++) {
50 200 : i = II/n; j = II-i*n;
51 200 : w = 0.0;
52 200 : if (i>0) { PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES)); w=w+1.0; }
53 200 : if (i<m-1) { PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES)); w=w+1.0; }
54 200 : if (j>0) { PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES)); w=w+1.0; }
55 200 : if (j<n-1) { PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES)); w=w+1.0; }
56 200 : PetscCall(MatSetValue(A,II,II,w,INSERT_VALUES));
57 : }
58 :
59 2 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
60 2 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
61 :
62 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63 : Create the eigensolver and set various options
64 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65 :
66 : /*
67 : Create eigensolver context
68 : */
69 2 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
70 :
71 : /*
72 : Set operators. In this case, it is a standard eigenvalue problem
73 : */
74 2 : PetscCall(EPSSetOperators(eps,A,NULL));
75 2 : PetscCall(EPSSetProblemType(eps,EPS_HEP));
76 :
77 : /*
78 : Select portion of spectrum
79 : */
80 2 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
81 :
82 : /*
83 : Set solver parameters at runtime
84 : */
85 2 : PetscCall(EPSSetFromOptions(eps));
86 :
87 : /*
88 : Attach deflation space: in this case, the matrix has a constant
89 : nullspace, [1 1 ... 1]^T is the eigenvector of the zero eigenvalue
90 : */
91 2 : PetscCall(MatCreateVecs(A,&x,NULL));
92 2 : PetscCall(VecSet(x,1.0));
93 2 : PetscCall(EPSSetDeflationSpace(eps,1,&x));
94 2 : PetscCall(VecDestroy(&x));
95 :
96 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97 : Solve the eigensystem
98 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99 :
100 2 : PetscCall(EPSSolve(eps));
101 :
102 : /*
103 : Optional: Get some information from the solver and display it
104 : */
105 2 : PetscCall(EPSGetType(eps,&type));
106 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
107 2 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
108 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
109 :
110 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111 : Display solution and clean up
112 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113 :
114 : /* show detailed info unless -terse option is given by user */
115 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
116 2 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
117 : else {
118 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
119 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
120 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
121 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
122 : }
123 2 : PetscCall(EPSDestroy(&eps));
124 2 : PetscCall(MatDestroy(&A));
125 2 : PetscCall(SlepcFinalize());
126 : return 0;
127 : }
128 :
129 : /*TEST
130 :
131 : testset:
132 : args: -eps_nev 4 -terse
133 : output_file: output/ex11_1.out
134 : test:
135 : suffix: 1
136 : args: -eps_krylovschur_restart .2
137 : test:
138 : suffix: 2
139 : args: -eps_ncv 20 -eps_target 0 -st_type sinvert -st_ksp_type cg -st_pc_type jacobi
140 :
141 : TEST*/
|
