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> */
5
6 #include <petscmat.h>
7
main(int argc,char ** args)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;
15
16 PetscFunctionBeginUser;
17 PetscCall(PetscInitialize(&argc, &args, NULL, 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.");
21
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));
26
27 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A1));
28 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
29
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));
34
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));
39
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));
44
45 /* dstMat += A */
46 PetscCall(MatAXPY(dstMat, 1.0, A, SUBSET_NONZERO_PATTERN));
47
48 PetscCall(MatDestroy(&A));
49 PetscCall(MatDestroy(&A1));
50 PetscCall(MatDestroy(&A2));
51 PetscCall(MatDestroy(&dstMat));
52 PetscCall(PetscFinalize());
53 return 0;
54 }
55