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*/