Actual source code: ex7.c
slepc-3.20.0 2023-09-29
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: PetscFunctionBeginUser;
37: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
39: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
40: Load the matrices that define the eigensystem, Ax=kBx
41: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
43: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized eigenproblem stored in file.\n\n"));
44: PetscCall(PetscOptionsGetString(NULL,NULL,"-f1",filename,sizeof(filename),&flg));
45: PetscCheck(flg,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: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n"));
49: #else
50: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n"));
51: #endif
52: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
53: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
54: PetscCall(MatSetFromOptions(A));
55: PetscCall(MatLoad(A,viewer));
56: PetscCall(PetscViewerDestroy(&viewer));
58: PetscCall(PetscOptionsGetString(NULL,NULL,"-f2",filename,sizeof(filename),&flg));
59: if (flg) {
60: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
61: PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
62: PetscCall(MatSetFromOptions(B));
63: PetscCall(MatLoad(B,viewer));
64: PetscCall(PetscViewerDestroy(&viewer));
65: } else {
66: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n"));
67: B = NULL;
68: }
70: PetscCall(MatCreateVecs(A,NULL,&xr));
71: PetscCall(MatCreateVecs(A,NULL,&xi));
73: /*
74: Read user constraints if available
75: */
76: PetscCall(PetscOptionsGetInt(NULL,NULL,"-nconstr",&ncon,&flg));
77: if (flg) {
78: PetscCheck(ncon>0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"The number of constraints must be >0");
79: PetscCall(PetscOptionsGetString(NULL,NULL,"-fconstr",filename,sizeof(filename),&flg));
80: PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must specify the name of the file storing the constraints");
81: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
82: PetscCall(VecDuplicateVecs(xr,ncon,&Cv));
83: for (i=0;i<ncon;i++) PetscCall(VecLoad(Cv[i],viewer));
84: PetscCall(PetscViewerDestroy(&viewer));
85: }
87: /*
88: Read initial guesses if available
89: */
90: PetscCall(PetscOptionsGetInt(NULL,NULL,"-ninitial",&nini,&flg));
91: if (flg) {
92: PetscCheck(nini>0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"The number of initial vectors must be >0");
93: PetscCall(PetscOptionsGetString(NULL,NULL,"-finitial",filename,sizeof(filename),&flg));
94: PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must specify the name of the file containing the initial vectors");
95: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
96: PetscCall(VecDuplicateVecs(xr,nini,&Iv));
97: for (i=0;i<nini;i++) PetscCall(VecLoad(Iv[i],viewer));
98: PetscCall(PetscViewerDestroy(&viewer));
99: }
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Create the eigensolver and set various options
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: /*
106: Create eigensolver context
107: */
108: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
110: /*
111: Set operators. In this case, it is a generalized eigenvalue problem
112: */
113: PetscCall(EPSSetOperators(eps,A,B));
115: /*
116: If the user provided initial guesses or constraints, pass them here
117: */
118: PetscCall(EPSSetInitialSpace(eps,nini,Iv));
119: PetscCall(EPSSetDeflationSpace(eps,ncon,Cv));
121: /*
122: Set solver parameters at runtime
123: */
124: PetscCall(EPSSetFromOptions(eps));
126: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: Solve the eigensystem
128: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: PetscCall(EPSSolve(eps));
132: /*
133: Optional: Get some information from the solver and display it
134: */
135: PetscCall(EPSGetIterationNumber(eps,&its));
136: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",its));
137: PetscCall(EPSGetST(eps,&st));
138: PetscCall(STGetKSP(st,&ksp));
139: PetscCall(KSPGetTotalIterations(ksp,&lits));
140: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations of the method: %" PetscInt_FMT "\n",lits));
141: PetscCall(EPSGetType(eps,&type));
142: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
143: PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
144: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
145: PetscCall(EPSGetTolerances(eps,&tol,&maxit));
146: PetscCall(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: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
156: if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
157: else {
158: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
159: PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
160: PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
161: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
162: }
164: /*
165: Save eigenvectors, if requested
166: */
167: PetscCall(PetscOptionsGetString(NULL,NULL,"-evecs",filename,sizeof(filename),&evecs));
168: PetscCall(EPSGetConverged(eps,&nconv));
169: if (nconv>0 && evecs) {
170: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer));
171: PetscCall(EPSIsHermitian(eps,&ishermitian));
172: for (i=0;i<nconv;i++) {
173: PetscCall(EPSGetEigenvector(eps,i,xr,xi));
174: PetscCall(VecView(xr,viewer));
175: #if !defined(PETSC_USE_COMPLEX)
176: if (!ishermitian) PetscCall(VecView(xi,viewer));
177: #endif
178: }
179: PetscCall(PetscViewerDestroy(&viewer));
180: }
182: /*
183: Free work space
184: */
185: PetscCall(EPSDestroy(&eps));
186: PetscCall(MatDestroy(&A));
187: PetscCall(MatDestroy(&B));
188: PetscCall(VecDestroy(&xr));
189: PetscCall(VecDestroy(&xi));
190: if (nini > 0) PetscCall(VecDestroyVecs(nini,&Iv));
191: if (ncon > 0) PetscCall(VecDestroyVecs(ncon,&Cv));
192: PetscCall(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*/