Actual source code: ex56.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[] = "Eigenvalue problem with Hamiltonian structure.\n\n"
12: "The command line options are:\n"
13: " -n <n>, where <n> = dimension of the blocks.\n\n";
15: #include <slepceps.h>
17: /*
18: This example computes eigenvalues of a matrix
20: H = [ A B
21: C -A^* ],
23: where B and C are Hermitian or real symmetric. In particular, A, B and C have
24: the following Toeplitz structure:
26: A = pentadiag{a,b,c,d,e}
27: B = tridiag{a,c,conj(a)}
28: C = tridiag{b,e,conj(b)}
30: where a,b,d are complex scalars, and c is real.
31: */
33: int main(int argc,char **argv)
34: {
35: Mat H,A,B,C; /* problem matrices */
36: EPS eps; /* eigenproblem solver context */
37: PetscScalar a,b,c,d,e;
38: PetscInt n=24,Istart,Iend,i;
39: PetscBool terse;
41: PetscFunctionBeginUser;
42: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
44: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
45: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nHamiltonian eigenproblem, n=%" PetscInt_FMT "\n\n",n));
47: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: Compute the problem matrices R and C
49: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51: #if defined(PETSC_USE_COMPLEX)
52: a = PetscCMPLX(-0.1,0.2);
53: b = PetscCMPLX(1.0,0.5);
54: d = PetscCMPLX(2.0,0.2);
55: #else
56: a = -0.1;
57: b = 1.0;
58: d = 2.0;
59: #endif
60: c = 4.5;
61: e = -2.5;
63: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
64: PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
65: PetscCall(MatSetFromOptions(A));
67: PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
68: PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n));
69: PetscCall(MatSetFromOptions(B));
71: PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
72: PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
73: PetscCall(MatSetFromOptions(C));
75: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
76: for (i=Istart;i<Iend;i++) {
77: if (i>1) PetscCall(MatSetValue(A,i,i-2,a,INSERT_VALUES));
78: if (i>0) PetscCall(MatSetValue(A,i,i-1,b,INSERT_VALUES));
79: PetscCall(MatSetValue(A,i,i,c,INSERT_VALUES));
80: if (i<n-1) PetscCall(MatSetValue(A,i,i+1,d,INSERT_VALUES));
81: if (i<n-2) PetscCall(MatSetValue(A,i,i+2,e,INSERT_VALUES));
82: }
84: PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
85: for (i=Istart;i<Iend;i++) {
86: if (i>0) PetscCall(MatSetValue(B,i,i-1,a,INSERT_VALUES));
87: PetscCall(MatSetValue(B,i,i,c,INSERT_VALUES));
88: if (i<n-1) PetscCall(MatSetValue(B,i,i+1,PetscConj(a),INSERT_VALUES));
89: }
91: PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
92: for (i=Istart;i<Iend;i++) {
93: if (i>0) PetscCall(MatSetValue(C,i,i-1,b,INSERT_VALUES));
94: PetscCall(MatSetValue(C,i,i,e,INSERT_VALUES));
95: if (i<n-1) PetscCall(MatSetValue(C,i,i+1,PetscConj(b),INSERT_VALUES));
96: }
98: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
99: PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
100: PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
101: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
102: PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
103: PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
105: PetscCall(MatCreateHamiltonian(A,B,C,&H));
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Create the eigensolver and set various options
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
112: PetscCall(EPSSetOperators(eps,H,NULL));
113: PetscCall(EPSSetProblemType(eps,EPS_HAMILT));
114: PetscCall(EPSSetFromOptions(eps));
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Solve the eigensystem and display solution
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: PetscCall(EPSSolve(eps));
122: /* show detailed info unless -terse option is given by user */
123: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
124: if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
125: else {
126: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
127: PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
128: PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
129: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
130: }
132: PetscCall(EPSDestroy(&eps));
133: PetscCall(MatDestroy(&A));
134: PetscCall(MatDestroy(&B));
135: PetscCall(MatDestroy(&C));
136: PetscCall(MatDestroy(&H));
137: PetscCall(SlepcFinalize());
138: return 0;
139: }
141: /*TEST
143: testset:
144: args: -eps_nev 8 -eps_ncv 28 -terse
145: nsize: {{1 2}}
146: requires: double !complex
147: output_file: output/ex56_1.out
148: test:
149: suffix: 1
150: test:
151: args: -eps_non_hermitian
152: suffix: 1_nhep
154: testset:
155: args: -eps_nev 4 -eps_ncv 16 -terse
156: nsize: {{1 2}}
157: requires: double complex
158: output_file: output/ex56_1_complex.out
159: test:
160: TODO: no support for complex scalars yet
161: suffix: 1_complex
162: test:
163: args: -eps_non_hermitian
164: suffix: 1_complex_nhep
166: TEST*/