Actual source code: ex7.c
slepc-3.18.0 2022-10-01
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, 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);
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);
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) {
79: PetscOptionsGetString(NULL,NULL,"-fconstr",filename,sizeof(filename),&flg);
81: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
82: VecDuplicateVecs(xr,ncon,&Cv);
83: for (i=0;i<ncon;i++) VecLoad(Cv[i],viewer);
84: PetscViewerDestroy(&viewer);
85: }
87: /*
88: Read initial guesses if available
89: */
90: PetscOptionsGetInt(NULL,NULL,"-ninitial",&nini,&flg);
91: if (flg) {
93: PetscOptionsGetString(NULL,NULL,"-finitial",filename,sizeof(filename),&flg);
95: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
96: VecDuplicateVecs(xr,nini,&Iv);
97: for (i=0;i<nini;i++) VecLoad(Iv[i],viewer);
98: PetscViewerDestroy(&viewer);
99: }
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Create the eigensolver and set various options
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: /*
106: Create eigensolver context
107: */
108: EPSCreate(PETSC_COMM_WORLD,&eps);
110: /*
111: Set operators. In this case, it is a generalized eigenvalue problem
112: */
113: EPSSetOperators(eps,A,B);
115: /*
116: If the user provided initial guesses or constraints, pass them here
117: */
118: EPSSetInitialSpace(eps,nini,Iv);
119: EPSSetDeflationSpace(eps,ncon,Cv);
121: /*
122: Set solver parameters at runtime
123: */
124: EPSSetFromOptions(eps);
126: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: Solve the eigensystem
128: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: EPSSolve(eps);
132: /*
133: Optional: Get some information from the solver and display it
134: */
135: EPSGetIterationNumber(eps,&its);
136: PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",its);
137: EPSGetST(eps,&st);
138: STGetKSP(st,&ksp);
139: KSPGetTotalIterations(ksp,&lits);
140: PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations of the method: %" PetscInt_FMT "\n",lits);
141: EPSGetType(eps,&type);
142: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
143: EPSGetDimensions(eps,&nev,NULL,NULL);
144: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev);
145: EPSGetTolerances(eps,&tol,&maxit);
146: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit);
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Display solution and clean up
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: /*
153: Show detailed info unless -terse option is given by user
154: */
155: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
156: if (terse) EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
157: else {
158: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
159: EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
160: EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
161: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
162: }
164: /*
165: Save eigenvectors, if requested
166: */
167: PetscOptionsGetString(NULL,NULL,"-evecs",filename,sizeof(filename),&evecs);
168: EPSGetConverged(eps,&nconv);
169: if (nconv>0 && evecs) {
170: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);
171: EPSIsHermitian(eps,&ishermitian);
172: for (i=0;i<nconv;i++) {
173: EPSGetEigenvector(eps,i,xr,xi);
174: VecView(xr,viewer);
175: #if !defined(PETSC_USE_COMPLEX)
176: if (!ishermitian) VecView(xi,viewer);
177: #endif
178: }
179: PetscViewerDestroy(&viewer);
180: }
182: /*
183: Free work space
184: */
185: EPSDestroy(&eps);
186: MatDestroy(&A);
187: MatDestroy(&B);
188: VecDestroy(&xr);
189: VecDestroy(&xi);
190: if (nini > 0) VecDestroyVecs(nini,&Iv);
191: if (ncon > 0) VecDestroyVecs(ncon,&Cv);
192: SlepcFinalize();
193: return 0;
194: }
196: /*TEST
198: test:
199: suffix: 1
200: args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62b.petsc -eps_nev 4 -terse
201: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
203: test:
204: suffix: ciss_1
205: 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
206: requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
207: timeoutfactor: 2
209: test:
210: suffix: 3 # test problem (A,A)
211: args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -eps_nev 4 -terse
212: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
214: TEST*/