Actual source code: test25.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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 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";
16: #include <slepceps.h>
18: int main(int argc,char **argv)
19: {
20: EPS eps;
21: RG rg;
22: Mat A,B;
23: PetscBool flg;
24: EPSCISSExtraction extr;
25: EPSCISSQuadRule quad;
26: char filename[PETSC_MAX_PATH_LEN];
27: PetscViewer viewer;
28: PetscErrorCode ierr;
30: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
31: PetscPrintf(PETSC_COMM_WORLD,"\nGNHEP problem with contour integral\n\n");
33: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
34: Load the matrices that define the eigensystem, Ax=kBx
35: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
37: PetscOptionsGetString(NULL,NULL,"-f1",filename,PETSC_MAX_PATH_LEN,&flg);
38: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for matrix A with the -f1 option");
40: #if defined(PETSC_USE_COMPLEX)
41: PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n");
42: #else
43: PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n");
44: #endif
45: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
46: MatCreate(PETSC_COMM_WORLD,&A);
47: MatSetFromOptions(A);
48: MatLoad(A,viewer);
49: PetscViewerDestroy(&viewer);
51: PetscOptionsGetString(NULL,NULL,"-f2",filename,PETSC_MAX_PATH_LEN,&flg);
52: if (flg) {
53: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
54: MatCreate(PETSC_COMM_WORLD,&B);
55: MatSetFromOptions(B);
56: MatLoad(B,viewer);
57: PetscViewerDestroy(&viewer);
58: } else {
59: PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n");
60: B = NULL;
61: }
63: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: Create the eigensolver and solve the problem
65: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67: EPSCreate(PETSC_COMM_WORLD,&eps);
68: EPSSetOperators(eps,A,B);
69: EPSSetProblemType(eps,EPS_GNHEP);
70: EPSSetTolerances(eps,1e-9,PETSC_DEFAULT);
72: /* set CISS solver with various options */
73: EPSSetType(eps,EPSCISS);
74: EPSCISSSetExtraction(eps,EPS_CISS_EXTRACTION_HANKEL);
75: EPSCISSSetQuadRule(eps,EPS_CISS_QUADRULE_CHEBYSHEV);
76: EPSCISSSetUseST(eps,PETSC_TRUE);
77: EPSGetRG(eps,&rg);
78: RGSetType(rg,RGINTERVAL);
79: RGIntervalSetEndpoints(rg,-3000.0,0.0,0.0,0.0);
81: EPSSetFromOptions(eps);
83: EPSSolve(eps);
84: PetscObjectTypeCompare((PetscObject)eps,EPSCISS,&flg);
85: if (flg) {
86: EPSCISSGetExtraction(eps,&extr);
87: EPSCISSGetQuadRule(eps,&quad);
88: PetscPrintf(PETSC_COMM_WORLD," Solved with CISS using %s extraction with %s quadrature rule\n\n",EPSCISSExtractions[extr],EPSCISSQuadRules[quad]);
89: }
91: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: Display solution and clean up
93: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95: EPSErrorView(eps,EPS_ERROR_BACKWARD,NULL);
96: EPSDestroy(&eps);
97: MatDestroy(&A);
98: MatDestroy(&B);
99: SlepcFinalize();
100: return ierr;
101: }
103: /*TEST
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 !define(PETSC_USE_64BIT_INDICES)
111: test:
112: suffix: 2_cuda
113: args: -mat_type aijcusparse
114: requires: cuda double !complex !define(PETSC_USE_64BIT_INDICES)
116: TEST*/