Actual source code: sleeper.c
slepc-3.18.2 2023-01-26
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;
36: SlepcInitialize(&argc,&argv,(char*)0,help);
38: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
39: PetscPrintf(PETSC_COMM_WORLD,"\nRailtrack resting on sleepers, n=%" PetscInt_FMT "\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++) MatSetValue(M,i,i,1.0,INSERT_VALUES);
108: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
109: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
111: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: Create the eigensolver and solve the problem
113: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: PEPCreate(PETSC_COMM_WORLD,&pep);
116: A[0] = K; A[1] = C; A[2] = M;
117: PEPSetOperators(pep,3,A);
118: PEPSetFromOptions(pep);
119: PEPSolve(pep);
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Display solution and clean up
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: /* show detailed info unless -terse option is given by user */
126: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
127: if (terse) PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
128: else {
129: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
130: PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
131: PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD);
132: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
133: }
134: PEPDestroy(&pep);
135: MatDestroy(&M);
136: MatDestroy(&C);
137: MatDestroy(&K);
138: SlepcFinalize();
139: return 0;
140: }
142: /*TEST
144: testset:
145: args: -n 100 -pep_nev 4 -pep_ncv 24 -st_type sinvert -terse
146: output_file: output/sleeper_1.out
147: requires: !single
148: filter: sed -e "s/[+-]0\.0*i//g"
149: test:
150: suffix: 1
151: args: -pep_type {{toar linear}} -pep_ncv 20
152: test:
153: suffix: 1_qarnoldi
154: args: -pep_type qarnoldi -pep_qarnoldi_restart 0.4
156: testset:
157: args: -n 24 -pep_nev 4 -pep_ncv 9 -pep_target -.62 -terse
158: output_file: output/sleeper_2.out
159: test:
160: suffix: 2_toar
161: args: -pep_type toar -pep_toar_restart .3 -st_type sinvert
162: test:
163: suffix: 2_jd
164: args: -pep_type jd -pep_jd_restart .3 -pep_jd_projection orthogonal
166: test:
167: suffix: 3
168: args: -n 275 -pep_type stoar -pep_hermitian -st_type sinvert -pep_nev 2 -pep_target -.89 -terse
169: requires: !single
171: test:
172: suffix: 4
173: args: -n 270 -pep_type stoar -pep_hermitian -pep_interval -3,-2.51 -st_type sinvert -st_pc_type cholesky -terse
174: requires: !single
176: TEST*/