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 = 10 * 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*/