Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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.
17 :
18 : The sleeper problem is a proportionally damped QEP describing the
19 : oscillations of a rail track resting on sleepers.
20 : */
21 :
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";
25 :
26 : #include <slepcpep.h>
27 :
28 7 : int main(int argc,char **argv)
29 : {
30 7 : Mat M,C,K,A[3]; /* problem matrices */
31 7 : PEP pep; /* polynomial eigenproblem solver context */
32 7 : PetscInt n=10,Istart,Iend,i;
33 7 : PetscBool terse;
34 :
35 7 : PetscFunctionBeginUser;
36 7 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
37 :
38 7 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
39 7 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nRailtrack resting on sleepers, n=%" PetscInt_FMT "\n\n",n));
40 :
41 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42 : Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
43 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
44 :
45 : /* K is a pentadiagonal */
46 7 : PetscCall(MatCreate(PETSC_COMM_WORLD,&K));
47 7 : PetscCall(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n));
48 7 : PetscCall(MatSetFromOptions(K));
49 :
50 7 : PetscCall(MatGetOwnershipRange(K,&Istart,&Iend));
51 900 : for (i=Istart;i<Iend;i++) {
52 893 : if (i==0) {
53 7 : PetscCall(MatSetValue(K,i,n-1,-3.0,INSERT_VALUES));
54 7 : PetscCall(MatSetValue(K,i,n-2,1.0,INSERT_VALUES));
55 : }
56 893 : if (i==1) PetscCall(MatSetValue(K,i,n-1,1.0,INSERT_VALUES));
57 893 : if (i>0) PetscCall(MatSetValue(K,i,i-1,-3.0,INSERT_VALUES));
58 886 : if (i>1) PetscCall(MatSetValue(K,i,i-2,1.0,INSERT_VALUES));
59 893 : PetscCall(MatSetValue(K,i,i,5.0,INSERT_VALUES));
60 893 : if (i==n-1) {
61 7 : PetscCall(MatSetValue(K,i,0,-3.0,INSERT_VALUES));
62 7 : PetscCall(MatSetValue(K,i,1,1.0,INSERT_VALUES));
63 : }
64 893 : if (i==n-2) PetscCall(MatSetValue(K,i,0,1.0,INSERT_VALUES));
65 893 : if (i<n-1) PetscCall(MatSetValue(K,i,i+1,-3.0,INSERT_VALUES));
66 893 : if (i<n-2) PetscCall(MatSetValue(K,i,i+2,1.0,INSERT_VALUES));
67 : }
68 :
69 7 : PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
70 7 : PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));
71 :
72 : /* C is a circulant matrix */
73 7 : PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
74 7 : PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
75 7 : PetscCall(MatSetFromOptions(C));
76 :
77 7 : PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
78 900 : for (i=Istart;i<Iend;i++) {
79 893 : if (i==0) {
80 7 : PetscCall(MatSetValue(C,i,n-1,-4.0,INSERT_VALUES));
81 7 : PetscCall(MatSetValue(C,i,n-2,1.0,INSERT_VALUES));
82 : }
83 893 : if (i==1) PetscCall(MatSetValue(C,i,n-1,1.0,INSERT_VALUES));
84 893 : if (i>0) PetscCall(MatSetValue(C,i,i-1,-4.0,INSERT_VALUES));
85 886 : if (i>1) PetscCall(MatSetValue(C,i,i-2,1.0,INSERT_VALUES));
86 893 : PetscCall(MatSetValue(C,i,i,7.0,INSERT_VALUES));
87 893 : if (i==n-1) {
88 7 : PetscCall(MatSetValue(C,i,0,-4.0,INSERT_VALUES));
89 7 : PetscCall(MatSetValue(C,i,1,1.0,INSERT_VALUES));
90 : }
91 893 : if (i==n-2) PetscCall(MatSetValue(C,i,0,1.0,INSERT_VALUES));
92 893 : if (i<n-1) PetscCall(MatSetValue(C,i,i+1,-4.0,INSERT_VALUES));
93 893 : if (i<n-2) PetscCall(MatSetValue(C,i,i+2,1.0,INSERT_VALUES));
94 : }
95 :
96 7 : PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
97 7 : PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
98 :
99 : /* M is the identity matrix */
100 7 : PetscCall(MatCreate(PETSC_COMM_WORLD,&M));
101 7 : PetscCall(MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n));
102 7 : PetscCall(MatSetFromOptions(M));
103 7 : PetscCall(MatGetOwnershipRange(M,&Istart,&Iend));
104 900 : for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(M,i,i,1.0,INSERT_VALUES));
105 7 : PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
106 7 : PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));
107 :
108 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109 : Create the eigensolver and solve the problem
110 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111 :
112 7 : PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
113 7 : A[0] = K; A[1] = C; A[2] = M;
114 7 : PetscCall(PEPSetOperators(pep,3,A));
115 7 : PetscCall(PEPSetFromOptions(pep));
116 7 : PetscCall(PEPSolve(pep));
117 :
118 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119 : Display solution and clean up
120 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121 :
122 : /* show detailed info unless -terse option is given by user */
123 7 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
124 7 : if (terse) PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
125 : else {
126 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
127 0 : PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
128 0 : PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD));
129 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
130 : }
131 7 : PetscCall(PEPDestroy(&pep));
132 7 : PetscCall(MatDestroy(&M));
133 7 : PetscCall(MatDestroy(&C));
134 7 : PetscCall(MatDestroy(&K));
135 7 : PetscCall(SlepcFinalize());
136 : return 0;
137 : }
138 :
139 : /*TEST
140 :
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
152 :
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
162 :
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
167 :
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
172 :
173 : TEST*/
|
