Actual source code: sleeper.c
slepc-3.17.2 2022-08-09
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: 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;
35: SlepcInitialize(&argc,&argv,(char*)0,help);
37: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
38: PetscPrintf(PETSC_COMM_WORLD,"\nRailtrack resting on sleepers, n=%" PetscInt_FMT "\n\n",n);
40: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41: Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
42: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
44: /* K is a pentadiagonal */
45: MatCreate(PETSC_COMM_WORLD,&K);
46: MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n);
47: MatSetFromOptions(K);
48: MatSetUp(K);
50: MatGetOwnershipRange(K,&Istart,&Iend);
51: for (i=Istart;i<Iend;i++) {
52: if (i==0) {
53: MatSetValue(K,i,n-1,-3.0,INSERT_VALUES);
54: MatSetValue(K,i,n-2,1.0,INSERT_VALUES);
55: }
56: if (i==1) MatSetValue(K,i,n-1,1.0,INSERT_VALUES);
57: if (i>0) MatSetValue(K,i,i-1,-3.0,INSERT_VALUES);
58: if (i>1) MatSetValue(K,i,i-2,1.0,INSERT_VALUES);
59: MatSetValue(K,i,i,5.0,INSERT_VALUES);
60: if (i==n-1) {
61: MatSetValue(K,i,0,-3.0,INSERT_VALUES);
62: MatSetValue(K,i,1,1.0,INSERT_VALUES);
63: }
64: if (i==n-2) MatSetValue(K,i,0,1.0,INSERT_VALUES);
65: if (i<n-1) MatSetValue(K,i,i+1,-3.0,INSERT_VALUES);
66: if (i<n-2) MatSetValue(K,i,i+2,1.0,INSERT_VALUES);
67: }
69: MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
70: MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);
72: /* C is a circulant matrix */
73: MatCreate(PETSC_COMM_WORLD,&C);
74: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);
75: MatSetFromOptions(C);
76: MatSetUp(C);
78: MatGetOwnershipRange(C,&Istart,&Iend);
79: for (i=Istart;i<Iend;i++) {
80: if (i==0) {
81: MatSetValue(C,i,n-1,-4.0,INSERT_VALUES);
82: MatSetValue(C,i,n-2,1.0,INSERT_VALUES);
83: }
84: if (i==1) MatSetValue(C,i,n-1,1.0,INSERT_VALUES);
85: if (i>0) MatSetValue(C,i,i-1,-4.0,INSERT_VALUES);
86: if (i>1) MatSetValue(C,i,i-2,1.0,INSERT_VALUES);
87: MatSetValue(C,i,i,7.0,INSERT_VALUES);
88: if (i==n-1) {
89: MatSetValue(C,i,0,-4.0,INSERT_VALUES);
90: MatSetValue(C,i,1,1.0,INSERT_VALUES);
91: }
92: if (i==n-2) MatSetValue(C,i,0,1.0,INSERT_VALUES);
93: if (i<n-1) MatSetValue(C,i,i+1,-4.0,INSERT_VALUES);
94: if (i<n-2) MatSetValue(C,i,i+2,1.0,INSERT_VALUES);
95: }
97: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
98: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
100: /* M is the identity matrix */
101: MatCreate(PETSC_COMM_WORLD,&M);
102: MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n);
103: MatSetFromOptions(M);
104: MatSetUp(M);
105: MatGetOwnershipRange(M,&Istart,&Iend);
106: for (i=Istart;i<Iend;i++) MatSetValue(M,i,i,1.0,INSERT_VALUES);
107: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
108: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Create the eigensolver and solve the problem
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114: PEPCreate(PETSC_COMM_WORLD,&pep);
115: A[0] = K; A[1] = C; A[2] = M;
116: PEPSetOperators(pep,3,A);
117: PEPSetFromOptions(pep);
118: PEPSolve(pep);
120: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: Display solution and clean up
122: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: /* show detailed info unless -terse option is given by user */
125: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
126: if (terse) PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
127: else {
128: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
129: PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
130: PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD);
131: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
132: }
133: PEPDestroy(&pep);
134: MatDestroy(&M);
135: MatDestroy(&C);
136: MatDestroy(&K);
137: SlepcFinalize();
138: return 0;
139: }
141: /*TEST
143: testset:
144: args: -n 100 -pep_nev 4 -pep_ncv 24 -st_type sinvert -terse
145: output_file: output/sleeper_1.out
146: requires: !single
147: filter: sed -e "s/[+-]0\.0*i//g"
148: test:
149: suffix: 1
150: args: -pep_type {{toar linear}} -pep_ncv 20
151: test:
152: suffix: 1_qarnoldi
153: args: -pep_type qarnoldi -pep_qarnoldi_restart 0.4
155: testset:
156: args: -n 24 -pep_nev 4 -pep_ncv 9 -pep_target -.62 -terse
157: output_file: output/sleeper_2.out
158: test:
159: suffix: 2_toar
160: args: -pep_type toar -pep_toar_restart .3 -st_type sinvert
161: test:
162: suffix: 2_jd
163: args: -pep_type jd -pep_jd_restart .3 -pep_jd_projection orthogonal
165: test:
166: suffix: 3
167: args: -n 275 -pep_type stoar -pep_hermitian -st_type sinvert -pep_nev 2 -pep_target -.89 -terse
168: requires: !single
170: test:
171: suffix: 4
172: args: -n 270 -pep_type stoar -pep_hermitian -pep_interval -3,-2.51 -st_type sinvert -st_pc_type cholesky -terse
173: requires: !single
175: TEST*/