1 static char help[] = "Test MatGetInertia().\n\n";
2
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 */
8
9 #include <petscksp.h>
main(int argc,char ** args)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;
22
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));
34
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 }
44
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));
54
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)); */
83
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));
91
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)); */
100
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 }
108
109 /* Test MatGetInertia() */
110 /* if A is symmetric, set its flag -- required by MatGetInertia() */
111 PetscCall(MatIsSymmetric(A, 0.0, &flag));
112
113 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
114 PetscCall(KSPSetType(ksp, KSPPREONLY));
115 PetscCall(KSPSetOperators(ksp, A, A));
116
117 PetscCall(KSPGetPC(ksp, &pc));
118 PetscCall(PCSetType(pc, PCCHOLESKY));
119 PetscCall(PCSetFromOptions(pc));
120
121 PetscCall(PCSetUp(pc));
122 PetscCall(PCFactorGetMatrix(pc, &F));
123 PetscCall(MatGetInertia(F, &nneg, &nzero, &npos));
124
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));
127
128 /* Destroy */
129 PetscCall(KSPDestroy(&ksp));
130 PetscCall(MatDestroy(&A));
131 PetscCall(MatDestroy(&B));
132 PetscCall(PetscFinalize());
133 return 0;
134 }
135
136 /*TEST
137
138 test:
139 args: -sigma 2.0 -mat_type {{aij baij}}
140 requires: !complex
141 output_file: output/ex33.out
142
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
148
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
155
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
161
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
167
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
174
175 TEST*/
176