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(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL));

 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"));

 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;

 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));

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

 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));

 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:   }

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

 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));

 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: }

110: /*TEST

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}}

117:     test:
118:       suffix: mumps_1
119:       args: -B ${DATAFILESPATH}/matrices/factorSchur/B1.dat

121:     test:
122:       suffix: mumps_2
123:       args: -B ${DATAFILESPATH}/matrices/factorSchur/B2.dat

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

130:     test:
131:       suffix: mumps_1_single
132:       args: -B ${DATAFILESPATH}/matrices/factorSchur/B1.dat

134:     test:
135:       suffix: mumps_2_single
136:       args: -B ${DATAFILESPATH}/matrices/factorSchur/B2.dat

138: TEST*/