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, (char *)0, 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(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")); 32 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; 44 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)); 51 52 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ns", &ns, NULL)); 53 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)); 58 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 } 74 75 PetscCall(MatFactorCreateSchurComplement(F, &S, NULL)); 76 PetscCall(MatViewFromOptions(S, NULL, "-S_view")); 77 PetscCall(MatDestroy(&S)); 78 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)); 85 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 } 107 108 /*TEST 109 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}} 115 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}} 121 122 TEST*/ 123