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 feeding exact eigenvectors as initial vectors of a second solve.\n\n"
12 : "Based on ex2.\n"
13 : "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 5 : int main(int argc,char **argv)
20 : {
21 5 : Mat A;
22 5 : EPS eps;
23 5 : PetscInt N,n=10,m,Istart,Iend,II,nev,nconv,i,j,neigs=5;
24 5 : PetscBool flag,terse;
25 5 : Vec v,*X;
26 :
27 5 : PetscFunctionBeginUser;
28 5 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
29 5 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
30 5 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
31 5 : if (!flag) m=n;
32 5 : N = n*m;
33 5 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-neigs",&neigs,NULL));
34 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n2-D Laplacian Eigenproblem, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid), neigs=%" PetscInt_FMT "\n\n",N,n,m,neigs));
35 :
36 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
37 : Create the 2-D Laplacian
38 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
39 :
40 5 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
41 5 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
42 5 : PetscCall(MatSetFromOptions(A));
43 5 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
44 505 : for (II=Istart;II<Iend;II++) {
45 500 : i = II/n; j = II-i*n;
46 500 : if (i>0) PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES));
47 500 : if (i<m-1) PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES));
48 500 : if (j>0) PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES));
49 500 : if (j<n-1) PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES));
50 500 : PetscCall(MatSetValue(A,II,II,4.0,INSERT_VALUES));
51 : }
52 5 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
53 5 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
54 5 : PetscCall(MatCreateVecs(A,&v,NULL));
55 :
56 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57 : Create the eigensolver and set various options
58 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59 :
60 5 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
61 5 : PetscCall(EPSSetOperators(eps,A,NULL));
62 5 : PetscCall(EPSSetProblemType(eps,EPS_HEP));
63 5 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
64 5 : PetscCall(EPSSetDimensions(eps,neigs,PETSC_DEFAULT,PETSC_DEFAULT));
65 5 : PetscCall(EPSSetFromOptions(eps));
66 :
67 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68 : Solve the eigensystem for the first time
69 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70 :
71 5 : PetscCall(EPSSolve(eps));
72 5 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
73 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
74 :
75 5 : PetscCall(EPSGetConverged(eps,&nconv));
76 5 : PetscCheck(nconv>=neigs,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Only %" PetscInt_FMT " eigenvalues have converged, %" PetscInt_FMT " requested",nconv,neigs);
77 5 : PetscCall(VecDuplicateVecs(v,neigs,&X));
78 30 : for (i=0;i<neigs;i++) PetscCall(EPSGetEigenvector(eps,i,X[i],NULL));
79 :
80 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81 : Display solution of first solve
82 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
83 :
84 5 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
85 5 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
86 : else {
87 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
88 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
89 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
90 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
91 : }
92 :
93 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94 : Solve the eigensystem again, feeding the initial vectors
95 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96 :
97 5 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solving again with eigenvectors as initial guesses\n"));
98 5 : PetscCall(EPSSetInitialSpace(eps,neigs,X));
99 5 : PetscCall(EPSSolve(eps));
100 :
101 5 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
102 : else {
103 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
104 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
105 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
106 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
107 : }
108 :
109 5 : PetscCall(VecDestroy(&v));
110 5 : PetscCall(VecDestroyVecs(neigs,&X));
111 5 : PetscCall(EPSDestroy(&eps));
112 5 : PetscCall(MatDestroy(&A));
113 5 : PetscCall(SlepcFinalize());
114 : return 0;
115 : }
116 :
117 : /*TEST
118 :
119 : test:
120 : suffix: 1
121 : args: -eps_type {{gd jd rqcg lobpcg}} -terse
122 :
123 : testset:
124 : args: -eps_interval .17,1.3 -terse
125 : filter: grep -v "requested"
126 : output_file: output/test27_2.out
127 : test:
128 : suffix: 2
129 : args: -st_type filter -st_filter_degree 150 -eps_nev 1
130 : requires: !single
131 : test:
132 : suffix: 2_evsl
133 : nsize: {{1 2}}
134 : args: -eps_type evsl
135 : requires: evsl
136 :
137 : TEST*/
|
