Actual source code: ex165.c


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

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

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

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

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

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