Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
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";
19 :
20 : #include <slepceps.h>
21 :
22 2 : int main(int argc,char **argv)
23 : {
24 2 : Mat A,B; /* matrices */
25 2 : EPS eps; /* eigenproblem solver context */
26 2 : ST st;
27 2 : KSP ksp;
28 2 : EPSType type;
29 2 : PetscReal tol;
30 2 : Vec xr,xi,*Iv,*Cv;
31 2 : PetscInt nev,maxit,i,its,lits,nconv,nini=0,ncon=0;
32 2 : char filename[PETSC_MAX_PATH_LEN];
33 2 : PetscViewer viewer;
34 2 : PetscBool flg,evecs,ishermitian,terse;
35 :
36 2 : PetscFunctionBeginUser;
37 2 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
38 :
39 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
40 : Load the matrices that define the eigensystem, Ax=kBx
41 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
42 :
43 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized eigenproblem stored in file.\n\n"));
44 2 : PetscCall(PetscOptionsGetString(NULL,NULL,"-f1",filename,sizeof(filename),&flg));
45 2 : PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -f1 option");
46 :
47 : #if defined(PETSC_USE_COMPLEX)
48 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n"));
49 : #else
50 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n"));
51 : #endif
52 2 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
53 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
54 2 : PetscCall(MatSetFromOptions(A));
55 2 : PetscCall(MatLoad(A,viewer));
56 2 : PetscCall(PetscViewerDestroy(&viewer));
57 :
58 2 : PetscCall(PetscOptionsGetString(NULL,NULL,"-f2",filename,sizeof(filename),&flg));
59 2 : if (flg) {
60 2 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
61 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
62 2 : PetscCall(MatSetFromOptions(B));
63 2 : PetscCall(MatLoad(B,viewer));
64 2 : PetscCall(PetscViewerDestroy(&viewer));
65 : } else {
66 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n"));
67 0 : B = NULL;
68 : }
69 :
70 2 : PetscCall(MatCreateVecs(A,NULL,&xr));
71 2 : PetscCall(MatCreateVecs(A,NULL,&xi));
72 :
73 : /*
74 : Read user constraints if available
75 : */
76 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-nconstr",&ncon,&flg));
77 2 : if (flg) {
78 0 : PetscCheck(ncon>0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"The number of constraints must be >0");
79 0 : PetscCall(PetscOptionsGetString(NULL,NULL,"-fconstr",filename,sizeof(filename),&flg));
80 0 : PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must specify the name of the file storing the constraints");
81 0 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
82 0 : PetscCall(VecDuplicateVecs(xr,ncon,&Cv));
83 0 : for (i=0;i<ncon;i++) PetscCall(VecLoad(Cv[i],viewer));
84 0 : PetscCall(PetscViewerDestroy(&viewer));
85 : }
86 :
87 : /*
88 : Read initial guesses if available
89 : */
90 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-ninitial",&nini,&flg));
91 2 : if (flg) {
92 0 : PetscCheck(nini>0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"The number of initial vectors must be >0");
93 0 : PetscCall(PetscOptionsGetString(NULL,NULL,"-finitial",filename,sizeof(filename),&flg));
94 0 : PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must specify the name of the file containing the initial vectors");
95 0 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
96 0 : PetscCall(VecDuplicateVecs(xr,nini,&Iv));
97 0 : for (i=0;i<nini;i++) PetscCall(VecLoad(Iv[i],viewer));
98 0 : PetscCall(PetscViewerDestroy(&viewer));
99 : }
100 :
101 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102 : Create the eigensolver and set various options
103 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104 :
105 : /*
106 : Create eigensolver context
107 : */
108 2 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
109 :
110 : /*
111 : Set operators. In this case, it is a generalized eigenvalue problem
112 : */
113 2 : PetscCall(EPSSetOperators(eps,A,B));
114 :
115 : /*
116 : If the user provided initial guesses or constraints, pass them here
117 : */
118 2 : PetscCall(EPSSetInitialSpace(eps,nini,Iv));
119 2 : PetscCall(EPSSetDeflationSpace(eps,ncon,Cv));
120 :
121 : /*
122 : Set solver parameters at runtime
123 : */
124 2 : PetscCall(EPSSetFromOptions(eps));
125 :
126 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127 : Solve the eigensystem
128 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129 :
130 2 : PetscCall(EPSSolve(eps));
131 :
132 : /*
133 : Optional: Get some information from the solver and display it
134 : */
135 2 : PetscCall(EPSGetIterationNumber(eps,&its));
136 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",its));
137 2 : PetscCall(EPSGetST(eps,&st));
138 2 : PetscCall(STGetKSP(st,&ksp));
139 2 : PetscCall(KSPGetTotalIterations(ksp,&lits));
140 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations of the method: %" PetscInt_FMT "\n",lits));
141 2 : PetscCall(EPSGetType(eps,&type));
142 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
143 2 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
144 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
145 2 : PetscCall(EPSGetTolerances(eps,&tol,&maxit));
146 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit));
147 :
148 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149 : Display solution and clean up
150 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151 :
152 : /*
153 : Show detailed info unless -terse option is given by user
154 : */
155 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
156 2 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
157 : else {
158 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
159 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
160 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
161 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
162 : }
163 :
164 : /*
165 : Save eigenvectors, if requested
166 : */
167 2 : PetscCall(PetscOptionsGetString(NULL,NULL,"-evecs",filename,sizeof(filename),&evecs));
168 2 : PetscCall(EPSGetConverged(eps,&nconv));
169 2 : if (nconv>0 && evecs) {
170 0 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer));
171 0 : PetscCall(EPSIsHermitian(eps,&ishermitian));
172 0 : for (i=0;i<nconv;i++) {
173 0 : PetscCall(EPSGetEigenvector(eps,i,xr,xi));
174 0 : PetscCall(VecView(xr,viewer));
175 : #if !defined(PETSC_USE_COMPLEX)
176 0 : if (!ishermitian) PetscCall(VecView(xi,viewer));
177 : #endif
178 : }
179 0 : PetscCall(PetscViewerDestroy(&viewer));
180 : }
181 :
182 : /*
183 : Free work space
184 : */
185 2 : PetscCall(EPSDestroy(&eps));
186 2 : PetscCall(MatDestroy(&A));
187 2 : PetscCall(MatDestroy(&B));
188 2 : PetscCall(VecDestroy(&xr));
189 2 : PetscCall(VecDestroy(&xi));
190 2 : if (nini > 0) PetscCall(VecDestroyVecs(nini,&Iv));
191 2 : if (ncon > 0) PetscCall(VecDestroyVecs(ncon,&Cv));
192 2 : PetscCall(SlepcFinalize());
193 : return 0;
194 : }
195 :
196 : /*TEST
197 :
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)
202 :
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
208 :
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)
213 :
214 : TEST*/
|
