Actual source code: ex113.c

  1: static char help[] = "Tests sequential and parallel MatMatMult() and MatAXPY(...,SUBSET_NONZERO_PATTERN) \n\
  2: Input arguments are:\n\
  3:   -f <input_file>  : file to load\n\n";
  4: /* e.g., mpiexec -n 3 ./ex113 -f <file> */

  6: #include <petscmat.h>

  8: int main(int argc, char **args)
  9: {
 10:   Mat         A, A1, A2, Mtmp, dstMat;
 11:   PetscViewer viewer;
 12:   PetscReal   fill = 4.0;
 13:   char        file[128];
 14:   PetscBool   flg;

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 18:   /*  Load the matrix A */
 19:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
 20:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate a file name for matrix A with the -f option.");

 22:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &viewer));
 23:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 24:   PetscCall(MatLoad(A, viewer));
 25:   PetscCall(PetscViewerDestroy(&viewer));

 27:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A1));
 28:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));

 30:   /* dstMat = A*A1*A2 */
 31:   PetscCall(MatMatMult(A1, A2, MAT_INITIAL_MATRIX, fill, &Mtmp));
 32:   PetscCall(MatMatMult(A, Mtmp, MAT_INITIAL_MATRIX, fill, &dstMat));
 33:   PetscCall(MatDestroy(&Mtmp));

 35:   /* dstMat += A1*A2 */
 36:   PetscCall(MatMatMult(A1, A2, MAT_INITIAL_MATRIX, fill, &Mtmp));
 37:   PetscCall(MatAXPY(dstMat, 1.0, Mtmp, SUBSET_NONZERO_PATTERN));
 38:   PetscCall(MatDestroy(&Mtmp));

 40:   /* dstMat += A*A1 */
 41:   PetscCall(MatMatMult(A, A1, MAT_INITIAL_MATRIX, fill, &Mtmp));
 42:   PetscCall(MatAXPY(dstMat, 1.0, Mtmp, SUBSET_NONZERO_PATTERN));
 43:   PetscCall(MatDestroy(&Mtmp));

 45:   /* dstMat += A */
 46:   PetscCall(MatAXPY(dstMat, 1.0, A, SUBSET_NONZERO_PATTERN));

 48:   PetscCall(MatDestroy(&A));
 49:   PetscCall(MatDestroy(&A1));
 50:   PetscCall(MatDestroy(&A2));
 51:   PetscCall(MatDestroy(&dstMat));
 52:   PetscCall(PetscFinalize());
 53:   return 0;
 54: }