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;
18: PetscInitialize(&argc, &argv, 0, help);
20: PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), NULL);
21: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd);
22: MatCreate(PETSC_COMM_WORLD, &A);
23: MatLoad(A, fd);
24: PetscViewerDestroy(&fd);
25: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
27: MatCreateVecs(A, &b, &x);
28: VecSetRandom(b, NULL);
29: VecDuplicateVecs(b, N, &S);
30: PetscMalloc2(N, &tetar, N, &tetai);
32: KSPCreate(PETSC_COMM_WORLD, &ksp);
33: KSPSetType(ksp, KSPGMRES);
34: KSPSetTolerances(ksp, 10000 * PETSC_MACHINE_EPSILON, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
35: KSPGetPC(ksp, &pc);
36: PCSetType(pc, PCNONE);
37: KSPSetComputeRitz(ksp, PETSC_TRUE);
38: KSPSetOperators(ksp, A, A);
39: KSPSetFromOptions(ksp);
40: KSPSolve(ksp, b, x);
41: PetscOptionsHasName(NULL, NULL, "-harmonic", &harmonic);
42: KSPComputeRitz(ksp, harmonic ? PETSC_FALSE : PETSC_TRUE, PETSC_TRUE, &Na, S, tetar, tetai);
44: PetscPrintf(PETSC_COMM_WORLD, "%% Number of Ritz pairs %" PetscInt_FMT "\n", Na);
45: for (i = 0; i < Na; i++) {
46: PetscPrintf(PETSC_COMM_WORLD, "%% Eigenvalue(s) %g ", (double)tetar[i]);
47: if (tetai[i]) {
48: PetscPrintf(PETSC_COMM_WORLD, "%+gi", (double)tetai[i]);
49: #if !defined(PETSC_USE_COMPLEX)
50: PetscPrintf(PETSC_COMM_WORLD, " %g ", (double)tetar[i]);
51: PetscPrintf(PETSC_COMM_WORLD, "%+gi", (double)-tetai[i]);
52: #endif
53: }
54: PetscPrintf(PETSC_COMM_WORLD, "\n");
55: PetscPrintf(PETSC_COMM_WORLD, "%% Eigenvector\n");
56: VecView(S[i], PETSC_VIEWER_STDOUT_WORLD);
57: #if !defined(PETSC_USE_COMPLEX)
58: if (tetai[i]) {
59: PetscPrintf(PETSC_COMM_WORLD, "%% Imaginary part of Eigenvector\n");
60: VecView(S[i + 1], PETSC_VIEWER_STDOUT_WORLD);
61: i++;
62: }
63: #endif
64: }
66: PetscFree2(tetar, tetai);
67: VecDestroy(&x);
68: VecDestroy(&b);
69: VecDestroyVecs(N, &S);
70: KSPDestroy(&ksp);
71: MatDestroy(&A);
72: 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*/