Actual source code: ex33.c
1: static char help[] = "Test MatGetInertia().\n\n";
3: /*
4: Examples of command line options:
5: ./ex33 -sigma 2.0 -pc_factor_mat_solver_type mumps -mat_mumps_icntl_13 1 -mat_mumps_icntl_24 1
6: ./ex33 -sigma <shift> -fA <matrix_file>
7: */
9: #include <petscksp.h>
10: int main(int argc, char **args)
11: {
12: Mat A, B, F;
13: KSP ksp;
14: PC pc;
15: PetscInt N, n = 10, m, Istart, Iend, II, J, i, j;
16: PetscInt nneg, nzero, npos;
17: PetscScalar v, sigma;
18: PetscBool flag, loadA = PETSC_FALSE, loadB = PETSC_FALSE;
19: char file[2][PETSC_MAX_PATH_LEN];
20: PetscViewer viewer;
21: PetscMPIInt rank;
24: PetscInitialize(&argc, &args, (char *)0, help);
25: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
26: Compute the matrices that define the eigensystem, Ax=kBx
27: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
28: PetscOptionsGetString(NULL, NULL, "-fA", file[0], sizeof(file[0]), &loadA);
29: if (loadA) {
30: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[0], FILE_MODE_READ, &viewer);
31: MatCreate(PETSC_COMM_WORLD, &A);
32: MatLoad(A, viewer);
33: PetscViewerDestroy(&viewer);
35: PetscOptionsGetString(NULL, NULL, "-fB", file[1], sizeof(file[1]), &loadB);
36: if (loadB) {
37: /* load B to get A = A + sigma*B */
38: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[1], FILE_MODE_READ, &viewer);
39: MatCreate(PETSC_COMM_WORLD, &B);
40: MatLoad(B, viewer);
41: PetscViewerDestroy(&viewer);
42: }
43: }
45: if (!loadA) { /* Matrix A is copied from slepc-3.0.0-p6/src/examples/ex13.c. */
46: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
47: PetscOptionsGetInt(NULL, NULL, "-m", &m, &flag);
48: if (!flag) m = n;
49: N = n * m;
50: MatCreate(PETSC_COMM_WORLD, &A);
51: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N);
52: MatSetFromOptions(A);
53: MatSetUp(A);
55: MatGetOwnershipRange(A, &Istart, &Iend);
56: for (II = Istart; II < Iend; II++) {
57: v = -1.0;
58: i = II / n;
59: j = II - i * n;
60: if (i > 0) {
61: J = II - n;
62: MatSetValues(A, 1, &II, 1, &J, &v, INSERT_VALUES);
63: }
64: if (i < m - 1) {
65: J = II + n;
66: MatSetValues(A, 1, &II, 1, &J, &v, INSERT_VALUES);
67: }
68: if (j > 0) {
69: J = II - 1;
70: MatSetValues(A, 1, &II, 1, &J, &v, INSERT_VALUES);
71: }
72: if (j < n - 1) {
73: J = II + 1;
74: MatSetValues(A, 1, &II, 1, &J, &v, INSERT_VALUES);
75: }
76: v = 4.0;
77: MatSetValues(A, 1, &II, 1, &II, &v, INSERT_VALUES);
78: }
79: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
81: }
82: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
84: if (!loadB) {
85: MatGetLocalSize(A, &m, &n);
86: MatCreate(PETSC_COMM_WORLD, &B);
87: MatSetSizes(B, m, n, PETSC_DECIDE, PETSC_DECIDE);
88: MatSetFromOptions(B);
89: MatSetUp(B);
90: MatGetOwnershipRange(A, &Istart, &Iend);
92: for (II = Istart; II < Iend; II++) {
93: v = 1.0;
94: MatSetValues(B, 1, &II, 1, &II, &v, INSERT_VALUES);
95: }
96: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
97: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
98: }
99: /* MatView(B,PETSC_VIEWER_STDOUT_WORLD); */
101: /* Set a shift: A = A - sigma*B */
102: PetscOptionsGetScalar(NULL, NULL, "-sigma", &sigma, &flag);
103: if (flag) {
104: sigma = -1.0 * sigma;
105: MatAXPY(A, sigma, B, DIFFERENT_NONZERO_PATTERN); /* A <- A - sigma*B */
106: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
107: }
109: /* Test MatGetInertia() */
110: /* if A is symmetric, set its flag -- required by MatGetInertia() */
111: MatIsSymmetric(A, 0.0, &flag);
113: KSPCreate(PETSC_COMM_WORLD, &ksp);
114: KSPSetType(ksp, KSPPREONLY);
115: KSPSetOperators(ksp, A, A);
117: KSPGetPC(ksp, &pc);
118: PCSetType(pc, PCCHOLESKY);
119: PCSetFromOptions(pc);
121: PCSetUp(pc);
122: PCFactorGetMatrix(pc, &F);
123: MatGetInertia(F, &nneg, &nzero, &npos);
125: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
126: if (rank == 0) PetscPrintf(PETSC_COMM_SELF, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos);
128: /* Destroy */
129: KSPDestroy(&ksp);
130: MatDestroy(&A);
131: MatDestroy(&B);
132: PetscFinalize();
133: return 0;
134: }
136: /*TEST
138: test:
139: args: -sigma 2.0
140: requires: !complex
141: output_file: output/ex33.out
143: test:
144: suffix: mumps
145: args: -sigma 2.0 -pc_factor_mat_solver_type mumps -mat_mumps_icntl_13 1 -mat_mumps_icntl_24 1
146: requires: mumps !complex
147: output_file: output/ex33.out
149: test:
150: suffix: mumps_2
151: args: -sigma 2.0 -pc_factor_mat_solver_type mumps -mat_mumps_icntl_13 1 -mat_mumps_icntl_24 1
152: requires: mumps !complex
153: nsize: 3
154: output_file: output/ex33.out
156: test:
157: suffix: mkl_pardiso
158: args: -sigma 2.0 -pc_factor_mat_solver_type mkl_pardiso -mat_type sbaij
159: requires: mkl_pardiso !complex
160: output_file: output/ex33.out
162: test:
163: suffix: superlu_dist
164: args: -sigma 2.0 -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
165: requires: superlu_dist !complex
166: output_file: output/ex33.out
168: test:
169: suffix: superlu_dist_2
170: args: -sigma 2.0 -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
171: requires: superlu_dist !complex
172: nsize: 3
173: output_file: output/ex33.out
175: TEST*/