Actual source code: ex16.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7:
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[] = "Quadratic eigenproblem for testing the QEP object.\n\n"
23: "The command line options are:\n"
24: " -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
25: " -m <m>, where <m> = number of grid subdivisions in y dimension.\n\n";
27: #include slepcqep.h
31: int main( int argc, char **argv )
32: {
33: Mat M, C, K; /* problem matrices */
34: QEP qep; /* quadratic eigenproblem solver context */
35: const QEPType type;
36: PetscReal error, tol, re, im;
37: PetscScalar kr, ki;
39: PetscInt N, n=10, m, Istart, Iend, II, nev, maxit, i, j, its, nconv;
40: PetscTruth flag;
42: SlepcInitialize(&argc,&argv,(char*)0,help);
44: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
45: PetscOptionsGetInt(PETSC_NULL,"-m",&m,&flag);
46: if(!flag) m=n;
47: N = n*m;
48: PetscPrintf(PETSC_COMM_WORLD,"\nQuadratic Eigenproblem, N=%d (%dx%d grid)\n\n",N,n,m);
50: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: Compute the matrices that define the eigensystem, (k^2*K+k*X+M)x=0
52: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: /* K is the 2-D Laplacian */
55: MatCreate(PETSC_COMM_WORLD,&K);
56: MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,N,N);
57: MatSetFromOptions(K);
58:
59: MatGetOwnershipRange(K,&Istart,&Iend);
60: for( II=Istart; II<Iend; II++ ) {
61: i = II/n; j = II-i*n;
62: if(i>0) { MatSetValue(K,II,II-n,-1.0,INSERT_VALUES); }
63: if(i<m-1) { MatSetValue(K,II,II+n,-1.0,INSERT_VALUES); }
64: if(j>0) { MatSetValue(K,II,II-1,-1.0,INSERT_VALUES); }
65: if(j<n-1) { MatSetValue(K,II,II+1,-1.0,INSERT_VALUES); }
66: MatSetValue(K,II,II,4.0,INSERT_VALUES);
67: }
69: MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
70: MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);
72: /* C is the zero matrix */
73: MatCreate(PETSC_COMM_WORLD,&C);
74: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N);
75: MatSetFromOptions(C);
76: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
77: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
78:
79: /* M is the identity matrix */
80: MatCreate(PETSC_COMM_WORLD,&M);
81: MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,N,N);
82: MatSetFromOptions(M);
83: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
84: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
85: MatShift(M,1.0);
86:
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Create the eigensolver and set various options
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91: /*
92: Create eigensolver context
93: */
94: QEPCreate(PETSC_COMM_WORLD,&qep);
96: /*
97: Set matrices and problem type
98: */
99: QEPSetOperators(qep,M,C,K);
100: QEPSetProblemType(qep,QEP_GENERAL);
102: /*
103: Set solver parameters at runtime
104: */
105: QEPSetFromOptions(qep);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Solve the eigensystem
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: QEPSolve(qep);
113: /*
114: Optional: Get some information from the solver and display it
115: */
116: QEPGetIterationNumber(qep, &its);
117: PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);
118: QEPGetType(qep,&type);
119: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
120: QEPGetDimensions(qep,&nev,PETSC_NULL,PETSC_NULL);
121: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %d\n",nev);
122: QEPGetTolerances(qep,&tol,&maxit);
123: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Display solution and clean up
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129: /*
130: Get number of converged approximate eigenpairs
131: */
132: QEPGetConverged(qep,&nconv);
133: PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %d\n\n",nconv);
134:
136: if (nconv>0) {
137: /*
138: Display eigenvalues and relative errors
139: */
140: PetscPrintf(PETSC_COMM_WORLD,
141: " k ||(k^2M+Ck+K)x||/||kx||\n"
142: " ----------------- -------------------------\n" );
144: for( i=0; i<nconv; i++ ) {
145: /*
146: Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
147: ki (imaginary part)
148: */
149: QEPGetEigenpair(qep,i,&kr,&ki,PETSC_NULL,PETSC_NULL);
150: /*
151: Compute the relative error associated to each eigenpair
152: */
153: QEPComputeRelativeError(qep,i,&error);
155: #ifdef PETSC_USE_COMPLEX
156: re = PetscRealPart(kr);
157: im = PetscImaginaryPart(kr);
158: #else
159: re = kr;
160: im = ki;
161: #endif
162: if (im!=0.0) {
163: PetscPrintf(PETSC_COMM_WORLD," %9f%+9f j %12g\n",re,im,error);
164: } else {
165: PetscPrintf(PETSC_COMM_WORLD," %12f %12g\n",re,error);
166: }
167: }
168: PetscPrintf(PETSC_COMM_WORLD,"\n" );
169: }
170:
171: /*
172: Free work space
173: */
174: QEPDestroy(qep);
175: MatDestroy(M);
176: MatDestroy(C);
177: MatDestroy(K);
178: SlepcFinalize();
179: return 0;
180: }