Actual source code: ex1.c
slepc-3.7.4 2017-05-17
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2016, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: static char help[] = "Standard symmetric eigenproblem corresponding to the Laplacian operator in 1 dimension.\n\n"
23: "The command line options are:\n"
24: " -n <n>, where <n> = number of grid subdivisions = matrix dimension.\n\n";
26: #include <slepceps.h>
30: int main(int argc,char **argv)
31: {
32: Mat A; /* problem matrix */
33: EPS eps; /* eigenproblem solver context */
34: EPSType type;
35: PetscReal error,tol,re,im;
36: PetscScalar kr,ki;
37: Vec xr,xi;
38: PetscInt n=30,i,Istart,Iend,nev,maxit,its,nconv;
41: SlepcInitialize(&argc,&argv,(char*)0,help);
43: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
44: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem, n=%D\n\n",n);
46: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: Compute the operator matrix that defines the eigensystem, Ax=kx
48: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50: MatCreate(PETSC_COMM_WORLD,&A);
51: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
52: MatSetFromOptions(A);
53: MatSetUp(A);
55: MatGetOwnershipRange(A,&Istart,&Iend);
56: for (i=Istart;i<Iend;i++) {
57: if (i>0) { MatSetValue(A,i,i-1,-1.0,INSERT_VALUES); }
58: if (i<n-1) { MatSetValue(A,i,i+1,-1.0,INSERT_VALUES); }
59: MatSetValue(A,i,i,2.0,INSERT_VALUES);
60: }
61: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
62: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
64: MatCreateVecs(A,NULL,&xr);
65: MatCreateVecs(A,NULL,&xi);
67: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: Create the eigensolver and set various options
69: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: /*
71: Create eigensolver context
72: */
73: EPSCreate(PETSC_COMM_WORLD,&eps);
75: /*
76: Set operators. In this case, it is a standard eigenvalue problem
77: */
78: EPSSetOperators(eps,A,NULL);
79: EPSSetProblemType(eps,EPS_HEP);
81: /*
82: Set solver parameters at runtime
83: */
84: EPSSetFromOptions(eps);
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Solve the eigensystem
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: EPSSolve(eps);
91: /*
92: Optional: Get some information from the solver and display it
93: */
94: EPSGetIterationNumber(eps,&its);
95: PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
96: EPSGetType(eps,&type);
97: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
98: EPSGetDimensions(eps,&nev,NULL,NULL);
99: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
100: EPSGetTolerances(eps,&tol,&maxit);
101: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Display solution and clean up
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: /*
107: Get number of converged approximate eigenpairs
108: */
109: EPSGetConverged(eps,&nconv);
110: PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %D\n\n",nconv);
112: if (nconv>0) {
113: /*
114: Display eigenvalues and relative errors
115: */
116: PetscPrintf(PETSC_COMM_WORLD,
117: " k ||Ax-kx||/||kx||\n"
118: " ----------------- ------------------\n");
120: for (i=0;i<nconv;i++) {
121: /*
122: Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
123: ki (imaginary part)
124: */
125: EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);
126: /*
127: Compute the relative error associated to each eigenpair
128: */
129: EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);
131: #if defined(PETSC_USE_COMPLEX)
132: re = PetscRealPart(kr);
133: im = PetscImaginaryPart(kr);
134: #else
135: re = kr;
136: im = ki;
137: #endif
138: if (im!=0.0) {
139: PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12g\n",(double)re,(double)im,(double)error);
140: } else {
141: PetscPrintf(PETSC_COMM_WORLD," %12f %12g\n",(double)re,(double)error);
142: }
143: }
144: PetscPrintf(PETSC_COMM_WORLD,"\n");
145: }
147: /*
148: Free work space
149: */
150: EPSDestroy(&eps);
151: MatDestroy(&A);
152: VecDestroy(&xr);
153: VecDestroy(&xi);
154: SlepcFinalize();
155: return ierr;
156: }