Actual source code: ex134.c
1: static const char help[] = "Test parallel assembly of SBAIJ matrices\n\n";
3: #include <petscmat.h>
5: PetscErrorCode Assemble(MPI_Comm comm, PetscInt bs, MatType mtype)
6: {
7: const PetscInt rc[] = {0, 1, 2, 3};
8: const PetscScalar vals[] = {100, 2, 3, 4, 5, 600, 7, 8, 9, 100, 11, 1200, 13, 14, 15, 1600, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 2800, 29, 30, 31, 32,
9: 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 49, 60, 61, 62, 63, 64};
10: Mat A;
11: #if defined(PETSC_HAVE_MUMPS) || defined(PETSC_HAVE_MKL_CPARDISO)
12: Mat F;
13: MatSolverType stype = MATSOLVERPETSC;
14: PetscRandom rdm;
15: Vec b, x, y;
16: PetscInt i, j;
17: PetscReal norm2, tol = 100 * PETSC_SMALL;
18: PetscBool issbaij;
19: #endif
20: PetscViewer viewer;
22: PetscFunctionBegin;
23: PetscCall(MatCreate(comm, &A));
24: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 4 * bs, 4 * bs));
25: PetscCall(MatSetType(A, mtype));
26: PetscCall(MatMPIBAIJSetPreallocation(A, bs, 2, NULL, 2, NULL));
27: PetscCall(MatMPISBAIJSetPreallocation(A, bs, 2, NULL, 2, NULL));
28: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
29: /* All processes contribute a global matrix */
30: PetscCall(MatSetValuesBlocked(A, 4, rc, 4, rc, vals, ADD_VALUES));
31: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
32: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
33: PetscCall(PetscPrintf(comm, "Matrix %s(%" PetscInt_FMT ")\n", mtype, bs));
34: PetscCall(PetscViewerASCIIGetStdout(comm, &viewer));
35: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO_DETAIL));
36: PetscCall(MatView(A, viewer));
37: PetscCall(PetscViewerPopFormat(viewer));
38: PetscCall(MatView(A, viewer));
39: #if defined(PETSC_HAVE_MUMPS) || defined(PETSC_HAVE_MKL_CPARDISO)
40: PetscCall(PetscStrcmp(mtype, MATMPISBAIJ, &issbaij));
41: if (!issbaij) PetscCall(MatShift(A, 10));
42: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
43: PetscCall(PetscRandomSetFromOptions(rdm));
44: PetscCall(MatCreateVecs(A, &x, &y));
45: PetscCall(VecDuplicate(x, &b));
46: for (j = 0; j < 2; j++) {
47: #if defined(PETSC_HAVE_MUMPS)
48: if (j == 0) stype = MATSOLVERMUMPS;
49: #else
50: if (j == 0) continue;
51: #endif
52: #if defined(PETSC_HAVE_MKL_CPARDISO)
53: if (j == 1) stype = MATSOLVERMKL_CPARDISO;
54: #else
55: if (j == 1) continue;
56: #endif
57: if (issbaij) {
58: PetscCall(MatGetFactor(A, stype, MAT_FACTOR_CHOLESKY, &F));
59: PetscCall(MatCholeskyFactorSymbolic(F, A, NULL, NULL));
60: PetscCall(MatCholeskyFactorNumeric(F, A, NULL));
61: } else {
62: PetscCall(MatGetFactor(A, stype, MAT_FACTOR_LU, &F));
63: PetscCall(MatLUFactorSymbolic(F, A, NULL, NULL, NULL));
64: PetscCall(MatLUFactorNumeric(F, A, NULL));
65: }
66: for (i = 0; i < 10; i++) {
67: PetscCall(VecSetRandom(b, rdm));
68: PetscCall(MatSolve(F, b, y));
69: /* Check the error */
70: PetscCall(MatMult(A, y, x));
71: PetscCall(VecAXPY(x, -1.0, b));
72: PetscCall(VecNorm(x, NORM_2, &norm2));
73: if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error:MatSolve(), norm2: %g\n", (double)norm2));
74: }
75: PetscCall(MatDestroy(&F));
76: }
77: PetscCall(VecDestroy(&x));
78: PetscCall(VecDestroy(&y));
79: PetscCall(VecDestroy(&b));
80: PetscCall(PetscRandomDestroy(&rdm));
81: #endif
82: PetscCall(MatDestroy(&A));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: int main(int argc, char *argv[])
87: {
88: MPI_Comm comm;
89: PetscMPIInt size;
91: PetscFunctionBeginUser;
92: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
93: comm = PETSC_COMM_WORLD;
94: PetscCallMPI(MPI_Comm_size(comm, &size));
95: PetscCheck(size == 2, comm, PETSC_ERR_USER, "This example must be run with exactly two processes");
96: PetscCall(Assemble(comm, 2, MATMPIBAIJ));
97: PetscCall(Assemble(comm, 2, MATMPISBAIJ));
98: PetscCall(Assemble(comm, 1, MATMPIBAIJ));
99: PetscCall(Assemble(comm, 1, MATMPISBAIJ));
100: PetscCall(PetscFinalize());
101: return 0;
102: }
104: /*TEST
106: test:
107: nsize: 2
108: args: -mat_ignore_lower_triangular
109: filter: sed -e "s~mem [0-9]*~mem~g"
111: TEST*/