Actual source code: test1.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2013, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
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(NULL,"-n",&n,NULL);
40: PetscOptionsGetInt(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: MatSetUp(A);
54: MatCreate(PETSC_COMM_WORLD,&B);
55: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,N,N);
56: MatSetFromOptions(B);
57: MatSetUp(B);
59: MatGetOwnershipRange(A,&Istart,&Iend);
60: for (II=Istart;II<Iend;II++) {
61: i = II/n; j = II-i*n;
62: if (i>0) { MatSetValue(A,II,II-n,-1.0,INSERT_VALUES); }
63: if (i<m-1) { MatSetValue(A,II,II+n,-1.0,INSERT_VALUES); }
64: if (j>0) { MatSetValue(A,II,II-1,-1.0,INSERT_VALUES); }
65: if (j<n-1) { MatSetValue(A,II,II+1,-1.0,INSERT_VALUES); }
66: MatSetValue(A,II,II,4.0,INSERT_VALUES);
67: MatSetValue(B,II,II,2.0/PetscLogScalar(II+2),INSERT_VALUES);
68: }
70: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
71: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
72: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
73: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
74: MatGetVecs(B,&v,NULL);
76: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77: Create the eigensolver and set various options
78: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80: EPSCreate(PETSC_COMM_WORLD,&eps);
81: EPSSetOperators(eps,A,B);
82: EPSSetProblemType(eps,EPS_GHEP);
83: EPSSetTolerances(eps,tol,PETSC_DECIDE);
84: EPSSetFromOptions(eps);
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Solve the eigensystem
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: EPSSolve(eps);
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Display solution and clean up
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: EPSPrintSolution(eps,NULL);
97: EPSGetConverged(eps,&nconv);
98: if (nconv>0) {
99: VecDuplicateVecs(v,nconv,&X);
100: for (i=0;i<nconv;i++) {
101: EPSGetEigenvector(eps,i,X[i],NULL);
102: }
103: SlepcCheckOrthogonality(X,nconv,NULL,nconv,B,NULL,&lev);
104: if (lev<10*tol) {
105: PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality below the tolerance\n");
106: } else {
107: PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %G\n",lev);
108: }
109: VecDestroyVecs(nconv,&X);
110: }
112: EPSDestroy(&eps);
113: MatDestroy(&A);
114: MatDestroy(&B);
115: VecDestroy(&v);
116: SlepcFinalize();
117: return 0;
118: }