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_SQRT_MACHINE_EPSILON;
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) {
54: if (issbaij && PetscDefined(USE_COMPLEX)) continue;
55: stype = MATSOLVERMKL_CPARDISO;
56: }
57: #else
58: if (j == 1) continue;
59: #endif
60: if (issbaij) {
61: PetscCall(MatGetFactor(A, stype, MAT_FACTOR_CHOLESKY, &F));
62: PetscCall(MatCholeskyFactorSymbolic(F, A, NULL, NULL));
63: PetscCall(MatCholeskyFactorNumeric(F, A, NULL));
64: } else {
65: PetscCall(MatGetFactor(A, stype, MAT_FACTOR_LU, &F));
66: PetscCall(MatLUFactorSymbolic(F, A, NULL, NULL, NULL));
67: PetscCall(MatLUFactorNumeric(F, A, NULL));
68: }
69: for (i = 0; i < 10; i++) {
70: PetscCall(VecSetRandom(b, rdm));
71: PetscCall(MatSolve(F, b, y));
72: /* Check the error */
73: PetscCall(MatMult(A, y, x));
74: PetscCall(VecAXPY(x, -1.0, b));
75: PetscCall(VecNorm(x, NORM_2, &norm2));
76: if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error:MatSolve(), norm2: %g\n", (double)norm2));
77: }
78: PetscCall(MatDestroy(&F));
79: }
80: PetscCall(VecDestroy(&x));
81: PetscCall(VecDestroy(&y));
82: PetscCall(VecDestroy(&b));
83: PetscCall(PetscRandomDestroy(&rdm));
84: #endif
85: PetscCall(MatDestroy(&A));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: int main(int argc, char *argv[])
90: {
91: MPI_Comm comm;
92: PetscMPIInt size;
94: PetscFunctionBeginUser;
95: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
96: comm = PETSC_COMM_WORLD;
97: PetscCallMPI(MPI_Comm_size(comm, &size));
98: PetscCheck(size == 2, comm, PETSC_ERR_USER, "This example must be run with exactly two processes");
99: PetscCall(Assemble(comm, 2, MATMPIBAIJ));
100: PetscCall(Assemble(comm, 2, MATMPISBAIJ));
101: PetscCall(Assemble(comm, 1, MATMPIBAIJ));
102: PetscCall(Assemble(comm, 1, MATMPISBAIJ));
103: PetscCall(PetscFinalize());
104: return 0;
105: }
107: /*TEST
109: test:
110: nsize: 2
111: args: -mat_ignore_lower_triangular
112: filter: sed -e "s~mem [0-9]*~mem~g"
114: TEST*/