Actual source code: test6.c
slepc-3.18.1 2022-11-02
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: Example based on spring problem in NLEVP collection [1]. See the parameters
12: meaning at Example 2 in [2].
14: [1] T. Betcke, N. J. Higham, V. Mehrmann, C. Schroder, and F. Tisseur,
15: NLEVP: A Collection of Nonlinear Eigenvalue Problems, MIMS EPrint
16: 2010.98, November 2010.
17: [2] F. Tisseur, Backward error and condition of polynomial eigenvalue
18: problems, Linear Algebra and its Applications, 309 (2000), pp. 339--361,
19: April 2000.
20: */
22: static char help[] = "Tests multiple calls to PEPSolve with different matrix of different size.\n\n"
23: "This is based on the spring problem from NLEVP collection.\n\n"
24: "The command line options are:\n"
25: " -n <n> ... number of grid subdivisions.\n"
26: " -mu <value> ... mass (default 1).\n"
27: " -tau <value> ... damping constant of the dampers (default 10).\n"
28: " -kappa <value> ... damping constant of the springs (default 5).\n"
29: " -initv ... set an initial vector.\n\n";
31: #include <slepcpep.h>
33: int main(int argc,char **argv)
34: {
35: Mat M,C,K,A[3]; /* problem matrices */
36: PEP pep; /* polynomial eigenproblem solver context */
37: PetscInt n=30,Istart,Iend,i,nev;
38: PetscReal mu=1.0,tau=10.0,kappa=5.0;
39: PetscBool terse=PETSC_FALSE;
42: SlepcInitialize(&argc,&argv,(char*)0,help);
44: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
45: PetscOptionsGetReal(NULL,NULL,"-mu",&mu,NULL);
46: PetscOptionsGetReal(NULL,NULL,"-tau",&tau,NULL);
47: PetscOptionsGetReal(NULL,NULL,"-kappa",&kappa,NULL);
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53: /* K is a tridiagonal */
54: MatCreate(PETSC_COMM_WORLD,&K);
55: MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n);
56: MatSetFromOptions(K);
57: MatSetUp(K);
59: MatGetOwnershipRange(K,&Istart,&Iend);
60: for (i=Istart;i<Iend;i++) {
61: if (i>0) MatSetValue(K,i,i-1,-kappa,INSERT_VALUES);
62: MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES);
63: if (i<n-1) MatSetValue(K,i,i+1,-kappa,INSERT_VALUES);
64: }
66: MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
67: MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);
69: /* C is a tridiagonal */
70: MatCreate(PETSC_COMM_WORLD,&C);
71: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);
72: MatSetFromOptions(C);
73: MatSetUp(C);
75: MatGetOwnershipRange(C,&Istart,&Iend);
76: for (i=Istart;i<Iend;i++) {
77: if (i>0) MatSetValue(C,i,i-1,-tau,INSERT_VALUES);
78: MatSetValue(C,i,i,tau*3.0,INSERT_VALUES);
79: if (i<n-1) MatSetValue(C,i,i+1,-tau,INSERT_VALUES);
80: }
82: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
83: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
85: /* M is a diagonal matrix */
86: MatCreate(PETSC_COMM_WORLD,&M);
87: MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n);
88: MatSetFromOptions(M);
89: MatSetUp(M);
90: MatGetOwnershipRange(M,&Istart,&Iend);
91: for (i=Istart;i<Iend;i++) MatSetValue(M,i,i,mu,INSERT_VALUES);
92: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
93: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Create the eigensolver and set various options
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: PEPCreate(PETSC_COMM_WORLD,&pep);
100: A[0] = K; A[1] = C; A[2] = M;
101: PEPSetOperators(pep,3,A);
102: PEPSetProblemType(pep,PEP_GENERAL);
103: PEPSetTolerances(pep,PETSC_SMALL,PETSC_DEFAULT);
104: PEPSetFromOptions(pep);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Solve the eigensystem
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: PEPSolve(pep);
111: PEPGetDimensions(pep,&nev,NULL,NULL);
112: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev);
114: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: Display solution of first solve
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
118: if (terse) PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
119: else {
120: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
121: PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
122: PEPErrorView(pep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
123: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
124: }
125: MatDestroy(&M);
126: MatDestroy(&C);
127: MatDestroy(&K);
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Compute the eigensystem, (k^2*M+k*C+K)x=0 for bigger n
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: n *= 2;
134: /* K is a tridiagonal */
135: MatCreate(PETSC_COMM_WORLD,&K);
136: MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n);
137: MatSetFromOptions(K);
138: MatSetUp(K);
140: MatGetOwnershipRange(K,&Istart,&Iend);
141: for (i=Istart;i<Iend;i++) {
142: if (i>0) MatSetValue(K,i,i-1,-kappa,INSERT_VALUES);
143: MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES);
144: if (i<n-1) MatSetValue(K,i,i+1,-kappa,INSERT_VALUES);
145: }
147: MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
148: MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);
150: /* C is a tridiagonal */
151: MatCreate(PETSC_COMM_WORLD,&C);
152: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);
153: MatSetFromOptions(C);
154: MatSetUp(C);
156: MatGetOwnershipRange(C,&Istart,&Iend);
157: for (i=Istart;i<Iend;i++) {
158: if (i>0) MatSetValue(C,i,i-1,-tau,INSERT_VALUES);
159: MatSetValue(C,i,i,tau*3.0,INSERT_VALUES);
160: if (i<n-1) MatSetValue(C,i,i+1,-tau,INSERT_VALUES);
161: }
163: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
164: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
166: /* M is a diagonal matrix */
167: MatCreate(PETSC_COMM_WORLD,&M);
168: MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n);
169: MatSetFromOptions(M);
170: MatSetUp(M);
171: MatGetOwnershipRange(M,&Istart,&Iend);
172: for (i=Istart;i<Iend;i++) MatSetValue(M,i,i,mu,INSERT_VALUES);
173: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
174: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
176: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177: Solve again, calling PEPReset() since matrix size has changed
178: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
179: /* PEPReset(pep); */ /* not required, will be called in PEPSetOperators() */
180: A[0] = K; A[1] = C; A[2] = M;
181: PEPSetOperators(pep,3,A);
182: PEPSolve(pep);
184: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: Display solution and clean up
186: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
187: if (terse) PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
188: else {
189: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
190: PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
191: PEPErrorView(pep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
192: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
193: }
194: PEPDestroy(&pep);
195: MatDestroy(&M);
196: MatDestroy(&C);
197: MatDestroy(&K);
198: SlepcFinalize();
199: return 0;
200: }
202: /*TEST
204: test:
205: suffix: 1
206: args: -pep_type {{toar qarnoldi linear}} -pep_nev 4 -terse
207: requires: double
209: test:
210: suffix: 2
211: args: -pep_type stoar -pep_hermitian -pep_nev 4 -terse
212: requires: !single
214: TEST*/