Actual source code: ex7.c
slepc-3.16.1 2021-11-17
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: static char help[] = "Solves a generalized eigensystem Ax=kBx with matrices loaded from a file.\n"
12: "The command line options are:\n"
13: " -f1 <filename> -f2 <filename>, PETSc binary files containing A and B.\n"
14: " -evecs <filename>, output file to save computed eigenvectors.\n"
15: " -ninitial <nini>, number of user-provided initial guesses.\n"
16: " -finitial <filename>, binary file containing <nini> vectors.\n"
17: " -nconstr <ncon>, number of user-provided constraints.\n"
18: " -fconstr <filename>, binary file containing <ncon> vectors.\n\n";
20: #include <slepceps.h>
22: int main(int argc,char **argv)
23: {
24: Mat A,B; /* matrices */
25: EPS eps; /* eigenproblem solver context */
26: ST st;
27: KSP ksp;
28: EPSType type;
29: PetscReal tol;
30: Vec xr,xi,*Iv,*Cv;
31: PetscInt nev,maxit,i,its,lits,nconv,nini=0,ncon=0;
32: char filename[PETSC_MAX_PATH_LEN];
33: PetscViewer viewer;
34: PetscBool flg,evecs,ishermitian,terse;
37: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
39: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
40: Load the matrices that define the eigensystem, Ax=kBx
41: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
43: PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized eigenproblem stored in file.\n\n");
44: PetscOptionsGetString(NULL,NULL,"-f1",filename,sizeof(filename),&flg);
45: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -f1 option");
47: #if defined(PETSC_USE_COMPLEX)
48: PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n");
49: #else
50: PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n");
51: #endif
52: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
53: MatCreate(PETSC_COMM_WORLD,&A);
54: MatSetFromOptions(A);
55: MatLoad(A,viewer);
56: PetscViewerDestroy(&viewer);
58: PetscOptionsGetString(NULL,NULL,"-f2",filename,sizeof(filename),&flg);
59: if (flg) {
60: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
61: MatCreate(PETSC_COMM_WORLD,&B);
62: MatSetFromOptions(B);
63: MatLoad(B,viewer);
64: PetscViewerDestroy(&viewer);
65: } else {
66: PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n");
67: B = NULL;
68: }
70: MatCreateVecs(A,NULL,&xr);
71: MatCreateVecs(A,NULL,&xi);
73: /*
74: Read user constraints if available
75: */
76: PetscOptionsGetInt(NULL,NULL,"-nconstr",&ncon,&flg);
77: if (flg) {
78: if (ncon<=0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"The number of constraints must be >0");
79: PetscOptionsGetString(NULL,NULL,"-fconstr",filename,sizeof(filename),&flg);
80: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must specify the name of the file storing the constraints");
81: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
82: VecDuplicateVecs(xr,ncon,&Cv);
83: for (i=0;i<ncon;i++) {
84: VecLoad(Cv[i],viewer);
85: }
86: PetscViewerDestroy(&viewer);
87: }
89: /*
90: Read initial guesses if available
91: */
92: PetscOptionsGetInt(NULL,NULL,"-ninitial",&nini,&flg);
93: if (flg) {
94: if (nini<=0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"The number of initial vectors must be >0");
95: PetscOptionsGetString(NULL,NULL,"-finitial",filename,sizeof(filename),&flg);
96: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must specify the name of the file containing the initial vectors");
97: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
98: VecDuplicateVecs(xr,nini,&Iv);
99: for (i=0;i<nini;i++) {
100: VecLoad(Iv[i],viewer);
101: }
102: PetscViewerDestroy(&viewer);
103: }
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Create the eigensolver and set various options
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: /*
110: Create eigensolver context
111: */
112: EPSCreate(PETSC_COMM_WORLD,&eps);
114: /*
115: Set operators. In this case, it is a generalized eigenvalue problem
116: */
117: EPSSetOperators(eps,A,B);
119: /*
120: If the user provided initial guesses or constraints, pass them here
121: */
122: EPSSetInitialSpace(eps,nini,Iv);
123: EPSSetDeflationSpace(eps,ncon,Cv);
125: /*
126: Set solver parameters at runtime
127: */
128: EPSSetFromOptions(eps);
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Solve the eigensystem
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: EPSSolve(eps);
136: /*
137: Optional: Get some information from the solver and display it
138: */
139: EPSGetIterationNumber(eps,&its);
140: PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
141: EPSGetST(eps,&st);
142: STGetKSP(st,&ksp);
143: KSPGetTotalIterations(ksp,&lits);
144: PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations of the method: %D\n",lits);
145: EPSGetType(eps,&type);
146: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
147: EPSGetDimensions(eps,&nev,NULL,NULL);
148: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
149: EPSGetTolerances(eps,&tol,&maxit);
150: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);
152: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: Display solution and clean up
154: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
156: /*
157: Show detailed info unless -terse option is given by user
158: */
159: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
160: if (terse) {
161: EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
162: } else {
163: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
164: EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
165: EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
166: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
167: }
169: /*
170: Save eigenvectors, if requested
171: */
172: PetscOptionsGetString(NULL,NULL,"-evecs",filename,sizeof(filename),&evecs);
173: EPSGetConverged(eps,&nconv);
174: if (nconv>0 && evecs) {
175: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);
176: EPSIsHermitian(eps,&ishermitian);
177: for (i=0;i<nconv;i++) {
178: EPSGetEigenvector(eps,i,xr,xi);
179: VecView(xr,viewer);
180: #if !defined(PETSC_USE_COMPLEX)
181: if (!ishermitian) { VecView(xi,viewer); }
182: #endif
183: }
184: PetscViewerDestroy(&viewer);
185: }
187: /*
188: Free work space
189: */
190: EPSDestroy(&eps);
191: MatDestroy(&A);
192: MatDestroy(&B);
193: VecDestroy(&xr);
194: VecDestroy(&xi);
195: if (nini > 0) {
196: VecDestroyVecs(nini,&Iv);
197: }
198: if (ncon > 0) {
199: VecDestroyVecs(ncon,&Cv);
200: }
201: SlepcFinalize();
202: return ierr;
203: }
205: /*TEST
207: test:
208: suffix: 1
209: args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62b.petsc -eps_nev 4 -terse
210: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
212: test:
213: suffix: ciss_1
214: args: -f1 ${DATAFILESPATH}/matrices/complex/mhd1280a.petsc -f2 ${DATAFILESPATH}/matrices/complex/mhd1280b.petsc -eps_type ciss -eps_ciss_usest 0 -eps_ciss_quadrule chebyshev -rg_type ring -rg_ring_center 0 -rg_ring_radius .49 -rg_ring_width 0.2 -rg_ring_startangle .25 -rg_ring_endangle .49 -terse -eps_max_it 1
215: requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
216: timeoutfactor: 2
218: test:
219: suffix: 3 # test problem (A,A)
220: args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -eps_nev 4 -terse
221: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
223: TEST*/