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