Actual source code: ex7.c
slepc-3.17.1 2022-04-11
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;
36: SlepcInitialize(&argc,&argv,(char*)0,help);
38: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39: Load the matrices that define the eigensystem, Ax=kBx
40: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
42: PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized eigenproblem stored in file.\n\n");
43: PetscOptionsGetString(NULL,NULL,"-f1",filename,sizeof(filename),&flg);
46: #if defined(PETSC_USE_COMPLEX)
47: PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n");
48: #else
49: PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n");
50: #endif
51: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
52: MatCreate(PETSC_COMM_WORLD,&A);
53: MatSetFromOptions(A);
54: MatLoad(A,viewer);
55: PetscViewerDestroy(&viewer);
57: PetscOptionsGetString(NULL,NULL,"-f2",filename,sizeof(filename),&flg);
58: if (flg) {
59: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
60: MatCreate(PETSC_COMM_WORLD,&B);
61: MatSetFromOptions(B);
62: MatLoad(B,viewer);
63: PetscViewerDestroy(&viewer);
64: } else {
65: PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n");
66: B = NULL;
67: }
69: MatCreateVecs(A,NULL,&xr);
70: MatCreateVecs(A,NULL,&xi);
72: /*
73: Read user constraints if available
74: */
75: PetscOptionsGetInt(NULL,NULL,"-nconstr",&ncon,&flg);
76: if (flg) {
78: PetscOptionsGetString(NULL,NULL,"-fconstr",filename,sizeof(filename),&flg);
80: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
81: VecDuplicateVecs(xr,ncon,&Cv);
82: for (i=0;i<ncon;i++) VecLoad(Cv[i],viewer);
83: PetscViewerDestroy(&viewer);
84: }
86: /*
87: Read initial guesses if available
88: */
89: PetscOptionsGetInt(NULL,NULL,"-ninitial",&nini,&flg);
90: if (flg) {
92: PetscOptionsGetString(NULL,NULL,"-finitial",filename,sizeof(filename),&flg);
94: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
95: VecDuplicateVecs(xr,nini,&Iv);
96: for (i=0;i<nini;i++) VecLoad(Iv[i],viewer);
97: PetscViewerDestroy(&viewer);
98: }
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Create the eigensolver and set various options
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: /*
105: Create eigensolver context
106: */
107: EPSCreate(PETSC_COMM_WORLD,&eps);
109: /*
110: Set operators. In this case, it is a generalized eigenvalue problem
111: */
112: EPSSetOperators(eps,A,B);
114: /*
115: If the user provided initial guesses or constraints, pass them here
116: */
117: EPSSetInitialSpace(eps,nini,Iv);
118: EPSSetDeflationSpace(eps,ncon,Cv);
120: /*
121: Set solver parameters at runtime
122: */
123: EPSSetFromOptions(eps);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Solve the eigensystem
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129: EPSSolve(eps);
131: /*
132: Optional: Get some information from the solver and display it
133: */
134: EPSGetIterationNumber(eps,&its);
135: PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",its);
136: EPSGetST(eps,&st);
137: STGetKSP(st,&ksp);
138: KSPGetTotalIterations(ksp,&lits);
139: PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations of the method: %" PetscInt_FMT "\n",lits);
140: EPSGetType(eps,&type);
141: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
142: EPSGetDimensions(eps,&nev,NULL,NULL);
143: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev);
144: EPSGetTolerances(eps,&tol,&maxit);
145: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit);
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Display solution and clean up
149: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: /*
152: Show detailed info unless -terse option is given by user
153: */
154: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
155: if (terse) EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
156: else {
157: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
158: EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
159: EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
160: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
161: }
163: /*
164: Save eigenvectors, if requested
165: */
166: PetscOptionsGetString(NULL,NULL,"-evecs",filename,sizeof(filename),&evecs);
167: EPSGetConverged(eps,&nconv);
168: if (nconv>0 && evecs) {
169: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);
170: EPSIsHermitian(eps,&ishermitian);
171: for (i=0;i<nconv;i++) {
172: EPSGetEigenvector(eps,i,xr,xi);
173: VecView(xr,viewer);
174: #if !defined(PETSC_USE_COMPLEX)
175: if (!ishermitian) VecView(xi,viewer);
176: #endif
177: }
178: PetscViewerDestroy(&viewer);
179: }
181: /*
182: Free work space
183: */
184: EPSDestroy(&eps);
185: MatDestroy(&A);
186: MatDestroy(&B);
187: VecDestroy(&xr);
188: VecDestroy(&xi);
189: if (nini > 0) VecDestroyVecs(nini,&Iv);
190: if (ncon > 0) VecDestroyVecs(ncon,&Cv);
191: SlepcFinalize();
192: return 0;
193: }
195: /*TEST
197: test:
198: suffix: 1
199: args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62b.petsc -eps_nev 4 -terse
200: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
202: test:
203: suffix: ciss_1
204: 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
205: requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
206: timeoutfactor: 2
208: test:
209: suffix: 3 # test problem (A,A)
210: args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -eps_nev 4 -terse
211: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
213: TEST*/