Actual source code: sleeper.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: This example implements one of the problems found at
12: NLEVP: A Collection of Nonlinear Eigenvalue Problems,
13: The University of Manchester.
14: The details of the collection can be found at:
15: [1] T. Betcke et al., "NLEVP: A Collection of Nonlinear Eigenvalue
16: Problems", ACM Trans. Math. Software 39(2), Article 7, 2013.
18: The sleeper problem is a proportionally damped QEP describing the
19: oscillations of a rail track resting on sleepers.
20: */
22: static char help[] = "Oscillations of a rail track resting on sleepers.\n\n"
23: "The command line options are:\n"
24: " -n <n>, where <n> = dimension of the matrices.\n\n";
26: #include <slepcpep.h>
28: int main(int argc,char **argv)
29: {
30: Mat M,C,K,A[3]; /* problem matrices */
31: PEP pep; /* polynomial eigenproblem solver context */
32: PetscInt n=10,Istart,Iend,i;
33: PetscBool terse;
36: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
38: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
39: PetscPrintf(PETSC_COMM_WORLD,"\nRailtrack resting on sleepers, n=%D\n\n",n);
41: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
43: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
45: /* K is a pentadiagonal */
46: MatCreate(PETSC_COMM_WORLD,&K);
47: MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n);
48: MatSetFromOptions(K);
49: MatSetUp(K);
51: MatGetOwnershipRange(K,&Istart,&Iend);
52: for (i=Istart;i<Iend;i++) {
53: if (i==0) {
54: MatSetValue(K,i,n-1,-3.0,INSERT_VALUES);
55: MatSetValue(K,i,n-2,1.0,INSERT_VALUES);
56: }
57: if (i==1) { MatSetValue(K,i,n-1,1.0,INSERT_VALUES); }
58: if (i>0) { MatSetValue(K,i,i-1,-3.0,INSERT_VALUES); }
59: if (i>1) { MatSetValue(K,i,i-2,1.0,INSERT_VALUES); }
60: MatSetValue(K,i,i,5.0,INSERT_VALUES);
61: if (i==n-1) {
62: MatSetValue(K,i,0,-3.0,INSERT_VALUES);
63: MatSetValue(K,i,1,1.0,INSERT_VALUES);
64: }
65: if (i==n-2) { MatSetValue(K,i,0,1.0,INSERT_VALUES); }
66: if (i<n-1) { MatSetValue(K,i,i+1,-3.0,INSERT_VALUES); }
67: if (i<n-2) { MatSetValue(K,i,i+2,1.0,INSERT_VALUES); }
68: }
70: MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
71: MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);
73: /* C is a circulant matrix */
74: MatCreate(PETSC_COMM_WORLD,&C);
75: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);
76: MatSetFromOptions(C);
77: MatSetUp(C);
79: MatGetOwnershipRange(C,&Istart,&Iend);
80: for (i=Istart;i<Iend;i++) {
81: if (i==0) {
82: MatSetValue(C,i,n-1,-4.0,INSERT_VALUES);
83: MatSetValue(C,i,n-2,1.0,INSERT_VALUES);
84: }
85: if (i==1) { MatSetValue(C,i,n-1,1.0,INSERT_VALUES); }
86: if (i>0) { MatSetValue(C,i,i-1,-4.0,INSERT_VALUES); }
87: if (i>1) { MatSetValue(C,i,i-2,1.0,INSERT_VALUES); }
88: MatSetValue(C,i,i,7.0,INSERT_VALUES);
89: if (i==n-1) {
90: MatSetValue(C,i,0,-4.0,INSERT_VALUES);
91: MatSetValue(C,i,1,1.0,INSERT_VALUES);
92: }
93: if (i==n-2) { MatSetValue(C,i,0,1.0,INSERT_VALUES); }
94: if (i<n-1) { MatSetValue(C,i,i+1,-4.0,INSERT_VALUES); }
95: if (i<n-2) { MatSetValue(C,i,i+2,1.0,INSERT_VALUES); }
96: }
98: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
99: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
101: /* M is the identity matrix */
102: MatCreate(PETSC_COMM_WORLD,&M);
103: MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n);
104: MatSetFromOptions(M);
105: MatSetUp(M);
106: MatGetOwnershipRange(M,&Istart,&Iend);
107: for (i=Istart;i<Iend;i++) {
108: MatSetValue(M,i,i,1.0,INSERT_VALUES);
109: }
110: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
111: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Create the eigensolver and solve the problem
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: PEPCreate(PETSC_COMM_WORLD,&pep);
118: A[0] = K; A[1] = C; A[2] = M;
119: PEPSetOperators(pep,3,A);
120: PEPSetFromOptions(pep);
121: PEPSolve(pep);
123: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: Display solution and clean up
125: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127: /* show detailed info unless -terse option is given by user */
128: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
129: if (terse) {
130: PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
131: } else {
132: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
133: PEPReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
134: PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD);
135: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
136: }
137: PEPDestroy(&pep);
138: MatDestroy(&M);
139: MatDestroy(&C);
140: MatDestroy(&K);
141: SlepcFinalize();
142: return ierr;
143: }
145: /*TEST
147: testset:
148: args: -n 100 -pep_nev 4 -pep_ncv 24 -st_type sinvert -terse
149: output_file: output/sleeper_1.out
150: test:
151: suffix: 1
152: args: -pep_type {{toar linear}} -pep_ncv 20
153: test:
154: suffix: 1_qarnoldi
155: args: -pep_type qarnoldi -pep_qarnoldi_restart 0.4
157: testset:
158: args: -n 24 -pep_nev 4 -pep_ncv 9 -pep_target -.62 -terse
159: output_file: output/sleeper_2.out
160: test:
161: suffix: 2_toar
162: args: -pep_type toar -pep_toar_restart .3 -st_type sinvert
163: requires: !single
164: test:
165: suffix: 2_jd
166: args: -pep_type jd -pep_jd_restart .3 -pep_jd_projection orthogonal
168: test:
169: suffix: 3
170: args: -n 275 -pep_type stoar -pep_hermitian -st_type sinvert -pep_nev 2 -pep_target -.89 -terse
171: requires: !single
173: test:
174: suffix: 4
175: args: -n 270 -pep_type stoar -pep_hermitian -pep_interval -3,-2.51 -st_type sinvert -st_pc_type cholesky -terse
176: requires: !complex !single
178: TEST*/