Actual source code: ex165.c

  1: static char help[] = "Tests C=A^T*B via MatTranspose() and MatMatMult(). \n\
  2:                      Contributed by Alexander Grayver, Jan. 2012 \n\n";
  3: /* Example:
  4:   mpiexec -n <np> ./ex165 -fA A.dat -fB B.dat -view_C
  5:  */

  7: #include <petscmat.h>
  8: int main(int argc, char **args)
  9: {
 10:   Mat         A, AT, B, C;
 11:   PetscViewer viewer;
 12:   PetscBool   flg;
 13:   char        file[PETSC_MAX_PATH_LEN];

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 17:   PetscCall(PetscOptionsGetString(NULL, NULL, "-fA", file, sizeof(file), &flg));
 18:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Input fileA not specified");
 19:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &viewer));
 20:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 21:   PetscCall(MatSetType(A, MATAIJ));
 22:   PetscCall(MatLoad(A, viewer));
 23:   PetscCall(PetscViewerDestroy(&viewer));

 25:   PetscCall(PetscOptionsGetString(NULL, NULL, "-fB", file, sizeof(file), &flg));
 26:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Input fileB not specified");
 27:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &viewer));
 28:   PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
 29:   PetscCall(MatSetType(B, MATDENSE));
 30:   PetscCall(MatLoad(B, viewer));
 31:   PetscCall(PetscViewerDestroy(&viewer));

 33:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
 34:   PetscCall(MatMatMult(AT, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C));

 36:   PetscCall(PetscOptionsHasName(NULL, NULL, "-view_C", &flg));
 37:   if (flg) {
 38:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "C.dat", FILE_MODE_WRITE, &viewer));
 39:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
 40:     PetscCall(MatView(C, viewer));
 41:     PetscCall(PetscViewerPopFormat(viewer));
 42:     PetscCall(PetscViewerDestroy(&viewer));
 43:   }
 44:   PetscCall(MatDestroy(&A));
 45:   PetscCall(MatDestroy(&B));
 46:   PetscCall(MatDestroy(&AT));
 47:   PetscCall(MatDestroy(&C));
 48:   PetscCall(PetscFinalize());
 49:   return 0;
 50: }