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