Actual source code: ex7.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[] = "Solves a generalized eigensystem Ax=kBx with matrices loaded from a file.\n"
23: "This example works for both real and complex numbers.\n\n"
24: "The command line options are:\n"
25: " -f1 <filename>, where <filename> = matrix (A) file in PETSc binary form.\n"
26: " -f2 <filename>, where <filename> = matrix (B) file in PETSc binary form.\n"
27: " -evecs <filename>, output file to save computed eigenvectors.\n"
28: " -ninitial <nini>, number of user-provided initial guesses.\n"
29: " -finitial <filename>, where <filename> contains <nini> vectors (binary).\n"
30: " -nconstr <ncon>, number of user-provided constraints.\n"
31: " -fconstr <filename>, where <filename> contains <ncon> vectors (binary).\n\n";
33: #include <slepceps.h>
37: int main(int argc,char **argv)
38: {
39: Mat A,B; /* matrices */
40: EPS eps; /* eigenproblem solver context */
41: EPSType type;
42: PetscReal tol;
43: Vec xr,xi,*Iv,*Cv;
44: PetscInt nev,maxit,i,its,lits,nconv,nini=0,ncon=0;
45: char filename[PETSC_MAX_PATH_LEN];
46: PetscViewer viewer;
47: PetscBool flg,evecs,ishermitian;
50: SlepcInitialize(&argc,&argv,(char*)0,help);
52: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: Load the matrices that define the eigensystem, Ax=kBx
54: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56: PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized eigenproblem stored in file.\n\n");
57: PetscOptionsGetString(NULL,"-f1",filename,PETSC_MAX_PATH_LEN,&flg);
58: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for matrix A with the -f1 option");
60: #if defined(PETSC_USE_COMPLEX)
61: PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n");
62: #else
63: PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n");
64: #endif
65: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
66: MatCreate(PETSC_COMM_WORLD,&A);
67: MatSetFromOptions(A);
68: MatLoad(A,viewer);
69: PetscViewerDestroy(&viewer);
71: PetscOptionsGetString(NULL,"-f2",filename,PETSC_MAX_PATH_LEN,&flg);
72: if (flg) {
73: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
74: MatCreate(PETSC_COMM_WORLD,&B);
75: MatSetFromOptions(B);
76: MatLoad(B,viewer);
77: PetscViewerDestroy(&viewer);
78: } else {
79: PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n");
80: B = NULL;
81: }
83: MatGetVecs(A,NULL,&xr);
84: MatGetVecs(A,NULL,&xi);
86: /*
87: Read user constraints if available
88: */
89: PetscOptionsGetInt(NULL,"-nconstr",&ncon,&flg);
90: if (flg) {
91: if (ncon<=0) SETERRQ(PETSC_COMM_WORLD,1,"The number of constraints must be >0");
92: PetscOptionsGetString(NULL,"-fconstr",filename,PETSC_MAX_PATH_LEN,&flg);
93: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must specify the name of the file storing the constraints");
94: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
95: VecDuplicateVecs(xr,ncon,&Cv);
96: for (i=0;i<ncon;i++) {
97: VecLoad(Cv[i],viewer);
98: }
99: PetscViewerDestroy(&viewer);
100: }
102: /*
103: Read initial guesses if available
104: */
105: PetscOptionsGetInt(NULL,"-ninitial",&nini,&flg);
106: if (flg) {
107: if (nini<=0) SETERRQ(PETSC_COMM_WORLD,1,"The number of initial vectors must be >0");
108: PetscOptionsGetString(NULL,"-finitial",filename,PETSC_MAX_PATH_LEN,&flg);
109: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must specify the name of the file containing the initial vectors");
110: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
111: VecDuplicateVecs(xr,nini,&Iv);
112: for (i=0;i<nini;i++) {
113: VecLoad(Iv[i],viewer);
114: }
115: PetscViewerDestroy(&viewer);
116: }
118: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: Create the eigensolver and set various options
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122: /*
123: Create eigensolver context
124: */
125: EPSCreate(PETSC_COMM_WORLD,&eps);
127: /*
128: Set operators. In this case, it is a generalized eigenvalue problem
129: */
130: EPSSetOperators(eps,A,B);
132: /*
133: If the user provided initial guesses or constraints, pass them here
134: */
135: EPSSetInitialSpace(eps,nini,Iv);
136: EPSSetDeflationSpace(eps,ncon,Cv);
138: /*
139: Set solver parameters at runtime
140: */
141: EPSSetFromOptions(eps);
143: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: Solve the eigensystem
145: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: EPSSolve(eps);
149: /*
150: Optional: Get some information from the solver and display it
151: */
152: EPSGetIterationNumber(eps,&its);
153: PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
154: EPSGetOperationCounters(eps,NULL,NULL,&lits);
155: PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations of the method: %D\n",lits);
156: EPSGetType(eps,&type);
157: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
158: EPSGetDimensions(eps,&nev,NULL,NULL);
159: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
160: EPSGetTolerances(eps,&tol,&maxit);
161: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4G, maxit=%D\n",tol,maxit);
163: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: Display solution and clean up
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167: EPSPrintSolution(eps,NULL);
168: /*
169: Save eigenvectors, if requested
170: */
171: PetscOptionsGetString(NULL,"-evecs",filename,PETSC_MAX_PATH_LEN,&evecs);
172: EPSGetConverged(eps,&nconv);
173: if (nconv>0 && evecs) {
174: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);
175: EPSIsHermitian(eps,&ishermitian);
176: for (i=0;i<nconv;i++) {
177: EPSGetEigenvector(eps,i,xr,xi);
178: VecView(xr,viewer);
179: if (!ishermitian) { VecView(xi,viewer); }
180: }
181: PetscViewerDestroy(&viewer);
182: }
184: /*
185: Free work space
186: */
187: EPSDestroy(&eps);
188: MatDestroy(&A);
189: MatDestroy(&B);
190: VecDestroy(&xr);
191: VecDestroy(&xi);
192: if (nini > 0) {
193: VecDestroyVecs(nini,&Iv);
194: }
195: if (ncon > 0) {
196: VecDestroyVecs(ncon,&Cv);
197: }
198: SlepcFinalize();
199: return 0;
200: }