1 static char help[] = "Tests MatMatSolve() in Schur complement mode.\n\n"; 2 3 #include <petscmat.h> 4 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; 17 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"); 22 23 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL)); 24 25 PetscCall(PetscOptionsGetString(NULL, NULL, "-A", file, sizeof(file), &flg)); 26 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Must provide a binary matrix with -A filename option"); 27 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd)); 28 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 29 PetscCall(MatLoad(A, fd)); 30 PetscCall(PetscViewerDestroy(&fd)); 31 PetscCall(MatGetSize(A, &m, &n)); 32 PetscCheck(m == n, PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n); 33 PetscCall(MatViewFromOptions(A, NULL, "-A_view")); 34 35 PetscCall(PetscOptionsGetString(NULL, NULL, "-B", file, sizeof(file), &flg)); 36 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Must provide a binary matrix with -B filename option"); 37 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd)); 38 PetscCall(MatCreate(PETSC_COMM_WORLD, &B)); 39 PetscCall(MatLoad(B, fd)); 40 PetscCall(PetscViewerDestroy(&fd)); 41 PetscCall(MatViewFromOptions(B, NULL, "-B_view")); 42 PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)B, &flg, MATSEQDENSE, MATMPIDENSE, NULL)); 43 if (!flg) PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &flg)); 44 if (!flg) { 45 Mat Bt; 46 47 PetscCall(MatCreateTranspose(B, &Bt)); 48 PetscCall(MatDestroy(&B)); 49 B = Bt; 50 } 51 PetscCall(PetscOptionsGetString(NULL, NULL, "-B_convert_type", converttype, sizeof(converttype), &flg)); 52 if (flg) PetscCall(MatConvert(B, converttype, MAT_INPLACE_MATRIX, &B)); 53 54 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ns", &ns, NULL)); 55 56 PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", solver, sizeof(solver), &flg)); 57 if (!flg) PetscCall(PetscStrncpy(solver, MATSOLVERMUMPS, sizeof(solver))); 58 PetscCall(PetscOptionsGetEnum(NULL, NULL, "-mat_factor_type", MatFactorTypes, (PetscEnum *)&factor, NULL)); 59 PetscCall(MatGetFactor(A, solver, factor, &F)); 60 61 PetscCall(ISCreateStride(PETSC_COMM_SELF, ns, m - ns, 1, &is_schur)); 62 PetscCall(MatFactorSetSchurIS(F, is_schur)); 63 PetscCall(ISDestroy(&is_schur)); 64 switch (factor) { 65 case MAT_FACTOR_LU: 66 PetscCall(MatLUFactorSymbolic(F, A, NULL, NULL, NULL)); 67 PetscCall(MatLUFactorNumeric(F, A, NULL)); 68 break; 69 case MAT_FACTOR_CHOLESKY: 70 PetscCall(MatCholeskyFactorSymbolic(F, A, NULL, NULL)); 71 PetscCall(MatCholeskyFactorNumeric(F, A, NULL)); 72 break; 73 default: 74 PetscCheck(PETSC_FALSE, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not coded for factor type %s", MatFactorTypes[factor]); 75 } 76 77 PetscCall(MatFactorCreateSchurComplement(F, &S, NULL)); 78 PetscCall(MatViewFromOptions(S, NULL, "-S_view")); 79 PetscCall(MatDestroy(&S)); 80 81 PetscCall(MatGetSize(B, NULL, &n)); 82 PetscCall(MatCreate(PETSC_COMM_WORLD, &X)); 83 PetscCall(MatSetSizes(X, m, PETSC_DECIDE, PETSC_DECIDE, n)); 84 PetscCall(MatSetType(X, MATDENSE)); 85 PetscCall(MatSetFromOptions(X)); 86 PetscCall(MatSetUp(X)); 87 88 PetscCall(MatMatSolve(F, B, X)); 89 PetscCall(MatViewFromOptions(X, NULL, "-X_view")); 90 PetscCall(MatMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Y)); 91 PetscCall(MatViewFromOptions(Y, NULL, "-Y_view")); 92 PetscCall(MatAXPY(Y, -1.0, B, SAME_NONZERO_PATTERN)); 93 PetscCall(MatViewFromOptions(Y, NULL, "-err_view")); 94 PetscCall(MatNorm(Y, NORM_FROBENIUS, &norm)); 95 if (norm > tol) { 96 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve: Norm of error %g\n", (double)norm)); 97 PetscCall(MatConvert(Y, MATAIJ, MAT_INPLACE_MATRIX, &Y)); 98 PetscCall(MatFilter(Y, PETSC_SMALL, PETSC_TRUE, PETSC_FALSE)); 99 PetscCall(MatViewFromOptions(Y, NULL, "-aij_err_view")); 100 } 101 PetscCall(MatDestroy(&A)); 102 PetscCall(MatDestroy(&X)); 103 PetscCall(MatDestroy(&F)); 104 PetscCall(MatDestroy(&B)); 105 PetscCall(MatDestroy(&Y)); 106 PetscCall(PetscFinalize()); 107 return 0; 108 } 109 110 /*TEST 111 112 testset: 113 output_file: output/empty.out 114 requires: datafilespath mumps double !complex !defined(PETSC_USE_64BIT_INDICES) 115 args: -A ${DATAFILESPATH}/matrices/factorSchur/A.dat -ns {{0 1}} 116 117 test: 118 suffix: mumps_1 119 args: -B ${DATAFILESPATH}/matrices/factorSchur/B1.dat 120 121 test: 122 suffix: mumps_2 123 args: -B ${DATAFILESPATH}/matrices/factorSchur/B2.dat 124 125 testset: 126 output_file: output/empty.out 127 requires: datafilespath mumps double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_MUMPS_MIXED_PRECISION) 128 args: -A ${DATAFILESPATH}/matrices/factorSchur/A.dat -ns {{0 1}} -pc_precision single -tol 1e-4 129 130 test: 131 suffix: mumps_1_single 132 args: -B ${DATAFILESPATH}/matrices/factorSchur/B1.dat 133 134 test: 135 suffix: mumps_2_single 136 args: -B ${DATAFILESPATH}/matrices/factorSchur/B2.dat 137 138 TEST*/ 139