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;
23: PetscFunctionBeginUser;
24: PetscCall(PetscInitialize(&argc, &args, NULL, help));
25: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
26: Compute the matrices that define the eigensystem, Ax=kBx
27: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
28: PetscCall(PetscOptionsGetString(NULL, NULL, "-fA", file[0], sizeof(file[0]), &loadA));
29: if (loadA) {
30: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[0], FILE_MODE_READ, &viewer));
31: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
32: PetscCall(MatLoad(A, viewer));
33: PetscCall(PetscViewerDestroy(&viewer));
35: PetscCall(PetscOptionsGetString(NULL, NULL, "-fB", file[1], sizeof(file[1]), &loadB));
36: if (loadB) {
37: /* load B to get A = A + sigma*B */
38: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[1], FILE_MODE_READ, &viewer));
39: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
40: PetscCall(MatLoad(B, viewer));
41: PetscCall(PetscViewerDestroy(&viewer));
42: }
43: }
45: if (!loadA) { /* Matrix A is copied from slepc-3.0.0-p6/src/examples/ex13.c. */
46: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
47: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, &flag));
48: if (!flag) m = n;
49: N = n * m;
50: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
51: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N));
52: PetscCall(MatSetFromOptions(A));
53: PetscCall(MatSetUp(A));
55: PetscCall(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: PetscCall(MatSetValues(A, 1, &II, 1, &J, &v, INSERT_VALUES));
63: }
64: if (i < m - 1) {
65: J = II + n;
66: PetscCall(MatSetValues(A, 1, &II, 1, &J, &v, INSERT_VALUES));
67: }
68: if (j > 0) {
69: J = II - 1;
70: PetscCall(MatSetValues(A, 1, &II, 1, &J, &v, INSERT_VALUES));
71: }
72: if (j < n - 1) {
73: J = II + 1;
74: PetscCall(MatSetValues(A, 1, &II, 1, &J, &v, INSERT_VALUES));
75: }
76: v = 4.0;
77: PetscCall(MatSetValues(A, 1, &II, 1, &II, &v, INSERT_VALUES));
78: }
79: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
80: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
81: }
82: /* PetscCall(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); */
84: if (!loadB) {
85: PetscCall(MatGetLocalSize(A, &m, &n));
86: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
87: PetscCall(MatSetSizes(B, m, n, PETSC_DECIDE, PETSC_DECIDE));
88: PetscCall(MatSetFromOptions(B));
89: PetscCall(MatSetUp(B));
90: PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
92: for (II = Istart; II < Iend; II++) {
93: v = 1.0;
94: PetscCall(MatSetValues(B, 1, &II, 1, &II, &v, INSERT_VALUES));
95: }
96: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
97: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
98: }
99: /* PetscCall(MatView(B,PETSC_VIEWER_STDOUT_WORLD)); */
101: /* Set a shift: A = A - sigma*B */
102: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigma", &sigma, &flag));
103: if (flag) {
104: sigma = -1.0 * sigma;
105: PetscCall(MatAXPY(A, sigma, B, DIFFERENT_NONZERO_PATTERN)); /* A <- A - sigma*B */
106: /* PetscCall(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); */
107: }
109: /* Test MatGetInertia() */
110: /* if A is symmetric, set its flag -- required by MatGetInertia() */
111: PetscCall(MatIsSymmetric(A, 0.0, &flag));
113: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
114: PetscCall(KSPSetType(ksp, KSPPREONLY));
115: PetscCall(KSPSetOperators(ksp, A, A));
117: PetscCall(KSPGetPC(ksp, &pc));
118: PetscCall(PCSetType(pc, PCCHOLESKY));
119: PetscCall(PCSetFromOptions(pc));
121: PetscCall(PCSetUp(pc));
122: PetscCall(PCFactorGetMatrix(pc, &F));
123: PetscCall(MatGetInertia(F, &nneg, &nzero, &npos));
125: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
126: if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos));
128: /* Destroy */
129: PetscCall(KSPDestroy(&ksp));
130: PetscCall(MatDestroy(&A));
131: PetscCall(MatDestroy(&B));
132: PetscCall(PetscFinalize());
133: return 0;
134: }
136: /*TEST
138: test:
139: args: -sigma 2.0 -mat_type {{aij baij}}
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*/