Actual source code: ex82.c
1: static const char help[] = "Uses KSPComputeRitz() on a matrix loaded from disk\n";
3: #include <petscksp.h>
5: int main(int argc, char **argv)
6: {
7: Mat A;
8: KSP ksp;
9: char file[PETSC_MAX_PATH_LEN];
10: PetscReal *tetar, *tetai;
11: Vec b, x, *S;
12: PetscInt i, N = 10, Na = N;
13: PetscViewer fd;
14: PC pc;
15: PetscBool harmonic = PETSC_FALSE;
17: PetscFunctionBeginUser;
18: PetscCall(PetscInitialize(&argc, &argv, 0, help));
20: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), NULL));
21: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
22: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
23: PetscCall(MatLoad(A, fd));
24: PetscCall(PetscViewerDestroy(&fd));
25: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
27: PetscCall(MatCreateVecs(A, &b, &x));
28: PetscCall(VecSetRandom(b, NULL));
29: PetscCall(VecDuplicateVecs(b, N, &S));
30: PetscCall(PetscMalloc2(N, &tetar, N, &tetai));
32: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
33: PetscCall(KSPSetType(ksp, KSPGMRES));
34: PetscCall(KSPSetTolerances(ksp, 10000 * PETSC_MACHINE_EPSILON, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
35: PetscCall(KSPGetPC(ksp, &pc));
36: PetscCall(PCSetType(pc, PCNONE));
37: PetscCall(KSPSetComputeRitz(ksp, PETSC_TRUE));
38: PetscCall(KSPSetOperators(ksp, A, A));
39: PetscCall(KSPSetFromOptions(ksp));
40: PetscCall(KSPSolve(ksp, b, x));
41: PetscCall(PetscOptionsHasName(NULL, NULL, "-harmonic", &harmonic));
42: PetscCall(KSPComputeRitz(ksp, harmonic ? PETSC_FALSE : PETSC_TRUE, PETSC_TRUE, &Na, S, tetar, tetai));
44: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%% Number of Ritz pairs %" PetscInt_FMT "\n", Na));
45: for (i = 0; i < Na; i++) {
46: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%% Eigenvalue(s) %g ", (double)tetar[i]));
47: if (tetai[i]) {
48: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%+gi", (double)tetai[i]));
49: #if !defined(PETSC_USE_COMPLEX)
50: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %g ", (double)tetar[i]));
51: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%+gi", (double)-tetai[i]));
52: #endif
53: }
54: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
55: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%% Eigenvector\n"));
56: PetscCall(VecView(S[i], PETSC_VIEWER_STDOUT_WORLD));
57: #if !defined(PETSC_USE_COMPLEX)
58: if (tetai[i]) {
59: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%% Imaginary part of Eigenvector\n"));
60: PetscCall(VecView(S[i + 1], PETSC_VIEWER_STDOUT_WORLD));
61: i++;
62: }
63: #endif
64: }
66: PetscCall(PetscFree2(tetar, tetai));
67: PetscCall(VecDestroy(&x));
68: PetscCall(VecDestroy(&b));
69: PetscCall(VecDestroyVecs(N, &S));
70: PetscCall(KSPDestroy(&ksp));
71: PetscCall(MatDestroy(&A));
72: PetscCall(PetscFinalize());
73: return 0;
74: }
76: /*TEST
78: test:
79: requires: !defined(PETSC_USE_64BIT_INDICES) !complex double
80: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/ritz2 -ksp_monitor
82: test:
83: suffix: 2
84: requires: !defined(PETSC_USE_64BIT_INDICES) !complex double
85: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/ritz5 -ksp_monitor
87: test:
88: suffix: harmonic
89: requires: !defined(PETSC_USE_64BIT_INDICES) !complex double
90: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/ritz5 -ksp_monitor -harmonic
92: TEST*/