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