Actual source code: ex154.c

  1: static char help[] = "Tests MatMatSolve() in Schur complement mode.\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **args)
  6: {
  7:   Mat           F, A, B, X, Y, S;
  8:   IS            is_schur;
  9:   PetscMPIInt   size;
 10:   PetscInt      ns = 0, m, n;
 11:   PetscReal     norm, tol = PETSC_SQRT_MACHINE_EPSILON;
 12:   MatFactorType factor = MAT_FACTOR_LU;
 13:   PetscViewer   fd;
 14:   char          solver[256], converttype[256];
 15:   char          file[PETSC_MAX_PATH_LEN];
 16:   PetscBool     flg;

 18:   PetscFunctionBeginUser;
 19:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 20:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 21:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor test");

 23:   PetscCall(PetscOptionsGetString(NULL, NULL, "-A", file, sizeof(file), &flg));
 24:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Must provide a binary matrix with -A filename option");
 25:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
 26:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 27:   PetscCall(MatLoad(A, fd));
 28:   PetscCall(PetscViewerDestroy(&fd));
 29:   PetscCall(MatGetSize(A, &m, &n));
 30:   PetscCheck(m == n, PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);
 31:   PetscCall(MatViewFromOptions(A, NULL, "-A_view"));

 33:   PetscCall(PetscOptionsGetString(NULL, NULL, "-B", file, sizeof(file), &flg));
 34:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Must provide a binary matrix with -B filename option");
 35:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
 36:   PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
 37:   PetscCall(MatLoad(B, fd));
 38:   PetscCall(PetscViewerDestroy(&fd));
 39:   PetscCall(MatViewFromOptions(B, NULL, "-B_view"));
 40:   PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)B, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
 41:   if (!flg) PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &flg));
 42:   if (!flg) {
 43:     Mat Bt;

 45:     PetscCall(MatCreateTranspose(B, &Bt));
 46:     PetscCall(MatDestroy(&B));
 47:     B = Bt;
 48:   }
 49:   PetscCall(PetscOptionsGetString(NULL, NULL, "-B_convert_type", converttype, sizeof(converttype), &flg));
 50:   if (flg) PetscCall(MatConvert(B, converttype, MAT_INPLACE_MATRIX, &B));

 52:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ns", &ns, NULL));

 54:   PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", solver, sizeof(solver), &flg));
 55:   if (!flg) PetscCall(PetscStrncpy(solver, MATSOLVERMUMPS, sizeof(solver)));
 56:   PetscCall(PetscOptionsGetEnum(NULL, NULL, "-mat_factor_type", MatFactorTypes, (PetscEnum *)&factor, NULL));
 57:   PetscCall(MatGetFactor(A, solver, factor, &F));

 59:   PetscCall(ISCreateStride(PETSC_COMM_SELF, ns, m - ns, 1, &is_schur));
 60:   PetscCall(MatFactorSetSchurIS(F, is_schur));
 61:   PetscCall(ISDestroy(&is_schur));
 62:   switch (factor) {
 63:   case MAT_FACTOR_LU:
 64:     PetscCall(MatLUFactorSymbolic(F, A, NULL, NULL, NULL));
 65:     PetscCall(MatLUFactorNumeric(F, A, NULL));
 66:     break;
 67:   case MAT_FACTOR_CHOLESKY:
 68:     PetscCall(MatCholeskyFactorSymbolic(F, A, NULL, NULL));
 69:     PetscCall(MatCholeskyFactorNumeric(F, A, NULL));
 70:     break;
 71:   default:
 72:     PetscCheck(PETSC_FALSE, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not coded for factor type %s", MatFactorTypes[factor]);
 73:   }

 75:   PetscCall(MatFactorCreateSchurComplement(F, &S, NULL));
 76:   PetscCall(MatViewFromOptions(S, NULL, "-S_view"));
 77:   PetscCall(MatDestroy(&S));

 79:   PetscCall(MatGetSize(B, NULL, &n));
 80:   PetscCall(MatCreate(PETSC_COMM_WORLD, &X));
 81:   PetscCall(MatSetSizes(X, m, PETSC_DECIDE, PETSC_DECIDE, n));
 82:   PetscCall(MatSetType(X, MATDENSE));
 83:   PetscCall(MatSetFromOptions(X));
 84:   PetscCall(MatSetUp(X));

 86:   PetscCall(MatMatSolve(F, B, X));
 87:   PetscCall(MatViewFromOptions(X, NULL, "-X_view"));
 88:   PetscCall(MatMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Y));
 89:   PetscCall(MatViewFromOptions(Y, NULL, "-Y_view"));
 90:   PetscCall(MatAXPY(Y, -1.0, B, SAME_NONZERO_PATTERN));
 91:   PetscCall(MatViewFromOptions(Y, NULL, "-err_view"));
 92:   PetscCall(MatNorm(Y, NORM_FROBENIUS, &norm));
 93:   if (norm > tol) {
 94:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve: Norm of error %g\n", (double)norm));
 95:     PetscCall(MatConvert(Y, MATAIJ, MAT_INPLACE_MATRIX, &Y));
 96:     PetscCall(MatFilter(Y, PETSC_SMALL, PETSC_TRUE, PETSC_FALSE));
 97:     PetscCall(MatViewFromOptions(Y, NULL, "-aij_err_view"));
 98:   }
 99:   PetscCall(MatDestroy(&A));
100:   PetscCall(MatDestroy(&X));
101:   PetscCall(MatDestroy(&F));
102:   PetscCall(MatDestroy(&B));
103:   PetscCall(MatDestroy(&Y));
104:   PetscCall(PetscFinalize());
105:   return 0;
106: }

108: /*TEST

110:   test:
111:     output_file: output/ex62_1.out
112:     suffix: mumps_1
113:     requires: datafilespath mumps double !complex !defined(PETSC_USE_64BIT_INDICES)
114:     args: -A ${DATAFILESPATH}/matrices/factorSchur/A.dat -B ${DATAFILESPATH}/matrices/factorSchur/B1.dat -ns {{0 1}}

116:   test:
117:     output_file: output/ex62_1.out
118:     suffix: mumps_2
119:     requires: datafilespath mumps double !complex !defined(PETSC_USE_64BIT_INDICES)
120:     args: -A ${DATAFILESPATH}/matrices/factorSchur/A.dat -B ${DATAFILESPATH}/matrices/factorSchur/B2.dat -ns {{0 1}}

122: TEST*/