1 #include <petscmat.h>
2
3 static char help[] = "Example of MatMat ops with MatDense in PETSc.\n";
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Mat A, P, PtAP, RARt, ABC, Pt;
8 PetscInt n = 4, m = 2; // Example dimensions
9 PetscMPIInt size;
10
11 PetscCall(PetscInitialize(&argc, &args, NULL, help));
12 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
13 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This example is for sequential runs only.");
14
15 // Create dense matrix P (n x m)
16 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n, m, NULL, &P));
17 PetscScalar P_values[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
18 PetscInt rows[] = {0, 1, 2, 3};
19 PetscInt cols[] = {0, 1};
20 PetscCall(MatSetValues(P, n, rows, m, cols, P_values, INSERT_VALUES));
21 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
22 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
23 PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &Pt));
24
25 // Create dense matrix A (n x n)
26 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n, n, NULL, &A));
27 PetscCall(MatSetBlockSize(A, n)); // Set block size for A
28 PetscScalar A_values[] = {4.0, 1.0, 2.0, 0.0, 1.0, 3.0, 0.0, 1.0, 2.0, 0.0, 5.0, 2.0, 0.0, 1.0, 2.0, 6.0};
29 PetscInt indices[] = {0, 1, 2, 3};
30 PetscCall(MatSetValues(A, n, indices, n, indices, A_values, INSERT_VALUES));
31 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
32 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
33
34 PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &PtAP));
35 PetscCall(MatRARt(A, Pt, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &RARt));
36 PetscCall(MatMatMatMult(Pt, A, P, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &ABC));
37
38 // View matrices
39 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix P:\n"));
40 PetscCall(MatView(P, PETSC_VIEWER_STDOUT_SELF));
41
42 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix A:\n"));
43 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF));
44
45 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix PtAP:\n"));
46 PetscCall(MatView(PtAP, PETSC_VIEWER_STDOUT_SELF));
47
48 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix RARt:\n"));
49 PetscCall(MatView(RARt, PETSC_VIEWER_STDOUT_SELF));
50
51 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix ABC:\n"));
52 PetscCall(MatView(ABC, PETSC_VIEWER_STDOUT_SELF));
53
54 // Clean up
55 PetscCall(MatDestroy(&P));
56 PetscCall(MatDestroy(&Pt));
57 PetscCall(MatDestroy(&A));
58 PetscCall(MatDestroy(&PtAP));
59 PetscCall(MatDestroy(&RARt));
60 PetscCall(MatDestroy(&ABC));
61
62 PetscCall(PetscFinalize());
63 return 0;
64 }
65
66 /*TEST
67
68 test:
69 diff_args: -j
70 suffix: 1
71
72 TEST*/
73