Actual source code: ex57.c
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: */
11: static char help[] = "Another eigenvalue problem with Hamiltonian structure.\n\n"
12: "Position and velocity control of a string of high-speed vehicles.\n"
13: "The command line options are:\n"
14: " -m <m>, where <n> = number of vehicles.\n\n";
16: #include <slepceps.h>
18: /*
19: This example computes eigenvalues of a matrix
21: H = [ A N
22: K -A^* ],
24: where A, N and K are n-by-n matrices, with n=2*m-1 and
26: N = diag( 1, 0, 1, 0,..., 0, 1)
27: K = diag( 0,10, 0,10,...,10, 0)
28: A = tridiag(b,d,c)
29: d = [-1, 0,-1, 0,...,-1, 0,-1]
30: b = [ 1, 0, 1, 0,..., 1, 0]
31: c = [ 0,-1, 0,-1,..., 0,-1]
33: References:
35: [1] W.R. Ferng, Wen-Wei Lin, Chern-Shuh Wang, The shift-inverted J-Lanczos algorithm
36: for the numerical solutions of large sparse algebraic Riccati equations,
37: Computers & Mathematics with Applications, Volume 33, Issue 10, 1997,
38: https://doi.org/10.1016/S0898-1221(97)00074-6.
40: */
42: static PetscErrorCode eigenCompare (PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *res,void *ctx)
43: {
44: PetscReal abs1, abs2, tol=1e-12, r_ar, r_ai, r_br, r_bi;
46: PetscFunctionBegin;
47: #if defined(PETSC_USE_COMPLEX)
48: r_ar = PetscRealPart(ar);
49: r_ai = PetscImaginaryPart(ar);
50: r_br = PetscRealPart(br);
51: r_bi = PetscImaginaryPart(br);
52: #else
53: r_ar = ar;
54: r_ai = ai;
55: r_br = br;
56: r_bi = bi;
57: #endif
58: if (PetscAbs(PetscAbs(r_ar)-PetscAbs(r_br)) < tol && PetscAbs(PetscAbs(r_ai)-PetscAbs(r_bi)) < tol) {
59: /* Negative real part first*/
60: if (r_ar<0 && r_br>0)
61: *res = -1;
62: else if (r_ar>0 && r_br<0)
63: *res = 1;
64: /* Positive imaginary part first*/
65: else if (r_ai>0 && r_bi<0)
66: *res = -1;
67: else if (r_ai<0 && r_bi>0)
68: *res = 1;
69: }
70: else {
71: abs1 = SlepcAbs(r_ar,r_ai);
72: abs2 = SlepcAbs(r_br,r_bi);
73: if (abs1 > abs2)
74: *res = -1;
75: else if (abs1 < abs2)
76: *res = 1;
77: else
78: *res = 0;
79: }
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: int main(int argc,char **argv)
84: {
85: Mat H,A,N,K; /* problem matrices */
86: EPS eps; /* eigenproblem solver context */
87: PetscInt m=12,n,Istart,Iend,i;
88: PetscBool terse, sort_hamilt;
90: PetscFunctionBeginUser;
91: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
93: PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
94: n = 2*m-1;
95: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nHamiltonian eigenproblem, n=%" PetscInt_FMT
96: " (m=%" PetscInt_FMT " vehicles)\n\n",n,m));
98: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: Compute the problem matrices A, N and K
100: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
103: PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
104: PetscCall(MatSetFromOptions(A));
106: PetscCall(MatCreate(PETSC_COMM_WORLD,&N));
107: PetscCall(MatSetSizes(N,PETSC_DECIDE,PETSC_DECIDE,n,n));
108: PetscCall(MatSetFromOptions(N));
110: PetscCall(MatCreate(PETSC_COMM_WORLD,&K));
111: PetscCall(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n));
112: PetscCall(MatSetFromOptions(K));
114: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
115: for (i=Istart;i<Iend;i++) {
116: if (i%2==0)
117: PetscCall(MatSetValue(A,i,i,-1.0,INSERT_VALUES));
118: else {
119: if (i>0) PetscCall(MatSetValue(A,i,i-1,1.0,INSERT_VALUES));
120: if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
121: }
122: }
124: PetscCall(MatGetOwnershipRange(N,&Istart,&Iend));
125: for (i=Istart+Istart%2;i<Iend;i+=2) {
126: PetscCall(MatSetValue(N,i,i,1.0,INSERT_VALUES));
127: }
129: PetscCall(MatGetOwnershipRange(K,&Istart,&Iend));
130: for (i=Istart+(Istart+1)%2;i<Iend;i+=2) {
131: PetscCall(MatSetValue(K,i,i,10.0,INSERT_VALUES));
132: }
134: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
135: PetscCall(MatAssemblyBegin(N,MAT_FINAL_ASSEMBLY));
136: PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
137: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
138: PetscCall(MatAssemblyEnd(N,MAT_FINAL_ASSEMBLY));
139: PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));
141: PetscCall(MatCreateHamiltonian(A,N,K,&H));
143: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: Create the eigensolver and set various options
145: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
148: PetscCall(EPSSetOperators(eps,H,NULL));
149: PetscCall(EPSSetProblemType(eps,EPS_HAMILT));
151: PetscCall(PetscOptionsHasName(NULL,NULL,"-sort_hamilt",&sort_hamilt));
152: if (sort_hamilt) {
153: /* Adjust ordering of non-hermitian solver so that it is the same as in as in EPS_HAMILT solver */
154: PetscCall(EPSSetEigenvalueComparison(eps,eigenCompare,NULL));
155: PetscCall(EPSSetWhichEigenpairs(eps,EPS_WHICH_USER));
156: }
157: PetscCall(EPSSetFromOptions(eps));
159: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: Solve the eigensystem and display solution
161: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163: PetscCall(EPSSolve(eps));
165: /* show detailed info unless -terse option is given by user */
166: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
167: if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
168: else {
169: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
170: PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
171: PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
172: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
173: }
175: PetscCall(EPSDestroy(&eps));
176: PetscCall(MatDestroy(&A));
177: PetscCall(MatDestroy(&N));
178: PetscCall(MatDestroy(&K));
179: PetscCall(MatDestroy(&H));
180: PetscCall(SlepcFinalize());
181: return 0;
182: }
184: /*TEST
186: testset:
187: args: -eps_nev 8 -eps_ncv 28 -terse
188: nsize: {{1 2}}
189: output_file: output/ex57_1.out
190: test:
191: requires: double !complex
192: suffix: 1
193: test:
194: requires: double !complex
195: args: -eps_non_hermitian -sort_hamilt
196: suffix: 1_nhep
197: test:
198: requires: double complex
199: TODO: no support for complex scalars yet
200: suffix: 1_complex
201: test:
202: requires: double complex
203: args: -eps_non_hermitian -sort_hamilt
204: suffix: 1_complex_nhep
206: testset:
207: args: -eps_nev 4 -eps_smallest_magnitude -eps_ncv 32 -terse
208: output_file: output/ex57_2.out
209: test:
210: requires: double !complex
211: suffix: 2
213: TEST*/