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 GNHEP problem with contour integral. "
12 : "Based on ex7.\n"
13 : "The command line options are:\n"
14 : " -f1 <filename> -f2 <filename>, PETSc binary files containing A and B.\n\n";
15 :
16 : #include <slepceps.h>
17 :
18 1 : int main(int argc,char **argv)
19 : {
20 1 : EPS eps;
21 1 : RG rg;
22 1 : Mat A,B;
23 1 : PetscBool flg;
24 1 : EPSCISSExtraction extr;
25 1 : EPSCISSQuadRule quad;
26 1 : char filename[PETSC_MAX_PATH_LEN];
27 1 : PetscViewer viewer;
28 :
29 1 : PetscFunctionBeginUser;
30 1 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
31 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGNHEP problem with contour integral\n\n"));
32 :
33 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
34 : Load the matrices that define the eigensystem, Ax=kBx
35 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
36 :
37 1 : PetscCall(PetscOptionsGetString(NULL,NULL,"-f1",filename,sizeof(filename),&flg));
38 1 : PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -f1 option");
39 :
40 : #if defined(PETSC_USE_COMPLEX)
41 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n"));
42 : #else
43 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n"));
44 : #endif
45 1 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
46 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
47 1 : PetscCall(MatSetFromOptions(A));
48 1 : PetscCall(MatLoad(A,viewer));
49 1 : PetscCall(PetscViewerDestroy(&viewer));
50 :
51 1 : PetscCall(PetscOptionsGetString(NULL,NULL,"-f2",filename,sizeof(filename),&flg));
52 1 : if (flg) {
53 0 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
54 0 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
55 0 : PetscCall(MatSetFromOptions(B));
56 0 : PetscCall(MatLoad(B,viewer));
57 0 : PetscCall(PetscViewerDestroy(&viewer));
58 : } else {
59 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n"));
60 1 : B = NULL;
61 : }
62 :
63 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64 : Create the eigensolver and solve the problem
65 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66 :
67 1 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
68 1 : PetscCall(EPSSetOperators(eps,A,B));
69 1 : PetscCall(EPSSetProblemType(eps,EPS_GNHEP));
70 1 : PetscCall(EPSSetTolerances(eps,1e-9,PETSC_DEFAULT));
71 :
72 : /* set CISS solver with various options */
73 1 : PetscCall(EPSSetType(eps,EPSCISS));
74 1 : PetscCall(EPSCISSSetExtraction(eps,EPS_CISS_EXTRACTION_HANKEL));
75 1 : PetscCall(EPSCISSSetQuadRule(eps,EPS_CISS_QUADRULE_CHEBYSHEV));
76 1 : PetscCall(EPSCISSSetUseST(eps,PETSC_TRUE));
77 1 : PetscCall(EPSGetRG(eps,&rg));
78 1 : PetscCall(RGSetType(rg,RGINTERVAL));
79 1 : PetscCall(RGIntervalSetEndpoints(rg,-3000.0,0.0,0.0,0.0));
80 :
81 1 : PetscCall(EPSSetFromOptions(eps));
82 :
83 1 : PetscCall(EPSSolve(eps));
84 1 : PetscCall(PetscObjectTypeCompare((PetscObject)eps,EPSCISS,&flg));
85 1 : if (flg) {
86 1 : PetscCall(EPSCISSGetExtraction(eps,&extr));
87 1 : PetscCall(EPSCISSGetQuadRule(eps,&quad));
88 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solved with CISS using %s extraction with %s quadrature rule\n\n",EPSCISSExtractions[extr],EPSCISSQuadRules[quad]));
89 : }
90 :
91 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92 : Display solution and clean up
93 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94 :
95 1 : PetscCall(EPSErrorView(eps,EPS_ERROR_BACKWARD,NULL));
96 1 : PetscCall(EPSDestroy(&eps));
97 1 : PetscCall(MatDestroy(&A));
98 1 : PetscCall(MatDestroy(&B));
99 1 : PetscCall(SlepcFinalize());
100 : return 0;
101 : }
102 :
103 : /*TEST
104 :
105 : testset:
106 : args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62b.petsc
107 : output_file: output/test25_1.out
108 : test:
109 : suffix: 1
110 : requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
111 : test:
112 : suffix: 1_cuda
113 : args: -mat_type aijcusparse -st_pc_factor_mat_solver_type cusparse
114 : requires: cuda double !complex !defined(PETSC_USE_64BIT_INDICES)
115 :
116 : testset:
117 : args: -f1 ${DATAFILESPATH}/matrices/complex/qc324.petsc -rg_type ellipse -rg_ellipse_center 1-0.09i -rg_ellipse_radius 0.15 -rg_ellipse_vscale 0.1
118 : output_file: output/test25_2.out
119 : test:
120 : suffix: 2
121 : requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
122 : test:
123 : suffix: 2_cuda
124 : args: -mat_type aijcusparse -st_pc_factor_mat_solver_type cusparse
125 : requires: cuda double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
126 :
127 : TEST*/
|
