Actual source code: ex7.c
slepc-3.7.4 2017-05-17
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2016, 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: "The command line options are:\n"
24: " -f1 <filename> -f2 <filename>, PETSc binary files containing A and B.\n"
25: " -evecs <filename>, output file to save computed eigenvectors.\n"
26: " -ninitial <nini>, number of user-provided initial guesses.\n"
27: " -finitial <filename>, binary file containing <nini> vectors.\n"
28: " -nconstr <ncon>, number of user-provided constraints.\n"
29: " -fconstr <filename>, binary file containing <ncon> vectors.\n\n";
31: #include <slepceps.h>
35: int main(int argc,char **argv)
36: {
37: Mat A,B; /* matrices */
38: EPS eps; /* eigenproblem solver context */
39: ST st;
40: KSP ksp;
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,terse;
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,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,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: MatCreateVecs(A,NULL,&xr);
84: MatCreateVecs(A,NULL,&xi);
86: /*
87: Read user constraints if available
88: */
89: PetscOptionsGetInt(NULL,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,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,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,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: EPSGetST(eps,&st);
155: STGetKSP(st,&ksp);
156: KSPGetTotalIterations(ksp,&lits);
157: PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations of the method: %D\n",lits);
158: EPSGetType(eps,&type);
159: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
160: EPSGetDimensions(eps,&nev,NULL,NULL);
161: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
162: EPSGetTolerances(eps,&tol,&maxit);
163: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);
165: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: Display solution and clean up
167: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: /*
170: Show detailed info unless -terse option is given by user
171: */
172: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
173: if (terse) {
174: EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
175: } else {
176: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
177: EPSReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
178: EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
179: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
180: }
182: /*
183: Save eigenvectors, if requested
184: */
185: PetscOptionsGetString(NULL,NULL,"-evecs",filename,PETSC_MAX_PATH_LEN,&evecs);
186: EPSGetConverged(eps,&nconv);
187: if (nconv>0 && evecs) {
188: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);
189: EPSIsHermitian(eps,&ishermitian);
190: for (i=0;i<nconv;i++) {
191: EPSGetEigenvector(eps,i,xr,xi);
192: VecView(xr,viewer);
193: #if !defined(PETSC_USE_COMPLEX)
194: if (!ishermitian) { VecView(xi,viewer); }
195: #endif
196: }
197: PetscViewerDestroy(&viewer);
198: }
200: /*
201: Free work space
202: */
203: EPSDestroy(&eps);
204: MatDestroy(&A);
205: MatDestroy(&B);
206: VecDestroy(&xr);
207: VecDestroy(&xi);
208: if (nini > 0) {
209: VecDestroyVecs(nini,&Iv);
210: }
211: if (ncon > 0) {
212: VecDestroyVecs(ncon,&Cv);
213: }
214: SlepcFinalize();
215: return ierr;
216: }