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, (char *)0, 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(PetscObjectTypeCompare((PetscObject)A, MATSEQBAIJ, &flag));
112: if (!flag) PetscCall(MatIsSymmetric(A, 0.0, &flag));
114: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
115: PetscCall(KSPSetType(ksp, KSPPREONLY));
116: PetscCall(KSPSetOperators(ksp, A, A));
118: PetscCall(KSPGetPC(ksp, &pc));
119: PetscCall(PCSetType(pc, PCCHOLESKY));
120: PetscCall(PCSetFromOptions(pc));
122: PetscCall(PCSetUp(pc));
123: PetscCall(PCFactorGetMatrix(pc, &F));
124: PetscCall(MatGetInertia(F, &nneg, &nzero, &npos));
126: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
127: if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos));
129: /* Destroy */
130: PetscCall(KSPDestroy(&ksp));
131: PetscCall(MatDestroy(&A));
132: PetscCall(MatDestroy(&B));
133: PetscCall(PetscFinalize());
134: return 0;
135: }
137: /*TEST
139: test:
140: args: -sigma 2.0 -mat_type {{aij baij}}
141: requires: !complex
142: output_file: output/ex33.out
144: test:
145: suffix: mumps
146: args: -sigma 2.0 -pc_factor_mat_solver_type mumps -mat_mumps_icntl_13 1 -mat_mumps_icntl_24 1
147: requires: mumps !complex
148: output_file: output/ex33.out
150: test:
151: suffix: mumps_2
152: args: -sigma 2.0 -pc_factor_mat_solver_type mumps -mat_mumps_icntl_13 1 -mat_mumps_icntl_24 1
153: requires: mumps !complex
154: nsize: 3
155: output_file: output/ex33.out
157: test:
158: suffix: mkl_pardiso
159: args: -sigma 2.0 -pc_factor_mat_solver_type mkl_pardiso -mat_type sbaij
160: requires: mkl_pardiso !complex
161: output_file: output/ex33.out
163: test:
164: suffix: superlu_dist
165: args: -sigma 2.0 -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
166: requires: superlu_dist !complex
167: output_file: output/ex33.out
169: test:
170: suffix: superlu_dist_2
171: args: -sigma 2.0 -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
172: requires: superlu_dist !complex
173: nsize: 3
174: output_file: output/ex33.out
176: TEST*/