Actual source code: test1.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7:
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: static char help[] = "Tests B-orthonormality of eigenvectors in a GHEP problem.\n\n";
24: #include <slepceps.h>
28: int main( int argc, char **argv )
29: {
30: Mat A,B; /* matrices */
31: EPS eps; /* eigenproblem solver context */
32: Vec *X,v;
33: PetscReal lev,tol=1000*PETSC_MACHINE_EPSILON;
34: PetscInt N,n=45,m,Istart,Iend,II,i,j,nconv;
35: PetscBool flag;
38: SlepcInitialize(&argc,&argv,(char*)0,help);
39: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
40: PetscOptionsGetInt(PETSC_NULL,"-m",&m,&flag);
41: if(!flag) m=n;
42: N = n*m;
43: PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized Symmetric Eigenproblem, N=%D (%Dx%D grid)\n\n",N,n,m);
45: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46: Compute the matrices that define the eigensystem, Ax=kBx
47: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
49: MatCreate(PETSC_COMM_WORLD,&A);
50: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
51: MatSetFromOptions(A);
52:
53: MatCreate(PETSC_COMM_WORLD,&B);
54: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,N,N);
55: MatSetFromOptions(B);
57: MatGetOwnershipRange(A,&Istart,&Iend);
58: for( II=Istart; II<Iend; II++ ) {
59: i = II/n; j = II-i*n;
60: if(i>0) { MatSetValue(A,II,II-n,-1.0,INSERT_VALUES); }
61: if(i<m-1) { MatSetValue(A,II,II+n,-1.0,INSERT_VALUES); }
62: if(j>0) { MatSetValue(A,II,II-1,-1.0,INSERT_VALUES); }
63: if(j<n-1) { MatSetValue(A,II,II+1,-1.0,INSERT_VALUES); }
64: MatSetValue(A,II,II,4.0,INSERT_VALUES);
65: MatSetValue(B,II,II,2.0/PetscLogScalar(II+2),INSERT_VALUES);
66: }
68: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
70: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
71: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
72: MatGetVecs(B,&v,PETSC_NULL);
74: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: Create the eigensolver and set various options
76: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78: EPSCreate(PETSC_COMM_WORLD,&eps);
79: EPSSetOperators(eps,A,B);
80: EPSSetProblemType(eps,EPS_GHEP);
81: EPSSetTolerances(eps,tol,PETSC_DECIDE);
82: EPSSetFromOptions(eps);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Solve the eigensystem
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: EPSSolve(eps);
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Display solution and clean up
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94: EPSPrintSolution(eps,PETSC_NULL);
95: EPSGetConverged(eps,&nconv);
96: if (nconv>0) {
97: VecDuplicateVecs(v,nconv,&X);
98: for( i=0; i<nconv; i++ ) {
99: EPSGetEigenvector(eps,i,X[i],PETSC_NULL);
100: }
101: SlepcCheckOrthogonality(X,nconv,PETSC_NULL,nconv,B,&lev);
102: if (lev<tol) {
103: PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality below the tolerance\n");
104: } else {
105: PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %G\n",lev);
106: }
107: VecDestroyVecs(nconv,&X);
108: }
109:
110: EPSDestroy(&eps);
111: MatDestroy(&A);
112: MatDestroy(&B);
113: VecDestroy(&v);
114: SlepcFinalize();
115: return 0;
116: }