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