Actual source code: test13.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[] = "Test EPSSetArbitrarySelection.\n\n";
24: #include <slepceps.h>
28: PetscErrorCode MyArbitrarySelection(PetscScalar eigr,PetscScalar eigi,Vec xr,Vec xi,PetscScalar *rr,PetscScalar *ri,void *ctx)
29: {
30: PetscErrorCode ierr;
31: Vec xref = *(Vec*)ctx;
34: VecDot(xr,xref,rr);
35: *rr = PetscAbsScalar(*rr);
36: if (ri) *ri = 0.0;
37: return(0);
38: }
42: int main(int argc,char **argv)
43: {
44: Mat A; /* problem matrices */
45: EPS eps; /* eigenproblem solver context */
46: PetscScalar seigr,seigi,value[3];
47: PetscReal tol=1000*PETSC_MACHINE_EPSILON;
48: Vec sxr,sxi;
49: PetscInt n=30,i,Istart,Iend,col[3],nconv;
50: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
53: SlepcInitialize(&argc,&argv,(char*)0,help);
55: PetscOptionsGetInt(NULL,"-n",&n,NULL);
56: PetscPrintf(PETSC_COMM_WORLD,"\nTridiagonal with zero diagonal, n=%D\n\n",n);
58: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: Create matrix tridiag([-1 0 -1])
60: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61: MatCreate(PETSC_COMM_WORLD,&A);
62: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
63: MatSetFromOptions(A);
64: MatSetUp(A);
66: MatGetOwnershipRange(A,&Istart,&Iend);
67: if (Istart==0) FirstBlock=PETSC_TRUE;
68: if (Iend==n) LastBlock=PETSC_TRUE;
69: value[0]=-1.0; value[1]=0.0; value[2]=-1.0;
70: for (i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++) {
71: col[0]=i-1; col[1]=i; col[2]=i+1;
72: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
73: }
74: if (LastBlock) {
75: i=n-1; col[0]=n-2; col[1]=n-1;
76: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
77: }
78: if (FirstBlock) {
79: i=0; col[0]=0; col[1]=1; value[0]=0.0; value[1]=-1.0;
80: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
81: }
83: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
84: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Create the eigensolver
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: EPSCreate(PETSC_COMM_WORLD,&eps);
90: EPSSetProblemType(eps,EPS_HEP);
91: EPSSetTolerances(eps,tol,PETSC_DECIDE);
92: EPSSetOperators(eps,A,NULL);
93: EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);
94: EPSSetFromOptions(eps);
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Solve eigenproblem and store some solution
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: EPSSolve(eps);
100: MatGetVecs(A,&sxr,NULL);
101: MatGetVecs(A,&sxi,NULL);
102: EPSGetConverged(eps,&nconv);
103: if (nconv>0) {
104: EPSGetEigenpair(eps,0,&seigr,&seigi,sxr,sxi);
105: EPSPrintSolution(eps,NULL);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Solve eigenproblem using an arbitrary selection
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: EPSSetArbitrarySelection(eps,MyArbitrarySelection,&sxr);
111: EPSSetWhichEigenpairs(eps,EPS_LARGEST_MAGNITUDE);
112: EPSSolve(eps);
113: EPSPrintSolution(eps,NULL);
114: } else {
115: PetscPrintf(PETSC_COMM_WORLD,"Problem: no eigenpairs converged.\n");
116: }
118: EPSDestroy(&eps);
119: VecDestroy(&sxr);
120: VecDestroy(&sxi);
121: MatDestroy(&A);
122: SlepcFinalize();
123: return 0;
124: }